Modelling to Generate Alternatives¶

In this example, we apply MGA ('modelling to generate alternatives') to a single node capacity expansion model in the style of model.energy.

The MGA algorithm, which can be called with n.optimize.optimize_mga(), tries to minimize or maximize investment or dispatch in (groups of) technologies within a set cost budget.

For instance, it can be used to minimize the amount of wind capacity while keeping costs within 5% of the cost-optimal solution in terms of system costs.

See also model.energy and Neumann and Brown (2021) which uses PyPSA for MGA-type analysis.

In [2]:

import pandas as pd

import pypsa

import pandas as pd import pypsa

Find cost-optimal solution¶

Running MGA requires knowledge of what the total system costs are in the optimum. So first, we need to solve for the cost-optimal solution.

In [3]:

n = pypsa.examples.model_energy()
n.optimize()
n.statistics.capex().sum() + n.statistics.opex().sum()

n = pypsa.examples.model_energy() n.optimize() n.statistics.capex().sum() + n.statistics.opex().sum()

/home/runner/work/PyPSA/PyPSA/pypsa/network/io.py:2082: FutureWarning: pandas infers the `str` dtype for string data since its version 3.0. PyPSA still converts it back to numpy object dtype on import, but will keep it from PyPSA 2.0 on. Set `pypsa.options.api.legacy_string_dtype` explicitly to suppress this warning.
  new_static = _coerce_string_dtypes(new_static)

INFO:pypsa.network.io:Imported network 'Model-Energy' has buses, carriers, generators, links, loads, storage_units, stores

/tmp/ipykernel_3559/2384500068.py:2: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize()

INFO:linopy.model: Solve problem using Gurobi solver

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

INFO:linopy.io:Writing objective.

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

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

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

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

INFO:linopy.io: Writing time: 0.14s

Set parameter WLSAccessID

Set parameter WLSSecret

Set parameter LicenseID to value 2537914

Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de

Read LP format model from file /tmp/linopy-problem-ieq0ctyv.lp

Reading time = 0.07 seconds

obj: 64246 rows, 29206 columns, 124144 nonzeros

Set parameter LogToConsole to value 0

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 29206 primals, 64246 duals
Objective: 8.08e+09
Solver: gurobi
Runtime: 1.26s
Dual bound: 8.08e+09
Solver model: available
Solver message: 2

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-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, Store-energy_balance were not assigned to the network.

Out[3]:

np.float64(8078135675.45119)

Extract cost-optimal results¶

Optimal total system cost by technology:

In [4]:

tsc = (
    pd.concat([n.statistics.capex(), n.statistics.opex()], axis=1).sum(axis=1).div(1e9)
)
optimal_cost = tsc.sum()
tsc

tsc = ( pd.concat([n.statistics.capex(), n.statistics.opex()], axis=1).sum(axis=1).div(1e9) ) optimal_cost = tsc.sum() tsc

Out[4]:

component    carrier         
Generator    solar               1.341015
             wind                3.300830
Link         electrolysis        0.570894
             turbine             1.172438
StorageUnit  battery storage     0.941196
Store        hydrogen storage    0.561618
Generator    load shedding       0.190144
dtype: float64

The optimised capacities in GW (GWh for Store component):

In [5]:

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

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

Out[5]:

component    carrier         
Generator    load shedding         10.901160
             solar                 26.116801
             wind                  32.474381
Link         electrolysis           3.025153
             turbine               10.073615
StorageUnit  battery storage       14.854330
Store        hydrogen storage    3786.558312
dtype: float64

Energy balances on electricity side (in TWh):

In [6]:

n.statistics.energy_balance(bus_carrier="electricity").sort_values().div(1e6)

n.statistics.energy_balance(bus_carrier="electricity").sort_values().div(1e6)

Out[6]:

component    carrier          bus_carrier
Load         -                electricity   -66.266089
Link         electrolysis     electricity   -15.295140
StorageUnit  battery storage  electricity    -0.821922
Generator    load shedding    electricity     0.095072
Link         turbine          electricity     3.898685
Generator    solar            electricity    30.603101
             wind             electricity    47.786293
dtype: float64

Energy balances plot as time series (in MW):

In [7]:

n.statistics.energy_balance.plot.area(linewidth=0, bus_carrier="electricity")

n.statistics.energy_balance.plot.area(linewidth=0, bus_carrier="electricity")

Out[7]:

(<Figure size 759x300 with 1 Axes>,
 <Axes: xlabel='snapshot', ylabel='Energy Balance []'>,
 <seaborn.axisgrid.FacetGrid at 0x7f56609d70e0>)

Find lowest wind capacity within 5% cost slack¶

The function n.optimize.optimize_mga takes three main arguments:

1. The slack for the allowed relative cost deviation from the cost-optimum (0.05 corresponds to 5%).
2. A dictionary of weights for defining the new objective function. The first level defines the component (e.g. "Generator"), the second level defines the optimisation variable (e.g. p_nom for investment), and the third level defines the component name (e.g. from n.generators.index).
3. The sense, noting whether to minimizes ("min") or maximize ("max") the new objective.

In [8]:

weights = {"Generator": {"p_nom": {"wind": 1}}}
n.optimize.optimize_mga(slack=0.05, weights=weights, sense="min")

weights = {"Generator": {"p_nom": {"wind": 1}}} n.optimize.optimize_mga(slack=0.05, weights=weights, sense="min")

/opt/hostedtoolcache/Python/3.13.15/x64/lib/python3.13/contextlib.py:85: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  return func(*args, **kwds)

INFO:linopy.model: Solve problem using Gurobi solver

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

INFO:linopy.io:Writing objective.

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

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

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

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

INFO:linopy.io: Writing time: 0.13s

Set parameter WLSAccessID

Set parameter WLSSecret

Set parameter LicenseID to value 2537914

Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de

Read LP format model from file /tmp/linopy-problem-i0p2_zar.lp

Reading time = 0.07 seconds

obj: 64247 rows, 29206 columns, 127070 nonzeros

Set parameter LogToConsole to value 0

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 29206 primals, 64247 duals
Objective: 2.13e+04
Solver: gurobi
Runtime: 1.63s
Dual bound: 2.13e+04
Solver model: available
Solver message: 2

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-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, Store-energy_balance, budget were not assigned to the network.

Out[8]:

('ok', 'optimal')

The breakdown of total system costs shifts from wind towards more solar.

In [9]:

tsc = (
    pd.concat([n.statistics.capex(), n.statistics.opex()], axis=1).sum(axis=1).div(1e9)
)
tsc

tsc = ( pd.concat([n.statistics.capex(), n.statistics.opex()], axis=1).sum(axis=1).div(1e9) ) tsc

Out[9]:

component    carrier         
Generator    solar               2.133254
             wind                2.168509
Link         electrolysis        0.539322
             turbine             1.191285
StorageUnit  battery storage     1.402346
Store        hydrogen storage    0.903488
Generator    load shedding       0.143837
dtype: float64

Up to numeric differences, it is 5% more expensive overall:

In [10]:

optimal_cost * 1.05

optimal_cost * 1.05

Out[10]:

np.float64(8.48204245922375)

In [11]:

tsc.sum()

tsc.sum()

Out[11]:

np.float64(8.482042459223749)

Overall, the wind capacity is cut by roughly a third:

In [12]:

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

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

Out[12]:

component    carrier         
Generator    load shedding         10.901160
             solar                 41.545979
             wind                  21.334328
Link         electrolysis           2.857854
             turbine               10.235553
StorageUnit  battery storage       22.132377
Store        hydrogen storage    6091.522243
dtype: float64

This is also evident in the energy balance:

In [13]:

n.statistics.energy_balance(bus_carrier="electricity").sort_values().div(1e6)

n.statistics.energy_balance(bus_carrier="electricity").sort_values().div(1e6)

Out[13]:

component    carrier          bus_carrier
Load         -                electricity   -66.266089
Link         electrolysis     electricity   -14.552445
StorageUnit  battery storage  electricity    -1.480815
Generator    load shedding    electricity     0.071919
Link         turbine          electricity     3.709375
Generator    wind             electricity    33.603394
             solar            electricity    44.914663
dtype: float64

And also recognizable in the energy balance time series plots:

In [14]:

n.statistics.energy_balance.plot.area(linewidth=0, bus_carrier="electricity")

n.statistics.energy_balance.plot.area(linewidth=0, bus_carrier="electricity")

Out[14]:

(<Figure size 759x300 with 1 Axes>,
 <Axes: xlabel='snapshot', ylabel='Energy Balance []'>,
 <seaborn.axisgrid.FacetGrid at 0x7f56007e3890>)