Myopic Pathway Planning¶

In this example, we demonstrate how to use PyPSA to plan a myopic pathway for emission reductions. That means that we optimize the system for each investment period sequentially subject to set emission reduction targets per period. In consequence, the model does not foresee future costs or benefits of investments. This is in contrast to a pathway optimisation with perfect foresight, where the model optimizes the system for all investment periods at once.

Here, we build a model from the single-node capacity expansion example and model a myopic decision process for the investment periods 2025, 2035, and 2045. The model is set up to reduce emissions from 30 MtCO2/a in 2025 to 10 MtCO2/a in 2035 and 3 MtCO2/a in 2045. Note that the goal here is not necessarily to showcase a realistic scenario, but rather to demonstrate how it is possible to model a myopic decision process in PyPSA.

In [2]:

import pandas as pd

import pypsa
from pypsa.costs import annuity

temp = pypsa.examples.model_energy()
temp.remove("Generator", "load shedding")
n = temp.copy()

INVESTMENT_PERIODS = [2025, 2035, 2045]

import pandas as pd import pypsa from pypsa.costs import annuity temp = pypsa.examples.model_energy() temp.remove("Generator", "load shedding") n = temp.copy() INVESTMENT_PERIODS = [2025, 2035, 2045]

INFO:pypsa.network.io:Imported network 'Model-Energy' has buses, carriers, generators, links, loads, storage_units, stores

First, we add a new component for each investment period and carrier to the network. These components are used to model the individual investment decisions for each period. For instance, for solar, there will be a Generator component for each investment period, e.g. solar-2025, solar-2035, and solar-2045 with the corresponding build_year attributes.

We can also keep pre-existing capacities as non-extendable components. Here, we keep 5 GW each of wind and solar capacities, assuming 2015 as build year.

In [3]:

# Assume 5 GW of wind and solar capacity each, built in 2015
n.generators.loc[:, ["build_year", "p_nom", "p_nom_extendable"]] = 2015, 5_000, False

# all other storage components need to be removed
n.remove("Link", n.links.index)
n.remove("StorageUnit", n.storage_units.index)
n.remove("Store", n.stores.index)

for year in INVESTMENT_PERIODS:
    n.add(
        "Generator",
        temp.generators.index + f"-{year}",
        build_year=year,
        p_max_pu=temp.generators_t.p_max_pu.rename(
            columns=lambda s: s + f"-{year}"
        ).loc[:, ::-1],
        **temp.generators.rename(index=lambda s: s + f"-{year}").drop(
            columns=["build_year", "p_max_pu"]
        ),
    )

    n.add(
        "Link",
        temp.links.index,
        suffix=f"-{year}",
        build_year=year,
        **temp.links.rename(index=lambda s: s + f"-{year}").drop(
            columns=["build_year"]
        ),
    )

    n.add(
        "StorageUnit",
        temp.storage_units.index,
        suffix=f"-{year}",
        build_year=year,
        **temp.storage_units.rename(index=lambda s: s + f"-{year}").drop(
            columns=["build_year"]
        ),
    )

    n.add(
        "Store",
        temp.stores.index,
        suffix=f"-{year}",
        build_year=year,
        **temp.stores.rename(index=lambda s: s + f"-{year}").drop(
            columns=["build_year"]
        ),
    )

# Assume 5 GW of wind and solar capacity each, built in 2015 n.generators.loc[:, ["build_year", "p_nom", "p_nom_extendable"]] = 2015, 5_000, False # all other storage components need to be removed n.remove("Link", n.links.index) n.remove("StorageUnit", n.storage_units.index) n.remove("Store", n.stores.index) for year in INVESTMENT_PERIODS: n.add( "Generator", temp.generators.index + f"-{year}", build_year=year, p_max_pu=temp.generators_t.p_max_pu.rename( columns=lambda s: s + f"-{year}" ).loc[:, ::-1], **temp.generators.rename(index=lambda s: s + f"-{year}").drop( columns=["build_year", "p_max_pu"] ), ) n.add( "Link", temp.links.index, suffix=f"-{year}", build_year=year, **temp.links.rename(index=lambda s: s + f"-{year}").drop( columns=["build_year"] ), ) n.add( "StorageUnit", temp.storage_units.index, suffix=f"-{year}", build_year=year, **temp.storage_units.rename(index=lambda s: s + f"-{year}").drop( columns=["build_year"] ), ) n.add( "Store", temp.stores.index, suffix=f"-{year}", build_year=year, **temp.stores.rename(index=lambda s: s + f"-{year}").drop( columns=["build_year"] ), )

In this way, the attributes of components can vary over the investment periods, e.g. the capital_cost of solar PV and battery storage can decrease over time, or the capacity factors (p_max_pu) of wind generators can increase over time.

In [4]:

n.generators.loc["solar-2035", "capital_cost"] *= 0.9
n.generators.loc["solar-2045", "capital_cost"] *= 0.85

n.storage_units.loc["battery storage-2035", "capital_cost"] *= 0.9
n.storage_units.loc["battery storage-2045", "capital_cost"] *= 0.85

n.generators_t.p_max_pu.loc[:, "wind-2035"] *= 1.1
n.generators_t.p_max_pu.loc[:, "wind-2045"] *= 1.15

n.generators.loc["solar-2035", "capital_cost"] *= 0.9 n.generators.loc["solar-2045", "capital_cost"] *= 0.85 n.storage_units.loc["battery storage-2035", "capital_cost"] *= 0.9 n.storage_units.loc["battery storage-2045", "capital_cost"] *= 0.85 n.generators_t.p_max_pu.loc[:, "wind-2035"] *= 1.1 n.generators_t.p_max_pu.loc[:, "wind-2045"] *= 1.15

Let us also add some pre-existing conventional generators to the network with emissions from burning fossil fuels.

We add some old coal generators, which will retire by 2035, and some gas generators, which will retire by 2045. We also let the model invest in new gas generators in 2035, with slightly higher cost but better efficiency.

In [5]:

n.add(
    "Carrier",
    ["coal", "gas"],
    color=["gray", "lightcoral"],
    co2_emissions=[0.336, 0.198],
)

n.add(
    "Generator",
    "coal",
    carrier="coal",
    bus="electricity",
    p_nom=4_000,
    p_nom_extendable=False,
    build_year=1990,
    lifetime=40,
    capital_cost=annuity(0.07, 40) * 4_000_000,
    efficiency=0.35,
    marginal_cost=30,
)

n.add(
    "Generator",
    "gas",
    carrier="gas",
    bus="electricity",
    p_nom=6_000,
    p_nom_extendable=False,
    build_year=2005,
    lifetime=35,
    capital_cost=annuity(0.07, 35) * 800_000,
    efficiency=0.5,
    marginal_cost=60,
)

n.add(
    "Generator",
    "gas-2035",
    carrier="gas",
    bus="electricity",
    p_nom_extendable=True,
    build_year=2035,
    lifetime=35,
    capital_cost=annuity(0.07, 35) * 900_000,
    efficiency=0.55,
    marginal_cost=55,
)

n.add( "Carrier", ["coal", "gas"], color=["gray", "lightcoral"], co2_emissions=[0.336, 0.198], ) n.add( "Generator", "coal", carrier="coal", bus="electricity", p_nom=4_000, p_nom_extendable=False, build_year=1990, lifetime=40, capital_cost=annuity(0.07, 40) * 4_000_000, efficiency=0.35, marginal_cost=30, ) n.add( "Generator", "gas", carrier="gas", bus="electricity", p_nom=6_000, p_nom_extendable=False, build_year=2005, lifetime=35, capital_cost=annuity(0.07, 35) * 800_000, efficiency=0.5, marginal_cost=60, ) n.add( "Generator", "gas-2035", carrier="gas", bus="electricity", p_nom_extendable=True, build_year=2035, lifetime=35, capital_cost=annuity(0.07, 35) * 900_000, efficiency=0.55, marginal_cost=55, )

In the next step, we define the investment periods and the emission reduction targets for each period. Note that each investment period represents a duration of 10 years, hence the annual targets are scaled to represent total emissions over the entire period.

In [6]:

n.set_investment_periods(INVESTMENT_PERIODS)
n.investment_period_weightings *= 10

REDUCTION_PATH = [300e6, 100e6, 20e6]

for i, year in enumerate(INVESTMENT_PERIODS):
    n.add(
        "GlobalConstraint",
        f"co2-limit-{year}",
        type="primary_energy",
        carrier_attribute="co2_emissions",
        investment_period=year,
        constant=REDUCTION_PATH[i],
        sense="<=",
    )

n.set_investment_periods(INVESTMENT_PERIODS) n.investment_period_weightings *= 10 REDUCTION_PATH = [300e6, 100e6, 20e6] for i, year in enumerate(INVESTMENT_PERIODS): n.add( "GlobalConstraint", f"co2-limit-{year}", type="primary_energy", carrier_attribute="co2_emissions", investment_period=year, constant=REDUCTION_PATH[i], sense="<=", )

INFO:pypsa.network.index:Repeating time-series for each investment period and converting snapshots to a pandas.MultiIndex.

To optimize only a single investment period at a time, we have to make a selection of the relevant snapshots (i.e. all that fall into the investment period). Passing this selection to n.optimize() will exclude investment variables of any inactive components (where the build year falls into future investment periods), and also only optimise the dispatch variables for the time steps in the investment period.

In [7]:

def optimize_period(n: pypsa.Network, period: int):
    snapshots = n.snapshots[n.snapshots.get_level_values("period") == period]
    return n.optimize(
        multi_investment_periods=True,
        snapshots=snapshots,
    )


optimize_period(n, 2025)

display(n.generators.p_nom_opt)

display(n.global_constraints)

def optimize_period(n: pypsa.Network, period: int): snapshots = n.snapshots[n.snapshots.get_level_values("period") == period] return n.optimize( multi_investment_periods=True, snapshots=snapshots, ) optimize_period(n, 2025) display(n.generators.p_nom_opt) display(n.global_constraints)

/tmp/ipykernel_5294/36310589.py:3: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  return n.optimize(

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io:Writing objective.

Writing constraints.:   0%|          | 0/22 [00:00<?, ?it/s]

Writing constraints.:  77%|███████▋  | 17/22 [00:00<00:00, 164.33it/s]

Writing constraints.: 100%|██████████| 22/22 [00:00<00:00, 147.55it/s]

Writing continuous variables.:   0%|          | 0/11 [00:00<?, ?it/s]

Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 318.23it/s]

INFO:linopy.io: Writing time: 0.2s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 37979 primals, 81780 duals
Objective: 2.05e+10
Solver model: available
Solver message: Optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.

name
wind          5000.000000
solar         5000.000000
wind-2025     4361.538462
solar-2025    5351.989786
wind-2035       -0.000000
solar-2035      -0.000000
wind-2045       -0.000000
solar-2045      -0.000000
coal          4000.000000
gas           6000.000000
gas-2035        -0.000000
Name: p_nom_opt, dtype: float64

	type	investment_period	bus	carrier_attribute	sense	constant	mu
name
co2-limit-2025	primary_energy	2025.0		co2_emissions	<=	300000000.0	-20.011164
co2-limit-2035	primary_energy	2035.0		co2_emissions	<=	100000000.0	0.000000
co2-limit-2045	primary_energy	2045.0		co2_emissions	<=	20000000.0	0.000000

In the results we can see that any components after 2025 are not touched in the first investment period. The CO2 constraint is mildly binding at 20 €/tCO2, which is the price of CO2 in the first period. The energy mix is still dominated by coal and gas, with some solar and wind added to the system.

In [8]:

n.statistics.energy_balance().div(1e6).round(2)

n.statistics.energy_balance().div(1e6).round(2)

Out[8]:

		period	2025	2035	2045
component	carrier	bus_carrier
Generator	coal	electricity	27.89	NaN	NaN
	gas	electricity	8.14	NaN	NaN
	solar	electricity	12.10	NaN	NaN
	wind	electricity	18.15	NaN	NaN
Link	electrolysis	electricity	NaN	NaN	NaN
	turbine	hydrogen	NaN	NaN	NaN
	electrolysis	hydrogen	NaN	NaN	NaN
	turbine	electricity	NaN	NaN	NaN
Load	-	electricity	-66.27	NaN	NaN
StorageUnit	battery storage	electricity	-0.01	NaN	NaN
Store	hydrogen storage	hydrogen	NaN	NaN	NaN

For the second investment period, we need to freeze the capacities that were built in the first period, so that they are not changed in the second period (see freeze_period() function below). Then we can again select the relevant snapshots and optimize the system for the second investment period -- and so on for any further investment periods.

In [9]:

def freeze_period(n: pypsa.Network, period: int):
    for c in n.components[["Generator", "Link", "StorageUnit", "Store"]]:
        attr = "e_nom" if c.name == "Store" else "p_nom"
        c.static[attr] = c.static[attr + "_opt"]
        c.static.loc[c.static.build_year == period, attr + "_extendable"] = False


freeze_period(n, 2025)
optimize_period(n, 2035)

freeze_period(n, 2035)
optimize_period(n, 2045)

def freeze_period(n: pypsa.Network, period: int): for c in n.components[["Generator", "Link", "StorageUnit", "Store"]]: attr = "e_nom" if c.name == "Store" else "p_nom" c.static[attr] = c.static[attr + "_opt"] c.static.loc[c.static.build_year == period, attr + "_extendable"] = False freeze_period(n, 2025) optimize_period(n, 2035) freeze_period(n, 2035) optimize_period(n, 2045)

/tmp/ipykernel_5294/36310589.py:3: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  return n.optimize(

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io:Writing objective.

Writing constraints.:   0%|          | 0/32 [00:00<?, ?it/s]

Writing constraints.:  56%|█████▋    | 18/32 [00:00<00:00, 172.39it/s]

Writing constraints.: 100%|██████████| 32/32 [00:00<00:00, 149.05it/s]

Writing continuous variables.:   0%|          | 0/11 [00:00<?, ?it/s]

Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 288.93it/s]

INFO:linopy.io: Writing time: 0.26s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 64253 primals, 134334 duals
Objective: 2.32e+10
Solver model: available
Solver message: Optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Link-fix-p-lower, Link-fix-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-fix-e-lower, Store-fix-e-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.

/tmp/ipykernel_5294/36310589.py:3: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  return n.optimize(

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io:Writing objective.

Writing constraints.:   0%|          | 0/32 [00:00<?, ?it/s]

Writing constraints.:  50%|█████     | 16/32 [00:00<00:00, 151.14it/s]

Writing constraints.: 100%|██████████| 32/32 [00:00<00:00, 128.00it/s]

Writing constraints.: 100%|██████████| 32/32 [00:00<00:00, 130.53it/s]

Writing continuous variables.:   0%|          | 0/11 [00:00<?, ?it/s]

Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 261.58it/s]

INFO:linopy.io: Writing time: 0.3s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 87606 primals, 181047 duals
Objective: 3.33e+10
Solver model: available
Solver message: Optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Link-fix-p-lower, Link-fix-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-fix-e-lower, Store-fix-e-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.

Out[9]:

('ok', 'optimal')

After the loop has finished, we can inspect the results of all investment periods. For instance, we can see how the endogenous CO2 price increases with each investment period.

In [10]:

n.global_constraints

n.global_constraints

Out[10]:

	type	investment_period	bus	carrier_attribute	sense	constant	mu
name
co2-limit-2025	primary_energy	2025.0		co2_emissions	<=	300000000.0	-20.011164
co2-limit-2035	primary_energy	2035.0		co2_emissions	<=	100000000.0	-32.857033
co2-limit-2045	primary_energy	2045.0		co2_emissions	<=	20000000.0	-162.180346

The energy balance and CAPEX statistics also tell us the story that the coal power plants that are phased out in 2035 are largely replaced by a large amount of new gas power plants. Only in 2045, the model starts to double down on renewables and storage investments to meet more ambitious emission reduction targets. By then, the new gas power plants are barely used, while their capital costs are still to be paid off.

In [11]:

display(n.statistics.energy_balance().div(1e6).round(1))

display(n.statistics.capex().div(1e6).round(1))

display(n.statistics.energy_balance().div(1e6).round(1)) display(n.statistics.capex().div(1e6).round(1))

		period	2025	2035	2045
component	carrier	bus_carrier
Generator	coal	electricity	27.9	NaN	NaN
	gas	electricity	8.1	27.2	5.6
	solar	electricity	12.1	13.7	28.4
	wind	electricity	18.1	25.4	35.2
Link	electrolysis	electricity	NaN	NaN	-2.5
	turbine	hydrogen	NaN	NaN	-1.6
	electrolysis	hydrogen	NaN	NaN	1.6
	turbine	electricity	NaN	NaN	0.6
Load	-	electricity	-66.3	-66.3	-66.3
StorageUnit	battery storage	electricity	-0.0	-0.0	-0.9
Store	hydrogen storage	hydrogen	NaN	NaN	NaN

		2025	2035	2045
component	carrier
Generator	coal	12001.5	NaN	NaN
	gas	3707.2	6329.8	2622.6
	solar	5315.4	6259.5	13244.2
	wind	9515.5	13717.8	20120.2
Link	electrolysis	NaN	NaN	885.2
	turbine	NaN	NaN	3572.3
StorageUnit	battery storage	83.8	271.9	12049.0
Store	hydrogen storage	NaN	NaN	592.8

Total annual system costs increase as emissions must be reduced. The data shown does not account for a social discount rate.

In [12]:

(n.statistics.capex().sum() + n.statistics.opex().sum()).div(1e6).round(1)

(n.statistics.capex().sum() + n.statistics.opex().sum()).div(1e6).round(1)

Out[12]:

2025    43873.4
2035    41828.4
2045    56141.8
dtype: float64

Some data wrangling is needed to plot the progression of operational generation capacities over time.

In [13]:

df = (
    pd.concat(
        {
            year: n.generators.query(
                f"build_year <= {year} and build_year + lifetime > {year}"
            )
            .groupby("carrier")
            .p_nom_opt.sum()
            for year in INVESTMENT_PERIODS
        },
        axis=1,
    )
    .fillna(0)
    .T.div(1e3)
)

df.plot.area(
    stacked=True,
    linewidth=0,
    ylabel="GW",
    xlabel="Year",
    title="Operational Capacity",
    color=df.columns.map(n.carriers.color),
)

df = ( pd.concat( { year: n.generators.query( f"build_year <= {year} and build_year + lifetime > {year}" ) .groupby("carrier") .p_nom_opt.sum() for year in INVESTMENT_PERIODS }, axis=1, ) .fillna(0) .T.div(1e3) ) df.plot.area( stacked=True, linewidth=0, ylabel="GW", xlabel="Year", title="Operational Capacity", color=df.columns.map(n.carriers.color), )

Out[13]:

<Axes: title={'center': 'Operational Capacity'}, xlabel='Year', ylabel='GW'>

Where to go from here?

A few ideas to explore further in this notebook:

• Try different exogenous emission reduction paths. How do they affect total costs and cumulative emissions?
• Decrease the capital cost of components over time in different ways. How does this affect the optimal build-out?
• Decrease the capital cost of components over time as a function of the cumulative installed capacity to simulate learning effects.
• Increase the capacity factors of wind and solar components over time to varying extents. How does this affect the optimal build-out?
• Compare the myopic pathway with a full pathway optimization. How do they differ in terms of costs, cumulative emissions, stranded assets? A code snippet is provided blow.

In [14]:

# for c in n.iterate_components({"Generator", "Link", "StorageUnit", "Store"}):
#     attr = "e_nom" if c.name == "Store" else "p_nom"
#     c.df.loc[c.df.build_year >= 2025, attr + "_extendable"] = True

# n.optimize(
#     multi_investment_periods=True,
#     log_to_console=True,
# )

# for c in n.iterate_components({"Generator", "Link", "StorageUnit", "Store"}): # attr = "e_nom" if c.name == "Store" else "p_nom" # c.df.loc[c.df.build_year >= 2025, attr + "_extendable"] = True # n.optimize( # multi_investment_periods=True, # log_to_console=True, # )