diff --git a/.gitignore b/.gitignore index d5a86f8..023a632 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,5 @@ dist/ docs/build docs/source/_autosummary poetry.lock +uv.lock +sodym.egg-info diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 3825429..f4985a8 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,4 +1,11 @@ repos: + - repo: https://github.com/mwouts/jupytext + rev: v1.16.4 + hooks: + - id: jupytext + args: [--sync, --pipe, black] + additional_dependencies: + - black==22.3.0 # Matches hook - repo: https://github.com/psf/black rev: 22.3.0 hooks: diff --git a/examples/example1.ipynb b/examples/example1.ipynb index f2e42c1..2ecd111 100644 --- a/examples/example1.ipynb +++ b/examples/example1.ipynb @@ -2061,7 +2061,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": ".venv", "language": "python", "name": "python3" }, diff --git a/examples/example1.py b/examples/example1.py new file mode 100644 index 0000000..06a04d5 --- /dev/null +++ b/examples/example1.py @@ -0,0 +1,122 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.16.4 +# kernelspec: +# display_name: .venv +# language: python +# name: python3 +# --- + +# %% [markdown] +# # Example 1. System with two processes, two parameters, one material. +# +# A simple MFA system with one material, a time horizon of 30 years [1980-2010], two processes, and a time-dependent parameter is analysed. +# +# Diagram of the MFA system +# +# The model equations are as follows: +# * $a(t) = D(t)$ (exogenous input flow) +# * $d(t) = \alpha (t)\cdot b(t)$ (recovery efficiency parameter) +# * $a(t) + d(t) = b(t) $ (mass balance process 1) +# * $b(t) = c(t) + d(t) $ (mass balance process 2) +# +# From these equations the system solution follows: +# * $c(t) = a(t) = D(t) $ +# * $b(t) = \frac{1}{1-\alpha (t)}\cdot D(t) $ +# * $c(t) = \frac{\alpha}{1-\alpha (t)}\cdot D(t) $ +# + +# %% [markdown] +# ## 1. Load sodym and other useful packages + +# %% +import numpy as np +import plotly.express as px + +from sodym import ( + Dimension, DimensionSet, Parameter, Process, FlowDefinition, Flow, MFASystem, +) +from sodym.flow_helper import make_empty_flows + +# %% [markdown] +# ## 2. Load data +# Normally data would be loaded from a file / database, but since this is just a small example, the values are input directly into the code below. + +# %% +time = Dimension(name='Time', letter='t', items=list(range(1980,2011))) +elements = Dimension(name='Elements', letter='e', items=['single material', ]) +dimensions = DimensionSet(dim_list=[time, elements]) + +parameters = { + 'D': Parameter(name='inflow', dims=dimensions, values=np.arange(0, 31).reshape(31, 1)), + 'alpha': Parameter(name='recovery rate', dims=dimensions, values=np.arange(2, 33).reshape(31, 1)/34) +} + +processes = { + 'sysenv': Process(name='sysenv', id=0), + 'process 1': Process(name='process 1', id=1), + 'process 2': Process(name='process 2', id=2) +} + +# %% [markdown] +# ## 3. Define flows and initialise them with zero values. +# By defining them with the shape specified by the relevant dimensions, we can later ensure that when we update the values, these have the correct dimensions. +# Note that flows are automatically asigned their names based on the names of the processes they are connecting. + +# %% +flow_definitions = [ + FlowDefinition(from_process_name='sysenv', to_process_name='process 1', dim_letters=('t', 'e')), # input + FlowDefinition(from_process_name='process 1', to_process_name='process 2', dim_letters=('t', 'e')), # consumption + FlowDefinition(from_process_name='process 2', to_process_name='sysenv', dim_letters=('t', 'e')), # output + FlowDefinition(from_process_name='process 2', to_process_name='process 1', dim_letters=('t', 'e')), # recovered material +] +flows = make_empty_flows(processes=processes, flow_definitions=flow_definitions, dims=dimensions) + + +# %% [markdown] +# ## 4. Define the MFA System equations +# We define a class with our system equations in the compute method. Afterwards we create an instance of this class, using the input data defined above. The class (system equations) can then easily be reused with different input data. +# +# We just need to define the compute method with our system equations, as all the other things we need are inherited from the MFASystem class. + +# %% +class SimpleMFA(MFASystem): + + def compute(self): + self.flows['sysenv => process 1'][...] = self.parameters['D'] # the elipsis slice [...] ensures the dimensionality of the flow is not changed + self.flows['process 1 => process 2'][...] = 1 / (1 - self.parameters['alpha']) * self.parameters['D'] + self.flows['process 2 => sysenv'][...] = self.parameters['D'] + self.flows['process 2 => process 1'][...] = self.parameters['alpha'] / (1 - self.parameters['alpha']) * self.parameters['D'] + + + +# %% +mfa_example = SimpleMFA(dims=dimensions, processes=processes, parameters=parameters, flows=flows, stocks={}) +mfa_example.compute() + +# %% +flow_a = mfa_example.flows['sysenv => process 1'] +fig = px.line( + x=flow_a.dims['t'].items, + y=flow_a['single material'].values, + title=flow_a.name, + labels={'x': 'Year', 'y': 'Mt/yr'}, +) +fig.show() + +# %% +flow_b = mfa_example.flows['process 1 => process 2'] +fig = px.line( + x=flow_b.dims['t'].items, + y=flow_b['single material'].values, + title=flow_b.name, + labels={'x': 'Year', 'y': 'Mt/yr'}, +) +fig.show() + +# %% diff --git a/examples/example2.ipynb b/examples/example2.ipynb index 25b80be..af7e2e6 100644 --- a/examples/example2.ipynb +++ b/examples/example2.ipynb @@ -5161,7 +5161,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": ".venv", "language": "python", "name": "python3" }, diff --git a/examples/example2.py b/examples/example2.py new file mode 100644 index 0000000..a618106 --- /dev/null +++ b/examples/example2.py @@ -0,0 +1,280 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.16.4 +# kernelspec: +# display_name: .venv +# language: python +# name: python3 +# --- + +# %% [markdown] +# # Example 2. Alloying elements in recycling. +# +# A recycling system with two end-of-life (EoL) products, two scrap types, one secondary material, and several types of losses are studied. Three chemical elements are considered: iron, copper, and manganese. A time horizon of 30 years [1980-2010], five processes, and time-dependent parameters are analysed. The processes have element-specific yield factors, meaning that the loss rates depend on the chemical element considered. These values are given below. +# +# The research questions are: +# +# * How much copper accumulates in the secondary steel assuming that all available scrap is remelted? +# * How much manganese is lost in the remelting process assuming that all available scrap is remelted? +# * What is more effective in reducing the copper concentraction of secondary steel: A reduction of the shredding yield factor for copper from EoL machines into steel scrap of 25% or an increase in the EoL buildings flow by 25%? (All other variables and parameters remaining equal) +# +# Simple MFA system +# +# The model equations are as follows: +# * $F_{1\_3}(t,e) = \Gamma_1(e) \cdot F_{0\_1}(t,e) $ (shredder yield factor) +# * $F_{1\_0}(t,e) = (1 - \Gamma_1(e)) \cdot F_{0\_1}(t,e) $ (mass balance) +# +# * $F_{2\_3}(t,e) = \Gamma_2(e) \cdot F_{0\_2}(t,e) $ (demolition yield factor) +# * $F_{2\_4}(t,e) = (1 - \Gamma_2(e)) \cdot F_{0\_2}(t,e) $ (mass balance) +# +# * $F_{3\_0}(t,e) = \Gamma_3(e) \cdot (F_{1\_3}(t,e)+F_{2\_3}(t,e)) $ (remelting yield factor) +# * $F_{3\_5}(t,e) = (1 - \Gamma_3(e)) \cdot (F_{1\_3}(t,e)+F_{2\_3}(t,e)) $ (mass balance) +# +# Here the index letters t denote the model time and e the chemical element. + +# %% [markdown] +# ## 1. Load sodym and other useful packages + +# %% +from copy import deepcopy +import numpy as np +import pandas as pd + +from sodym.data_reader import DataReader +from sodym import ( + DimensionDefinition, Dimension, DimensionSet, ParameterDefinition, Parameter, Process, FlowDefinition, Flow, StockDefinition, MFASystem, NamedDimArray +) +from sodym.flow_helper import make_empty_flows +from sodym.stock_helper import make_empty_stocks +from sodym.flow_naming import process_names_with_arrow +from sodym.export.array_plotter import PlotlyArrayPlotter + +# %% [markdown] +# ## 2. Define the data requirements, flows, stocks and MFA system equations +# +# We define the dimensions that are relevant for our system and the model parameters, processes, stocks and flows. We further define a class with our system equations in the compute method. + +# %% +dimension_definitions = [ + DimensionDefinition(letter='t', name='Time', dtype=int), + DimensionDefinition(letter='e', name='Material', dtype=str) +] + +parameter_definitions = [ + ParameterDefinition(name='eol machines', dim_letters=('t', )), + ParameterDefinition(name='eol buildings', dim_letters=('t', )), + ParameterDefinition(name='composition eol machines', dim_letters=('e', )), + ParameterDefinition(name='composition eol buildings', dim_letters=('e', )), + ParameterDefinition(name='shredder yield', dim_letters=('e', )), + ParameterDefinition(name='demolition yield', dim_letters=('e', )), + ParameterDefinition(name='remelting yield', dim_letters=('e', )), +] + +# %% +process_names = ['sysenv', 'shredder', 'demolition', 'remelting', 'landfills', 'slag piles'] +processes = {name: Process(name=name, id=index) for index, name in enumerate(process_names)} + +# %% +flow_definitions = [ + FlowDefinition(from_process_name='sysenv', to_process_name='shredder', dim_letters=('t', 'e')), + FlowDefinition(from_process_name='sysenv', to_process_name='demolition', dim_letters=('t', 'e')), + FlowDefinition(from_process_name='shredder', to_process_name='remelting', dim_letters=('t', 'e')), # scrap type 1 + FlowDefinition(from_process_name='shredder', to_process_name='sysenv', dim_letters=('t', 'e')), # shredder residue + FlowDefinition(from_process_name='demolition', to_process_name='remelting', dim_letters=('t', 'e')), # scrap type 2 + FlowDefinition(from_process_name='demolition', to_process_name='landfills', dim_letters=('t', 'e')), # loss + FlowDefinition(from_process_name='remelting', to_process_name='slag piles', dim_letters=('t', 'e')), # secondary steel + FlowDefinition(from_process_name='remelting', to_process_name='sysenv', dim_letters=('t', 'e')), # slag +] + +# %% +stock_definitions = [ + StockDefinition(name='landfills', process='landfills', dim_letters=('t', 'e'),), + StockDefinition(name='slag piles', process='slag piles', dim_letters=('t', 'e'),), +] + + +# %% [markdown] +# +# We just need to define the compute method with our system equations, as all the other things we need are inherited from the MFASystem class. The flow names are generated from the processes each flow connects, in this case with the naming function `process_names_with_arrow`, which is passed to the flow initialization below. + +# %% +class SimpleMFA(MFASystem): + def compute(self): + self.flows['sysenv => shredder'][...] = self.parameters['eol machines'] * self.parameters['composition eol machines'] + self.flows['sysenv => demolition'][...] = self.parameters['eol buildings'] * self.parameters['composition eol buildings'] + self.flows['shredder => remelting'][...] = self.flows['sysenv => shredder'] * self.parameters['shredder yield'] + self.flows['shredder => sysenv'][...] = self.flows['sysenv => shredder'] * (1 - self.parameters['shredder yield']) + self.flows['demolition => remelting'][...] = self.flows['sysenv => demolition'] * self.parameters['demolition yield'] + self.flows['demolition => landfills'][...] = self.flows['sysenv => demolition'] * (1 - self.parameters['demolition yield']) + self.flows['remelting => sysenv'][...] = ( + self.flows['shredder => remelting'] + self.flows['demolition => remelting']) * self.parameters['remelting yield'] + self.flows['remelting => slag piles'][...] = ( + self.flows['shredder => remelting'] + self.flows['demolition => remelting']) * (1 - self.parameters['remelting yield']) + self.stocks['landfills'].inflow[...] = self.flows['demolition => landfills'] + self.stocks['landfills'].compute() + self.stocks['slag piles'].inflow[...] = self.flows['shredder => remelting'] + self.stocks['slag piles'].compute() + + + +# %% [markdown] +# ## 3. Load data and initialise flows and stocks +# Now that we have defined the MFA system and know what data we need, we can load the data. +# Even though this is only a small system, we will load the data from excel files, as an example for more complex systems with larger datasets. To do this data loading, we define a DataReader class. Such a class can be reused with different datasets of the same format by passing attributes, e.g. data paths, in the init function. +# +# The methods `read_dimensions` and `read_parameters` are already defined in the parent DataReader class, and loop over the methods `read_dimension` and `read_parameter_values` that we specify for our usecase here. + +# %% +class CustomDataReader(DataReader): + + def __init__(self, path_to_time_parameters, path_to_element_parameters): + self.time_parameters = pd.read_excel(path_to_time_parameters, index_col=0) + self.element_parameters = pd.read_excel(path_to_element_parameters, index_col=0) + + def read_dimension(self, dimension_definition: DimensionDefinition) -> Dimension: + if dimension_definition.letter == 't': + data = list(self.time_parameters.index) + elif dimension_definition.letter == 'e': + data = list(self.element_parameters.index) + return Dimension(name=dimension_definition.name, letter=dimension_definition.letter, items=data) + + def read_parameter_values(self, parameter: str, dims: DimensionSet) -> Parameter: + if dims.letters == ('t', ): + data = self.time_parameters[parameter].values + elif dims.letters == ('e', ): + data = self.element_parameters[parameter].values + return Parameter(dims=dims, values=data) + + +# %% +data_reader = CustomDataReader( + path_to_time_parameters='example2_temporal_parameters.xlsx', + path_to_element_parameters='example2_material_parameters.xlsx' +) +dimensions = data_reader.read_dimensions(dimension_definitions) +parameters = data_reader.read_parameters(parameter_definitions, dims=dimensions) + +# %% +flows = make_empty_flows(processes=processes, flow_definitions=flow_definitions, dims=dimensions, + naming=process_names_with_arrow) + +# %% +stocks = make_empty_stocks(stock_definitions, dims=dimensions, processes=processes) # flow-driven stock objects + +# %% [markdown] +# ## 4. Put the pieces together +# Create a SimpleMFA instance by passing the loaded dimension and parameter data, as well as the initialised flow and stock objects. Solve the system equations by running the `compute` method. + +# %% +mfa_example = SimpleMFA(dims=dimensions, processes=processes, parameters=parameters, flows=flows, stocks=stocks) +mfa_example.compute() + +# %% [markdown] +# ## 5. Results +# Here we answer the research questions from the beginning of the notebook. +# +# **How much copper accumulates in the secondary steel assuming that all available scrap is remelted?** +# +# Clicking on the `Fe` entry of the plot legend hides it and adjusts the y-axis to better display the trace elements `Mn` and `Cu`. + +# %% +remelted = mfa_example.flows['remelting => sysenv'] + +plotter = PlotlyArrayPlotter( + array=remelted, + intra_line_dim='Time', + linecolor_dim='Material', + title="GDP-per-capita") +fig = plotter.plot(do_show=True) + + +# %% +remelted_shares = NamedDimArray(dims=remelted.dims) +remelted_shares[...] = remelted / remelted.sum_nda_over(('e',)) + +plotter = PlotlyArrayPlotter( + array=remelted_shares, + intra_line_dim='Time', + linecolor_dim='Material', + title="Share of copper and manganese in secondary steel") +fig = plotter.plot(do_show=True) + + +# %% [markdown] +# The copper flow in the secondary steel increases linearly from 0.34 kt/yr in 1980 to 0.78 kt/yr in 2010. The concentration of copper declines in a hyperbolic curve from 0.294% in 1980 to 0.233% in 2010. +# +# That concentration is below 0.4% at all times, the latter being the treshold for construction grade steel, but above 0.04%, which is the threshold for automotive steel. + +# %% [markdown] +# **How much manganese is lost in the remelting process assuming that all available scrap is remelted?** + +# %% +manganese_to_slag = mfa_example.flows['remelting => slag piles']['Mn'] + +plotter = PlotlyArrayPlotter( + array=manganese_to_slag, + intra_line_dim='Time', + ylabel='kt/yr', + title='Manganese lost in the remelting process') +fig = plotter.plot(do_show=True) + +# %% [markdown] +# **What is more effective in reducing the copper concentraction of secondary steel: A reduction of the shredding yield factor for copper from EoL machines into steel scrap of 25% or an increase in the EoL buildings flow by 25%? (All other variables and parameters remaining equal)** +# +# To answer this we change the parameter values and recalculate the entire system. +# In case a, we update the shredder yield, and in case b we increase the EoL buildings flow. +# We could load new datasets for the parameters, but since we are only changing one value, we will just update that value. +# + +# %% +parameters_a = deepcopy(parameters) +parameters_a['shredder yield'].set_values(np.array([0.92, 0.075, 0.92])) + +mfa_example_a = SimpleMFA(dims=dimensions, processes=processes, parameters=parameters_a, flows=deepcopy(flows), stocks=deepcopy(stocks)) +mfa_example_a.compute() + +# %% +parameters_b = deepcopy(parameters) +parameters_b['eol buildings'][...] = 1.25 * parameters_b['eol buildings'] + +mfa_example_b = SimpleMFA(dims=dimensions, processes=processes, parameters=parameters_b, flows=deepcopy(flows), stocks=deepcopy(stocks)) +mfa_example_b.compute() + +# %% +flow_a = mfa_example_a.flows['remelting => sysenv'] +shares_shredder = flow_a / flow_a.sum_nda_over(('e')) + +flow_b = mfa_example_b.flows['remelting => sysenv'] +shares_demolition = flow_b / flow_b.sum_nda_over(('e')) + +plotter = PlotlyArrayPlotter( + array=remelted_shares, + intra_line_dim='Time', + subplot_dim='Material', + line_label='Standard', + title='Material concentration in secondary steel') +fig = plotter.plot() +plotter = PlotlyArrayPlotter( + array=shares_shredder, + intra_line_dim='Time', + subplot_dim='Material', + line_label='Updated shredder yield', + fig=fig) +fig = plotter.plot() +plotter = PlotlyArrayPlotter( + array=shares_demolition, + intra_line_dim='Time', + subplot_dim='Material', + line_label='Increased buildings demolition', + fig=fig) +fig = plotter.plot(do_show=True) + +# %% [markdown] +# We can see that both measures reduce the copper concentration in the secondary steel. For the first year, the copper concentration is reduced from 0.294% to 0.244% if the Cu-yield into steel scrap of the shredder is reduced and to 0.259% if the EoL building flow treated is increased by 25%. The yield measure thus has a slightly higher impact on the copper contentration than the increase of a copper-poor scrap flow for dilution. In both cases the impact is not high enough to bring the copper concentration to values below 0.04%, which is necessary for automotive applications. + +# %% diff --git a/examples/example3.ipynb b/examples/example3.ipynb index 9559466..e92d30f 100644 --- a/examples/example3.ipynb +++ b/examples/example3.ipynb @@ -6389,7 +6389,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": ".venv", "language": "python", "name": "python3" }, @@ -6403,7 +6403,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.7" + "version": "3.10.10" } }, "nbformat": 4, diff --git a/examples/example3.py b/examples/example3.py new file mode 100644 index 0000000..beb8433 --- /dev/null +++ b/examples/example3.py @@ -0,0 +1,166 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.16.4 +# kernelspec: +# display_name: .venv +# language: python +# name: python3 +# --- + +# %% [markdown] +# # Example 3. Dynamic Stock modelling +# +# sodym defines the class DynamicStockModel for handling the inflow-driven and stock driven model of in-use stocks, see methods section 3 of the [uni-freiburg industrial ecology course](http://www.teaching.industrialecology.uni-freiburg.de/). In this notebook, we show how the dynamic stock model is used in the sodym framework. Other methods of the dynamic_stock_modelling class can be used in a similar way. +# +# The research question is: +# * How large are in-use stocks of steel in selected countries? +# * What is the ration between steel in EoL (end-of-life) products to final steel consumption in selected countries? +# To answer that question the system definition is chosen as in the figure below. +# +# Simple MFA system +# +# Stocks S and outflow O are calculated from apparent final consumption i(t), which is obtained from statistics, cf. DOI 10.1016/j.resconrec.2012.11.008 +# The model equations are as follows: +# +# First, we compute the outflow o_c(t,c) of each historic inflow/age-cohort i(c) in year t as +# $ o\_c(t,c) = i(c) \cdot sf(t,c) $ +# where sf is the survival function of the age cohort, which is 1-cdf, see the [wikipedia page on the survival function](https://en.wikipedia.org/wiki/Survival_function). +# The total outflow o(t) in a given year is then +# $ o(t) = \sum_{c\leq t} o\_c(t,c) $ +# The mass balance leads to the stock change $dS$: +# $ dS(t) = i(t) - o(t)$ +# And the stock finally is computed as +# $ S(t) = \sum_{t'\leq t} ds(t') $ + +# %% [markdown] +# ## 1. Load sodym and useful packages + +# %% +import numpy as np +import pandas as pd +import plotly.express as px + +from sodym import DimensionDefinition, ParameterDefinition, Dimension, DimensionSet, Parameter, Process, StockArray +from sodym.data_reader import DataReader +from sodym.survival_functions import NormalSurvival +from sodym.stocks import InflowDrivenDSM + +# %% [markdown] +# ## 2. Define system dimensions and load data +# +# First, we specify the dimensions that are relevant to our system. These will get passed to our data reader class and thereby we can ensure that the data we are reading has the correct shape. +# +# Even though this is only a small system, we will load data from an excel file, as an example for more complex systems with larger datasets. As mentioned above, we define a data reader class to do read the data and put it into the desired python objects. Such a class can be reused with different datasets of the same format by passing attributes, e.g. the file path, in the init function. + +# %% +dimension_definitions = [ + DimensionDefinition(letter='t', name='Time', dtype=int), + DimensionDefinition(letter='r', name='Region', dtype=str), +] + +parameter_definitions = [ + ParameterDefinition(name='inflow', dim_letters=('t', 'r', )), + ParameterDefinition(name='tau', dim_letters=('r', )), + ParameterDefinition(name='sigma', dim_letters=('r', )), +] + + +# %% +class LittleDataReader(DataReader): + + def __init__(self, country_lifetimes, steel_consumption_file): + self.country_lifetimes = country_lifetimes + self.steel_consumption = self.prepare_steel_consumption_data(steel_consumption_file) + + def prepare_steel_consumption_data(self, steel_consumption_file): + steel_consumption = pd.read_excel(steel_consumption_file) + steel_consumption = steel_consumption[['CS', 'T', 'V']] + return steel_consumption.rename(columns={'CS': 'r', 'T': 't', 'V': 'inflow'}) + + def read_dimension(self, dimension_definition: DimensionDefinition) -> Dimension: + if dimension_definition.letter == 't': + data = list(self.steel_consumption['t'].unique()) + elif dimension_definition.letter == 'r': + data = list(self.country_lifetimes.keys()) + return Dimension(name=dimension_definition.name, letter=dimension_definition.letter, items=data) + + def read_parameter_values(self, parameter: str, dims: DimensionSet) -> Parameter: + if parameter == 'tau': + data = np.array(list(country_lifetimes.values())) + elif parameter == 'sigma': + data = np.array([0.3 * lifetime for lifetime in country_lifetimes.values()]) + elif parameter == 'inflow': + multiindex = self.steel_consumption.set_index(['t', 'r']) + data = multiindex.unstack().values[:, :] + return Parameter(dims=dims, values=data) + +country_lifetimes = { + 'Argentina': 45, + 'Brazil': 25, + 'Canada': 35, + 'Denmark': 55, + 'Ethiopia': 70, + 'France': 45, + 'Greece': 70, + 'Hungary': 30, + 'Indonesia': 30, +} +data_reader = LittleDataReader(country_lifetimes=country_lifetimes, steel_consumption_file='example3_steel_consumption.xlsx') +dimensions = data_reader.read_dimensions(dimension_definitions) +parameters = data_reader.read_parameters(parameter_definitions, dimensions) + +# %% [markdown] +# ## 3. Perform dynamic stock modelling +# +# In this example, we do not need to build a whole MFA system, since we are only considering one dynamic stock. To make a dynamic stock in sodym, we first need to define a survival model; in this case we assume a normal distribution of lifetimes. Then, we can initiate the dynamic stock model. Here we choose an inflow driven stock model, because we have data that specifies the inflow and from this and the survival model we want to calculate the stock and the outflow. + +# %% +normal_survival_model = NormalSurvival( + dims=dimensions, + lifetime_mean=parameters['tau'], + lifetime_std=parameters['sigma'], +) + +inflow_stock = StockArray(dims=dimensions, values=parameters['inflow'].values) + +dynamic_stock = InflowDrivenDSM( + name='steel', + process=Process(name='in use', id=1), + survival_model=normal_survival_model, + inflow=inflow_stock, +) +dynamic_stock.compute() + +# %% [markdown] +# ## 4. Take a look at the results +# First, we plot the steel stock in the different countries over time. + +# %% +stock_df = dynamic_stock.stock.to_df(index=False) +stock_df = stock_df.pivot(index='Time', columns='Region')['value'] + +fig = px.line(stock_df, title='In-use stocks of steel') +fig.show() + +# %% [markdown] +# We then plot the ratio of outflow over inflow, which is a measure of the stationarity of a stock, and can be interpreted as one indicator for a circular economy. + +# %% +inflow = dynamic_stock.inflow +outflow = dynamic_stock.outflow +with np.errstate(divide='ignore'): + ratio_df = (outflow/inflow).to_df(index=False) +ratio_df = ratio_df.pivot(index='Time', columns='Region')['value'] + +fig = px.line(ratio_df, title='Ratio outflow:inflow') +fig.show() + +# %% [markdown] +# We see that for the rich countries France and Canada the share has been steadily growing since WW2. Upheavals such as wars and major economic crises can also be seen, in particular for Hungary. + +# %% diff --git a/examples/example5.ipynb b/examples/example5.ipynb index c4a7550..bb47f9e 100644 --- a/examples/example5.ipynb +++ b/examples/example5.ipynb @@ -4512,7 +4512,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": ".venv", "language": "python", "name": "python3" }, diff --git a/examples/example5.py b/examples/example5.py new file mode 100644 index 0000000..9d92c18 --- /dev/null +++ b/examples/example5.py @@ -0,0 +1,286 @@ +# %% [markdown] +# # Example 5. Estimating the material content of the global vehicle fleet +# +# This example shows a fully-fledged MFA system, designed to estimate the material composition of the global passenger vehicle fleet in 2017, covering 130 countries, 25 age-cohorts, and 25 materials. +# +# The research questions asked are: How big is the material stock currently embodied in the global passenger vehicle fleet, and when will this material become available for recycling? +# +# To answer these questions a dynamic material flow analysis of the global passenger vehicle fleet and the waste management industries is performed. +# +# The dynamic MFA model has the following indices: +# +# * t: time (1990-2017) +# * c: age-cohort (1990-2017) +# * r: region (130 countries accounting for most of the global vehicle fleet) +# * g: good (passenger vehicle) +# * p: process (vehicle markets, use phase, waste management industries, scrap markets) +# * m: engineering materials (25) +# * e: chemical elements (all) +# * w: waste types (steel, Al, Cu scrap, plastics, glass, and other waste) +# +# The system definition of the model is given in the figure below. The data availability limits the temporal scope to 2017. The figure also shows the aspects of the different system variables. The total registration of vehicles, for example, is broken down into individual materials, whereas the flow of deregistered vehicles is broken down into regions, age-cohorts, and materials. +# +# MFA system +# +# The model equations are as follows: +# +# 1) inflow-driven dynamic stock model, where $F_{1-2}$ is the historic inflow, $Sf$ is the survival function of the age-cohort (1-sum(pdf of discard)), and $S_2$ is the stock: +# $S_2(t,c,r,g) = F_{1-2}(c,r,g)\cdot Sf(t,c,r,g)$ +# 2) Calculation of difference between inflow-driven stock (covering only the age-cohorts 2005-2017 due to data availability) and the 2015 reported stock and distribution of the difference to the years 1990-2005 (constant inflow assumed for these years) +# 3) Calculation of material composition of the fleet $S_2$ with +# $S_2(t,c,r,g,m) = \mu(c,r,g,m)\cdot S_2(t,c,r,g)$ +# 4) Estimation of available future end-of-life vehicle scrap $F_{3-4}$ with +# $F_{3-4}(r,g,w,m) = \sum_{t,c}EoL_eff(r,g,m,w)\cdot M(t,c,r,g,m)$ +# +# The remaining system variables are calculated by mass balance. + +# %% [markdown] +# ## 1. Load sodym and other useful packages + +# %% +from os.path import join + +import numpy as np +import pandas as pd +import plotly.express as px + +from sodym.data_reader import DataReader +from sodym import ( + DimensionDefinition, Dimension, DimensionSet, ParameterDefinition, Parameter, Process, FlowDefinition, StockDefinition, MFASystem, +) +from sodym.stocks import InflowDrivenDSM +from sodym.survival_functions import NormalSurvival +from sodym.flow_helper import make_empty_flows +from sodym.stock_helper import make_empty_stocks + +# %% [markdown] +# ## 2. Define the data requirements, flows, stocks and MFA system equations +# +# We define the dimensions that are relevant for our system and the model parameters, processes, stocks and flows. We further define a class with our system equations in the compute method. + +# %% +dimension_definitions = [ + DimensionDefinition(letter='t', name='time', dtype=int), + DimensionDefinition(letter='r', name='region', dtype=str), + DimensionDefinition(letter='m', name='material', dtype=str), + DimensionDefinition(letter='w', name='waste', dtype=str), +] + +parameter_definitions = [ + ParameterDefinition(name='vehicle lifetime', dim_letters=('r', )), + ParameterDefinition(name='vehicle material content', dim_letters=('m', )), + ParameterDefinition(name='vehicle new registration', dim_letters=('t', 'r')), + ParameterDefinition(name='vehicle stock', dim_letters=('r', )), + ParameterDefinition(name='eol recovery rate', dim_letters=('m', 'w')) +] + +# %% +process_names = ['sysenv', 'market', 'use', 'waste', 'scrap'] +processes = {name: Process(name=name, id=index) for index, name in enumerate(process_names)} + +# %% +flow_definitions = [ + FlowDefinition(from_process_name='sysenv', to_process_name='market', dim_letters=('t', 'r')), + FlowDefinition(from_process_name='market', to_process_name='use', dim_letters=('t', 'm', 'r')), + FlowDefinition(from_process_name='use', to_process_name='waste', dim_letters=('t', 'm', 'r')), + FlowDefinition(from_process_name='waste', to_process_name='scrap', dim_letters=('t', 'w', 'm')), + FlowDefinition(from_process_name='waste', to_process_name='sysenv', dim_letters=('t', 'm')), + FlowDefinition(from_process_name='scrap', to_process_name='sysenv', dim_letters=('t', 'w', 'm')), +] +stock_definitions = [ + StockDefinition(name='in use', process='use', dim_letters=('t', 'r'),) +] + +# %% +class VehicleMFA(MFASystem): + """We just need to define the compute method with our system equations, + as all the other things we need are inherited from the MFASystem class.""" + def compute(self): + stock_diff = self.compute_stock() + self.compute_flows() + return stock_diff + + def compute_stock(self): + self.flows['sysenv => market'][...] = self.parameters['vehicle new registration'] + self.stocks['in use'].inflow[...] = self.flows['sysenv => market'] + survival_model = NormalSurvival( + dims=self.stocks['in use'].inflow.dims, + lifetime_mean=self.parameters['vehicle lifetime'].values, + lifetime_std=self.parameters['vehicle lifetime'].values*0.3 + ) + if not isinstance(self.stocks['in use'], InflowDrivenDSM): + self.stocks['in use'] = self.stocks['in use'].to_stock_type(InflowDrivenDSM, survival_model=survival_model) + else: + self.stocks['in use'].survival_model = survival_model + self.stocks['in use'].compute() + stock_diff = self.get_new_array(dim_letters=('r')) + stock_diff[...] = 1000 * self.parameters['vehicle stock'] - self.stocks['in use'].stock[{'t': 2015}] + return (stock_diff * 1e-6).to_df() # in millions + + def compute_flows(self): + self.flows['market => use'][...] = self.parameters['vehicle material content'] * self.parameters['vehicle new registration'] * 1e-9 + self.flows['use => waste'][...] = self.parameters['vehicle material content'] * self.stocks['in use'].outflow * 1e-9 + self.flows['waste => scrap'][...] = self.parameters['eol recovery rate'] * self.flows['use => waste'] + self.flows['scrap => sysenv'][...] = self.flows['waste => scrap'] + self.flows['waste => sysenv'][...] = self.flows['use => waste'] - self.flows['waste => scrap'] + +# %% [markdown] +# ## 3. Define our data reader +# Now that we have defined the MFA system and know what data we need, we can load the data. +# To do the data loading, we define a DataReader class. Such a class can be reused with different datasets of the same format by passing attributes, e.g. the directory where the data is stored, in the init function. In this example, we will also build upon this data reader in a following step. + +# %% +class CustomDataReader(DataReader): + """The methods `read_dimensions` and `read_parameters` are already defined in the parent + DataReader class, and loop over the methods `read_dimension` and `read_parameter_values` + that we specify for our usecase here. + """ + def __init__(self, data_directory): + self.data_directory = data_directory + + @property + def years(self): + return list(range(2012, 2018)) + + def read_dimension(self, dimension_definition: DimensionDefinition) -> Dimension: + if (dim_name := dimension_definition.name) == 'region': + data = pd.read_excel(join(self.data_directory, 'vehicle_lifetime.xlsx'), 'Data') + other_data = pd.read_excel(join(self.data_directory, 'vehicle_stock.xlsx'), 'Data') + data = (set(data[dim_name].unique())).intersection(set(other_data[dim_name].unique())) + data = list(data) + data.sort() + elif (dim_name := dimension_definition.name) in ['waste', 'material']: + data = pd.read_excel(join(self.data_directory, 'eol_recovery_rate.xlsx'), 'Data') + data.columns = [x.lower() for x in data.columns] + data = list(data[dim_name].unique()) + data.sort() + elif dimension_definition.name == 'time': + data = self.years + return Dimension(name=dimension_definition.name, letter=dimension_definition.letter, items=data) + + def read_parameter_values(self, parameter: str, dims: DimensionSet) -> Parameter: + data = pd.read_excel(join(self.data_directory, (parameter.replace(' ', '_') + '.xlsx')), 'Data') + data = data.fillna(0) + if 'r' in dims.letters: # remove unwanted regions + data = data[data['region'].isin(dims['r'].items)] + + if parameter == 'vehicle new registration': + return self.vehicle_new_registration(data, dims) + + data.columns = [x.lower() for x in data.columns] + data = data[[dim.name for dim in dims] + ['value']] # select only relevant information + data.sort_values([dim.name for dim in dims], inplace=True) # sort to ensure correct order + if len(dims.dim_list) == 1: + return Parameter(dims=dims, values=data['value'].values) + elif len(dims.dim_list) == 2: + multiindex = data.set_index([dim.name for dim in dims]) + data = multiindex.unstack().values[:, :] + return Parameter(dims=dims, values=data) + + def vehicle_new_registration(self, data, dims): + data.sort_values('region', inplace=True) + data = data[self.years] + return Parameter(dims=dims, values=data.values.T) + +# %% [markdown] +# ## 4. Put everything together +# We make an instance of our `CustomDataReader`, read in the data and use it to create an instance of our `VehicleMFA` class. Then we can run the calculations, and check what our estimate of vehicle stocks looks like compared to the data for 2015 in the `vehicle_stock.xlsx` file. + +# %% +data_reader = CustomDataReader(data_directory='example5_data') +dimensions = data_reader.read_dimensions(dimension_definitions) +parameters = data_reader.read_parameters(parameter_definitions, dims=dimensions) + +stocks = make_empty_stocks(stock_definitions=stock_definitions, processes=processes, dims=dimensions) +flows = make_empty_flows(processes=processes, dims=dimensions, flow_definitions=flow_definitions) + +vehicle_mfa = VehicleMFA(dims=dimensions, parameters=parameters, flows=flows, stocks=stocks, processes=processes) +stock_diff = vehicle_mfa.compute_stock() # compute 1st stock estimate and difference to data + +# %% +stock_diff + +# %% [markdown] +# ## 5. Improve the model +# +# We can see that the results from our inflow-driven dynamic stock model differ significantly from the data about vehicle stocks. One reason for this is the missing inflows for years before our model starts. In a next step, we extend the years back to 1990 and copy the inflow data (given by the `new_vehicle_registration` parameter) by region for each of these earlier years; assuming that the inflow stayed constant between 1990 and 2004. +# +# In addition to this and to answer our research questions, we will extend the timeseries out to 2050, with no more inflow afer 2017. +# +# To make these changes in the model dimensions and parameters, we build on the `CustomDataReader` class defined above, extending the time dimension to earlier years, as well as extending the data for the `vehicle_new_registration` parameter, which is our only time-dependent parameter. + +# %% +class AnotherCustomDataReader(CustomDataReader): + @property + def years(self): + return self.historical_years + self.future_years + + @property + def historical_years(self): + return list(range(1990, 2018)) + + @property + def future_years(self): + return list(range(2018, 2051)) + + def vehicle_new_registration(self, data, dims): + data.sort_values('region', inplace=True) + data.set_index('region', inplace=True) + data_years = [k for k in data.keys() if isinstance(k, int)] + data = data[data_years] + + repeated_values = np.tile(data[min(data_years)].values, (len(self.historical_years), 1)) + historical_data = pd.DataFrame(index=data.index, columns=self.historical_years, data=repeated_values.T) + future_data = pd.DataFrame(index=data.index, columns=self.future_years, data=0) + data = data.combine_first(historical_data).combine_first(future_data) + return Parameter(dims=dims, values=data.values.T) + +# %% +data_reader_2 = AnotherCustomDataReader(data_directory='example5_data') +dimensions_2 = data_reader_2.read_dimensions(dimension_definitions) +parameters_2 = data_reader_2.read_parameters(parameter_definitions, dims=dimensions_2) + +stocks_2 = make_empty_stocks(stock_definitions=stock_definitions, processes=processes, dims=dimensions_2) +flows_2 = make_empty_flows(processes=processes, dims=dimensions_2, flow_definitions=flow_definitions) + +vehicle_mfa_2 = VehicleMFA(dims=dimensions_2, parameters=parameters_2, flows=flows_2, stocks=stocks_2, processes=processes) +stock_diff_2 = vehicle_mfa_2.compute() # compute 2nd stock estimate and difference to data, as well as the MFA flows + +# %% +stock_diff_2[stock_diff_2 <= stock_diff] + +# %% [markdown] +# The stock values in 2015 have improved as compared to the data, but are still quite far away from the dataset. + +# %% [markdown] +# ## 6. Answer research questions +# Approximately how big was the material stock embodied in the global passenger vehicle fleet in 2017? +# And when will this material become available for recycling? + +# %% +stock_by_material_type = vehicle_mfa_2.stocks['in use'].stock * vehicle_mfa_2.parameters['vehicle material content'] * 1e-9 +global_stock_by_material_type = stock_by_material_type.sum_nda_over(sum_over_dims=('r')) +global_stock_by_material_type_in_2017 = global_stock_by_material_type[{'t': 2017}] + +stock_df = global_stock_by_material_type_in_2017.to_df(index=False) +fig = px.bar(stock_df, x='material', y='value') +fig.show() + +# %% +np.nan_to_num(vehicle_mfa_2.flows["scrap => sysenv"].values, copy=False) +scrap_outflow = vehicle_mfa_2.flows["scrap => sysenv"].sum_nda_over(sum_over_dims=('r', 'm')) +outflow_df = pd.DataFrame( + data=scrap_outflow.values, + columns=scrap_outflow.dims['w'].items, + index=scrap_outflow.dims['t'].items, +) +outflow_df = outflow_df[outflow_df.index > 2017] +fig = px.line(outflow_df, title='Scrap outflow') +fig.show() +fig.update_yaxes(type="log") +fig.show() + +# %% + + + diff --git a/examples/jupytext.toml b/examples/jupytext.toml new file mode 100644 index 0000000..cdc7434 --- /dev/null +++ b/examples/jupytext.toml @@ -0,0 +1 @@ +formats = "ipynb,py:percent" \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 9867202..257d2f5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,7 @@ dependencies = [ "docformatter>=1.7.5", "openpyxl>=3.1.5", "nbformat>=5.10.4", + "jupytext>=1.16.4", ] [project.optional-dependencies]