Capacity Expansion Planning

Capacity Expansion Planning#

Note

If you have not yet set up Python on your computer, you can execute this tutorial in your browser via Google Colab. Click on the rocket in the top right corner and launch “Colab”. If that doesn’t work download the .ipynb file and import it in Google Colab.

Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.

!pip install pypsa pandas highspy "plotly<6"

In this tutorial, we want to build a replica of model.energy. This tool calculates the cost of meeting a constant electricity demand from a combination of wind power, solar power and storage for different regions of the world. We deviate from model.energy by including electricity demand profiles rather than a constant electricity demand.

Note

From electricity market modelling to capacity expansion planning#

Review the problem formulation of the electricity market model from the previous tutorial. 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{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*}\]

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.

Note

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

Package Imports#

import pypsa
import pandas as pd
import plotly.express as px
import plotly.io as pio
import plotly.offline as py
from plotly.subplots import make_subplots
pd.options.plotting.backend = "plotly"

Techology Data and Costs#

At TU Berlin, 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 use for our research.

Reading this data into a useable pandas.DataFrame requires some pre-processing (e.g. converting units, setting defaults, re-arranging dimensions):

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

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.

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

annuity(0.07, 20)

0.09439292574325567

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

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

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

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

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

Capacity Factor and Load Time Series#

We are also going to need some time series for wind, solar and load.

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

ts.head(3)

	load_mw	pv_pu	wind_pu
timestamp
2019-01-01 00:00:00	5719.26	0.0	0.1846
2019-01-01 01:00:00	5677.73	0.0	0.2293
2019-01-01 02:00:00	5622.20	0.0	0.2718

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

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

	load_mw	pv_pu	wind_pu
timestamp
2019-01-01 00:00:00	5719.26	0.000	0.1846
2019-01-01 03:00:00	5474.74	0.000	0.3146
2019-01-01 06:00:00	5413.39	0.000	0.4957
2019-01-01 09:00:00	5891.23	0.059	0.6824
2019-01-01 12:00:00	6662.42	0.063	0.7519
...	...	...	...
2019-12-31 09:00:00	7029.73	0.195	0.3049
2019-12-31 12:00:00	7342.83	0.135	0.3944
2019-12-31 15:00:00	7267.22	0.000	0.3788
2019-12-31 18:00:00	6821.12	0.000	0.3157
2019-12-31 21:00:00	6191.86	0.000	0.2636

2920 rows × 3 columns

Building the Model#

Model Initialisation#

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

n = pypsa.Network()

Then, we add a single bus…

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

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

n.set_snapshots(ts.index)

n.snapshots[:5]

DatetimeIndex(['2019-01-01 00:00:00', '2019-01-01 03:00:00',
               '2019-01-01 06:00:00', '2019-01-01 09:00:00',
               '2019-01-01 12:00:00'],
              dtype='datetime64[ns]', name='snapshot', freq='3h')

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.

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

n.snapshot_weightings.head(3)

	objective	stores	generators
snapshot
2019-01-01 00:00:00	3.0	3.0	3.0
2019-01-01 03:00:00	3.0	3.0	3.0
2019-01-01 06:00:00	3.0	3.0	3.0

Adding Components#

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

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

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

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

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

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

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

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

n.loads_t.p_set.plot()

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}\).

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

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.

n.add(
    "Generator",
    "wind",
    bus="electricity",
    carrier="wind",
    p_max_pu=ts.wind_pu,
    capital_cost=costs.at["onwind", "capital_cost"],
    marginal_cost=costs.at["onwind", "marginal_cost"],
    p_nom_extendable=True,
)

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

n.add(
    "Generator",
    "solar",
    bus="electricity",
    carrier="solar",
    p_max_pu=ts.pv_pu,
    capital_cost=costs.at["solar", "capital_cost"],
    marginal_cost=costs.at["solar", "marginal_cost"],
    p_nom_extendable=True,
)

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

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

n.generators_t.p_max_pu.loc["2019-03"].plot()

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.

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 generators have carriers which are not defined:
Index(['wind'], dtype='object', name='Generator')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io: Writing time: 0.13s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 8763 primals, 20443 duals
Objective: 4.31e+09
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.

Running HiGHS 1.9.0 (git hash: fa40bdf): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [2e-04, 1e+00]
  Cost   [3e-02, 1e+05]
  Bound  [0e+00, 0e+00]
  RHS    [5e+03, 1e+04]
Presolving model
8836 rows, 5919 cols, 19170 nonzeros  0s
8836 rows, 5919 cols, 19170 nonzeros  0s
Presolve : Reductions: rows 8836(-11607); columns 5919(-2844); elements 19170(-14451)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Ph1: 0(0) 0s
       5937     4.3107252364e+09 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-fv1154d3
Model status        : Optimal
Simplex   iterations: 5937
Objective value     :  4.3107252364e+09
Relative P-D gap    :  1.7698633322e-15
HiGHS run time      :          0.12
Writing the solution to /tmp/linopy-solve-f3puia0t.sol

('ok', 'optimal')

Model Evaluation#

The total system cost in billion Euros per year:

n.objective / 1e9

4.310725236366917

The optimised capacities in GW:

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

Generator
OCGT     10.118517
wind     10.692584
solar    11.493762
Name: p_nom_opt, dtype: float64

The energy balance by component in TWh:

n.statistics.energy_balance().sort_values().div(1e6)  # TWh

component  carrier  bus_carrier
Load       -        electricity   -66.266089
Generator  solar    electricity    13.517459
           wind     electricity    19.939127
           OCGT     electricity    32.809502
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:

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

component  carrier
Generator  OCGT       2605.090402
           solar       590.311502
           wind       1115.323332
dtype: float64

Possibly, we are also interested in the total emissions:

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

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

15.844588885480375

Plotting Optimal Dispatch#

We want to plot the supply and withdrawal as a stacked area chart for electricity feed-in and storage charging.

def plot_dispatch(n):
    p = (
        n.statistics.energy_balance(aggregate_time=False)
        .groupby("carrier")
        .sum()
        .div(1e3)
        .T
    )

    supply = (
        p.where(p > 0, 0)
        .stack()
        .reset_index()
        .rename(columns={0: "GW"})
    )

    withdrawal = (
        p.where(p < 0, 0)
        .stack()
        .reset_index()
        .rename(columns={0: "GW"})
    )

    fig = make_subplots(rows=2, cols=1, shared_xaxes=True, vertical_spacing=0)

    for data, row, yaxis_title in [
        (supply, 1, "Supply (GW)"),
        (withdrawal, 2, "Consumption (GW)"),
    ]:
        fig_data = px.area(
            data,
            x="snapshot",
            color="carrier",
            y="GW",
            line_group="carrier",
            height=400,
        )["data"]
        for trace in fig_data:
            trace.update(line=dict(width=0))
            fig.add_trace(trace, row=row, col=1)
        fig.update_yaxes(title_text=yaxis_title, row=row, col=1)

    return fig

Let’s test it:

plot_dispatch(n)

Adding Storage Units#

Alright, but there are a few important components missing for a system with high shares of renewables? What about short-term storage options (e.g. batteries) and long-term storage options (e.g. hydrogen storage)? Let’s add them, too.

First, the battery storage. We are going to assume a fixed energy-to-power ratio of 4 hours, i.e. if fully charged, the battery can discharge at full capacity for 4 hours.

For the capital cost, we have to factor in both the capacity and energy cost of the storage. We are also going to enforce a cyclic state-of-charge condition, i.e. the state of charge at the beginning of the optimisation period must equal the final state of charge.

n.add(
    "StorageUnit",
    "battery storage",
    bus="electricity",
    carrier="battery storage",
    max_hours=4,
    capital_cost=costs.at["battery inverter", "capital_cost"]
    + 4 * costs.at["battery storage", "capital_cost"],
    efficiency_store=costs.at["battery inverter", "efficiency"],
    efficiency_dispatch=costs.at["battery inverter", "efficiency"],
    p_nom_extendable=True,
    cyclic_state_of_charge=True,
)

Index(['battery storage'], dtype='object')

Second, the hydrogen storage. This one is composed of an electrolysis to convert electricity to hydrogen, a fuel cell to re-convert hydrogen to electricity and underground storage (e.g. in salt caverns). We assume an energy-to-power ratio of 336 hours, such that this type of storage can be used for weekly balancing.

capital_costs = (
    costs.at["electrolysis", "capital_cost"]
    + costs.at["fuel cell", "capital_cost"]
    + 336 * costs.at["hydrogen storage underground", "capital_cost"]
)

n.add(
    "StorageUnit",
    "hydrogen storage underground",
    bus="electricity",
    carrier="hydrogen storage underground",
    max_hours=336,
    capital_cost=capital_costs,
    efficiency_store=costs.at["electrolysis", "efficiency"],
    efficiency_dispatch=costs.at["fuel cell", "efficiency"],
    p_nom_extendable=True,
    cyclic_state_of_charge=True,
)

Index(['hydrogen storage underground'], dtype='object')

Ok, lets run the again, now with storage, and see what’s changed.

n.optimize(solver_name="highs")

Show code cell output Hide code cell output

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')

WARNING:pypsa.consistency:The following generators have carriers which are not defined:
Index(['wind'], dtype='object', name='Generator')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io:Writing objective.

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

Writing constraints.:  50%|█████     | 7/14 [00:00<00:00, 67.50it/s]

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

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

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

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

INFO:linopy.io: Writing time: 0.42s

Running HiGHS 1.9.0 (git hash: fa40bdf): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [2e-04, 3e+02]
  Cost   [3e-02, 5e+05]
  Bound  [0e+00, 0e+00]
  RHS    [5e+03, 1e+04]
Presolving model
33618 rows, 24863 cols, 92094 nonzeros  0s
33618 rows, 24863 cols, 92094 nonzeros  0s
Presolve : Reductions: rows 33618(-27707); columns 24863(-1422); elements 92094(-29129)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2920(1.5003e+09) 0s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 26285 primals, 61325 duals
Objective: 4.31e+09
Solver model: available
Solver message: optimal

      37407     4.3098930942e+09 Pr: 0(0); Du: 0(2.25719e-13) 4s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-7rrab0ou
Model status        : Optimal
Simplex   iterations: 37407
Objective value     :  4.3098930942e+09
Relative P-D gap    :  1.3276537894e-15
HiGHS run time      :          3.93
Writing the solution to /tmp/linopy-solve-4_oc23su.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.

('ok', 'optimal')

n.statistics.optimal_capacity().div(1e3)  # GW

component    carrier        
StorageUnit  battery storage     0.214420
Generator    OCGT                9.902907
             solar              11.699476
             wind               10.861762
dtype: float64

n.statistics.energy_balance().sort_values().div(1e6)  # TWh

component    carrier          bus_carrier
Load         -                electricity   -66.266089
StorageUnit  battery storage  electricity    -0.009501
Generator    solar            electricity    13.763493
             wind             electricity    20.245123
             OCGT             electricity    32.266974
dtype: float64

pd.concat({
    "capex": n.statistics.capex(),
    "opex": n.statistics.opex(),
}, axis=1).div(1e9).round(2) # bn€/a

		capex	opex
component	carrier
StorageUnit	battery storage	0.02	NaN
Generator	OCGT	0.47	2.09
	solar	0.60	0.00
	wind	1.10	0.03

n.buses_t.marginal_price.sort_values(by="electricity", ascending=False).reset_index(
    drop=True
).plot(title="price duration curve [€/MWh]")

Adding emission limits#

Now, let’s model a 100% renewable electricity system by adding a CO\(_2\) emission limit as global constraint:

n.add(
    "GlobalConstraint",
    "CO2Limit",
    carrier_attribute="co2_emissions",
    sense="<=",
    constant=0,
)

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

When we run the model now…

n.optimize(solver_name="highs")

Show code cell output Hide code cell output

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')

WARNING:pypsa.consistency:The following generators have carriers which are not defined:
Index(['wind'], dtype='object', name='Generator')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io:Writing objective.

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

Writing constraints.:  53%|█████▎    | 8/15 [00:00<00:00, 62.51it/s]

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

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

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

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

INFO:linopy.io: Writing time: 0.4s

Running HiGHS 1.9.0 (git hash: fa40bdf): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [2e-04, 3e+02]
  Cost   [3e-02, 5e+05]
  Bound  [0e+00, 0e+00]
  RHS    [5e+03, 1e+04]
Presolving model
30698 rows, 21942 cols, 83334 nonzeros  0s
30698 rows, 21942 cols, 83334 nonzeros  0s
Presolve : Reductions: rows 30698(-30628); columns 21942(-4343); elements 83334(-40809)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2920(1.44316e+09) 0s

      13695     4.1437919525e+09 Pr: 8410(1.18244e+12) 5s

      20302     8.0566543065e+09 Pr: 1130(1.36011e+09); Du: 0(6.55046e-07) 10s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 26285 primals, 61326 duals
Objective: 8.10e+09
Solver model: available
Solver message: optimal

      21314     8.0989693901e+09 Pr: 0(0) 11s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-dq4yrqob
Model status        : Optimal
Simplex   iterations: 21314
Objective value     :  8.0989693901e+09
Relative P-D gap    :  3.2970714635e-15
HiGHS run time      :         10.99
Writing the solution to /tmp/linopy-solve-fh9ikhoj.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.

('ok', 'optimal')

n.statistics.optimal_capacity().div(1e3)  # GW

component    carrier                     
StorageUnit  battery storage                 12.342147
             hydrogen storage underground     5.446858
Generator    solar                           18.677068
             wind                            35.426277
dtype: float64

n.statistics.energy_balance().sort_values().div(1e6)  # TWh

component    carrier                       bus_carrier
Load         -                             electricity   -66.266089
StorageUnit  hydrogen storage underground  electricity   -10.304619
             battery storage               electricity    -0.495359
Generator    solar                         electricity    22.055171
             wind                          electricity    55.010896
dtype: float64

pd.concat({
    "capex": n.statistics.capex(),
    "opex": n.statistics.opex(),
}, axis=1).div(1e9).round(2) # bn€/a

		capex	opex
component	carrier
StorageUnit	battery storage	0.94	NaN
StorageUnit	hydrogen storage underground	2.52	NaN
Generator	solar	0.96	0.00
Generator	wind	3.60	0.08

n.storage_units.p_nom_opt.div(1e3) * n.storage_units.max_hours  # GWh

StorageUnit
battery storage                   49.368588
hydrogen storage underground    1830.144210
dtype: float64

plot_dispatch(n)

n.buses_t.marginal_price.sort_values(by="electricity", ascending=False).reset_index(
    drop=True
).plot(title="price duration curve [€/MWh]")

Finally, we will export this network so that we can build on it when adding further sectors (e.g. electric vehicles and heat pumps) in the next tutorial:

n.export_to_netcdf("electricity-network.nc");

INFO:pypsa.io:Exported network 'electricity-network.nc' contains: global_constraints, buses, storage_units, carriers, loads, generators

Exercises#

Explore how the model reacts to changing assumptions and available technologies. Here are a few inspirations, but choose in any order according to your interests:

What if the model were rerun with assumptions for 2050?
What if either hydrogen or battery storage cannot be expanded?
What if you could either only build solar or only build wind?
Vary the energy-to-power ratio of the hydrogen storage. What ratio leads to lowest costs?
On model.energy, you can download capacity factors for onshore wind and solar for any region in the world. What changes?
Add nuclear as another dispatchable low-emission generation technology (similar to OCGT). Perform a sensitivity analysis trying to answer how low the capital cost of a nuclear plant would need to be to be chosen.
How inaccurate is the 3-hourly resolution used for demonstration? How does it compare to hourly resolution? How much longer does it take to solve?

Capacity Expansion Planning

Contents

Capacity Expansion Planning#

From electricity market modelling to capacity expansion planning#

Package Imports#

Techology Data and Costs#

Capacity Factor and Load Time Series#

Building the Model#

Model Initialisation#

Adding Components#

Model Run#

Model Evaluation#

Plotting Optimal Dispatch#

Adding Storage Units#

Adding emission limits#

Exercises#