Problem 9.1

Integrated Energy Grids

Problem 9.2. Join capacity and dispatch optimization with storage and CO2 emissions limits.

Use the model built in PyPSA described in Problem 8.2 and assume that methane gas emits 0.198 tCO2 per MWh of thermal energy contained in the gas. Limit the maximum CO2 emissions to 5 MtCO2/year.

a) Calculate the optimal installed capacities and plot the hourly generation and demand during January.

b) What is the CO2 tax required to attain such CO2 emissions limit?

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 pandas and numpy by executing the following command in a Jupyter cell at the top of the notebook.

!pip install pandas pypsa

Note

See also https://model.energy.

In this exercise, 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.

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

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

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

Now we can for example read the capital and marginal cost of onshore wind and solar, or the emissions factors of the carrier gas used in and OCGT

costs.at["onwind", "capital_cost"] #EUR/MW/a

np.float64(101644.12332388277)

costs.at["solar", "capital_cost"] #EUR/MW/a

np.float64(51346.82981964593)

costs.at["CCGT", "CO2 intensity"] #tCO2/MWh_th

np.float64(0.198)

Retrieving time series data

In this example, wind data from https://zenodo.org/record/3253876#.XSiVOEdS8l0 and solar PV data from https://zenodo.org/record/2613651#.X0kbhDVS-uV is used. The data is downloaded in csv format and saved in the ‘data’ folder. The Pandas package is used as a convenient way of managing the datasets.

For convenience, the column including date information is converted into Datetime and set as index

data_solar = pd.read_csv('data/pv_optimal.csv',sep=';')
data_solar.index = pd.DatetimeIndex(data_solar['utc_time'])

data_wind = pd.read_csv('data/onshore_wind_1979-2017.csv',sep=';')
data_wind.index = pd.DatetimeIndex(data_wind['utc_time'])

data_el = pd.read_csv('data/electricity_demand.csv',sep=';')
data_el.index = pd.DatetimeIndex(data_el['utc_time'])

The data format can now be analyzed using the .head() function to show the first lines of the data set

data_solar.head()

	utc_time	AUT	BEL	BGR	BIH	CHE	CYP	CZE	DEU	DNK	...	MLT	NLD	NOR	POL	PRT	ROU	SRB	SVK	SVN	SWE
utc_time
1979-01-01 00:00:00+00:00	1979-01-01T00:00:00Z	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
1979-01-01 01:00:00+00:00	1979-01-01T01:00:00Z	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
1979-01-01 02:00:00+00:00	1979-01-01T02:00:00Z	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
1979-01-01 03:00:00+00:00	1979-01-01T03:00:00Z	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
1979-01-01 04:00:00+00:00	1979-01-01T04:00:00Z	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

5 rows × 33 columns

We will use timeseries for Portugal in this excercise

country = 'PRT'

Join capacity and dispatch optimization

For building the model, we start again by initialising an empty network, adding the snapshots, and the electricity bus.

n = pypsa.Network()
hours_in_2015 = pd.date_range('2015-01-01 00:00Z',
                              '2015-12-31 23:00Z',
                              freq='h')

n.set_snapshots(hours_in_2015.values)

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

n.snapshots

DatetimeIndex(['2015-01-01 00:00:00', '2015-01-01 01:00:00',
               '2015-01-01 02:00:00', '2015-01-01 03:00:00',
               '2015-01-01 04:00:00', '2015-01-01 05:00:00',
               '2015-01-01 06:00:00', '2015-01-01 07:00:00',
               '2015-01-01 08:00:00', '2015-01-01 09:00:00',
               ...
               '2015-12-31 14:00:00', '2015-12-31 15:00:00',
               '2015-12-31 16:00:00', '2015-12-31 17:00:00',
               '2015-12-31 18:00:00', '2015-12-31 19:00:00',
               '2015-12-31 20:00:00', '2015-12-31 21:00:00',
               '2015-12-31 22:00:00', '2015-12-31 23:00:00'],
              dtype='datetime64[ns]', name='snapshot', length=8760, freq=None)

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

carriers = [
    "onwind",
    "solar",
    "OCGT",
    "CCGT",
    "battery storage",
]

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

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

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

# add load to the bus
n.add("Load",
      "demand",
      bus="electricity",
      p_set=data_el[country].values)

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

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

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

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

We add now the generators and set up their capacities to be extendable so that they can be optimized together with the dispatch time series. For the wind and solar generator, we need to indicate the capacity factor or maximum power per unit ‘p_max_pu’

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

CF_wind = data_wind[country][[hour.strftime("%Y-%m-%dT%H:%M:%SZ") for hour in n.snapshots]]
n.add(
        "Generator",
        "onwind",
        bus="electricity",
        carrier="onwind",
        p_max_pu=CF_wind.values,
        capital_cost=costs.at["onwind", "capital_cost"],
        marginal_cost=costs.at["onwind", "marginal_cost"],
        efficiency=costs.at["onwind", "efficiency"],
        p_nom_extendable=True,
    )

CF_solar = data_solar[country][[hour.strftime("%Y-%m-%dT%H:%M:%SZ") for hour in n.snapshots]]
n.add(
        "Generator",
        "solar",
        bus="electricity",
        carrier="solar",
        p_max_pu= CF_solar.values,
        capital_cost=costs.at["solar", "capital_cost"],
        marginal_cost=costs.at["solar", "marginal_cost"],
        efficiency=costs.at["solar", "efficiency"],
        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["2015-01"].plot(figsize=(6, 2), ylabel="CF")

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

We add the battery storage, assuming a fixed energy-to-power ratio of 2 hours, i.e. if fully charged, the battery can discharge at full capacity for 2 hours.

For the capital cost, we have to factor in both the capacity and energy cost of the storage.

We include the charging and discharging efficiencies we 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=2,
    capital_cost=costs.at["battery inverter", "capital_cost"]
    + 2 * 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')

We add the Combined Cycle Gas Turbine (CCGT). In this case, its capacity is not extendable but fixed to 1 GW.

n.add(
    "Generator",
    "CCGT",
    bus="electricity",
    carrier="CCGT",
    capital_cost=costs.at["CCGT", "capital_cost"],
    marginal_cost=costs.at["CCGT", "marginal_cost"],
    efficiency=costs.at["CCGT", "efficiency"],
    p_nom=6000, #6 Gw
)

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

Model Run

Before solving the problem we add a CO\(_2\) emission limit as global constraint:

n.add(
    "GlobalConstraint",
    "CO2Limit",
    carrier_attribute="co2_emissions",
    sense="<=",
    constant=5000000, #5MtCO2
)

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

We can already solved the model using the open-solver “highs” or the commercial solver “gurobi” with the academic license

n.optimize(solver_name="highs")

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

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io:Writing objective.

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

Writing constraints.:  41%|████      | 7/17 [00:00<00:00, 50.24it/s]

Writing constraints.:  76%|███████▋  | 13/17 [00:00<00:00, 30.18it/s]

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

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

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

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

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

INFO:linopy.io: Writing time: 0.8s

Running HiGHS 1.10.0 (git hash: fd86653): Copyright (c) 2025 HiGHS under MIT licence terms
LP   linopy-problem-jzpsxx8v has 140165 rows; 61324 cols; 275962 nonzeros
Coefficient ranges:
  Matrix [1e-03, 2e+00]
  Cost   [1e-02, 1e+05]
  Bound  [0e+00, 0e+00]
  RHS    [3e+03, 5e+06]
Presolving model
65719 rows, 56962 cols, 197154 nonzeros  0s
Dependent equations search running on 17520 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
65719 rows, 56962 cols, 197154 nonzeros  0s
Presolve : Reductions: rows 65719(-74446); columns 56962(-4362); elements 197154(-78808)
Solving the presolved LP

Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(1.53237e+09) 0s

      21775     1.2576741977e+09 Pr: 5885(2.87659e+09) 6s

      29815     1.5787401769e+09 Pr: 5945(1.14792e+10); Du: 0(1.14654e-07) 11s

      38888     2.0183783593e+09 Pr: 5960(2.16594e+09); Du: 0(3.12448e-08) 16s

      44292     2.2997610910e+09 Pr: 4928(5.84017e+08); Du: 0(1.07187e-07) 21s

      49008     2.5299859408e+09 Pr: 3962(1.03815e+09); Du: 0(7.59424e-08) 27s

Now, we can look at the results and evaluate the total system cost (in billion Euros per year)

n.objective / 1e9

2.7599916471178565

The optimised capacities in GW:

n.generators.p_nom_opt.div(1e3)  # MW -> GW

Generator
OCGT      -0.000000
onwind    12.476224
solar     10.342911
CCGT       6.000000
Name: p_nom_opt, dtype: float64

n.storage_units.p_nom_opt.div(1e3)  # MW -> GW

StorageUnit
battery storage    4.915387
Name: p_nom_opt, dtype: float64

The total energy generation by technology in TWh:

n.generators_t.p.sum().div(1e6)  # MWh -> TWh

Generator
OCGT       0.000000
onwind    18.934473
solar     15.583091
CCGT      14.646465
dtype: float64

We can plot the dispatch of every generator thoughout January

n.generators_t.p.loc["2015-01"].plot.area(figsize=(6, 2), ylabel="dispatch")

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

We can also plot the charging and discharging of the battery

n.storage_units_t.p.loc["2015-01"].plot(figsize=(6, 2), ylabel="battery")

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

b) What is the CO2 tax required to attain such CO2 emissions limit?

We can read the Lagrange/KKT multiplier associated to the maximum CO2 limit constraint

n.global_constraints

	type	investment_period	carrier_attribute	sense	constant	mu
GlobalConstraint
CO2Limit	primary_energy	NaN	co2_emissions	<=	5000000.0	-189.568204

In this case, we get a negative multiplier indicating that when the constant in the constraint equation (allowed CO2 emisssions) increases, the objective function (total system costs) decreases.