Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Implementing spinning reserve constraints

Objective

When using a linear programming approach in an energy system, the optimization usually results in all generators either running at full capacity or not running at all, depending on their position in the merit order. In the real world, generators often run at partial load. Unfortunately, most of the reasons why generators run at partial load are difficult to account for in modeling. One reason for generators to run at partial load is to be able to act as spinning reserves, i.e. to be able to ramp up or down quickly when needed to maintain system stability.

In this example we will implement spinning reserve constraints in PyPSA in a very simplified way. The objective is to force some generators to provide reserve power by running below their rated capacity, but still maintain a linear programming problem formulation.

To do this, we need to implement additional variables and additional constraints in the model.

Methodology

We follow the approach presented by Andreas Hösl et al here: https://www.youtube.com/watch?v=fmwDxNpSMM4&t=8043s

The basic idea is that each generator must provide reserve power symmetrically. This means that it must be able to increase and decrease its output by the same amount in order to contribute to meeting reserve requirements. This ensures that generators must operate at partial load to provide reserve power.

The following changes need to be made to the linopy model in a PyPSA network:

• a new variable \(p_{\text{reserve}}(g,t)\) representing the reserve power provided by generator \(g\) at time step \(t\).

• a constraint that ensures that for each time step \(t\), the sum of all reserve power provided is greater than or equal to the required reserves.

  \[\forall t: \sum_{g} p_{\text{reserve}}(g,t) \geq \text{reserve requirement}\]

• A constraint to ensure that the reserve power of a generator is less than or equal to the difference between its power \(p\) and its nominal power \(p_\text{nom}\), multiplied by a scalar coefficient \(a\). This coefficient can take any value between 0 and 1 and represents the technical availability of a generator to provide reserve power.

  \[\forall g, t: p_\text{reserve}(g, t) \leq a(g) p_\text{nom}(g) - p(g,t)\]

• a constraint to ensure that the reserve power of a generator is less than or equal to its actual power \(p\) multiplied by a scalar coefficient \(b\). This coefficient can take any value between 0 and 1 and represents the technical availability of a generator to provide reserve power.

\[\forall g, t: p_\text{reserve}(g, t) \leq b(g) p(g,t)\]

The relationships between the variables \(a\), \(b\), \(p_\text{nom}\), \(p\) and \(p_\text{reserve}\) are shown in the following schematic diagram.

Limitations and other approaches

Note that this is an oversimplified approach that has significant limitations:

• It does not distinguish between different categories of reserves, such as primary or secondary reserves.

• Reserves are provided symmetrically; there is no distinction between positive and negative reserves.

• The approach only considers the provision of reserve power, not the actual delivery. The additional constraints simply force some generators to run at partial load so that they could ramp up or down when reserves are required.

All of these issues can be addressed in a MIP unit commitment model, albeit at a much higher numerical cost.

Implementation

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import pypsa

Basic model

Our toy model consists of a single bus with generators that have different marginal costs. We use a sine function for the load profile.

n_basic = pypsa.Network()
n_basic.add("Carrier", name="carrier1")
n_basic.add("Bus", name="bus1", carrier="carrier1")

# add generators with increasing marginal cost
n_basic.add("Generator", name="gen1", bus="bus1", p_nom=10, marginal_cost=1)
n_basic.add("Generator", name="gen2", bus="bus1", p_nom=10, marginal_cost=2)
n_basic.add("Generator", name="gen3", bus="bus1", p_nom=10, marginal_cost=3)
n_basic.add("Generator", name="gen4", bus="bus1", p_nom=10, marginal_cost=4)

# create 48 snapshots
snapshots = np.arange(1, 49)
n_basic.set_snapshots(snapshots)

# create load
load_max = 30
load_profile = np.sin(snapshots / 12 * np.pi) + 3.5
load_profile = load_profile / load_profile.max() * load_max
n_basic.add("Load", name="load1", bus="bus1", p_set=load_profile)

Index(['load1'], dtype='object')

We make a copy of the basic model, which we later modify to add reserve power constraints.

n_reserve = n_basic.copy()

As a reference point, we solve the model without any additional constraints.

n_basic.optimize()

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io: Writing time: 0.02s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 192 primals, 432 duals
Objective: 1.95e+03
Solver model: available
Solver message: optimal

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

Running HiGHS 1.10.0 (git hash: fd86653): Copyright (c) 2025 HiGHS under MIT licence terms
LP   linopy-problem-97uzb6l9 has 432 rows; 192 cols; 576 nonzeros
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [1e+00, 4e+00]
  Bound  [0e+00, 0e+00]
  RHS    [1e+01, 3e+01]
Presolving model
48 rows, 192 cols, 192 nonzeros  0s
Dependent equations search running on 46 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
46 rows, 184 cols, 184 nonzeros  0s
Presolve : Reductions: rows 46(-386); columns 184(-8); elements 184(-392)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     1.2000000000e+02 Pr: 46(1060) 0s
         46     1.9543748803e+03 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
         46     1.8343748803e+03 Pr: 2(60) 0s
         52     1.9543748803e+03 Pr: 0(0) 0s
Required 6 simplex iterations after postsolve
Model name          : linopy-problem-97uzb6l9
Model status        : Optimal
Simplex   iterations: 52
Objective value     :  1.9543748803e+03
Relative P-D gap    :  1.1634087080e-16
HiGHS run time      :          0.00
Writing the solution to /tmp/linopy-solve-w716fyiq.sol

('ok', 'optimal')

We plot the dispatch over time. As expected, the generators are dispatched strictly according to their marginal cost, each one running at nominal capacity until demand is met.

n_basic.generators_t["p"].plot.area(lw=0).legend(
    loc="upper left", bbox_to_anchor=(1.0, 1.0)
)
plt.show()

Modified model with custum variable and constraints

Now let’s modify the model by adding some additional constraints. We create a new network and create a model instance attached to it. Now we can inspect the model instance to get a list of variables and constraints:

n_reserve.optimize.create_model()
n_reserve.model

Linopy LP model
===============

Variables:
----------
 * Generator-p (snapshot, Generator)

Constraints:
------------
 * Generator-fix-p-lower (snapshot, Generator-fix)
 * Generator-fix-p-upper (snapshot, Generator-fix)
 * Bus-nodal_balance (Bus, snapshot)

Status:
-------
initialized

We now add a new variable p_reserve which represents the reserve power. It has a lower bound of zero, is defined for all dispatchable generators and has a time index.

v_rp = n_reserve.model.add_variables(
    lower=0,
    coords=[n_reserve.snapshots, n_reserve.generators.index],
    name="Generator-p_reserve",
)
v_rp

Variable (snapshot: 48, Generator: 4)
-------------------------------------
[1, gen1]: Generator-p_reserve[1, gen1] ∈ [0, inf]
[1, gen2]: Generator-p_reserve[1, gen2] ∈ [0, inf]
[1, gen3]: Generator-p_reserve[1, gen3] ∈ [0, inf]
[1, gen4]: Generator-p_reserve[1, gen4] ∈ [0, inf]
[2, gen1]: Generator-p_reserve[2, gen1] ∈ [0, inf]
[2, gen2]: Generator-p_reserve[2, gen2] ∈ [0, inf]
[2, gen3]: Generator-p_reserve[2, gen3] ∈ [0, inf]
                ...
[47, gen2]: Generator-p_reserve[47, gen2] ∈ [0, inf]
[47, gen3]: Generator-p_reserve[47, gen3] ∈ [0, inf]
[47, gen4]: Generator-p_reserve[47, gen4] ∈ [0, inf]
[48, gen1]: Generator-p_reserve[48, gen1] ∈ [0, inf]
[48, gen2]: Generator-p_reserve[48, gen2] ∈ [0, inf]
[48, gen3]: Generator-p_reserve[48, gen3] ∈ [0, inf]
[48, gen4]: Generator-p_reserve[48, gen4] ∈ [0, inf]

Next, we define a new constraint which ensures that for each snapshot the total reserve requirement is satisfied by the sum of the reserve power provided by all generators.

reserve_req = 10

c_sum = n_reserve.model.add_constraints(
    v_rp.sum("Generator") >= reserve_req, name="GlobalConstraint-sum_of_reserves"
)
c_sum

Constraint `GlobalConstraint-sum_of_reserves` [snapshot: 48]:
-------------------------------------------------------------
[1]: +1 Generator-p_reserve[1, gen1] + 1 Generator-p_reserve[1, gen2] + 1 Generator-p_reserve[1, gen3] + 1 Generator-p_reserve[1, gen4]      ≥ 10.0
[2]: +1 Generator-p_reserve[2, gen1] + 1 Generator-p_reserve[2, gen2] + 1 Generator-p_reserve[2, gen3] + 1 Generator-p_reserve[2, gen4]      ≥ 10.0
[3]: +1 Generator-p_reserve[3, gen1] + 1 Generator-p_reserve[3, gen2] + 1 Generator-p_reserve[3, gen3] + 1 Generator-p_reserve[3, gen4]      ≥ 10.0
[4]: +1 Generator-p_reserve[4, gen1] + 1 Generator-p_reserve[4, gen2] + 1 Generator-p_reserve[4, gen3] + 1 Generator-p_reserve[4, gen4]      ≥ 10.0
[5]: +1 Generator-p_reserve[5, gen1] + 1 Generator-p_reserve[5, gen2] + 1 Generator-p_reserve[5, gen3] + 1 Generator-p_reserve[5, gen4]      ≥ 10.0
[6]: +1 Generator-p_reserve[6, gen1] + 1 Generator-p_reserve[6, gen2] + 1 Generator-p_reserve[6, gen3] + 1 Generator-p_reserve[6, gen4]      ≥ 10.0
[7]: +1 Generator-p_reserve[7, gen1] + 1 Generator-p_reserve[7, gen2] + 1 Generator-p_reserve[7, gen3] + 1 Generator-p_reserve[7, gen4]      ≥ 10.0
                ...
[42]: +1 Generator-p_reserve[42, gen1] + 1 Generator-p_reserve[42, gen2] + 1 Generator-p_reserve[42, gen3] + 1 Generator-p_reserve[42, gen4] ≥ 10.0
[43]: +1 Generator-p_reserve[43, gen1] + 1 Generator-p_reserve[43, gen2] + 1 Generator-p_reserve[43, gen3] + 1 Generator-p_reserve[43, gen4] ≥ 10.0
[44]: +1 Generator-p_reserve[44, gen1] + 1 Generator-p_reserve[44, gen2] + 1 Generator-p_reserve[44, gen3] + 1 Generator-p_reserve[44, gen4] ≥ 10.0
[45]: +1 Generator-p_reserve[45, gen1] + 1 Generator-p_reserve[45, gen2] + 1 Generator-p_reserve[45, gen3] + 1 Generator-p_reserve[45, gen4] ≥ 10.0
[46]: +1 Generator-p_reserve[46, gen1] + 1 Generator-p_reserve[46, gen2] + 1 Generator-p_reserve[46, gen3] + 1 Generator-p_reserve[46, gen4] ≥ 10.0
[47]: +1 Generator-p_reserve[47, gen1] + 1 Generator-p_reserve[47, gen2] + 1 Generator-p_reserve[47, gen3] + 1 Generator-p_reserve[47, gen4] ≥ 10.0
[48]: +1 Generator-p_reserve[48, gen1] + 1 Generator-p_reserve[48, gen2] + 1 Generator-p_reserve[48, gen3] + 1 Generator-p_reserve[48, gen4] ≥ 10.0

Now we need to limit the amount of reserve power that each generator can provide. The following constraint ensures that the reserve power provided by a generator must be less than or equal to the difference between its power p and its nominal power p_nom:

a = 1

c_rpos = n_reserve.model.add_constraints(
    v_rp
    <= -n_reserve.model.variables["Generator-p"] + a * n_reserve.generators["p_nom"],
    name="Generator-reserve_upper_limit",
)
c_rpos

Constraint `Generator-reserve_upper_limit` [snapshot: 48, Generator: 4]:
------------------------------------------------------------------------
[1, gen1]: +1 Generator-p_reserve[1, gen1] + 1 Generator-p[1, gen1]    ≤ 10.0
[1, gen2]: +1 Generator-p_reserve[1, gen2] + 1 Generator-p[1, gen2]    ≤ 10.0
[1, gen3]: +1 Generator-p_reserve[1, gen3] + 1 Generator-p[1, gen3]    ≤ 10.0
[1, gen4]: +1 Generator-p_reserve[1, gen4] + 1 Generator-p[1, gen4]    ≤ 10.0
[2, gen1]: +1 Generator-p_reserve[2, gen1] + 1 Generator-p[2, gen1]    ≤ 10.0
[2, gen2]: +1 Generator-p_reserve[2, gen2] + 1 Generator-p[2, gen2]    ≤ 10.0
[2, gen3]: +1 Generator-p_reserve[2, gen3] + 1 Generator-p[2, gen3]    ≤ 10.0
                ...
[47, gen2]: +1 Generator-p_reserve[47, gen2] + 1 Generator-p[47, gen2] ≤ 10.0
[47, gen3]: +1 Generator-p_reserve[47, gen3] + 1 Generator-p[47, gen3] ≤ 10.0
[47, gen4]: +1 Generator-p_reserve[47, gen4] + 1 Generator-p[47, gen4] ≤ 10.0
[48, gen1]: +1 Generator-p_reserve[48, gen1] + 1 Generator-p[48, gen1] ≤ 10.0
[48, gen2]: +1 Generator-p_reserve[48, gen2] + 1 Generator-p[48, gen2] ≤ 10.0
[48, gen3]: +1 Generator-p_reserve[48, gen3] + 1 Generator-p[48, gen3] ≤ 10.0
[48, gen4]: +1 Generator-p_reserve[48, gen4] + 1 Generator-p[48, gen4] ≤ 10.0

Finally, we add a constraint to ensure that the reserve power provided by a generator must be less than or equal to its actual power p multiplied by a scalar coefficient b. This coefficient can take any value between 0 and 1 and represents the technical availability of a generator to provide reserve power.

b = 0.7

c_rneg = n_reserve.model.add_constraints(
    v_rp <= b * n_reserve.model.variables["Generator-p"],
    name="Generator-reserve_lower_limit",
)
c_rneg

Constraint `Generator-reserve_lower_limit` [snapshot: 48, Generator: 4]:
------------------------------------------------------------------------
[1, gen1]: +1 Generator-p_reserve[1, gen1] - 0.7 Generator-p[1, gen1]    ≤ -0.0
[1, gen2]: +1 Generator-p_reserve[1, gen2] - 0.7 Generator-p[1, gen2]    ≤ -0.0
[1, gen3]: +1 Generator-p_reserve[1, gen3] - 0.7 Generator-p[1, gen3]    ≤ -0.0
[1, gen4]: +1 Generator-p_reserve[1, gen4] - 0.7 Generator-p[1, gen4]    ≤ -0.0
[2, gen1]: +1 Generator-p_reserve[2, gen1] - 0.7 Generator-p[2, gen1]    ≤ -0.0
[2, gen2]: +1 Generator-p_reserve[2, gen2] - 0.7 Generator-p[2, gen2]    ≤ -0.0
[2, gen3]: +1 Generator-p_reserve[2, gen3] - 0.7 Generator-p[2, gen3]    ≤ -0.0
                ...
[47, gen2]: +1 Generator-p_reserve[47, gen2] - 0.7 Generator-p[47, gen2] ≤ -0.0
[47, gen3]: +1 Generator-p_reserve[47, gen3] - 0.7 Generator-p[47, gen3] ≤ -0.0
[47, gen4]: +1 Generator-p_reserve[47, gen4] - 0.7 Generator-p[47, gen4] ≤ -0.0
[48, gen1]: +1 Generator-p_reserve[48, gen1] - 0.7 Generator-p[48, gen1] ≤ -0.0
[48, gen2]: +1 Generator-p_reserve[48, gen2] - 0.7 Generator-p[48, gen2] ≤ -0.0
[48, gen3]: +1 Generator-p_reserve[48, gen3] - 0.7 Generator-p[48, gen3] ≤ -0.0
[48, gen4]: +1 Generator-p_reserve[48, gen4] - 0.7 Generator-p[48, gen4] ≤ -0.0

We can now inspect the model formulation. We can see that our new variables and constraints have been successfully added:

n_reserve.model

Linopy LP model
===============

Variables:
----------
 * Generator-p (snapshot, Generator)
 * Generator-p_reserve (snapshot, Generator)

Constraints:
------------
 * Generator-fix-p-lower (snapshot, Generator-fix)
 * Generator-fix-p-upper (snapshot, Generator-fix)
 * Bus-nodal_balance (Bus, snapshot)
 * GlobalConstraint-sum_of_reserves (snapshot)
 * Generator-reserve_upper_limit (snapshot, Generator)
 * Generator-reserve_lower_limit (snapshot, Generator)

Status:
-------
initialized

We can now solve the modified model:

n_reserve.optimize.solve_model()

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io: Writing time: 0.03s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 384 primals, 864 duals
Objective: 2.30e+03
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, GlobalConstraint-sum_of_reserves, Generator-reserve_upper_limit, Generator-reserve_lower_limit were not assigned to the network.

Running HiGHS 1.10.0 (git hash: fd86653): Copyright (c) 2025 HiGHS under MIT licence terms
LP   linopy-problem-536bw98k has 864 rows; 384 cols; 1536 nonzeros
Coefficient ranges:
  Matrix [7e-01, 1e+00]
  Cost   [1e+00, 4e+00]
  Bound  [0e+00, 0e+00]
  RHS    [1e+01, 3e+01]
Presolving model
480 rows, 384 cols, 1152 nonzeros  0s
432 rows, 336 cols, 1152 nonzeros  0s
Presolve : Reductions: rows 432(-432); columns 336(-48); elements 1152(-384)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Ph1: 0(0) 0s
        296     2.3027070110e+03 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-536bw98k
Model status        : Optimal
Simplex   iterations: 296
Objective value     :  2.3027070110e+03
Relative P-D gap    :  6.3194818789e-15
HiGHS run time      :          0.00
Writing the solution to /tmp/linopy-solve-9p3brx5m.sol

('ok', 'optimal')

Examine the results

We create a plot to examine the results of the modified model. On the left subplot we plot the active power generation p over time for each generator. On the right subplot we plot the reserve power p_reserve over time for each generator. The following observations can be made:

• The reserve requirement of 20 MW is met in every time step.

• In order to provide reserves, some generators must always run below their nominal power.

• Among all running generators, those with the highest marginal costs provide as much reserve capacity as possible.

fig, axs = plt.subplots(1, 2, sharey=True, figsize=(10, 5))
n_reserve.generators_t["p"].plot.area(
    ax=axs[0], title="p", legend=False, ylabel="p [MW]"
)
n_reserve.generators_t["p_reserve"].plot.area(ax=axs[1], title="p_reserve")
plt.tight_layout()
plt.show()

Looking at the average reserve power provided by a generator, we can see that the cheapest and most expensive generators provide less reserve power on average than the other two generators.

n_reserve.generators_t["p_reserve"].mean().plot(
    kind="bar", ylabel="mean(p_reserve) [MW]"
)
plt.show()

Comparison of the two model versions

To visually compare the base model with the modified model, we plot the active power generation p over time for both models. On the left we plot the base model, on the right we plot the model with additional constraints.

fig, axs = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(10, 4))
for i, (n, r) in enumerate([(n_basic, 0), (n_reserve, reserve_req)]):
    n.generators_t["p"].plot.area(
        ax=axs[i], ylabel="p [MW]", title=f"{r} MW reserve required", legend=False, lw=0
    )

plt.tight_layout()
plt.show()

We can also compare average power and reserve power over time. The graph shows that adding reserve constraints reduces the average generation of the cheaper generators and increases the average generation of the more expensive generators.

data = pd.concat(
    [n.generators_t.get("p").mean() for n in [n_basic, n_reserve]],
    axis=1,
    keys=["0 MW", f"{reserve_req} MW"],
)
data.plot(kind="bar", ylabel="mean(p) [MW]")
plt.show()

And that’s it. Feel free to change the values of reserve_req, a and b and see how this affects the results. However, be aware that reserve requirements that are too high can make the model infeasible.

You can also try out an interactive dashboard to play around with a small example model where you can try out different reserve settings at https://pypsa-reserves-dashboard.streamlit.app/