Stationary example of pathway optimisation¶

This example solves a bounded stationary pathway optimisation problem. It is based on the Italy model developed in https://github.com/FLomb/Calliope-Italy, with lower spatial resolution.

In [1]:

import calliope
import calliope_pathways
import plotly.express as px

OUTPUT_PATH = "outputs/"

calliope.set_log_verbosity("INFO", include_solver_output=False)

import calliope import calliope_pathways import plotly.express as px OUTPUT_PATH = "outputs/" calliope.set_log_verbosity("INFO", include_solver_output=False)

Model input¶

In [2]:

# Initialise
model = calliope_pathways.models.italy()

# Initialise model = calliope_pathways.models.italy()

[2024-07-19 11:27:28] INFO     Model: initialising

[2024-07-19 11:27:28] INFO     Model: preprocessing stage 1 (model_run)

[2024-07-19 11:27:34] INFO     Model: preprocessing stage 2 (model_data)

[2024-07-19 11:27:34] INFO     Model: preprocessing complete

Assessing the input data¶

The model can be configured by running the script in src/calliope_pathways/models/italy/pre_processing/parse_lombardi.py. It automatically assigns random decommission rates for each technology, accounting for their lifetimes. The number of vintagesteps and investsteps can be altered when calling the model (e.g., :::python calliope_pathways.models.italy(first_year=2030, investstep_resolution=10)).

In [3]:

model.inputs.investsteps

model.inputs.investsteps

Out[3]:

<xarray.DataArray 'investsteps' (investsteps: 6)> Size: 48B
array(['2025-01-01T00:00:00.000000000', '2030-01-01T00:00:00.000000000',
       '2035-01-01T00:00:00.000000000', '2040-01-01T00:00:00.000000000',
       '2045-01-01T00:00:00.000000000', '2050-01-01T00:00:00.000000000'],
      dtype='datetime64[ns]')
Coordinates:
  * investsteps  (investsteps) datetime64[ns] 48B 2025-01-01 ... 2050-01-01

xarray.DataArray

'investsteps'

• investsteps: 6

• 2025-01-01 2030-01-01 2035-01-01 2040-01-01 2045-01-01 2050-01-01

  array(['2025-01-01T00:00:00.000000000', '2030-01-01T00:00:00.000000000',
         '2035-01-01T00:00:00.000000000', '2040-01-01T00:00:00.000000000',
         '2045-01-01T00:00:00.000000000', '2050-01-01T00:00:00.000000000'],
        dtype='datetime64[ns]')

• Coordinates: (1)
  • investsteps
    (investsteps)
    datetime64[ns]
    2025-01-01 ... 2050-01-01

    array(['2025-01-01T00:00:00.000000000', '2030-01-01T00:00:00.000000000',
           '2035-01-01T00:00:00.000000000', '2040-01-01T00:00:00.000000000',
           '2045-01-01T00:00:00.000000000', '2050-01-01T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Indexes: (1)
  • investsteps
    PandasIndex

    PandasIndex(DatetimeIndex(['2025-01-01', '2030-01-01', '2035-01-01', '2040-01-01',
                   '2045-01-01', '2050-01-01'],
                  dtype='datetime64[ns]', name='investsteps', freq=None))

• Attributes: (0)

In [4]:

model.inputs.vintagesteps

model.inputs.vintagesteps

Out[4]:

<xarray.DataArray 'vintagesteps' (vintagesteps: 6)> Size: 48B
array(['2025-01-01T00:00:00.000000000', '2030-01-01T00:00:00.000000000',
       '2035-01-01T00:00:00.000000000', '2040-01-01T00:00:00.000000000',
       '2045-01-01T00:00:00.000000000', '2050-01-01T00:00:00.000000000'],
      dtype='datetime64[ns]')
Coordinates:
  * vintagesteps  (vintagesteps) datetime64[ns] 48B 2025-01-01 ... 2050-01-01

xarray.DataArray

'vintagesteps'

• vintagesteps: 6

• 2025-01-01 2030-01-01 2035-01-01 2040-01-01 2045-01-01 2050-01-01

  array(['2025-01-01T00:00:00.000000000', '2030-01-01T00:00:00.000000000',
         '2035-01-01T00:00:00.000000000', '2040-01-01T00:00:00.000000000',
         '2045-01-01T00:00:00.000000000', '2050-01-01T00:00:00.000000000'],
        dtype='datetime64[ns]')

• Coordinates: (1)
  • vintagesteps
    (vintagesteps)
    datetime64[ns]
    2025-01-01 ... 2050-01-01

    array(['2025-01-01T00:00:00.000000000', '2030-01-01T00:00:00.000000000',
           '2035-01-01T00:00:00.000000000', '2040-01-01T00:00:00.000000000',
           '2045-01-01T00:00:00.000000000', '2050-01-01T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Indexes: (1)
  • vintagesteps
    PandasIndex

    PandasIndex(DatetimeIndex(['2025-01-01', '2030-01-01', '2035-01-01', '2040-01-01',
                   '2045-01-01', '2050-01-01'],
                  dtype='datetime64[ns]', name='vintagesteps', freq=None))

• Attributes: (0)

Stationary means that both costs and demand remain constant across the years:

In [5]:

model.inputs.cost_flow_cap.to_series().dropna()

model.inputs.cost_flow_cap.to_series().dropna()

Out[5]:

techs               costs   
ac_CNOR_to_CSUD     monetary     450.0
ac_CNOR_to_SARD     monetary    1660.0
ac_CSUD_SUD         monetary     450.0
ac_CSUD_to_SARD     monetary     750.0
ac_NORD_to_CNOR     monetary     450.0
ac_SUD_to_SICI      monetary     636.0
bioenergy           monetary    4137.0
ccgt                monetary     709.0
coal                monetary    1584.0
export_electricity  monetary     450.0
geothermal          monetary    4062.0
hydropower          monetary    4022.0
import_electricity  monetary     450.0
oil                 monetary     490.0
pv                  monetary    1322.0
waste               monetary    3743.0
wind                monetary    1451.0
Name: cost_flow_cap, dtype: float64

Initial capacity¶

The model is initialized with the capacity of Italy in 2015.

In [6]:

model.inputs.flow_cap_initial.to_series().dropna()

model.inputs.flow_cap_initial.to_series().dropna()

Out[6]:

techs            nodes
ac_CNOR_to_CSUD  CNOR     1300000.0
                 CSUD     1300000.0
ac_CNOR_to_SARD  CNOR      300000.0
                 SARD      300000.0
ac_CSUD_SUD      CSUD     4600000.0
                            ...    
wind             CSUD     1632100.0
                 NORD      115600.0
                 SARD     1005500.0
                 SICI     1757600.0
                 SUD      4517600.0
Name: flow_cap_initial, Length: 74, dtype: float64

This technology capacity can then be phased out as we step through to the end of our time horizon:

In [7]:

# This is a fraction of the initial capacity that remains available in investment steps
model.inputs.available_initial_cap.to_series().dropna()

# This is a fraction of the initial capacity that remains available in investment steps model.inputs.available_initial_cap.to_series().dropna()

Out[7]:

techs            nodes  investsteps
ac_CNOR_to_CSUD  CNOR   2025-01-01     1.000000
                        2030-01-01     1.000000
                        2035-01-01     1.000000
                        2040-01-01     1.000000
                        2045-01-01     1.000000
                                         ...   
wind             SUD    2030-01-01     0.995501
                        2035-01-01     0.767926
                        2040-01-01     0.057528
                        2045-01-01     0.000000
                        2050-01-01     0.000000
Name: available_initial_cap, Length: 396, dtype: float64

End-of-life decommissioning is tracked with a similar matrix.

Note how vintages are never available in investsteps that are in their future.

In [8]:

model.inputs.available_vintages.to_series().dropna().unstack("investsteps")

model.inputs.available_vintages.to_series().dropna().unstack("investsteps")

Out[8]:

	investsteps	2025-01-01	2030-01-01	2035-01-01	2040-01-01	2045-01-01	2050-01-01
techs	vintagesteps
ac_CNOR_to_CSUD	2025-01-01	1.0	1.0	1.0	1.0	1.0	1.0
	2030-01-01	NaN	1.0	1.0	1.0	1.0	1.0
	2035-01-01	NaN	NaN	1.0	1.0	1.0	1.0
	2040-01-01	NaN	NaN	NaN	1.0	1.0	1.0
	2045-01-01	NaN	NaN	NaN	NaN	1.0	1.0
...	...	...	...	...	...	...	...
wind	2030-01-01	NaN	1.0	1.0	1.0	1.0	0.0
	2035-01-01	NaN	NaN	1.0	1.0	1.0	1.0
	2040-01-01	NaN	NaN	NaN	1.0	1.0	1.0
	2045-01-01	NaN	NaN	NaN	NaN	1.0	1.0
	2050-01-01	NaN	NaN	NaN	NaN	NaN	1.0

114 rows × 6 columns

Building a pathways optimisation problem¶

We have created a math YAML file with updates to all the pre-defined math to handle the existence of investsteps and vintagesteps. Tracking new capacity in each investment period and linking it to technology vintages requires new variables and constraints

Variables¶

In [9]:

# Note the "investsteps" dimension added to this pre-defined variable
model.math["variables"]["flow_cap"]

# Note the "investsteps" dimension added to this pre-defined variable model.math["variables"]["flow_cap"]

Out[9]:

{'description': "A technology's flow capacity, also known as its nominal or nameplate capacity.",
 'default': 0,
 'unit': 'power',
 'foreach': ['nodes', 'techs', 'carriers', 'investsteps'],
 'bounds': {'min': 'flow_cap_min', 'max': 'flow_cap_max'}}

In [10]:

# This new variable tracks the amount of each technology vintage that exists
model.math["variables"]["flow_cap_new"]

# This new variable tracks the amount of each technology vintage that exists model.math["variables"]["flow_cap_new"]

Out[10]:

{'bounds': {'max': 'flow_cap_new_max', 'min': 0},
 'description': 'Additional flow capacity commissioned in an investstep.',
 'foreach': ['nodes', 'techs', 'carriers', 'vintagesteps'],
 'unit': 'power',
 'where': 'flow_cap'}

Constraints¶

In [11]:

# All existing constraints have the "investsteps" dimension added.
# This allows the dispatch decisions to be optimised individually for each investment step.
model.math["constraints"]["system_balance"]

# All existing constraints have the "investsteps" dimension added. # This allows the dispatch decisions to be optimised individually for each investment step. model.math["constraints"]["system_balance"]

Out[11]:

{'description': 'Set the global carrier balance of the optimisation problem by fixing the total production of a given carrier to equal the total consumption of that carrier at every node in every timestep.',
 'foreach': ['nodes', 'carriers', 'timesteps', 'investsteps'],
 'equations': [{'expression': 'sum(flow_out, over=techs) - sum(flow_in, over=techs) - $flow_export + $unmet_demand_and_unused_supply == 0'}],
 'sub_expressions': {'flow_export': [{'where': 'any(carrier_export, over=techs)',
    'expression': 'sum(flow_export, over=techs)'},
   {'where': 'NOT any(carrier_export, over=techs)', 'expression': '0'}],
  'unmet_demand_and_unused_supply': [{'where': 'config.ensure_feasibility=True',
    'expression': 'unmet_demand + unused_supply'},
   {'where': 'NOT config.ensure_feasibility=True', 'expression': '0'}]}}

In [12]:

# In each investment period, capacities are a combination of all available vintages.
model.math["constraints"]["flow_cap_bounding"]

# In each investment period, capacities are a combination of all available vintages. model.math["constraints"]["flow_cap_bounding"]

Out[12]:

{'description': 'Flow capacity in a given investment period is the capacity in the previous period plus the new capacity installed minus old capacity that has been decommissioned between the previous period and this one.',
 'equations': [{'expression': 'flow_cap == sum(flow_cap_new * available_vintages, over=vintagesteps) + flow_cap_initial * available_initial_cap'}],
 'foreach': ['nodes', 'techs', 'carriers', 'investsteps'],
 'where': 'flow_cap'}

In [13]:

# Similarly, note that some technologies cannot go beyond their initial installed capacity.
model.inputs.flow_cap_max.to_series().dropna()

# Similarly, note that some technologies cannot go beyond their initial installed capacity. model.inputs.flow_cap_max.to_series().dropna()

Out[13]:

techs        nodes
battery_phs  CSUD      1681600.0
             NORD      5064300.0
             SARD       236900.0
             SICI       573700.0
geothermal   CNOR       768000.0
hydropower   CNOR      1100900.0
             CSUD      1053400.0
             NORD     11191600.0
             SARD       223800.0
             SICI       145800.0
             SUD        946700.0
waste        CNOR        23000.0
             CSUD       127700.0
             NORD       384700.0
             SARD         3900.0
             SUD         32700.0
Name: flow_cap_max, dtype: float64

Building¶

In [14]:

model.build()

model.build()

[2024-07-19 11:27:34] INFO     Optimisation Model | parameters | Generated.

[2024-07-19 11:27:35] INFO     Optimisation Model | variables | Generated.

[2024-07-19 11:27:36] INFO     Optimisation Model | global_expressions | Generated.

[2024-07-19 11:27:39] INFO     Optimisation Model | constraints | Generated.

[2024-07-19 11:27:39] INFO     Optimisation Model | objectives | Generated.

Analyse results¶

In [15]:

model.solve()

model.solve()

[2024-07-19 11:27:39] INFO     Optimisation model | starting model in plan mode.

[2024-07-19 11:27:41] INFO     Backend: solver finished running. Time since start of solving optimisation problem: 0:00:02.583646

[2024-07-19 11:27:41] INFO     Postprocessing: started

[2024-07-19 11:27:41] INFO     Postprocessing: All values < 1e-10 set to 0 in flow_cap, flow_out, flow_in, source_use, flow_cap_new, flow_out_inc_eff, flow_in_inc_eff, cost_var, cost_investment, cost, capacity_factor, systemwide_capacity_factor

[2024-07-19 11:27:41] INFO     Postprocessing: ended. Time since start of solving optimisation problem: 0:00:02.709833

[2024-07-19 11:27:41] INFO     Model: loaded model_data

In [16]:

df_capacity = (
    model.results.flow_cap.where(model.results.techs != "demand_electricity")
    .sel(carriers="electricity")
    .sum("nodes")
    .to_series()
    .where(lambda x: x != 0)
    .dropna()
    .to_frame("Flow capacity (kW)")
    .reset_index()
)

print(df_capacity.head())

fig = px.bar(
    df_capacity,
    x="investsteps",
    y="Flow capacity (kW)",
    color="techs",
    color_discrete_map=model.inputs.color.to_series().to_dict(),
)
fig.show()

df_capacity = ( model.results.flow_cap.where(model.results.techs != "demand_electricity") .sel(carriers="electricity") .sum("nodes") .to_series() .where(lambda x: x != 0) .dropna() .to_frame("Flow capacity (kW)") .reset_index() ) print(df_capacity.head()) fig = px.bar( df_capacity, x="investsteps", y="Flow capacity (kW)", color="techs", color_discrete_map=model.inputs.color.to_series().to_dict(), ) fig.show()

             techs investsteps  Flow capacity (kW)
0  ac_CNOR_to_CSUD  2025-01-01           2600000.0
1  ac_CNOR_to_CSUD  2030-01-01           2600000.0
2  ac_CNOR_to_CSUD  2035-01-01           2600000.0
3  ac_CNOR_to_CSUD  2040-01-01           2600000.0
4  ac_CNOR_to_CSUD  2045-01-01           2600000.0

In [17]:

df_capacity = (
    model.results.storage_cap.sum("nodes")
    .to_series()
    .where(lambda x: x != 0)
    .dropna()
    .to_frame("Storage capacity (kWh)")
    .reset_index()
)

print(df_capacity.head())

fig = px.bar(
    df_capacity,
    x="investsteps",
    y="Storage capacity (kWh)",
    color="techs",
    color_discrete_map=model.inputs.color.to_series().to_dict(),
)
fig.show()

df_capacity = ( model.results.storage_cap.sum("nodes") .to_series() .where(lambda x: x != 0) .dropna() .to_frame("Storage capacity (kWh)") .reset_index() ) print(df_capacity.head()) fig = px.bar( df_capacity, x="investsteps", y="Storage capacity (kWh)", color="techs", color_discrete_map=model.inputs.color.to_series().to_dict(), ) fig.show()

         techs investsteps  Storage capacity (kWh)
0  battery_phs  2025-01-01             699868500.0
1  battery_phs  2030-01-01             689017872.0
2  battery_phs  2035-01-01             641971765.0
3  battery_phs  2040-01-01             624877169.0
4  battery_phs  2045-01-01             581134595.0

In [18]:

df_capacity = (
    model.results.flow_cap_new.where(model.results.techs != "demand_electricity")
    .sel(carriers="electricity")
    .sum("nodes")
    .to_series()
    .where(lambda x: x != 0)
    .dropna()
    .to_frame("New flow capacity (kW)")
    .reset_index()
)

print(df_capacity.head())

fig = px.bar(
    df_capacity,
    x="vintagesteps",
    y="New flow capacity (kW)",
    color="techs",
    color_discrete_map=model.inputs.color.to_series().to_dict(),
)
fig.show()

df_capacity = ( model.results.flow_cap_new.where(model.results.techs != "demand_electricity") .sel(carriers="electricity") .sum("nodes") .to_series() .where(lambda x: x != 0) .dropna() .to_frame("New flow capacity (kW)") .reset_index() ) print(df_capacity.head()) fig = px.bar( df_capacity, x="vintagesteps", y="New flow capacity (kW)", color="techs", color_discrete_map=model.inputs.color.to_series().to_dict(), ) fig.show()

                techs vintagesteps  New flow capacity (kW)
0  export_electricity   2025-01-01            1.327492e+07
1  export_electricity   2030-01-01            3.260194e-09
2  export_electricity   2050-01-01           -1.223068e-09
3  import_electricity   2025-01-01            1.327492e+07
4                 oil   2035-01-01            1.077884e+06

In [19]:

df_outflow = (
    (model.results.flow_out.fillna(0) - model.results.flow_in.fillna(0))
    .sel(carriers="electricity")
    .sum(["nodes", "timesteps"], min_count=1)
    .to_series()
    .where(lambda x: x > 1)
    .dropna()
    .to_frame("Annual outflow (kWh)")
    .reset_index()
)

print(df_capacity.head())

fig = px.bar(
    df_outflow,
    x="investsteps",
    y="Annual outflow (kWh)",
    color="techs",
    color_discrete_map=model.inputs.color.to_series().to_dict(),
)
df_demand = (
    model.results.flow_in.sel(techs="demand_electricity", carriers="electricity")
    .sum(["nodes", "timesteps"])
    .to_series()
    .reset_index()
)
fig.add_scatter(
    x=df_demand.investsteps, y=df_demand.flow_in, line={"color": "black"}, name="Demand"
)
fig.show()

df_outflow = ( (model.results.flow_out.fillna(0) - model.results.flow_in.fillna(0)) .sel(carriers="electricity") .sum(["nodes", "timesteps"], min_count=1) .to_series() .where(lambda x: x > 1) .dropna() .to_frame("Annual outflow (kWh)") .reset_index() ) print(df_capacity.head()) fig = px.bar( df_outflow, x="investsteps", y="Annual outflow (kWh)", color="techs", color_discrete_map=model.inputs.color.to_series().to_dict(), ) df_demand = ( model.results.flow_in.sel(techs="demand_electricity", carriers="electricity") .sum(["nodes", "timesteps"]) .to_series() .reset_index() ) fig.add_scatter( x=df_demand.investsteps, y=df_demand.flow_in, line={"color": "black"}, name="Demand" ) fig.show()

                techs vintagesteps  New flow capacity (kW)
0  export_electricity   2025-01-01            1.327492e+07
1  export_electricity   2030-01-01            3.260194e-09
2  export_electricity   2050-01-01           -1.223068e-09
3  import_electricity   2025-01-01            1.327492e+07
4                 oil   2035-01-01            1.077884e+06

In [20]:

df_electricity = (
    (model.results.flow_out.fillna(0) - model.results.flow_in.fillna(0))
    .sel(carriers="electricity")
    .sum("nodes")
    .to_series()
    .where(lambda x: x != 0)
    .dropna()
    .to_frame("Flow in/out (kWh)")
    .reset_index()
)
df_electricity_demand = df_electricity[df_electricity.techs == "demand_electricity"]
df_electricity_other = df_electricity[df_electricity.techs != "demand_electricity"]

invest_order = sorted(df_electricity.investsteps.unique())

fig = px.bar(
    df_electricity_other,
    x="timesteps",
    y="Flow in/out (kWh)",
    facet_row="investsteps",
    color="techs",
    height=1000,
    category_orders={"investsteps": invest_order},
    color_discrete_map=model.inputs.color.to_series().to_dict(),
)

showlegend = True
# we reverse the investment year order (`[::-1]`) because the rows are numbered from bottom to top.
for row, year in enumerate(invest_order[::-1]):
    demand_ = df_electricity_demand.loc[(df_electricity_demand.investsteps == year)]
    fig.add_scatter(
        x=demand_["timesteps"],
        y=-1 * demand_["Flow in/out (kWh)"],
        row=row + 1,
        col="all",
        marker_color="black",
        name="Demand",
        legendgroup="demand",
        showlegend=showlegend,
    )
    showlegend = False
fig.update_yaxes(matches=None)
fig.show()

df_electricity = ( (model.results.flow_out.fillna(0) - model.results.flow_in.fillna(0)) .sel(carriers="electricity") .sum("nodes") .to_series() .where(lambda x: x != 0) .dropna() .to_frame("Flow in/out (kWh)") .reset_index() ) df_electricity_demand = df_electricity[df_electricity.techs == "demand_electricity"] df_electricity_other = df_electricity[df_electricity.techs != "demand_electricity"] invest_order = sorted(df_electricity.investsteps.unique()) fig = px.bar( df_electricity_other, x="timesteps", y="Flow in/out (kWh)", facet_row="investsteps", color="techs", height=1000, category_orders={"investsteps": invest_order}, color_discrete_map=model.inputs.color.to_series().to_dict(), ) showlegend = True # we reverse the investment year order (`[::-1]`) because the rows are numbered from bottom to top. for row, year in enumerate(invest_order[::-1]): demand_ = df_electricity_demand.loc[(df_electricity_demand.investsteps == year)] fig.add_scatter( x=demand_["timesteps"], y=-1 * demand_["Flow in/out (kWh)"], row=row + 1, col="all", marker_color="black", name="Demand", legendgroup="demand", showlegend=showlegend, ) showlegend = False fig.update_yaxes(matches=None) fig.show()