Security of Supply Analysis - PyPSA Workshop Birmingham

Note

On the workshop JupyterHub, everything is already installed - just log in and run the cells below, nothing to set up.

No hub access? You can also run this notebook in your browser on Google Colab (a Google account is required). On Colab, first run this in a cell at the top:

!pip install pypsa "plotly<6"

In this tutorial, we will analyse the security of supply of different scenarios using a pre-built 3-node model of the Philippine power system (Luzon, Visayas, Mindanao). In this setting, we will focus on dispatch optimisation under different stress events and evaluate the performance of the system using different metrics.

We start by importing the necessary libraries and loading the base network. This network is amended by a load shedding generator at each node and we also resample the time series to 4h resolution to speed up the optimisation.

import pandas as pd
import pypsa

pypsa.options.params.optimize.include_objective_constant = False
pypsa.options.params.optimize.log_to_console = False

pd.options.plotting.backend = "plotly"

# Backup: https://cloud.pypsalabs.org/s/KPD4AcTtJx2b9Hk/download/phl-network.nc
url = "https://raw.githubusercontent.com/pypsalabs/workshop-202609/main/bham/data/phl-network.nc"

base = pypsa.Network(url)

INFO:pypsa.network.io:Retrieving network data from https://raw.githubusercontent.com/pypsalabs/workshop-202609/main/bham/data/phl-network.nc.

WARNING:pypsa.network.io:Importing network from PyPSA version v1.1.2 while current version is v1.2.4. Read the release notes at `https://go.pypsa.org/release-notes` to prepare your network for import.

INFO:pypsa.network.io:New version 1.3.0 available! (Current: 1.2.4)

INFO:pypsa.network.io:Imported network 'Unnamed Network' has buses, carriers, generators, links, loads, storage_units

base.add(
    "Generator",
    base.buses.index + " load-shedding",
    bus=base.buses.index,
    carrier="load-shedding",
    marginal_cost=2000,
    p_nom=base.loads_t.p_set.max().max(),
)
base.add("Carrier", "load-shedding", color="crimson")

base = base.cluster.temporal.resample("8h")

n0 = base.copy()
n0.name = "base"

Events¶

Power Plant Decommissioning¶

We decommission 2 GW of coal generation capacity in Luzon and all coal generation capacity in Mindanao.

n1 = base.copy()
n1.name = "coal-decommissioning"
n1.generators.loc["Luzon coal-subcritical", "p_nom"] -= 2000
n1.generators.loc["Mindanao coal-subcritical", "p_nom"] -= 836

Power Plant Outage¶

We simulate the outage of Luzon’s supercritical coal power plants in all of February.

n2 = base.copy()
n2.name = "coal-generator-outage"

n2.generators_t.p_max_pu.loc["2025-02", "Luzon coal-supercritical"] = 0
n2.generators_t.p_max_pu.fillna(1, inplace=True)

Interconnector Outage¶

We simulate the outage of the interconnectors between all three regions throughout the full year (i.e. no power can be transmitted between the regions).

n3 = base.copy()
n3.name = "interconnector-outage"

n3.links.p_max_pu = 0
n3.links.p_min_pu = 0

Volcano Eruption¶

We simulate a volcano eruption that reduces solar photovoltaics generation by 80% from February to June.

n4 = base.copy()
n4.name = "volcano-eruption"

pv_gens = n4.generators.carrier == "photovoltaic"

n4.generators_t.p_max_pu.loc[
    "2025-02":"2025-06",
    pv_gens,
] *= 0.2

Price Shock in Oil and Gas¶

We simulate a price shock in oil and gas (e.g. due to a geopolitical crisis) that doubles the fuel costs of oil and gas generators for the full year.

n5 = base.copy()
n5.name = "price-shock-oil-gas"

selection = n5.generators.carrier.str.contains("oil|gas")
n5.generators.loc[selection, "marginal_cost"] *= 2

Demand Increase¶

We simulate a demand increase of 20% across all three regions for the full year (e.g. due to economic growth).

n6 = base.copy()
n6.name = "demand-increase"

n6.loads_t.p_set *= 1.2

Drought¶

We simulate a drought event that reduces hydro generation by 50% for the full year.

n7 = base.copy()
n7.name = "drought"

# drought reduces hydro generation by 50%
hydro = n7.generators.carrier == "hydro-reservoir-and-run-of-river"
n7.generators_t.p_max_pu.loc[:, hydro] *= 0.5

Coal Import Limitation¶

We simulate a coal usage restriction that limits the total coal consumption across the system to 80 TWh thermal for the full year (e.g. due to a coal import limitation).

Note, that we need to factor in the efficiency of the coal generators to calculate the thermal energy consumption from the electrical energy generation.

This is a constraint that is not directly supported by PyPSA, so we need to add it manually to the model as a custom constraint after the model has been created but before it is optimised.

n8 = base.copy()
n8.name = "coal-import-limit"

n8.optimize.create_model()
coal_i = n8.generators.loc[n8.generators.carrier.str.contains("coal")].index
lhs = (
    n8.model["Generator-p"]
    .loc[:, coal_i]
    .mul(n8.snapshot_weightings.generators)
    .sum(dim="snapshot")
    .div(n8.generators.loc[coal_i, "efficiency"])
    .sum()
)
n8.model.add_constraints(lhs <= 80e6, name="Generator-coal_limit");

Reserve Constraints¶

In a similar way, we can also add constraints to ensure that there is sufficient reserve capacity available in the system to cover for unexpected events. For example, we can require that there is always at least 20% of the hour’s demand or at least 3 GW of reserve capacity available.

n9 = base.copy()
n9.name = "reserve-constraint"

n9.optimize.create_model()

regex = "coal|gas|geothermal|oil|bioenergy"
reserve_i = n9.generators.loc[n9.generators.carrier.str.contains(regex)].index

We use the dispatch of firm generators (i.e. generators that are not renewable):

p_reserve = n9.model["Generator-p"].loc[:, reserve_i].sum(dim="name")
p_reserve

LinearExpression [snapshot: 1095]:
----------------------------------
[2025-01-01 08:00:00]: +1 Generator-p[2025-01-01 08:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-01-01 08:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-01-01 08:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-01-01 08:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-01-01 08:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-01-01 08:00:00, Visayas gas-ccs]
[2025-01-01 16:00:00]: +1 Generator-p[2025-01-01 16:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-01-01 16:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-01-01 16:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-01-01 16:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-01-01 16:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-01-01 16:00:00, Visayas gas-ccs]
[2025-01-02 00:00:00]: +1 Generator-p[2025-01-02 00:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-01-02 00:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-01-02 00:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-01-02 00:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-01-02 00:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-01-02 00:00:00, Visayas gas-ccs]
[2025-01-02 08:00:00]: +1 Generator-p[2025-01-02 08:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-01-02 08:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-01-02 08:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-01-02 08:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-01-02 08:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-01-02 08:00:00, Visayas gas-ccs]
[2025-01-02 16:00:00]: +1 Generator-p[2025-01-02 16:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-01-02 16:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-01-02 16:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-01-02 16:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-01-02 16:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-01-02 16:00:00, Visayas gas-ccs]
[2025-01-03 00:00:00]: +1 Generator-p[2025-01-03 00:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-01-03 00:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-01-03 00:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-01-03 00:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-01-03 00:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-01-03 00:00:00, Visayas gas-ccs]
[2025-01-03 08:00:00]: +1 Generator-p[2025-01-03 08:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-01-03 08:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-01-03 08:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-01-03 08:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-01-03 08:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-01-03 08:00:00, Visayas gas-ccs]
		...
[2025-12-30 00:00:00]: +1 Generator-p[2025-12-30 00:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-12-30 00:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-12-30 00:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-12-30 00:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-12-30 00:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-12-30 00:00:00, Visayas gas-ccs]
[2025-12-30 08:00:00]: +1 Generator-p[2025-12-30 08:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-12-30 08:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-12-30 08:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-12-30 08:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-12-30 08:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-12-30 08:00:00, Visayas gas-ccs]
[2025-12-30 16:00:00]: +1 Generator-p[2025-12-30 16:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-12-30 16:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-12-30 16:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-12-30 16:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-12-30 16:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-12-30 16:00:00, Visayas gas-ccs]
[2025-12-31 00:00:00]: +1 Generator-p[2025-12-31 00:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-12-31 00:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-12-31 00:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-12-31 00:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-12-31 00:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-12-31 00:00:00, Visayas gas-ccs]
[2025-12-31 08:00:00]: +1 Generator-p[2025-12-31 08:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-12-31 08:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-12-31 08:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-12-31 08:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-12-31 08:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-12-31 08:00:00, Visayas gas-ccs]
[2025-12-31 16:00:00]: +1 Generator-p[2025-12-31 16:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2025-12-31 16:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2025-12-31 16:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2025-12-31 16:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2025-12-31 16:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2025-12-31 16:00:00, Visayas gas-ccs]
[2026-01-01 00:00:00]: +1 Generator-p[2026-01-01 00:00:00, Mindanao coal-ultrasupercritical] + 1 Generator-p[2026-01-01 00:00:00, Mindanao geothermal-unspecified] + 1 Generator-p[2026-01-01 00:00:00, Mindanao bioenergy-unspecified] ... +1 Generator-p[2026-01-01 00:00:00, Visayas coal-circulating-fluidized-bed] + 1 Generator-p[2026-01-01 00:00:00, Visayas gas-combined-cycle] + 1 Generator-p[2026-01-01 00:00:00, Visayas gas-ccs]

as well as the capacity of firm generators to calculate the available reserve capacity in each hour and add a constraint to ensure that this is always above the required reserve margin:

p_nom_reserve = n9.generators.p_nom.loc[reserve_i].sum()
p_nom_reserve

np.float64(22957.079999999998)

The reserve requirement we calculate as the maximum of 20% of the demand or 3 GW:

reserve_required = n9.loads_t.p_set.sum(axis=1).mul(0.2).clip(lower=3000)
reserve_required.plot()

This is added as a constraint to the model:

expr = p_nom_reserve - p_reserve >= reserve_required
n9.model.add_constraints(expr, name="reserve_constraints");

Compound Stress Events¶

You decide which events to combine.

n10 = base.copy()
n10.name = "compound-event"

Solving Scenarios¶

Now, we can solve the different scenarios. Note that the cases where we have added custom constraints do not require to rebuild the model, but just require it to be solved, so the syntax is slightly different for these cases.

for n in [n0, n1, n2, n3, n4, n5, n6, n7, n10]:
    n.optimize()

n8.optimize.solve_model()
n9.optimize.solve_model()

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io:Writing objective.

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

Writing constraints.:  92%|█████████▏| 11/12 [00:00<00:00, 84.47it/s]

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

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

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


INFO:linopy.io: Writing time: 0.2s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 77745 primals, 165345 duals
Objective: 4.35e+09
Solver: highs
Runtime: 0.49s
MIP gap: inf
Dual bound: 0.00e+00
Solver model: available
Solver message: Optimal

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

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io:Writing objective.

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

Writing constraints.:  92%|█████████▏| 11/12 [00:00<00:00, 91.71it/s]

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

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

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


INFO:linopy.io: Writing time: 0.18s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 77745 primals, 165345 duals
Objective: 4.45e+09
Solver: highs
Runtime: 0.82s
MIP gap: inf
Dual bound: 0.00e+00
Solver model: available
Solver message: Optimal

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

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io:Writing objective.

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

Writing constraints.:  92%|█████████▏| 11/12 [00:00<00:00, 101.45it/s]

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

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

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


INFO:linopy.io: Writing time: 0.16s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 77745 primals, 165345 duals
Objective: 4.36e+09
Solver: highs
Runtime: 0.72s
MIP gap: inf
Dual bound: 0.00e+00
Solver model: available
Solver message: Optimal

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

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io:Writing objective.

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

Writing constraints.:  92%|█████████▏| 11/12 [00:00<00:00, 89.18it/s]

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

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

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


INFO:linopy.io: Writing time: 0.19s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 77745 primals, 165345 duals
Objective: 4.39e+09
Solver: highs
Runtime: 0.54s
MIP gap: inf
Dual bound: 0.00e+00
Solver model: available
Solver message: Optimal

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

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io:Writing objective.

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

Writing constraints.:  92%|█████████▏| 11/12 [00:00<00:00, 87.67it/s]

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

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

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


INFO:linopy.io: Writing time: 0.19s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 77745 primals, 165345 duals
Objective: 4.45e+09
Solver: highs
Runtime: 0.68s
MIP gap: inf
Dual bound: 0.00e+00
Solver model: available
Solver message: Optimal

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

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io:Writing objective.

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

Writing constraints.:  92%|█████████▏| 11/12 [00:00<00:00, 90.32it/s]

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

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

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


INFO:linopy.io: Writing time: 0.18s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 77745 primals, 165345 duals
Objective: 4.36e+09
Solver: highs
Runtime: 0.73s
MIP gap: inf
Dual bound: 0.00e+00
Solver model: available
Solver message: Optimal

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

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io:Writing objective.

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

Writing constraints.:  83%|████████▎ | 10/12 [00:00<00:00, 99.87it/s]

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

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

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


INFO:linopy.io: Writing time: 0.2s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 77745 primals, 165345 duals
Objective: 5.77e+09
Solver: highs
Runtime: 0.54s
MIP gap: inf
Dual bound: 0.00e+00
Solver model: available
Solver message: Optimal

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

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io:Writing objective.

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

Writing constraints.:  92%|█████████▏| 11/12 [00:00<00:00, 102.11it/s]

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

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

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


INFO:linopy.io: Writing time: 0.16s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 77745 primals, 165345 duals
Objective: 4.68e+09
Solver: highs
Runtime: 0.51s
MIP gap: inf
Dual bound: 0.00e+00
Solver model: available
Solver message: Optimal

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

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io:Writing objective.

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

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

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

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

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


INFO:linopy.io: Writing time: 0.14s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 77745 primals, 165345 duals
Objective: 4.35e+09
Solver: highs
Runtime: 0.53s
MIP gap: inf
Dual bound: 0.00e+00
Solver model: available
Solver message: Optimal

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

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io:Writing objective.

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

Writing constraints.:  85%|████████▍ | 11/13 [00:00<00:00, 104.13it/s]

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

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

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


INFO:linopy.io: Writing time: 0.16s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 77745 primals, 165346 duals
Objective: 6.18e+09
Solver: highs
Runtime: 0.83s
MIP gap: inf
Dual bound: 0.00e+00
Solver model: available
Solver message: Optimal

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

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io:Writing objective.

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

Writing constraints.:  77%|███████▋  | 10/13 [00:00<00:00, 96.79it/s]

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

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

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


INFO:linopy.io: Writing time: 0.22s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 77745 primals, 166440 duals
Objective: 4.35e+09
Solver: highs
Runtime: 0.54s
MIP gap: inf
Dual bound: 0.00e+00
Solver model: available
Solver message: Optimal

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

('ok', 'optimal')

Network Collections¶

PyPSA has a pypsa.NetworkCollection class that provides a convenient way to manage and analyse multiple networks simultaneously, access their combined data, and generate combined statistics and plots. It uses the names of the networks as indices.

nc = pypsa.NetworkCollection([n0, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10])

Evaluation Metrics¶

Operational Cost¶

As absolute values by technology, as total, as relative difference, or as absolute difference to the base case.

opex = nc.statistics.opex().div(1e6).round(1).droplevel("component").unstack("network")

opex.sum()

network
base                     4348.6
coal-decommissioning     4446.3
coal-generator-outage    4359.5
coal-import-limit        6184.1
compound-event           4348.6
demand-increase          5773.0
drought                  4684.7
interconnector-outage    4388.0
price-shock-oil-gas      4357.4
reserve-constraint       4348.6
volcano-eruption         4450.3
dtype: float64

opex.sum() / opex.sum().min()

network
base                     1.000000
coal-decommissioning     1.022467
coal-generator-outage    1.002507
coal-import-limit        1.422090
compound-event           1.000000
demand-increase          1.327554
drought                  1.077289
interconnector-outage    1.009060
price-shock-oil-gas      1.002024
reserve-constraint       1.000000
volcano-eruption         1.023387
dtype: float64

opex.sum() - opex.sum().min()

network
base                        0.0
coal-decommissioning       97.7
coal-generator-outage      10.9
coal-import-limit        1835.5
compound-event              0.0
demand-increase          1424.4
drought                   336.1
interconnector-outage      39.4
price-shock-oil-gas         8.8
reserve-constraint          0.0
volcano-eruption          101.7
dtype: float64

Electricity Mix¶

energy = (
    nc.statistics.energy_balance()
    .div(1e6)
    .round(1)
    .groupby(["network", "carrier"])
    .sum()
    .unstack("network")
    .drop("AC")
)
energy.T.plot(kind="bar", barmode="stack", labels={"value": "Electricity Supply [TWh]"})

nc["demand-increase"].statistics.energy_balance.iplot()

Prices¶

We can look at average prices or just the frequency of hours with very high prices (e.g. above 100 $/MWh).

nc.statistics.prices.iplot()

(nc.buses_t.marginal_price > 100).sum().unstack("network").div(8760 / 8 / 100).round(1)

Expected Energy Not Served¶

The expected energy not served (EENS) is a common reliability metric that quantifies the expected amount of energy demand that cannot be met due to supply shortages. It is calculated as the energy used by the load shedding generators in our model.

# EENS (all scenarios share the same snapshot weightings, so we use the base network's)
(
    n0.snapshot_weightings.generators
    @ nc.generators_t.p.loc[:, nc.generators.carrier == "load-shedding"]
).unstack("network")

Loss of Load Probability¶

The loss of load probability (LOLP) is a reliability metric that quantifies the probability that the system will not be able to meet the demand in a given hour. It is calculated as the number of hours with load shedding divided by the total number of hours.

# LOLP
(
    nc.generators_t.p.loc[:, nc.generators.carrier == "load-shedding"] > 0.1
).sum().unstack("network") / (8760 / 4)

Maximum Loss of Load¶

The maximum loss of load (MLOL) is a reliability metric that quantifies the maximum amount of power demand that cannot be met in any given hour. It is calculated as the maximum power used by the load shedding generators in our model.

nc.generators_t.p.loc[:, nc.generators.carrier == "load-shedding"].max().unstack(
    "network"
)

Reserve Margin¶

The reserve margin is a static reliability metric that quantifies the amount of reserve capacity available in the system with the peak demand subtracted. Here, for the base network:

(
    n0.generators.filter(regex="coal|geothermal|oil|gas|bioenergy", axis=0)
    .groupby("bus")
    .p_nom.sum()
    - n0.loads_t.p_set.max()
)

bus
Luzon       3480.20
Mindanao     671.21
Visayas      772.67
dtype: float64

Exercises¶

Task 1: For scenario n10, create different combinations of the above events and analyse the results. For example, you could combine a drought with a coal import limitation, or a power plant outage with a price shock in oil and gas.

Task 2: Change the capacities of the system (e.g. increase solar capacity, add battery storage or increase gas capacity). Repeat the analysis for the different scenarios and compare the results. Can you identify effective strategies for improving the security of supply? Which changes cause more load shedding?