Exploring the Near-Optimal Feasible Space¶

Here, we use modelling to generate alternatives (MGA) to explore the near-opimal feasible space of a simple one-node model in the style of model.energy (the same as in the previous MGA example).

While we explored the use of n.optimize.optimize_mga() in the previous example, PyPSA includes a few additional functions making it convenient to work with near-optimal spaces in a geometric way. While the optimize_mga function solves a network with an alternative objective function (given by the weights argument), n.optimize_mga_in_direction is a simple alternative useful for exploring trade-offs between multiple alternative objectives. It takes two key arguments: a dictionary of dimensions, specifying multiple alternative objectives in the same format as the previously mentioned weights, and a direction representing a vector in the coordinate space defined by dimensions.

In this example, we explore the trade-offs between wind and solar expansion in a simple renewables-based system. Therefore, we define our dimensions (or alternative objectives) as total installed wind capacity and total installed solar capacity in MW, respectively. These dimensions can also be seen as defining a projection of near-optimal space of the model down to two dimensions.

:::{note} See also this research article for a background on dimension reduction and approximation of near-optimal spaces. :::

In [2]:

Copied!

import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from scipy.spatial import ConvexHull, convex_hull_plot_2d

import pypsa
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from scipy.spatial import ConvexHull, convex_hull_plot_2d

import pypsa

Before any near-optimal analysis can start, the model has to be solved to optimality first. This is in order to find the minimum system cost, which is later used to define a system cost bound for near-optimality (i.e. "slack", usually in percentage of minimum system cost).

In [3]:

Copied!





n = pypsa.examples.model_energy()
n.snapshots = n.snapshots[::3]  # Crude reduction of time resolution for speed
n.generators.at["load shedding", "marginal_cost"] = 20000
n.optimize(solver_name="gurobi")
n = pypsa.examples.model_energy()
n.snapshots = n.snapshots[::3]  # Crude reduction of time resolution for speed
n.generators.at["load shedding", "marginal_cost"] = 20000
n.optimize(solver_name="gurobi")

INFO:pypsa.network.io:Retrieving network data from https://github.com/PyPSA/PyPSA/raw/v1.1.2/examples/networks/model-energy/model-energy.nc.

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

/tmp/ipykernel_5184/2365194803.py:4: 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.

INFO:linopy.model: Solve problem using Gurobi solver

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

INFO:linopy.io: Writing time: 0.13s

Set parameter WLSAccessID

Set parameter WLSSecret

Set parameter LicenseID to value 2766203

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

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

Reading time = 0.02 seconds

obj: 21434 rows, 9746 columns, 41420 nonzeros

Set parameter LogToConsole to value 0

Warning: environment still referenced so free is deferred (Continue to use WLS)

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 9746 primals, 21434 duals
Objective: 6.90e+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]:

('ok', 'optimal')

Individual dimensions are specified just like the weights argument for the n.optimize.optimize_mga function: a nested dictionary of components, variables and finally a series or dictionary giving weights to individual components.

Here, we define two dimensions, and give them the user-defined names "wind" and "solar" (the keys in the outer dictionary). These names are arbitrary and don't have to match names of underlaying components.

In [4]:

Copied!





dimensions = {
    "wind": {"Generator": {"p_nom": {"wind": 1}}},
    "solar": {"Generator": {"p_nom": {"solar": 1}}},
}
dimensions = {
    "wind": {"Generator": {"p_nom": {"wind": 1}}},
    "solar": {"Generator": {"p_nom": {"solar": 1}}},
}

Suppose we are interested in minimizing wind capacity alone, and also interested in minimizing the total capacity of wind and solar. These correspond to the directions (-1, 0) and (-1, -1) in wind-solar coordinate space, respectively. Let's see what is possible within a 5% total system cost slack.

In [5]:

Copied!





_, _, min_wind = n.optimize.optimize_mga_in_direction(
    dimensions=dimensions,
    direction={
        "wind": -1,
        "solar": 0,
    },  # Use coordinate names matching those given in `dimensions`
    slack=0.05,
)
_, _, min_wind_solar = n.optimize.optimize_mga_in_direction(
    dimensions=dimensions,
    direction={"wind": -1, "solar": -1},
    slack=0.05,
)

display(min_wind)
display(min_wind_solar)
_, _, min_wind = n.optimize.optimize_mga_in_direction(
    dimensions=dimensions,
    direction={
        "wind": -1,
        "solar": 0,
    },  # Use coordinate names matching those given in `dimensions`
    slack=0.05,
)
_, _, min_wind_solar = n.optimize.optimize_mga_in_direction(
    dimensions=dimensions,
    direction={"wind": -1, "solar": -1},
    slack=0.05,
)

display(min_wind)
display(min_wind_solar)

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/checkouts/v1.1.2/pypsa/optimization/mga.py:502: 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.

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io: Writing time: 0.13s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 9746 primals, 21435 duals
Objective: 1.70e+04
Solver model: available
Solver message: Optimal

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.

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/checkouts/v1.1.2/pypsa/optimization/mga.py:502: 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.

INFO:linopy.model: Solve problem using Highs solver

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

INFO:linopy.io: Writing time: 0.13s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 9746 primals, 21435 duals
Objective: 4.75e+04
Solver model: available
Solver message: Optimal

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.

wind     17006.298742
solar    44921.817858
dtype: float64

wind     29563.944662
solar    17925.539753
dtype: float64

Results are returned as dictionaries represting coordinates in the wind-solar space. The dictionaries use the same user-given coordinate names that we passed in the dimensions argument. We see that minimizing wind alone leads to a large installation of solar, while minimizing both jointly leads to a more balanced solution.

In order to get a full picture of the near-optimal space (projected to the wind and solar dimensions), it is useful to optimize in many directions and consider the convex hull of the resulting points.

Of course, you can use optimize_mga_in_direction manually to do this. For convenience, PyPSA also provides the optimize_mga_in_multiple_directions function. This function takes a list of directions (or a pandas DataFrame with rows for directions), and uses built-in Python multiprocessing functionality to optimize in these directions in parallel. In the event of any individual optimisations failing, only the successful ones are returned.

While the user is free to select directions themselves, PyPSA also includes some simple functions (for instance, pypsa.optimization.mga.generate_directions_randomber of random directions in the format expected by optimize_mga_in_multiple_directions`.

In [6]:

Copied!





directions = pypsa.optimization.mga.generate_directions_random(
    dimensions.keys(), 10, seed=0
)
dirs, pts = n.optimize.optimize_mga_in_multiple_directions(
    dimensions=dimensions,
    directions=directions,
    slack=0.05,
    solver_name="gurobi",
    # Solve in up to 8 directions in parallel.
    max_parallel=8,
)
directions = pypsa.optimization.mga.generate_directions_random(
    dimensions.keys(), 10, seed=0
)
dirs, pts = n.optimize.optimize_mga_in_multiple_directions(
    dimensions=dimensions,
    directions=directions,
    slack=0.05,
    solver_name="gurobi",
    # Solve in up to 8 directions in parallel.
    max_parallel=8,
)

INFO:pypsa.network.io:Exported network 'Model-Energy' saved to '/tmp/tmpth2ijiew.nc contains: storage_units, links, sub_networks, buses, stores, carriers, generators, loads

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/checkouts/v1.1.2/pypsa/optimization/mga.py:502: 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.

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/checkouts/v1.1.2/pypsa/optimization/mga.py:502: 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.