Implement stochastic programming formulation in PyPSA

This is a draft.

This PR will:
-- add an option to cast a two-stage stochastic programming formulation with pypsa out-of-the-box.
-- add related documentation
-- add related notebook examples

This PR will close #575
Eventually, it should be possible to replicate this stochastic problem out-of-the-box with a few lines of code: https://pypsa.readthedocs.io/en/latest/examples/stochastic-problem.html

Update

there is an ongoing restructuring process in PyPSA with the introduction of the component class (#1075) and a component xarray view (#1076) and the coming update of the optimization module (all optimization parameters will be handled as xarray objects instead of dataframes. PR to be created soon).
This pipeline, once complete, will simplify the handling of stochastic programming in PyPSA by eliminating the need for case-specific equation handling, the StochasticNetwork class, and update of all relevant descriptors, as this PR suggests.
We complete the three PRs mentioned above and revise this PR accordingly.

work log

• introduce a StochasticNetwork class
• reindex snapshots with scenario names (SP: branches)
• update snapshot_weightings with scenario probabilities
• reindex time series parameters and variables (!) with scenario names (SP: branches)
• introduce scenario-dependent static components.
• add get_switchable_by_scenario in descriptors.py
• add get_switchable_as_dense_by_scenario in descriptors.py
• add scenario dim to operational variables (variables.py)
• broadcast static data w/o _s components
• make sure __repr__ and iterate_components are compatible
• fancy printout
• --------------------------------------------------------------------------------
• move new descriptors into StochNetwork class
• make sure get_switchable_as_dense is compatible for varying and static cases in SP
• make sure get_bounds_pu is compatible with SP
• make sure get_bounds_pu returns right xarray object for equations to work out-of-the-box
• make n.df(c) to be multi-index dataframe instead of a dict. Make everything above compatible.
• fix specific case of StorageUnit - state_of_charge - in get_bounds_pu for SP
• rework get_switchable_by_scenario for multi-index n.df(c)
• bring back explicit "scenario" dim to operational variables
• add explicit scenario dim to operational constraints
• compatible nodal balance constraints
• --------------------------------------------------------------------------------
• make objective compatible
• add dedicated scenario_weightings into objective
• optimize.create_model() correctly casts for a toy stochastic problem
• n.optimize() call works for a toy stochastic problem
• make required descriptors compatible (to be continued)
• --------------------------------------------------------------------------------
• consistency_check is compatible
• expand implementation to all types of variables
• expand implementation to all constraints
• comprehensive tests
• make sure it works with investment_periods
• correct writing of results back to Network
• all Network methods are be compatible
• --------------------------------------------------------------------------------
• n.optimize() call works for any StochasticNetwork
• REVERTED: remove explicit "scenario" dim from vars, broadcast it via snapshots
• CANCELLED: introduce optional argument stochastic_problem (True or False) into n.optimize() with default False

What it should do

The linear cost-minimization problem behind PyPSA's LOPF is deterministic and can be represented as:

$$ DP := \min_{x} c^T x $$

$$ s.t. Ax \geq b, \quad x \geq 0 $$

where x: a vector of variables; c and b: vectors of parameters; A: a matrix of parameters. The inequalities specify a convex polytope that constrains the objective function. The optimal deterministic solution is to find a vector x that minimises the objective function given a set of constraints.

This PR adds a possibility to cast a two-stage risk-neutral stochastic linear optimization problem that uses a set of scenarios provided by a modeller to describe the uncertain parameters. The optimal first-stage decisions must be made before the information about uncertain parameters is revealed, while the second-stage decisions are made after it is revealed:

$$ SP := \min_{x} c^T x + \mathbb{E}_s[Q(x,s)] $$

$$ \text{s.t.} \quad Ax \geq b, \quad x \geq 0 $$

$$ \text{where} \quad Q(x,s) := \min_{y} q_s^T y $$

$$ \text{s.t.} \quad T_s x + W_s y \geq d_s, \quad y \geq 0 $$

where x: first-stage variables; y: second-stage variables; s: vector of uncertain data (the discrete set of scenarios that include a finite number of scenarios). The first-stage parameters (c, A, b) are assumed to be known with certainty. The second-stage decisions are restricted by the first-stage decisions x; the parameters (q, T, W, d) are actual realizations of uncertain data (the solution must be feasible for every scenario s ∈ S). The first-stage problem seeks to minimize the sum of the costs of the first-stage decisions and the expected costs of the second-stage decisions. The second-stage problem seeks to minimize the second-stage costs.

How it works

Basic PyPSA network:

n = pypsa.Network(snapshots=pd.date_range("2024-01-01", periods=5, freq="H"))
n.add("Bus", "UA")
n.add("Generator", "nuclear", bus="UA", marginal_cost=10, p_nom=100)
n.add("Generator", "gas", bus="UA", marginal_cost=70, p_nom=50)
load_data = [80, 90, 110, 120, 100]
load_series = pd.Series(load_data, n.snapshots)
n.add("Load", "demand", bus="UA", p_set=load_series)

Network and Components

User specifies a dict with scenario names and probabilities:
scenarios = {"low": 0.33, "medium": 0.33, "high": 0.34}

This call converts the network, i.e., all static and time-dependent components are broadcasted over the provided list of scenarios.
n = StochasticNetwork(n, scenarios)

> PyPSA StochasticNetwork 🔀
> Components:
>  - Bus: 1
>  - Generator: 6
>  - Load: 3
> Snapshots: 15
> Scenarios: [low, medium, high]

Time-dependent components get new index "scenario" -- thus we use here a similar structure to the case when the network has "investment_periods".

> MultiIndex([(   'low', '2024-01-01 00:00:00'),
>             (   'low', '2024-01-01 01:00:00'),
>              ...
>             (  'high', '2024-01-01 03:00:00'),
>             (  'high', '2024-01-01 04:00:00')],
>            names=['scenario', 'snapshot'])

Static components are now a dictionary containing scenarios (keys) and component dataframes (values) . This will make possible to pass any static components' attribute as uncertainty vector.

A new function for retrieving scenario-dependent static component data:
get_switchable_by_scenario(n, "Generator", "p_nom")

In stochastic network, get_switchable_as_dense is redefined to work out-of-the-box. It can be called from network instance:
sn.get_switchable_as_dense("Generator", "p_max_pu")

Same applies for other descriptors, e.g., a call sn.get_bounds_pu("StorageUnit", sn.snapshots, None, "p_store") would return bounds per scenario.

Variables

Operational variables get new dimension "scenario". Thus for e.g. n.model.variables["Generator-p"]:

Variable (scenario: 3, snapshot: 5, Generator: 2)
-------------------------------------------------
[low, 2024-01-01 00:00:00, nuclear]: Generator-p[low, 2024-01-01 00:00:00, nuclear] ∈ [-inf, inf]
[low, 2024-01-01 00:00:00, gas]: Generator-p[low, 2024-01-01 00:00:00, gas] ∈ [-inf, inf]
		...
[high, 2024-01-01 04:00:00, nuclear]: Generator-p[high, 2024-01-01 04:00:00, nuclear] ∈ [-inf, inf]
[high, 2024-01-01 04:00:00, gas]: Generator-p[high, 2024-01-01 04:00:00, gas] ∈ [-inf, inf]

NB: We thought about keeping things aligned with "investment_periods", where periods are NOT explicitly added to variable dims, but rather broadcasted via snapshots, which contain information about period and timestep. After implementing it, we realized this approach is not a good idea for scenarios since we have combination of static and time-dependent attr. This could work -- but requires an update in every constraint. Explicit scenario dim keeps chaos under control.

Constraints

n.optimize()

PR Checklist

• Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in doc.
• Unit tests for new features were added (if applicable).
• Newly introduced dependencies are added to environment.yaml, environment_docs.yaml and setup.py (if applicable).
• A note for the release notes doc/release_notes.rst of the upcoming release is included.
• I consent to the release of this PR's code under the MIT license.