Multi Investment Optimization

Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Multi Investment Optimization#

In the following, we show how PyPSA can deal with multi-investment optimization, also known as multi-horizon optimization.

Here, the total set of snapshots is divided into investment periods. For the model, this translates into multi-indexed snapshots with the first level being the investment period and the second level the according time steps. In each investment period new asset may be added to the system. On the other hand assets may only operate as long as allowed by their lifetime.

In contrast to the ordinary optimisation, the following concepts have to be taken into account.

  1. investment_periods - pypsa.Network attribute. This is the set of periods which specify when new assets may be built. In the current implementation, these have to be the same as the first level values in the snapshots attribute.

  2. investment_period_weightings - pypsa.Network attribute. These specify the weighting of each period in the objective function.

  3. build_year - general component attribute. A single asset may only be built when the build year is smaller or equal to the current investment period. For example, assets with a build year 2029 are considered in the investment period 2030, but not in the period 2025.

  4. lifetime - general component attribute. An asset is only considered in an investment period if present at the beginning of an investment period. For example, an asset with build year 2029 and lifetime 30 is considered in the investment period 2055 but not in the period 2060.

In the following, we set up a three node network with generators, lines and storages and run a optimisation covering the time span from 2020 to 2050 and each decade is one investment period.

[1]:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import pypsa
ERROR 1: PROJ: proj_create_from_database: Open of /home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/share/proj failed

We set up the network with investment periods and snapshots.

[2]:
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 / float(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
/tmp/ipykernel_4098/981291686.py:7: FutureWarning: Non-integer 'periods' in pd.date_range, pd.timedelta_range, pd.period_range, and pd.interval_range are deprecated and will raise in a future version.
  period = pd.date_range(
/tmp/ipykernel_4098/981291686.py:7: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.
  period = pd.date_range(
[2]:
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

[3]:
n.investment_periods
[3]:
Index([2020, 2030, 2040, 2050], dtype='int64')

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.

[4]:
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
[4]:
objective years
2020 9.566018 10
2030 8.659991 10
2040 7.839777 10
2050 7.097248 10

Add the components

[5]:
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
[6]:
n.lines
[6]:
attribute 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
Line
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 × 30 columns

[7]:
# add some generators
p_nom_max = pd.Series(
    (np.random.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
[7]:
attribute bus control type p_nom p_nom_mod p_nom_extendable p_nom_min p_nom_max p_min_pu p_max_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
Generator
generator ext 0 2020 bus 0 PQ 50.0 0.0 True 0.0 inf 0.0 1.0 ... 0 0 1 0 NaN NaN 1.0 1.0 1.0 0.0
generator ext 0 2040 bus 0 PQ 50.0 0.0 True 0.0 inf 0.0 1.0 ... 0 0 1 0 NaN NaN 1.0 1.0 1.0 0.0
generator fix 1 2040 bus 1 PQ 50.0 0.0 False 0.0 inf 0.0 1.0 ... 0 0 1 0 NaN NaN 1.0 1.0 1.0 0.0

3 rows × 34 columns

[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
[8]:
attribute bus control type p_nom p_nom_mod p_nom_extendable p_nom_min p_nom_max p_min_pu p_max_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
StorageUnit
storageunit non-cyclic 2030 bus 2 PQ 0.0 0.0 False 0.0 inf -1.0 1.0 ... False NaN False True 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 -1.0 1.0 ... False NaN True True 1.0 1.0 1.0 0.0 0.0 0.0

2 rows × 30 columns

Add the load

[9]:
load_var = pd.Series(
    100 * np.random.rand(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)

Run the optimization

[10]:
n.loads_t.p_set
[10]:
Load load 2 load 1
period timestep
2020 2020-01-01 21.358629 75.0
2020-01-02 14.661023 75.0
2020-01-03 83.153735 75.0
2020-01-04 37.419201 75.0
2020-01-05 69.051851 75.0
... ... ... ...
2050 2050-12-27 9.331339 75.0
2050-12-28 68.268808 75.0
2050-12-29 56.862443 75.0
2050-12-30 60.494483 75.0
2050-12-31 39.427332 75.0

1460 rows × 2 columns

[11]:
n.optimize(multi_investment_periods=True)
WARNING:pypsa.components: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='Line')
WARNING:pypsa.components: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='Line')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 27/27 [00:00<00:00, 107.31it/s]
Writing continuous variables.: 100%|██████████| 9/9 [00:00<00:00, 258.38it/s]
INFO:linopy.io: Writing time: 0.3s
INFO:linopy.solvers:Log file at /tmp/highs.log.
Running HiGHS 1.7.0 (git hash: 27ccfaa): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [5e-03, 1e+01]
  Cost   [1e+00, 2e+02]
  Bound  [8e+03, 8e+03]
  RHS    [6e-02, 1e+02]
Presolving model
9855 rows, 6574 cols, 24455 nonzeros  0s
9125 rows, 5844 cols, 22995 nonzeros  0s
Presolve : Reductions: rows 9125(-23366); columns 5844(-6573); elements 22995(-40883)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     8.1156147079e+06 Pr: 1825(156775); Du: 0(4.54198e-10) 0s
       5615     1.8322505640e+07 Pr: 0(0); Du: 0(7.24754e-13) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 5615
Objective value     :  1.8322505640e+07
HiGHS run time      :          0.32
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 12417 primals, 32491 duals
Objective: 1.83e+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.
[11]:
('ok', 'optimal')
[12]:
c = "Generator"
df = pd.concat(
    {
        period: n.get_active_assets(c, period) * n.df(c).p_nom_opt
        for period in n.investment_periods
    },
    axis=1,
)
df.T.plot.bar(
    stacked=True,
    edgecolor="white",
    width=1,
    ylabel="Capacity",
    xlabel="Investment Period",
    rot=0,
    figsize=(10, 5),
)
plt.tight_layout()
../_images/examples_multi-investment-optimisation_17_0.png
[13]:
df = n.generators_t.p.sum(axis=0).T
df.T.plot.bar(
    stacked=True,
    edgecolor="white",
    width=1,
    ylabel="Generation",
    xlabel="Investment Period",
    rot=0,
    figsize=(10, 5),
)
[13]:
<Axes: xlabel='Investment Period', ylabel='Generation'>
../_images/examples_multi-investment-optimisation_18_1.png