Note

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

Single Node Capacity Expansion Planning#

From electricity market modelling to capacity expansion planning#

Review the problem formulation of the electricity market model. Below you can find an adapted version where the capacity limits have been promoted to decision variables with corresponding terms in the objective function and new constraints for their expansion limits (e.g. wind and solar potentials). This is known as capacity expansion problem.

\[\begin{equation*} \min_{g,e,f,G,E,F} \quad \sum_{i,s,t} w_t o_{s} g_{i,s,t} + \sum_{i,s} c_sG_{i,s} + c_{r,\text{dis/charge}}G_{i,r, \text{dis/charge}} + c_{r}E_{i,r} + c_\ell F_{\ell} \end{equation*}\]

such that

\[\begin{split}\begin{align*} d_{i,t} &= \sum_s g_{i,s,t} - \sum_\ell K_{i\ell} f_{\ell,t} & \text{energy balance} \\ 0 &\leq g_{i,s,t} \leq \hat{g}_{i,s,t} G_{i,s} & \text{generator limits}\\ 0 & \leq g_{i,r,t,\text{dis/charge}} \leq G_{i,r,\text{dis/charge}}& \text{storage dis/charge limits} \\ 0 & \leq e_{i,r,t} \leq E_{r} & \text{storage energy limits} \\ e_{i,r,t} &= \eta^0_{i,r,t} e_{i,r,t-1} + \eta^1_{r}g_{i,r,t,\text{charge}} - \frac{1}{\eta^2_{r}} g_{i,r,t,\text{discharge}} & \text{storage consistency} \\ -F_\ell &\leq f_{\ell,t} \leq F_\ell & \text{line limits} \\ 0 &= \sum_\ell C_{\ell c} x_\ell f_{\ell,t} & \text{KVL} \\ \underline{G}_{i,s} & \leq G_{i,s} \leq \overline{G}_{i,s} & \text{generator capacity expansion limits} \\ \underline{G}_{i,r, \text{dis/charge}} & \leq G_{i,r, \text{dis/charge}} \leq \overline{G}_{i,r, \text{dis/charge}} & \text{storage power capacity expansion limits} \\ \underline{E}_{i,r} & \leq E_{i,r} \leq \overline{E}_{i,r} & \text{storage energy expansion limits} \\ \underline{F}_{\ell} & \leq F_{\ell} \leq \overline{F}_{\ell} & \text{line capacity expansion limits} \end{align*}\end{split}\]

New decision variables for capacity expansion planning:

\(G_{i,s}\) is the generator capacity at bus \(i\), technology \(s\),
\(F_{\ell}\) is the transmission capacity of line \(\ell\),
\(G_{i,r,\text{dis-/charge}}\) denotes the charge and discharge capacities of storage unit \(r\) at bus \(i\),
\(E_{i,r}\) is the energy capacity of storage \(r\) at bus \(i\) and time step \(t\).

New parameters for capacity expansion planning:

\(c_{\star}\) is the capital cost of technology \(\star\) at bus \(i\)
\(w_t\) is the weighting of time step \(t\) (e.g. number of hours it represents)
\(\underline{G}_\star, \underline{F}_\star, \underline{E}_\star\) are the minimum capacities per technology and location/connection.
\(\underline{G}_\star, \underline{F}_\star, \underline{E}_\star\) are the maximum capacities per technology and location.

For a full reference to the optimisation problem description, see https://pypsa.readthedocs.io/en/latest/optimal_power_flow.html

First things first! We need a few packages for this tutorial:

[1]:

import matplotlib.pyplot as plt
import pandas as pd

import pypsa

plt.style.use("bmh")

Prerequisites: handling technology data and costs#

We maintain a database (PyPSA/technology-data) which collects assumptions and projections for energy system technologies (such as costs, efficiencies, lifetimes, etc.) for given years, which we can load into a pandas.DataFrame. This requires some pre-processing to load (e.g. converting units, setting defaults, re-arranging dimensions):

[2]:

year = 2030
url = f"https://raw.githubusercontent.com/PyPSA/technology-data/master/outputs/costs_{year}.csv"
costs = pd.read_csv(url, index_col=[0, 1])

[3]:

costs.loc[costs.unit.str.contains("/kW"), "value"] *= 1e3
costs.unit = costs.unit.str.replace("/kW", "/MW")

defaults = {
    "FOM": 0,
    "VOM": 0,
    "efficiency": 1,
    "fuel": 0,
    "investment": 0,
    "lifetime": 25,
    "CO2 intensity": 0,
    "discount rate": 0.07,
}
costs = costs.value.unstack().fillna(defaults)

costs.at["OCGT", "fuel"] = costs.at["gas", "fuel"]
costs.at["CCGT", "fuel"] = costs.at["gas", "fuel"]
costs.at["OCGT", "CO2 intensity"] = costs.at["gas", "CO2 intensity"]
costs.at["CCGT", "CO2 intensity"] = costs.at["gas", "CO2 intensity"]

Let’s also write a small utility function that calculates the annuity to annualise investment costs. The formula is

\[a(r, n) = \frac{r}{1-(1+r)^{-n}}\]

where \(r\) is the discount rate and \(n\) is the lifetime.

[4]:

def annuity(r, n):
    return r / (1.0 - 1.0 / (1.0 + r) ** n)

[5]:

annuity(0.07, 20)

[5]:

0.09439292574325567

Based on this, we can calculate the short-term marginal generation costs (STMGC, €/MWh):

[6]:

costs["marginal_cost"] = costs["VOM"] + costs["fuel"] / costs["efficiency"]

and the annualised investment costs (capital_cost in PyPSA terms, €/MW/a):

[7]:

annuity = costs.apply(lambda x: annuity(x["discount rate"], x["lifetime"]), axis=1)

[8]:

costs["capital_cost"] = (annuity + costs["FOM"] / 100) * costs["investment"]

Loading time series data#

We are also going to need some time series for wind, solar and load. For now, we are going to recycle the time series we used at the beginning of the course. They are given for Germany in the year 2015.

[9]:

url = (
    "https://tubcloud.tu-berlin.de/s/pKttFadrbTKSJKF/download/time-series-lecture-2.csv"
)
ts = pd.read_csv(url, index_col=0, parse_dates=True)

[10]:

ts.head(3)

[10]:

	load	onwind	offwind	prices
2015-01-01 00:00:00	41.151	0.1566	0.7030	NaN
2015-01-01 01:00:00	40.135	0.1659	0.6875	NaN
2015-01-01 02:00:00	39.106	0.1746	0.6535	NaN

Let’s convert the load time series from GW to MW, the base unit of PyPSA:

[11]:

ts.load *= 1e3

We are also going to adapt the temporal resolution of the time series, e.g. sample only every other hour, to save some time:

[12]:

resolution = 4
ts = ts.resample(f"{resolution}h").first()

Simple capacity expansion planning example#

Model Initialisation#

For building the model, we start again by initialising an empty network.

[13]:

n = pypsa.Network()

Then, we add a single bus…

[14]:

n.add("Bus", "electricity")

[14]:

Index(['electricity'], dtype='object')

…and tell the pypsa.Network object n what the snapshots of the model will be using the utility function n.set_snapshots().

[15]:

n.set_snapshots(ts.index)

[16]:

n.snapshots

[16]:

DatetimeIndex(['2015-01-01 00:00:00', '2015-01-01 04:00:00',
               '2015-01-01 08:00:00', '2015-01-01 12:00:00',
               '2015-01-01 16:00:00', '2015-01-01 20:00:00',
               '2015-01-02 00:00:00', '2015-01-02 04:00:00',
               '2015-01-02 08:00:00', '2015-01-02 12:00:00',
               ...
               '2015-12-30 08:00:00', '2015-12-30 12:00:00',
               '2015-12-30 16:00:00', '2015-12-30 20:00:00',
               '2015-12-31 00:00:00', '2015-12-31 04:00:00',
               '2015-12-31 08:00:00', '2015-12-31 12:00:00',
               '2015-12-31 16:00:00', '2015-12-31 20:00:00'],
              dtype='datetime64[ns]', name='snapshot', length=2190, freq='4h')

The weighting of the snapshots (e.g. how many hours they represent, see \(w_t\) in problem formulation above) can be set in n.snapshot_weightings.

[17]:

n.snapshot_weightings.head(3)

[17]:

	objective	stores	generators
snapshot
2015-01-01 00:00:00	1.0	1.0	1.0
2015-01-01 04:00:00	1.0	1.0	1.0
2015-01-01 08:00:00	1.0	1.0	1.0

[18]:

n.snapshot_weightings.loc[:, :] = resolution

[19]:

n.snapshot_weightings.head(3)

[19]:

	objective	stores	generators
snapshot
2015-01-01 00:00:00	4.0	4.0	4.0
2015-01-01 04:00:00	4.0	4.0	4.0
2015-01-01 08:00:00	4.0	4.0	4.0

Adding Components#

Then, we add all the technologies we are going to include as carriers.

[20]:

carriers = [
    "onwind",
    "offwind",
    "solar",
    "OCGT",
    "hydrogen storage underground",
    "battery storage",
]

n.add(
    "Carrier",
    carriers,
    color=["dodgerblue", "aquamarine", "gold", "indianred", "magenta", "yellowgreen"],
    co2_emissions=[costs.at[c, "CO2 intensity"] for c in carriers],
)

[20]:

Index(['onwind', 'offwind', 'solar', 'OCGT', 'hydrogen storage underground',
       'battery storage'],
      dtype='object')

Next, we add the demand time series to the model.

[21]:

n.add(
    "Load",
    "demand",
    bus="electricity",
    p_set=ts.load,
)

[21]:

Index(['demand'], dtype='object')

Let’s have a check whether the data was read-in correctly.

[22]:

n.loads_t.p_set.plot(figsize=(6, 2), ylabel="MW")

[22]:

<Axes: xlabel='snapshot', ylabel='MW'>

../_images/examples_capacity-expansion-planning-single-node_45_1.png

We are going to add one dispatchable generation technology to the model. This is an open-cycle gas turbine (OCGT) with CO\(_2\) emissions of 0.2 t/MWh\(_{th}\).

[23]:

n.add(
    "Generator",
    "OCGT",
    bus="electricity",
    carrier="OCGT",
    capital_cost=costs.at["OCGT", "capital_cost"],
    marginal_cost=costs.at["OCGT", "marginal_cost"],
    efficiency=costs.at["OCGT", "efficiency"],
    p_nom_extendable=True,
)

[23]:

Index(['OCGT'], dtype='object')

Adding the variable renewable generators works almost identically, but we also need to supply the capacity factors to the model via the attribute p_max_pu.

[24]:

for tech in ["onwind", "offwind", "solar"]:
    n.add(
        "Generator",
        tech,
        bus="electricity",
        carrier=tech,
        p_max_pu=ts[tech],
        capital_cost=costs.at[tech, "capital_cost"],
        marginal_cost=costs.at[tech, "marginal_cost"],
        efficiency=costs.at[tech, "efficiency"],
        p_nom_extendable=True,
    )

So let’s make sure the capacity factors are read-in correctly.

[25]:

n.generators_t.p_max_pu.loc["2015-03"].plot(figsize=(6, 2), ylabel="CF")

[25]:

<Axes: xlabel='snapshot', ylabel='CF'>

../_images/examples_capacity-expansion-planning-single-node_51_1.png

Model Run#

Then, we can already solve the model for the first time. At this stage, the model does not have any storage or emission limits implemented. It’s going to look for the least-cost combination of variable renewables and the gas turbine to supply demand.

[26]:

n.optimize(solver_name="highs")

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.12/site-packages/linopy/common.py:147: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io: Writing time: 0.19s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 4e+00]
  Cost   [4e-02, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [3e+04, 8e+04]
Presolving model
9900 rows, 7714 cols, 23130 nonzeros  0s
9900 rows, 7714 cols, 23130 nonzeros  0s
Presolve : Reductions: rows 9900(-9818); columns 7714(-1050); elements 23130(-19624)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2190(1.25606e+09) 0s
       7615     3.1624531941e+10 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 7615
Objective value     :  3.1624531941e+10
HiGHS run time      :          0.20
Writing the solution to /tmp/linopy-solve-9tiyatf6.sol

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 8764 primals, 19718 duals
Objective: 3.16e+10
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper were not assigned to the network.

[26]:

('ok', 'optimal')

Model Evaluation#

The total system cost in billion Euros per year:

[27]:

n.objective / 1e9

[27]:

31.62453194099228

The optimised capacities in GW:

[28]:

n.generators.p_nom_opt.div(1e3)  # GW

[28]:

Generator
OCGT       70.096357
onwind     -0.000000
offwind    49.326062
solar      85.967820
Name: p_nom_opt, dtype: float64

The total energy generation by technology in GW:

[29]:

n.snapshot_weightings.generators @ n.generators_t.p.div(1e6)  # TWh

[29]:

Generator
OCGT       235.778879
onwind       0.000000
offwind    150.379873
solar       92.766988
Name: generators, dtype: float64

While we get the objective value through n.objective, in many cases we want to know how the costs are distributed across the technologies. We can use the statistics module for this:

[30]:

(n.statistics.capex() + n.statistics.opex()).div(1e6)

[30]:

component  carrier
Generator  OCGT       18596.014468
           offwind     8613.359135
           solar       4415.158338
dtype: float64

Possibly, we are also interested in the total emissions:

[31]:

emissions = (
    n.generators_t.p
    / n.generators.efficiency
    * n.generators.carrier.map(n.carriers.co2_emissions)
)  # t/h

[32]:

n.snapshot_weightings.generators @ emissions.sum(axis=1).div(1e6)  # Mt

[32]:

113.86394645337052

Plotting Optimal Dispatch#

This function takes the network object n as an argument and, optionally, a time frame. We want to plot the load time series, and stacked area charts for electricity feed-in and storage charging. Technologies should be coloured by their color defined in n.carriers.

[33]:

def plot_dispatch(n, time="2015-07"):
    p_by_carrier = n.generators_t.p.groupby(n.generators.carrier, axis=1).sum().div(1e3)

    if not n.storage_units.empty:
        sto = n.storage_units_t.p.T.groupby(n.storage_units.carrier).sum().T.div(1e3)
        p_by_carrier = pd.concat([p_by_carrier, sto], axis=1)

    fig, ax = plt.subplots(figsize=(6, 3))

    color = p_by_carrier.columns.map(n.carriers.color)

    p_by_carrier.where(p_by_carrier > 0).loc[time].plot.area(
        ax=ax,
        linewidth=0,
        color=color,
    )

    charge = p_by_carrier.where(p_by_carrier < 0).dropna(how="all", axis=1).loc[time]

    if not charge.empty:
        charge.plot.area(
            ax=ax,
            linewidth=0,
            color=charge.columns.map(n.carriers.color),
        )

    n.loads_t.p_set.sum(axis=1).loc[time].div(1e3).plot(ax=ax, c="k")

    plt.legend(loc=(1.05, 0))
    ax.set_ylabel("GW")
    ax.set_ylim(-200, 200)

Oje, that was complicated. Let’s test it:

[34]:

plot_dispatch(n)

/tmp/ipykernel_1233/184362636.py:2: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  p_by_carrier = n.generators_t.p.groupby(n.generators.carrier, axis=1).sum().div(1e3)

../_images/examples_capacity-expansion-planning-single-node_71_1.png

Sensitivity Analysis#

Sensitivity analyses constitute a core activity of energy system modelling. Below, you can find sensitivity analyses regarding the

variation in allowed CO\(_2\) emissions
variation in solar overnight costs
variation in offshore wind potentials

[52]:

sensitivity = {}
for co2 in [150, 100, 50, 25, 0]:
    n.global_constraints.loc["CO2Limit", "constant"] = co2 * 1e6
    n.optimize(solver_name="highs")
    sensitivity[co2] = system_cost(n)

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.12/site-packages/linopy/common.py:147: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 17/17 [00:00<00:00, 45.72it/s]
Writing continuous variables.: 100%|██████████| 6/6 [00:00<00:00, 103.65it/s]
INFO:linopy.io: Writing time: 0.46s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 21906 primals, 50381 duals
Objective: 3.16e+10
Solver model: available
Solver message: optimal

  Matrix [1e-04, 2e+02]
  Cost   [4e-02, 4e+05]
  Bound  [0e+00, 0e+00]
  RHS    [3e+04, 2e+08]
Presolving model
27421 rows, 20856 cols, 77880 nonzeros  0s
27421 rows, 20856 cols, 77880 nonzeros  0s
Presolve : Reductions: rows 27421(-22960); columns 20856(-1050); elements 77880(-32766)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2190(1.67812e+09) 0s
      24532     3.1624531941e+10 Pr: 0(0); Du: 0(1.65201e-13) 5s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 24532
Objective value     :  3.1624531941e+10
HiGHS run time      :          4.63
Writing the solution to /tmp/linopy-solve-85yf53_y.sol

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.12/site-packages/linopy/common.py:147: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 17/17 [00:00<00:00, 45.70it/s]
Writing continuous variables.: 100%|██████████| 6/6 [00:00<00:00, 104.47it/s]
INFO:linopy.io: Writing time: 0.45s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 2e+02]
  Cost   [4e-02, 4e+05]
  Bound  [0e+00, 0e+00]
  RHS    [3e+04, 1e+08]
Presolving model
27421 rows, 20856 cols, 77880 nonzeros  0s
27421 rows, 20856 cols, 77880 nonzeros  0s
Presolve : Reductions: rows 27421(-22960); columns 20856(-1050); elements 77880(-32766)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2190(1.67812e+09) 0s
      19041     2.9245951705e+10 Pr: 3073(7.1811e+10); Du: 0(4.5186e-07) 5s
      25751     3.1764091630e+10 Pr: 0(0); Du: 0(3.11244e-11) 10s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 25751
Objective value     :  3.1764091630e+10
HiGHS run time      :          9.82
Writing the solution to /tmp/linopy-solve-8zp1eihh.sol

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 21906 primals, 50381 duals
Objective: 3.18e+10
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.12/site-packages/linopy/common.py:147: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 17/17 [00:00<00:00, 45.82it/s]
Writing continuous variables.: 100%|██████████| 6/6 [00:00<00:00, 103.75it/s]
INFO:linopy.io: Writing time: 0.45s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 21906 primals, 50381 duals
Objective: 3.58e+10
Solver model: available
Solver message: optimal

  Matrix [1e-04, 2e+02]
  Cost   [4e-02, 4e+05]
  Bound  [0e+00, 0e+00]
  RHS    [3e+04, 5e+07]
Presolving model
27421 rows, 20856 cols, 77880 nonzeros  0s
27421 rows, 20856 cols, 77880 nonzeros  0s
Presolve : Reductions: rows 27421(-22960); columns 20856(-1050); elements 77880(-32766)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2190(1.67812e+09) 0s
      15913     2.8565817800e+10 Pr: 10548(1.58651e+11); Du: 0(2.02512e-07) 5s
      25521     3.4479528678e+10 Pr: 1976(9.63063e+09); Du: 0(7.85692e-07) 10s
      30327     3.5768174048e+10 Pr: 0(0) 13s
      30327     3.5768174048e+10 Pr: 0(0) 13s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 30327
Objective value     :  3.5768174048e+10
HiGHS run time      :         13.42
Writing the solution to /tmp/linopy-solve-9_j7zqlu.sol

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.12/site-packages/linopy/common.py:147: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 17/17 [00:00<00:00, 45.79it/s]
Writing continuous variables.: 100%|██████████| 6/6 [00:00<00:00, 104.88it/s]
INFO:linopy.io: Writing time: 0.45s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 2e+02]
  Cost   [4e-02, 4e+05]
  Bound  [0e+00, 0e+00]
  RHS    [3e+04, 2e+07]
Presolving model
27421 rows, 20856 cols, 77880 nonzeros  0s
27421 rows, 20856 cols, 77880 nonzeros  0s
Presolve : Reductions: rows 27421(-22960); columns 20856(-1050); elements 77880(-32766)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2190(1.67812e+09) 0s
      11470     2.3500856251e+10 Pr: 11733(4.88208e+12); Du: 0(5.60196e-08) 5s
      17749     3.7543369478e+10 Pr: 11168(1.06863e+12); Du: 0(4.29846e-08) 10s
      22314     4.1469352890e+10 Pr: 0(0); Du: 0(8.992e-11) 14s
      22314     4.1469352890e+10 Pr: 0(0); Du: 0(8.992e-11) 14s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 22314
Objective value     :  4.1469352890e+10
HiGHS run time      :         14.08
Writing the solution to /tmp/linopy-solve-amu9x36r.sol

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 21906 primals, 50381 duals
Objective: 4.15e+10
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.12/site-packages/linopy/common.py:147: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 17/17 [00:00<00:00, 46.17it/s]
Writing continuous variables.: 100%|██████████| 6/6 [00:00<00:00, 103.41it/s]
INFO:linopy.io: Writing time: 0.45s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 2e+02]
  Cost   [4e-02, 4e+05]
  Bound  [0e+00, 0e+00]
  RHS    [3e+04, 8e+04]
Presolving model
25230 rows, 18665 cols, 69120 nonzeros  0s
25230 rows, 18665 cols, 69120 nonzeros  0s
Presolve : Reductions: rows 25230(-25151); columns 18665(-3241); elements 69120(-41526)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2190(1.52524e+09) 0s
      11432     2.6428138659e+10 Pr: 4637(1.5953e+13); Du: 0(2.8779e-08) 5s
      21606     7.9405597314e+10 Pr: 104(5.13604e+08); Du: 0(1.21127e-07) 11s
      21723     7.9406619875e+10 Pr: 0(0); Du: 0(1.42109e-14) 11s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 21723
Objective value     :  7.9406619875e+10
HiGHS run time      :         10.74
Writing the solution to /tmp/linopy-solve-v3qgdpg3.sol

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 21906 primals, 50381 duals
Objective: 7.94e+10
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.

[53]:

df = pd.DataFrame(sensitivity).T.div(1e3)  # billion Euro/a
df.plot.area(
    stacked=True,
    linewidth=0,
    color=df.columns.map(n.carriers.color),
    figsize=(4, 4),
    xlim=(0, 150),
    xlabel=r"CO$_2$ emissions [Mt/a]",
    ylabel="System cost [bn€/a]",
    ylim=(0, 100),
)
plt.legend(frameon=False, loc=(1.05, 0))

[53]:

<matplotlib.legend.Legend at 0x7fc04c6d5400>

../_images/examples_capacity-expansion-planning-single-node_107_1.png

[ ]:

Single Node Capacity Expansion Planning

Contents

Single Node Capacity Expansion Planning#

From electricity market modelling to capacity expansion planning#

Prerequisites: handling technology data and costs#

Loading time series data#

Simple capacity expansion planning example#

Model Initialisation#

Adding Components#

Model Run#

Model Evaluation#

Plotting Optimal Dispatch#

Adding Storage Units#

Adding emission limits#

Sensitivity Analysis#