Building a Simple Energy Model¶
The goal of this tutorial is to build a simple energy model using MESSAGEix with minimal features that can be expanded in future tutorials.
We will build the model component by component, focusing on both the how (code implementation) and why (mathematical formulation).
Follow Along at Home¶
The full model documentation is available online at https://messageix.iiasa.ac.at. Source code is at https://github.com/iiasa/message_ix.
And you can easily install MESSAGEix yourself and get all the tutorials:
$ conda install -c conda-forge message-ix
$ messageix-dl # install all tutorials to your current directory
MESSAGEix¶
- system of constraints
- optimizing an objective function, nominally minimizing total system cost
- decision VARIABLES ($x$) are capitalized
- parameters ($A$, $b$) are lower case
MESSAGEix: the Framework¶
A Brave New World: Westeros with Electricity!¶
 |  |
Setup¶
First, we import all the packages we need.
import pandas as pd
import ixmp as ix
import message_ix
from message_ix.utils import make_df
%matplotlib inline
The Platform is your connection to a database which will automatically store all data for you
mp = ix.Platform(dbtype='HSQLDB')
Once connected, we make our Scenario which we use to build our model.
scenario = message_ix.Scenario(mp, model='Westeros Electrified',
scen='baseline', version='new')
Model Structure¶
We start by defining basic characteristics of the model, including time, space, and the energy system structure.
The model horizon will span 3 decades. Let's assume that we're far in the future after the events of A Song of Ice and Fire (which occur ~300 years after Aegon the conqueror).
| Math Notation | Model Meaning |
|---|---|
| $y \in Y^H$ | time periods in history |
| $y \in Y^M$ | time periods in model horizon |
history = [690]
model_horizon = [700, 710, 720]
scenario.add_horizon({'year': history + model_horizon,
'firstmodelyear': model_horizon[0]})
Our model will have a single node, i.e., its spatial dimension.
| Math Notation | Model Meaning |
|---|---|
| $n \in N$ | node |
country = 'Westeros'
scenario.add_spatial_sets({'country': country})
And we fill in the energy system's commodities, levels, technologies, and mode (defining how certain technologies operate).
| Math Notation | Model Meaning |
|---|---|
| $c \in C$ | commodity |
| $l \in L$ | level |
| $t \in T$ | technology |
| $m \in M$ | mode |
scenario.add_set("commodity", ["electricity", "light"])
scenario.add_set("level", ["secondary", "final", "useful"])
scenario.add_set("technology", ['coal_ppl', 'wind_ppl', 'grid', 'bulb'])
scenario.add_set("mode", "standard")
Supply and Demand (or Balancing Commodities)¶
The fundamental premise of the model is to satisfy demand. To first order, demand for services like electricity track with economic productivity (GDP). We define a GDP profile similar to first-world GDP growth from 1900-1930:
gdp_profile = pd.Series([1., 1.5, 1.9], index=model_horizon)
gdp_profile.plot(title='Demand')
The COMMODITY_BALANCE equation ensures that demand for each commodity is met at each level in the energy system.
$\sum_{\substack{n^L,t,m \\ y^V \leq y}} output_{n^L,t,y^V,y,m,n,c,l} \cdot ACT_{n^L,t,y^V,y,m}$
$- \sum_{\substack{n^L,t,m, \\ y^V \leq y}} input_{n^L,t,y^V,y,m,n,c,l} \cdot ACT_{n^L,t,m,y}$
$\geq demand_{n,c,l,y} \quad \forall \ l \in L$
The COMMODITY_BALANCE equation is formulated as an inequality implying that demand must be met, but supply of a commodity can exceed demand. The formulation implicitly assumes "zero cost of disposal", as is common in economics. This implementation simplifies calibration and is in line with conventions in energy systems modelling.
First we establish demand. Let's assume
- 40 million people in 300 AC
- similar population growth to Earth in the same time frame (~factor of 12)
- a per capita demand for electricity of 1000 kWh
- and 8760 hours in a year (of course!)
Then we can add the demand parameter
demand_baseyear = 40e6 * 12 * 1000 * 1e6 / 8760
light_demand = pd.DataFrame({
'node': country,
'commodity': 'light',
'level': 'useful',
'year': model_horizon,
'time': 'year',
'value': demand_baseyear * gdp_profile,
'unit': 'GWa',
})
light_demand
scenario.add_par("demand", light_demand)
Working backwards along the Reference Energy System, we can add connections for the bulb
bulb_out = make_df(base_output, technology='bulb', commodity='light',
level='useful', value=1.0)
scenario.add_par('output', bulb_out)
bulb_in = make_df(base_input, technology='bulb', commodity='electricity',
level='final', value=1.0, unit='%')
scenario.add_par('input', bulb_in)
Next, the grid, with loses of 13%
grid_efficiency = 0.87
grid_out = make_df(base_output, technology='grid', commodity='electricity',
level='final', value=grid_efficiency)
scenario.add_par('output', grid_out)
grid_in = make_df(base_input, technology='grid', commodity='electricity',
level='secondary', value=1.0, unit='%')
scenario.add_par('input', grid_in)
And finally, our power plants. The model does not include the fossil resources used as input for coal plants; however, costs of coal extraction are included in the parameter $variable\_cost$.
coal_out = make_df(base_output, technology='coal_ppl', commodity='electricity',
level='secondary', value=1.)
scenario.add_par('output', coal_out)
wind_out = make_df(base_output, technology='wind_ppl', commodity='electricity',
level='secondary', value=1.)
scenario.add_par('output', wind_out)
Operational Constraints and Parameters¶
The model has a number of "reality" constraints, which relate built capacity to available power.
The CAPACITY_CONSTRAINT
$\sum_{m} ACT_{n,t,y^V,y,m,h} \leq duration\_time_{h} \cdot capacity\_factor_{n,t,y^V,y,h} \cdot CAP_{n,t,y^V,y} \quad t \ \in \ T^{INV}$
This requires us to provide capacity factors
capacity_factor = {
'coal_ppl': 0.85,
'wind_ppl': 0.2,
'bulb': 0.1,
}
for tec, val in capacity_factor.items():
df = make_df(base_capacity_factor, technology=tec, value=val)
scenario.add_par('capacity_factor', df)
The model can further be provided technical_lifetimes in order to properly manage deployed capacity and related costs via the CAPACITY_MAINTENENCE constraint:
$CAP_{n,t,y^V,y} \leq remaining\_capacity_{n,t,y^V,y} \cdot value \: \: \forall t \in T^{INV}$
Where value can take different forms depending on what time period is considered:
| Value | Condition |
|---|---|
| $\Delta_y historical\_new\_capacity_{n,t,y^V}$ | $y$ is first model period |
| $\Delta_y CAP\_NEW_{n,t,y^V}$ | $y = y^V$ |
| $CAP_{n,t,y^V,y-1}$ | if $y > y^V$ and $y - y^V < technical\_lifetime_{n,t,y^V}$ |
lifetime = {
'coal_ppl': 20,
'wind_ppl': 10,
'bulb': 1,
}
for tec, val in lifetime.items():
df = make_df(base_technical_lifetime, technology=tec, value=val)
scenario.add_par('technical_lifetime', df)
Technological Diffusion and Contraction¶
We know from historical precedent that energy systems can not be transformed instantaneously. Therefore, we use a family of constraints on activity (ACT) and capacity (CAP).
$\sum_{y^V \leq y,m} ACT_{n,t,y^V,y,m,h} \leq$
$initial\_activity\_up_{n,t,y,h}
\cdot \frac{ \Big( 1 + growth\_activity\_up_{n,t,y,h} \Big)^{|y|} - 1 }
{ growth\_activity\_up_{n,t,y,h} }+ \Big( 1 + growth\_activity\_up_{n,t,y,h} \Big)^{|y|} \cdot \Big( \sum_{y^V \leq y-1,m} ACT_{n,t,y^V,y-1,m,h} + \sum_{m} historical\_activity_{n,t,y-1,m,h}\Big)$
This example limits the ability for technologies to grow. To do so, we need to provide growth_activity_up values for each technology that we want to model as being diffusion constrained. Here, we set this constraint at 5% per year.
growth_technologies = [
"coal_ppl",
"wind_ppl",
]
for tec in growth_technologies:
df = make_df(base_growth, technology=tec, value=0.05)
scenario.add_par('growth_activity_up', df)
Defining an Energy Mix¶
To model the transition of an energy system, one must start with the existing system which are defined by the parameters historical_activity and historical_capacity. These parameters define the energy mix before the model horizon.
We begin by defining a few key values:
- how much useful energy was needed
- how much final energy was generated
- and the mix for different technologies
historic_demand = 0.85 * demand_baseyear
historic_generation = historic_demand / grid_efficiency
coal_fraction = 0.6
Then, we can define the activity and capacity in the historic period
old_activity = {
'coal_ppl': coal_fraction * historic_generation,
'wind_ppl': (1 - coal_fraction) * historic_generation,
}
for tec, val in old_activity.items():
df = make_df(base_activity, technology=tec, value=val)
scenario.add_par('historical_activity', df)
act_to_cap = {
'coal_ppl': 1 / 10 / capacity_factor['coal_ppl'] / 2, # 20 year lifetime
'wind_ppl': 1 / 10 / capacity_factor['wind_ppl'],
}
for tec in act_to_cap:
value = old_activity[tec] * act_to_cap[tec]
df = make_df(base_capacity, technology=tec, value=value)
scenario.add_par('historical_new_capacity', df)
Objective Function¶
The objective function drives the purpose of the optimization. Do we wish to seek maximum utility of the social planner, minimize carbon emissions, or something else? Energy-system focused IAMs seek to minimize total discounted system cost over space and time.
$\min \sum_{n,y \in Y^{M}} interestrate_{y} \cdot COST\_NODAL_{n,y}$
First, let's add the interest rate parameter.
rate = [0.05] * len(model_horizon)
unit = ['%'] * len(model_horizon)
scenario.add_par("interestrate", key=model_horizon, val=rate, unit=unit)
COST_NODAL is comprised of a variety of costs related to the use of different technologies.
Investment Costs¶
Capital, or investment, costs are invoked whenever a new plant or unit is built
$inv\_cost_{n,t,y} \cdot construction\_time\_factor_{n,t,y} \cdot CAP\_NEW_{n,t,y}$
# in $ / kW
costs = {
'coal_ppl': 1500,
'wind_ppl': 1100,
'bulb': 5,
}
for tec, val in costs.items():
df = make_df(base_inv_cost, technology=tec, value=val * 1e6)
scenario.add_par('inv_cost', df)
Fixed O&M Costs¶
Fixed cost are only relevant as long as the capacity is active. This formulation allows to include the potential cost savings from early retirement of installed capacity.
$\sum_{y^V \leq y} \ fix\_cost_{n,t,y^V,y} \cdot CAP_{n,t,y^V,y}$
# in $ / kW
costs = {
'coal_ppl': 40,
'wind_ppl': 30,
}
for tec, val in costs.items():
df = make_df(base_fix_cost, technology=tec, value=val * 1e6)
scenario.add_par('fix_cost', df)
Variable O&M Costs¶
Variable Operation and Maintence costs are associated with the costs of actively running the plant. Thus, they are not applied if a plant is on standby (i.e., constructed, but not currently in use).
$\sum_{\substack{y^V \leq y \\ m,h}} \ var\_cost_{n,t,y^V,y,m,h} \cdot ACT_{n,t,y^V,y,m,h} $
# in $ / MWh
costs = {
'coal_ppl': 24.4,
}
for tec, val in costs.items():
df = make_df(base_var_cost, technology=tec, value=val * 8760. * 1e3)
scenario.add_par('var_cost', df)
A full model will also have costs associated with
- resource extraction: $\sum_{c,g} \ resource\_cost_{n,c,g,y} \cdot EXT_{n,c,g,y} $
- emissions : $\sum_{e,t, \hat{e} \in E(e)} emission\_scaling_{e,\hat{e}} \cdot \ emission\_tax_{n,e,t,y} \cdot EMISS_{n,e,t,y}$
- land use (emulator): $\sum_{s} land\_cost_{n,s,y} \cdot LAND_{n,s,y}$
- exceedence on "hard" bounds of technology expansion and contraction: $\sum_{m,h} \ \Big( abs\_cost\_activity\_soft\_lo_{n,t,y,m,h} + level\_cost\_activity\_soft\_lo_{n,t,y,m,h} \cdot\ levelized\_cost_{n,t,y,m,h}\Big) \cdot ACT\_LO_{n,t,y,h}$
Time to Solve the Model¶
scenario.commit(comment='basic model of Westerosi electrification')
scenario.set_as_default()
scenario.solve()
scenario.var('OBJ')['lvl']
Plotting Results¶
from tools import Plots
p = Plots(scenario, country, firstyear=model_horizon[0])
Activity¶
How much energy is generated in each time period from the different energy supply technologies?
p.plot_activity(baseyear=True, subset=['coal_ppl', 'wind_ppl'])
Capacity¶
How many new plants are built?
p.plot_new_capacity(baseyear=True, subset=['coal_ppl', 'wind_ppl'])
Electricity Price¶
And how much does the electricity cost? These prices are given by the dual variables of the commodity balance constraint. Economists use the term shadow prices instead of dual variables. They reflect the marginal price of electricity, taken from the most expensive producer.
Note that prices are lower when the system does not depend on expensive technologies any longer because it sufficiently expanded capacity of cheap technologies.
p.plot_prices(subset=['light'], baseyear=True)
Fin¶
With that, you have built and run your very first MESSAGEix model. Welcome to the community!
Check us out on Github: https://github.com/iiasa/message_ix
Get in touch with us online: https://groups.google.com/forum/message-ix
And feel free to contact me with any further questions: gidden@iiasa.ac.at