Capacity Expansion

Capacity Expansion#

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 -q pypsa pandas highspy "plotly<6"

In this tutorial, we want to build a replica of model.energy with PyPSA. As previously shown in the slides, 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.

This type of problem follows the model of an investment planning optimisation problem as discussed in the slides.

Note

Package Imports#

import pypsa
import pandas as pd
import plotly.io as pio
import plotly.offline as py

pd.options.plotting.backend = "plotly"

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

Capacity Factor & Load Time Series#

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

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

	load_mw	pv_pu	wind_pu
timestamp
2019-01-01 00:00:00	5719.26	0.000	0.1846
2019-01-01 01:00:00	5677.73	0.000	0.2293
2019-01-01 02:00:00	5622.20	0.000	0.2718
2019-01-01 03:00:00	5474.74	0.000	0.3146
2019-01-01 04:00:00	5432.51	0.000	0.3552
2019-01-01 05:00:00	5446.35	0.000	0.4055
2019-01-01 06:00:00	5413.39	0.000	0.4957
2019-01-01 07:00:00	5450.17	0.004	0.5808
2019-01-01 08:00:00	5449.62	0.026	0.6450

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

	load_mw	pv_pu	wind_pu
timestamp
2019-01-01 00:00:00	5719.26	0.0	0.1846
2019-01-01 03:00:00	5474.74	0.0	0.3146
2019-01-01 06:00:00	5413.39	0.0	0.4957

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");

…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 = [
    "wind",
    "solar",
    "OCGT",
    "hydrogen storage underground",
    "battery storage",
    "electricity",
]

n.add(
    "Carrier",
    carriers,
    color=["dodgerblue", "gold", "indianred", "magenta", "yellowgreen", "black"],
    co2_emissions=[0, 0, 0.2, 0, 0, 0],
);

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

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

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=annuity(0.07, 30) * 500_000,
    marginal_cost=65,
    efficiency=0.4,
    p_nom_extendable=True,
);

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=annuity(0.07, 30) * 1_100_000,
    marginal_cost=1,
    p_nom_extendable=True,
);

n.add(
    "Generator",
    "solar",
    bus="electricity",
    carrier="solar",
    p_max_pu=ts.pv_pu,
    capital_cost=annuity(0.05, 25) * 500_000,
    marginal_cost=0.5,
    p_nom_extendable=True,
);

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(log_to_console=False)

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.12s

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

('ok', 'optimal')

Model Evaluation#

The total system cost in billion Euros per year:

n.objective / 1e9

3.88757918749457

The optimised capacities in GW:

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

Generator
OCGT     10.103965
wind     12.685117
solar    14.186638
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  electricity   -66.266089
Generator  solar        electricity    16.237620
           wind         electricity    21.988632
           OCGT         electricity    28.039837
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       2229.710534
           solar       511.407233
           wind       1146.461421
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

np.float64(14.019918660988745)

n.statistics.energy_balance.iplot()

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=annuity(0.07, 10) * 150_000 + 4 * annuity(0.07, 20) * 150_000,
    efficiency_store=0.95,
    efficiency_dispatch=0.95,
    p_nom_extendable=True,
    cyclic_state_of_charge=True,
);

Second, the hydrogen storage. This one is composed of an electrolysis to convert electricity to hydrogen, a fuel cell (or hydrogen turbine) 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 = (
    annuity(0.07, 25) * 1_500_000
    + annuity(0.07, 25) * 550_000
    + 336 * annuity(0.07, 100) * 2_000
)

n.add(
    "StorageUnit",
    "hydrogen storage underground",
    bus="electricity",
    carrier="hydrogen storage underground",
    max_hours=336,
    capital_cost=capital_costs,
    efficiency_store=0.65,
    efficiency_dispatch=0.45,
    p_nom_extendable=True,
    cyclic_state_of_charge=True,
);

Adding emission limits#

Let’s also constrain the model to 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,
);

Then, we can re-run the model with the new components and constraints:

n.optimize(log_to_console=False)

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

component    carrier                     
Generator    solar                           24.510751
             wind                            33.953472
StorageUnit  battery storage                  4.971099
             hydrogen storage underground     7.103009
dtype: float64

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

component    carrier                       bus_carrier
Load         electricity                   electricity   -66.266089
StorageUnit  hydrogen storage underground  electricity   -15.178400
             battery storage               electricity    -0.425633
Generator    solar                         electricity    28.943986
             wind                          electricity    52.926136
dtype: float64

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

		capex	opex
component	carrier
Generator	solar	0.87	0.01
Generator	wind	3.01	0.05
StorageUnit	battery storage	0.39	NaN
StorageUnit	hydrogen storage underground	1.58	NaN

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

StorageUnit
battery storage                   19.884394
hydrogen storage underground    2386.610887
dtype: float64

n.statistics.energy_balance.iplot()

n.buses_t.marginal_price.plot(title="prices [€/MWh]")

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

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 solar and batteries had 30% lower investment costs? You can alter the costs with n.storage_units.loc["StorageName", "capital_cost"] = new_value.
What if either hydrogen or battery storage cannot be expanded? You can remove components with n.remove("StorageUnit", "ComponentName").
What if you could either only build solar or only build wind? You can remove components with n.remove("Generator", "ComponentName").
Vary the energy-to-power ratio of the hydrogen storage. What ratio leads to lowest costs? You can change the ratio with n.storage_units.loc["StorageName", "max_hours"] = new_value.
On model.energy, you can download capacity factors for onshore wind and solar for any region in the world as CSV files. What changes if you select another region? You can read the CSV files from URL with pd.read_csv("URL", index_col=0, parse_dates=True).