Pathway Planning¶

In this example, we demonstrate how PyPSA can deal with optimisation problems spanning multiple investment periods, also known as pathway planning.

For models with multiple investment periods, the total set of snapshots is divided into investment periods, which translates into multi-indexed snapshots with investment periods as the first level time steps as the second level. In each investment period, new components may be added to the system. Additionally, any component may only operate as long as allowed by their lifetime.

In contrast to the models with a single investment period (overnight scenarios), the following concepts have to be taken into account.:

1. The network attribute n.investment_periods: This is the set of periods which specify when new components may be built. These have to be the same as the first level values in the n.snapshots index.
2. The network attribute n.investment_period_weightings: These specify the weighting of each period in the objective function and the global constraints.
3. The component attribute build_year: Any one component may only be built when the build year is equal to the current investment period, or larger than the previous investment period. That means, components with build_year=2029 are considered in the investment period 2030, but not in the period 2025.
4. The component attribute lifetime: Any one component is only considered for dispatch in an investment period if it is still active at the beginning of an investment period. That means, components with build_year=2029 and lifetime=30 are considered in the investment period 2055, but not in the period 2060.

In the following, we set up a three bus network with generators, lines and storage units and run a optimisation for the investment periods 2020, 2030, 2040 and 2050.

In [2]:

import numpy as np
import pandas as pd

import pypsa

rng = np.random.default_rng()  # Create a random number generator

import numpy as np import pandas as pd import pypsa rng = np.random.default_rng() # Create a random number generator

We set up the network with investment periods and snapshots.

In [3]:

n = pypsa.Network()
years = [2020, 2030, 2040, 2050]
freq = 24

snapshots = pd.DatetimeIndex([])
for year in years:
    period = pd.date_range(
        start=f"{year}-01-01 00:00",
        freq=f"{freq}h",
        periods=8760 // freq,
    )
    snapshots = snapshots.append(period)

# convert to multiindex and assign to network
n.snapshots = pd.MultiIndex.from_arrays([snapshots.year, snapshots])
n.investment_periods = years

n.snapshot_weightings

n = pypsa.Network() years = [2020, 2030, 2040, 2050] freq = 24 snapshots = pd.DatetimeIndex([]) for year in years: period = pd.date_range( start=f"{year}-01-01 00:00", freq=f"{freq}h", periods=8760 // freq, ) snapshots = snapshots.append(period) # convert to multiindex and assign to network n.snapshots = pd.MultiIndex.from_arrays([snapshots.year, snapshots]) n.investment_periods = years n.snapshot_weightings

Out[3]:

		objective	stores	generators
period	timestep
2020	2020-01-01	1.0	1.0	1.0
	2020-01-02	1.0	1.0	1.0
	2020-01-03	1.0	1.0	1.0
	2020-01-04	1.0	1.0	1.0
	2020-01-05	1.0	1.0	1.0
...	...	...	...	...
2050	2050-12-27	1.0	1.0	1.0
	2050-12-28	1.0	1.0	1.0
	2050-12-29	1.0	1.0	1.0
	2050-12-30	1.0	1.0	1.0
	2050-12-31	1.0	1.0	1.0

1460 rows × 3 columns

In [4]:

n.investment_periods

n.investment_periods

Out[4]:

Index([2020, 2030, 2040, 2050], dtype='int32', name='period')

Set the years and objective weighting per investment period. For the objective weighting, we consider a discount rate defined by $$ D(t) = \dfrac{1}{(1+r)^t} $$

where $r$ is the discount rate. For each period we sum up all discounts rates of the corresponding years which gives us the effective objective weighting.

In [5]:

n.investment_period_weightings["years"] = list(np.diff(years)) + [10]

r = 0.01
T = 0
for period, nyears in n.investment_period_weightings.years.items():
    discounts = [(1 / (1 + r) ** t) for t in range(T, T + nyears)]
    n.investment_period_weightings.at[period, "objective"] = sum(discounts)
    T += nyears
n.investment_period_weightings

n.investment_period_weightings["years"] = list(np.diff(years)) + [10] r = 0.01 T = 0 for period, nyears in n.investment_period_weightings.years.items(): discounts = [(1 / (1 + r) ** t) for t in range(T, T + nyears)] n.investment_period_weightings.at[period, "objective"] = sum(discounts) T += nyears n.investment_period_weightings

Out[5]:

	objective	years
period
2020	9.566018	10
2030	8.659991	10
2040	7.839777	10
2050	7.097248	10

Add the components

In [6]:

for i in range(3):
    n.add("Bus", f"bus {i}")

# add three lines in a ring
n.add(
    "Line",
    "line 0->1",
    bus0="bus 0",
    bus1="bus 1",
)

n.add(
    "Line",
    "line 1->2",
    bus0="bus 1",
    bus1="bus 2",
    capital_cost=10,
    build_year=2030,
)

n.add(
    "Line",
    "line 2->0",
    bus0="bus 2",
    bus1="bus 0",
)

n.lines["x"] = 0.0001
n.lines["s_nom_extendable"] = True

n.lines

for i in range(3): n.add("Bus", f"bus {i}") # add three lines in a ring n.add( "Line", "line 0->1", bus0="bus 0", bus1="bus 1", ) n.add( "Line", "line 1->2", bus0="bus 1", bus1="bus 2", capital_cost=10, build_year=2030, ) n.add( "Line", "line 2->0", bus0="bus 2", bus1="bus 0", ) n.lines["x"] = 0.0001 n.lines["s_nom_extendable"] = True n.lines

Out[6]:

	bus0	bus1	type	x	r	g	b	s_nom	s_nom_mod	s_nom_extendable	...	v_ang_min	v_ang_max	sub_network	x_pu	r_pu	g_pu	b_pu	x_pu_eff	r_pu_eff	s_nom_opt
name
line 0->1	bus 0	bus 1		0.0001	0.0	0.0	0.0	0.0	0.0	True	...	-inf	inf		0.0	0.0	0.0	0.0	0.0	0.0	0.0
line 1->2	bus 1	bus 2		0.0001	0.0	0.0	0.0	0.0	0.0	True	...	-inf	inf		0.0	0.0	0.0	0.0	0.0	0.0	0.0
line 2->0	bus 2	bus 0		0.0001	0.0	0.0	0.0	0.0	0.0	True	...	-inf	inf		0.0	0.0	0.0	0.0	0.0	0.0	0.0

3 rows × 35 columns

In [7]:

# add some generators
p_nom_max = pd.Series(
    (rng.uniform() for sn in range(len(n.snapshots))),
    index=n.snapshots,
    name="generator ext 2020",
)

# renewable (can operate 2020, 2030)
n.add(
    "Generator",
    "generator ext 0 2020",
    bus="bus 0",
    p_nom=50,
    build_year=2020,
    lifetime=20,
    marginal_cost=2,
    capital_cost=1,
    p_max_pu=p_nom_max,
    carrier="solar",
    p_nom_extendable=True,
)

# can operate 2040, 2050
n.add(
    "Generator",
    "generator ext 0 2040",
    bus="bus 0",
    p_nom=50,
    build_year=2040,
    lifetime=11,
    marginal_cost=25,
    capital_cost=10,
    carrier="OCGT",
    p_nom_extendable=True,
)

# can operate in 2040
n.add(
    "Generator",
    "generator fix 1 2040",
    bus="bus 1",
    p_nom=50,
    build_year=2040,
    lifetime=10,
    carrier="CCGT",
    marginal_cost=20,
    capital_cost=1,
)

n.generators

# add some generators p_nom_max = pd.Series( (rng.uniform() for sn in range(len(n.snapshots))), index=n.snapshots, name="generator ext 2020", ) # renewable (can operate 2020, 2030) n.add( "Generator", "generator ext 0 2020", bus="bus 0", p_nom=50, build_year=2020, lifetime=20, marginal_cost=2, capital_cost=1, p_max_pu=p_nom_max, carrier="solar", p_nom_extendable=True, ) # can operate 2040, 2050 n.add( "Generator", "generator ext 0 2040", bus="bus 0", p_nom=50, build_year=2040, lifetime=11, marginal_cost=25, capital_cost=10, carrier="OCGT", p_nom_extendable=True, ) # can operate in 2040 n.add( "Generator", "generator fix 1 2040", bus="bus 1", p_nom=50, build_year=2040, lifetime=10, carrier="CCGT", marginal_cost=20, capital_cost=1, ) n.generators

Out[7]:

	bus	control	type	p_nom	p_nom_mod	p_nom_extendable	p_nom_min	p_nom_max	p_nom_set	p_min_pu	...	min_up_time	min_down_time	up_time_before	down_time_before	ramp_limit_up	ramp_limit_down	ramp_limit_start_up	ramp_limit_shut_down	weight	p_nom_opt
name
generator ext 0 2020	bus 0	PQ		50.0	0.0	True	0.0	inf	NaN	0.0	...	0	0	1	0	NaN	NaN	NaN	NaN	1.0	0.0
generator ext 0 2040	bus 0	PQ		50.0	0.0	True	0.0	inf	NaN	0.0	...	0	0	1	0	NaN	NaN	NaN	NaN	1.0	0.0
generator fix 1 2040	bus 1	PQ		50.0	0.0	False	0.0	inf	NaN	0.0	...	0	0	1	0	NaN	NaN	NaN	NaN	1.0	0.0

3 rows × 42 columns

In [8]:

n.add(
    "StorageUnit",
    "storageunit non-cyclic 2030",
    bus="bus 2",
    p_nom=0,
    capital_cost=2,
    build_year=2030,
    lifetime=21,
    cyclic_state_of_charge=False,
    p_nom_extendable=False,
)

n.add(
    "StorageUnit",
    "storageunit periodic 2020",
    bus="bus 2",
    p_nom=0,
    capital_cost=1,
    build_year=2020,
    lifetime=21,
    cyclic_state_of_charge=True,
    cyclic_state_of_charge_per_period=True,
    p_nom_extendable=True,
)

n.storage_units

n.add( "StorageUnit", "storageunit non-cyclic 2030", bus="bus 2", p_nom=0, capital_cost=2, build_year=2030, lifetime=21, cyclic_state_of_charge=False, p_nom_extendable=False, ) n.add( "StorageUnit", "storageunit periodic 2020", bus="bus 2", p_nom=0, capital_cost=1, build_year=2020, lifetime=21, cyclic_state_of_charge=True, cyclic_state_of_charge_per_period=True, p_nom_extendable=True, ) n.storage_units

Out[8]:

	bus	control	type	p_nom	p_nom_mod	p_nom_extendable	p_nom_min	p_nom_max	p_nom_set	p_min_pu	...	state_of_charge_initial_per_period	state_of_charge_set	cyclic_state_of_charge	cyclic_state_of_charge_per_period	max_hours	efficiency_store	efficiency_dispatch	standing_loss	inflow	p_nom_opt
name
storageunit non-cyclic 2030	bus 2	PQ		0.0	0.0	False	0.0	inf	NaN	-1.0	...	False	NaN	False	False	1.0	1.0	1.0	0.0	0.0	0.0
storageunit periodic 2020	bus 2	PQ		0.0	0.0	True	0.0	inf	NaN	-1.0	...	False	NaN	True	True	1.0	1.0	1.0	0.0	0.0	0.0

2 rows × 39 columns

Add the load

In [9]:

load_var = pd.Series(
    100 * rng.random(size=len(n.snapshots)), index=n.snapshots, name="load"
)
n.add("Load", "load 2", bus="bus 2", p_set=load_var)

load_fix = pd.Series(75, index=n.snapshots, name="load")
n.add("Load", "load 1", bus="bus 1", p_set=load_fix)

n.loads_t.p_set.head()

load_var = pd.Series( 100 * rng.random(size=len(n.snapshots)), index=n.snapshots, name="load" ) n.add("Load", "load 2", bus="bus 2", p_set=load_var) load_fix = pd.Series(75, index=n.snapshots, name="load") n.add("Load", "load 1", bus="bus 1", p_set=load_fix) n.loads_t.p_set.head()

Out[9]:

	name	load 2	load 1
period	timestep
2020	2020-01-01	16.466597	75.0
	2020-01-02	72.529805	75.0
	2020-01-03	8.801329	75.0
	2020-01-04	63.074000	75.0
	2020-01-05	31.176991	75.0

Run the optimization

In [10]:

n.optimize(multi_investment_periods=True)

n.optimize(multi_investment_periods=True)

/tmp/ipykernel_5247/264390566.py:1: 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.
  n.optimize(multi_investment_periods=True)
WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['bus 0', 'bus 1', 'bus 2'], dtype='object', name='name')

WARNING:pypsa.consistency:The following generators have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['generator ext 0 2020', 'generator ext 0 2040', 'generator fix 1 2040'], dtype='object', name='name')

WARNING:pypsa.consistency:The following lines have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['line 0->1', 'line 1->2', 'line 2->0'], dtype='object', name='name')

WARNING:pypsa.consistency:The following lines have zero r, which could break the linear load flow:
Index(['line 0->1', 'line 1->2', 'line 2->0'], dtype='object', name='name')

WARNING:pypsa.optimization.constraints:StorageUnits ['storageunit periodic 2020']: Per-period cyclic (cyclic_state_of_charge_per_period=True) overrides global cyclic (cyclic_state_of_charge=True). Storage will cycle within each investment period, not across the entire horizon.

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/24 [00:00<?, ?it/s]

Writing constraints.:  88%|████████▊ | 21/24 [00:00<00:00, 203.33it/s]

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

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

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

INFO:linopy.io: Writing time: 0.17s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 12416 primals, 32491 duals
Objective: 1.81e+07
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, Line-ext-s-lower, Line-ext-s-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, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.

Out[10]:

('ok', 'optimal')

In [11]:

c = "Generator"
df = pd.concat(
    {
        period: n.components[c].get_active_assets(period)
        * n.components[c].static.p_nom_opt
        for period in n.investment_periods
    },
    axis=1,
)
df.T.plot.bar(
    stacked=True,
    edgecolor="white",
    width=1,
    ylabel="Capacity (MW)",
    xlabel="Investment Period",
    rot=0,
)

c = "Generator" df = pd.concat( { period: n.components[c].get_active_assets(period) * n.components[c].static.p_nom_opt for period in n.investment_periods }, axis=1, ) df.T.plot.bar( stacked=True, edgecolor="white", width=1, ylabel="Capacity (MW)", xlabel="Investment Period", rot=0, )

Out[11]:

<Axes: xlabel='Investment Period', ylabel='Capacity (MW)'>

In [12]:

df = n.generators_t.p.sum(axis=0).T.div(1e3)
df.T.plot.bar(
    stacked=True,
    edgecolor="white",
    width=1,
    ylabel="Generation (GWh)",
    xlabel="Investment Period",
    rot=0,
)

df = n.generators_t.p.sum(axis=0).T.div(1e3) df.T.plot.bar( stacked=True, edgecolor="white", width=1, ylabel="Generation (GWh)", xlabel="Investment Period", rot=0, )

Out[12]:

<Axes: xlabel='Investment Period', ylabel='Generation (GWh)'>