Problem 9.3

Problem 9.3#

Integrated Energy Grids

For the model built in PyPSA described in Problem 9.2, remove the CCGT with fix capacity, solve the optimization problem and:

a) Calculate the revenues collected by the OCGT plant throughout the year and show that their sum is equal to its costs.

b) Solve the problem for different CO\(_2\) values ranging from 5 MtCO2/year to zero. Plot the total system cost and the required CO\(_2\) prices as a function of the emissions allowance.

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

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["OCGT", "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'>

../_images/d40bdf91d7e9434172ab2ebd6dc441dc37a4d8cf3eac0000972359dac0943aa1.png

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'>

../_images/ed79caff22453fb39b57b258944527afb7c44d8d1041bd846450db5053dfabb6.png

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

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/15 [00:00<?, ?it/s]

Writing constraints.:  40%|████      | 6/15 [00:00<00:00, 31.27it/s]

Writing constraints.:  67%|██████▋   | 10/15 [00:00<00:00, 30.85it/s]

Writing constraints.:  93%|█████████▎| 14/15 [00:00<00:00, 24.33it/s]

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

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

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

INFO:linopy.io: Writing time: 0.7s

Running HiGHS 1.10.0 (git hash: fd86653): Copyright (c) 2025 HiGHS under MIT licence terms
LP   linopy-problem-95w649gm has 122645 rows; 52564 cols; 240922 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, 48202 cols, 179634 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, 48202 cols, 179634 nonzeros  0s
Presolve : Reductions: rows 65719(-56926); columns 48202(-4362); elements 179634(-61288)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(1.54638e+09) 0s

      22812     1.7882887540e+09 Pr: 5337(1.64127e+09); Du: 0(9.56358e-08) 6s

      31670     2.0918610969e+09 Pr: 7102(1.92092e+10); Du: 0(6.81106e-08) 11s

      41127     2.6199921359e+09 Pr: 7207(1.15154e+10); Du: 0(1.08686e-08) 16s

      47094     3.0134716311e+09 Pr: 32820(9.69131e+10); Du: 0(1.12399e-16) 21s

      52219     3.2791786966e+09 Pr: 3091(5.24729e+08); Du: 0(7.12469e-17) 26s

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

n.objective / 1e9
n.global_constraints

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

a) Calculate the revenues collected by the OCGT plant throughout the year and show that their sum is equal to its costs.

To calculate the revenues collected by every technology, we multiply the energy generated in every hour by the electricity price in that hour and sum for the entire year.

n.generators_t.p.multiply(n.buses_t.marginal_price.to_numpy()).sum().div(1e6) # EUR -> MEUR

Generator
OCGT      1415.170598
onwind    1340.130330
solar      675.780235
dtype: float64

To calculate the tax paid by the OCGT generator, we mutiply the sum of the energy generated by its emission factor, and the CO2 price.

(-n.global_constraints.mu)*n.generators_t.p.sum()['OCGT']*costs.at["OCGT", "CO2 intensity"]/costs.at["OCGT", "efficiency"]*0.000001 #tCO2/MtCO2

GlobalConstraint
CO2Limit    531.328616
Name: mu, dtype: float64

The market revenues, minus what is paid in CO2 tax,corresponds to the total cost for OCGT, which we can also read using the statistics module:

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

component    carrier        
Generator    OCGT                883.841982
             onwind             1340.130330
             solar               675.780235
StorageUnit  battery storage            NaN
dtype: float64