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 [ ]:

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 [ ]:

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

In [ ]:

n.investment_periods

n.investment_periods

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 [ ]:

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

Add the components

In [ ]:

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

In [ ]:

# 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

In [ ]:

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

Add the load

In [ ]:

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()

Run the optimization

In [ ]:

n.optimize(multi_investment_periods=True)

n.optimize(multi_investment_periods=True)

In [ ]:

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, )

In [ ]:

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, )