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 pypsa
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [1], line 1
----> 1 import pypsa
      2 import pandas as pd
      3 import numpy as np

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/__init__.py:10
      1 # -*- coding: utf-8 -*-
      4 """
      5 Python for Power Systems Analysis (PyPSA)
      6
      7 Grid calculation library.
      8 """
---> 10 from pypsa import (
     11     components,
     12     contingency,
     13     descriptors,
     14     examples,
     15     geo,
     16     io,
     17     linopf,
     18     linopt,
     19     networkclustering,
     20     opf,
     21     opt,
     22     optimization,
     23     pf,
     24     plot,
     25 )
     26 from pypsa.components import Network, SubNetwork
     28 __version__ = "0.21.2"

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/components.py:50
     37 from pypsa.io import (
     38     export_to_csv_folder,
     39     export_to_hdf5,
   (...)
     47     import_series_from_dataframe,
     48 )
     49 from pypsa.opf import network_lopf, network_opf
---> 50 from pypsa.optimization.optimize import OptimizationAccessor
     51 from pypsa.pf import (
     52     calculate_B_H,
     53     calculate_dependent_values,
   (...)
     62     sub_network_pf,
     63 )
     64 from pypsa.plot import iplot, plot

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/optimization/__init__.py:7
      1 #!/usr/bin/env python3
      2 # -*- coding: utf-8 -*-
      3 """
      4 Build optimisation problems from PyPSA networks with Linopy.
      5 """
----> 7 from pypsa.optimization import abstract, constraints, optimize, variables
      8 from pypsa.optimization.optimize import create_model

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/optimization/constraints.py:9
      6 import logging
      8 import pandas as pd
----> 9 from linopy.expressions import LinearExpression, merge
     10 from numpy import arange, cumsum, inf, nan, roll
     11 from scipy import sparse

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/__init__.py:9
      1 #!/usr/bin/env python3
      2 # -*- coding: utf-8 -*-
      3 """
      4 Created on Wed Mar 10 11:03:06 2021.
      5
      6 @author: fabulous
      7 """
----> 9 from linopy import model, remote
     10 from linopy.expressions import merge
     11 from linopy.io import read_netcdf

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/model.py:22
     20 from linopy import solvers
     21 from linopy.common import best_int, replace_by_map
---> 22 from linopy.constraints import (
     23     AnonymousConstraint,
     24     AnonymousScalarConstraint,
     25     Constraints,
     26 )
     27 from linopy.eval import Expr
     28 from linopy.expressions import LinearExpression, ScalarLinearExpression

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/constraints.py:21
     18 from scipy.sparse import coo_matrix
     19 from xarray import DataArray, Dataset
---> 21 from linopy import expressions, variables
     22 from linopy.common import (
     23     _merge_inplace,
     24     has_assigned_model,
   (...)
     27     replace_by_map,
     28 )
     31 class Constraint(DataArray):

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/expressions.py:23
     20 from xarray.core.dataarray import DataArrayCoordinates
     21 from xarray.core.groupby import _maybe_reorder, peek_at
---> 23 from linopy import constraints, variables
     24 from linopy.common import as_dataarray
     27 def exprwrap(method, *default_args, **new_default_kwargs):

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/variables.py:398
    393     roll = varwrap(DataArray.roll)
    395     rolling = varwrap(DataArray.rolling)
--> 398 @dataclass(repr=False)
    399 class Variables:
    400     """
    401     A variables container used for storing multiple variable arrays.
    402     """
    404     labels: Dataset = Dataset()

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:1211, in dataclass.<locals>.wrap(cls)
   1210 def wrap(cls):
-> 1211     return _process_class(cls, init, repr, eq, order, unsafe_hash,
   1212                           frozen, match_args, kw_only, slots,
   1213                           weakref_slot)

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:959, in _process_class(cls, init, repr, eq, order, unsafe_hash, frozen, match_args, kw_only, slots, weakref_slot)
    956         kw_only = True
    957     else:
    958         # Otherwise it's a field of some type.
--> 959         cls_fields.append(_get_field(cls, name, type, kw_only))
    961 for f in cls_fields:
    962     fields[f.name] = f

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:816, in _get_field(cls, a_name, a_type, default_kw_only)
    812 # For real fields, disallow mutable defaults.  Use unhashable as a proxy
    813 # indicator for mutability.  Read the __hash__ attribute from the class,
    814 # not the instance.
    815 if f._field_type is _FIELD and f.default.__class__.__hash__ is None:
--> 816     raise ValueError(f'mutable default {type(f.default)} for field '
    817                      f'{f.name} is not allowed: use default_factory')
    819 return f

ValueError: mutable default <class 'xarray.core.dataset.Dataset'> for field labels is not allowed: use default_factory

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="{}-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
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [2], line 1
----> 1 n = pypsa.Network()
      2 years = [2020, 2030, 2040, 2050]
      3 freq = "24"

NameError: name 'pypsa' is not defined
[3]:
n.investment_periods
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [3], line 1
----> 1 n.investment_periods

NameError: name 'n' is not defined

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
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [4], line 1
----> 1 n.investment_period_weightings["years"] = list(np.diff(years)) + [10]
      3 r = 0.01
      4 T = 0

NameError: name 'np' is not defined

Add the components

[5]:
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
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [5], line 2
      1 for i in range(3):
----> 2     n.add("Bus", "bus {}".format(i))
      4 # add three lines in a ring
      5 n.add(
      6     "Line",
      7     "line 0->1",
      8     bus0="bus 0",
      9     bus1="bus 1",
     10 )

NameError: name 'n' is not defined
[6]:
n.lines
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [6], line 1
----> 1 n.lines

NameError: name 'n' is not defined
[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
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [7], line 2
      1 # add some generators
----> 2 p_nom_max = pd.Series(
      3     (np.random.uniform() for sn in range(len(n.snapshots))),
      4     index=n.snapshots,
      5     name="generator ext 2020",
      6 )
      8 # renewable (can operate 2020, 2030)
      9 n.add(
     10     "Generator",
     11     "generator ext 0 2020",
   (...)
     20     p_nom_extendable=True,
     21 )

NameError: name 'pd' is not defined
[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
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [8], line 1
----> 1 n.add(
      2     "StorageUnit",
      3     "storageunit non-cyclic 2030",
      4     bus="bus 2",
      5     p_nom=0,
      6     capital_cost=2,
      7     build_year=2030,
      8     lifetime=21,
      9     cyclic_state_of_charge=False,
     10     p_nom_extendable=False,
     11 )
     13 n.add(
     14     "StorageUnit",
     15     "storageunit periodic 2020",
   (...)
     23     p_nom_extendable=True,
     24 )
     26 n.storage_units

NameError: name 'n' is not defined

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)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [9], line 1
----> 1 load_var = pd.Series(
      2     100 * np.random.rand(len(n.snapshots)), index=n.snapshots, name="load"
      3 )
      4 n.add("Load", "load 2", bus="bus 2", p_set=load_var)
      6 load_fix = pd.Series(75, index=n.snapshots, name="load")

NameError: name 'pd' is not defined

Run the optimization

[10]:
n.loads_t.p_set
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [10], line 1
----> 1 n.loads_t.p_set

NameError: name 'n' is not defined
[11]:
n.lopf(pyomo=False, multi_investment_periods=True)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [11], line 1
----> 1 n.lopf(pyomo=False, multi_investment_periods=True)

NameError: name 'n' is not defined
[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()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [12], line 2
      1 c = "Generator"
----> 2 df = pd.concat(
      3     {
      4         period: n.get_active_assets(c, period) * n.df(c).p_nom_opt
      5         for period in n.investment_periods
      6     },
      7     axis=1,
      8 )
      9 df.T.plot.bar(
     10     stacked=True,
     11     edgecolor="white",
   (...)
     16     figsize=(10, 5),
     17 )
     18 plt.tight_layout()

NameError: name 'pd' is not defined
[13]:
df = n.generators_t.p.sum(level=0).T
df.T.plot.bar(
    stacked=True,
    edgecolor="white",
    width=1,
    ylabel="Generation",
    xlabel="Investment Period",
    rot=0,
    figsize=(10, 5),
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [13], line 1
----> 1 df = n.generators_t.p.sum(level=0).T
      2 df.T.plot.bar(
      3     stacked=True,
      4     edgecolor="white",
   (...)
      9     figsize=(10, 5),
     10 )

NameError: name 'n' is not defined
[ ]: