We have been iteratively coming back to discussion that Stochastic Programming (SP) is one of the most wished features.
There is an example notebook with SP in PyPSA; however, that formulation is badly scalable. A stochastic investment planning feature in PyPSA itself would likely find many users.

Below is the a quick sketch how to implement SP with Linopy integration module in PyPSA. This would allow for modelling two-stage SP with investment decisions in the first stage, and operational decision in the second stage. There is a ton of energy-related literature with SP implementation, a really good example is here.

Eventually, a user should be able to define (i) list of scenarios to model, (ii) realization probabilities of defined scenarios, and (iii) a vector (or a matrix) of uncertain parameter(s). The latter could be a cost coefficient or a rhs parameter (like capacity constraint).

1. First, we add functionality to pypsa/optimization/optimize.py:

• create_model() receives new argument stochastic (boolean).
• If False, the function stays as is.
• If True:
  • all operational variables get new coordinate scen
  • all constraints get new coordinate scen
  • the global constraints likely stay same

2. define_objective() inherits argument stochastic. If True, the changes in function are so that line objective.append((operation * cost).sum()) can eventually be changed to objective.append((operation[scenario] * cost[scenario]* probability[scenario] )).sum()) The investment part stays same.

3. post_processing() inherits argument stochastic. If True, the p0 and p1 are defined per scen

4. assign_duals() inherits argument stochastic. if True, the duals are scaled back by 1/probability.

5. User has to pass three bits of information to model (via **kwargs ?)

• a list of scenario names -> this is read by create_model() to populate scen
• one vector (like pandas series) of scenario probabilities (or dictionary with scenario: probability pairs) -> this is passed to define_objective()
• a vector (e.g., CO2 prices) or matrix (e.g., time-series) of uncertain parameter
• Other parameters that are NOT changed by user-defined vector/matrix are populated with the same values across scenarios.

A more general form covering e.g. unit-commitment decisions under RES feed-in uncertainty is harder to generalise via pypsa/optimization/optimize.py; however, structurally the problem would be of the same form.