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.

import pypsa
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

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.

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

snapshots = pd.DatetimeIndex([])
for year in years:
    period = pd.date_range(
        start="{}-01-01 00:00".format(year),
        freq="{}H".format(freq),
        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_4141/3148991382.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_4141/3148991382.py:7: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.
  period = pd.date_range(

		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

n.investment_periods

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.

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

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

Add the components

for i in range(3):
    n.add("Bus", "bus {}".format(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

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

# 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

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

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

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

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

n.loads_t.p_set

	Load	load 2	load 1
period	timestep
2020	2020-01-01	31.333208	75.0
	2020-01-02	96.900782	75.0
	2020-01-03	60.763826	75.0
	2020-01-04	73.956441	75.0
	2020-01-05	66.200882	75.0
...	...	...	...
2050	2050-12-27	23.425833	75.0
	2050-12-28	16.482154	75.0
	2050-12-29	8.113703	75.0
	2050-12-30	14.347971	75.0
	2050-12-31	82.601711	75.0

1460 rows × 2 columns

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 Glpk solver

INFO:linopy.io:Writing objective.

Writing constraints.: 0%|<span style=”color: rgb(128,191,255)”> </span>| 0/27 [00:00&lt;?, ?it/s]

</pre>

Writing constraints.: 0%|deftcRGB{textcolor[RGB]}expandaftertcRGBexpandafter{detokenize{128,191,255}}{ }| 0/27 [00:00<?, ?it/s]

end{sphinxVerbatim}

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

Writing constraints.: 56%|<span style=”color: rgb(128,191,255)”>█████▌ </span>| 15/27 [00:00&lt;00:00, 144.29it/s]

</pre>

Writing constraints.: 56%|deftcRGB{textcolor[RGB]}expandaftertcRGBexpandafter{detokenize{128,191,255}}{█████▌ }| 15/27 [00:00<00:00, 144.29it/s]

end{sphinxVerbatim}

Writing constraints.: 56%|█████▌ | 15/27 [00:00<00:00, 144.29it/s]

Writing constraints.: 100%|<span style=”color: rgb(128,191,255)”>██████████</span>| 27/27 [00:00&lt;00:00, 118.49it/s]

</pre>

Writing constraints.: 100%|deftcRGB{textcolor[RGB]}expandaftertcRGBexpandafter{detokenize{128,191,255}}{██████████}| 27/27 [00:00<00:00, 118.49it/s]

end{sphinxVerbatim}

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

Writing continuous variables.: 0%|<span style=”color: rgb(128,191,255)”> </span>| 0/9 [00:00&lt;?, ?it/s]

</pre>

Writing continuous variables.: 0%|deftcRGB{textcolor[RGB]}expandaftertcRGBexpandafter{detokenize{128,191,255}}{ }| 0/9 [00:00<?, ?it/s]

end{sphinxVerbatim}

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

Writing continuous variables.: 100%|<span style=”color: rgb(128,191,255)”>██████████</span>| 9/9 [00:00&lt;00:00, 276.97it/s]

</pre>

Writing continuous variables.: 100%|deftcRGB{textcolor[RGB]}expandaftertcRGBexpandafter{detokenize{128,191,255}}{██████████}| 9/9 [00:00<00:00, 276.97it/s]

end{sphinxVerbatim}

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

INFO:linopy.io: Writing time: 0.27s

INFO:linopy.solvers:GLPSOL--GLPK LP/MIP Solver 5.0
Parameter(s) specified in the command line:
 --lp /tmp/linopy-problem-sgf5gp57.lp --output /tmp/linopy-solve-jng2qpoq.sol
Reading problem data from '/tmp/linopy-problem-sgf5gp57.lp'...
32491 rows, 12417 columns, 63878 non-zeros
175587 lines were read
GLPK Simplex Optimizer 5.0
32491 rows, 12417 columns, 63878 non-zeros
Preprocessing...
18250 rows, 8401 columns, 42705 non-zeros
Scaling...
 A: min|aij| =  2.272e-05  max|aij| =  1.000e+01  ratio =  4.401e+05
GM: min|aij| =  3.492e-01  max|aij| =  2.864e+00  ratio =  8.200e+00
EQ: min|aij| =  1.232e-01  max|aij| =  1.000e+00  ratio =  8.114e+00
Constructing initial basis...
Size of triangular part is 17882
      0: obj =   1.234429256e+07 inf =   2.454e+07 (5475)
   5478: obj =   1.834976421e+07 inf =   1.283e-12 (0) 52
*  6100: obj =   1.792097357e+07 inf =   8.171e-13 (0) 5
OPTIMAL LP SOLUTION FOUND
Time used:   3.3 secs
Memory used: 26.3 Mb (27549469 bytes)
Writing basic solution to '/tmp/linopy-solve-jng2qpoq.sol'...

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 12417 primals, 32491 duals
Objective: 1.79e+07
Solver model: not available
Solver message: optimal

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/pypsa/optimization/optimize.py:357: FutureWarning: A value is trying to be set on a copy of a DataFrame or Series through chained assignment using an inplace method.
The behavior will change in pandas 3.0. This inplace method will never work because the intermediate object on which we are setting values always behaves as a copy.

For example, when doing 'df[col].method(value, inplace=True)', try using 'df.method({col: value}, inplace=True)' or df[col] = df[col].method(value) instead, to perform the operation inplace on the original object.


  n.df(c)[attr + "_opt"].update(df)
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.

('ok', 'optimal')

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

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

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