Demand Elasticity¶
This example demonstrates how demand elasticity can be modelled in PyPSA, using a single node capacity expansion model in the style of model.energy.
See Brown, Neumann, Riepin (2025) for more details.
Preparations¶
We start by loading the required packages, an example network and defining an utility function to compute the price-duration curve.
The network used here (see pypsa.examples.model_energy for details) is a minimal renewable-based system. It consists of wind and solar generators, a hydrogen conversion chain via electrolysis and a turbine, and storage (battery and hydrogen), which allows energy to be shifted over time.
In the default configuration, the main components have extendable capacities, so the optimization jointly determines both capacity investments and dispatch subject to availability, conversion efficiencies, and storage dynamics.
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pypsa
plt.style.use("bmh")
n = pypsa.examples.model_energy()
n.remove("Generator", "load shedding")
INFO:pypsa.network.io:Imported network 'Model-Energy' has buses, carriers, generators, links, loads, storage_units, stores
def get_price_duration(n: pypsa.Network, bus: str = "electricity") -> pd.Series:
s = (
n.buses_t.marginal_price[bus]
.sort_values(ascending=False)
.reset_index(drop=True)
)
s.index = np.arange(0, 100, 100 / len(s.index))
return s
To save computation time, we will only sample the first day of each month. Each day will be considered at a 3-hourly resolution, so there will be eight snapshots per representative day.
For computational reasons, the correct snapshot weighting is not applied here. Using smaller weights keeps the objective coefficients smaller and improves numerical stability, which is necessary for reliable convergence with the default HiGHS solver. As a consequence, the resulting marginal costs should be interpreted in relative terms. When using a more robust solver such as Gurobi, the weighting can be re-enabled.
days = n.snapshots.normalize().to_series()
one_day_per_month = days.groupby(days.dt.to_period("M")).first().values
snapshots = n.snapshots[n.snapshots.normalize().isin(one_day_per_month)]
n.set_snapshots(snapshots)
n.snapshot_weightings[["objective", "generators"]] *= 5
Perfectly inelastic demand¶
In the original example, demand varies both seasonally and diurnally. This time-varying load is now replaced by a fixed electricity demand. Most commonly, capacity expansion models would prescribe a perfectly inelastic demand via the p_set attribute, e.g. 100 MW.
The utility drawn from this consumption is effectively infinite. The model has to find a way to satisfy it. Otherwise, the model is infeasible.
n.remove("Load", "demand")
n.add("Load", "demand", bus="electricity", p_set=100)
n.optimize(solver_name="gurobi")
/tmp/ipykernel_3687/1458849580.py:1: 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(solver_name="gurobi")
INFO:linopy.model: Solve problem using Gurobi solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.1s
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-tsnwscpd.lp
Reading time = 0.00 seconds
obj: 1926 rows, 870 columns, 3800 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: 870 primals, 1926 duals Objective: 6.97e+07 Solver model: available Solver message: 2
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints 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.
('ok', 'optimal')
Market clearing prices can spike to extreme values in few hours of the year, while remaining close to zero for a majority of time.
Although all operational marginal costs are set to zero (for generation, conversion and storage), the optimized network produces non-zero marginal prices. This is to be expected when capacities are chosen endogenously, since the objective still includes investment (capital) costs for extendable assets.
fig, ax = plt.subplots()
get_price_duration(n).plot(
ax=ax,
ylabel="Clearing Price (€/MWh)",
xlabel="Fraction of Time (%)",
label="default",
legend=True,
)
<Axes: xlabel='Fraction of Time (%)', ylabel='Clearing Price (€/MWh)'>
capacities = n.statistics.optimal_capacity(round=2).to_frame("inelastic")
capacities
| inelastic | ||
|---|---|---|
| component | carrier | |
| Generator | solar | 331.40 |
| wind | 278.24 | |
| Link | electrolysis | 9.61 |
| turbine | 30.40 | |
| StorageUnit | battery storage | 297.76 |
| Store | hydrogen storage | 912.02 |
Perfectly inelastic demand up to VOLL¶
One way to avoid the price spikes is to model demand as perfectly inelastic up to a predefined value of lost load (VOLL).
Effectively, this is defined by a utility function $U(d) = Vd$ with a constant value $V$ for consumption $d\in[0,D]$, for instance 1000 €/MWh.
In electricity market models, utility generally represents the willingness to pay that consumers obtain from consuming electricity. Formally, utility is a function $U(d)$ of demand $d$. Its key property is that its derivative equals the (inverse) demand curve:
$\frac{\partial U(d)}{\partial d} = p(d)$
where $p(d)$ is the maximum price consumers are willing to pay for the marginal unit of electricity.
In our example, the demand curve is a step function. It is perfectly inelastic up to a price of $V$ at which point it is perfectly elastic.
When we make the substitution $d=D-g$ (where $g$ denotes the unserved demand), we see that we can model the VOLL case with a load shedding generator with marginal costs of $V=1000$ €/MWh, omitting the term $VD$ as it is a constant and does not affect optimization.
$U(d) = Vd$
$U(d) = VD - Vg$
Note that the objective sense of PyPSA is to minimize costs in order to maximize utility, so any costs have a positive sign and utility gains have a negative sign in the objective.
Therefore, we add a load-shedding generator with marginal cost $V$ and capacity to meet the maximum demand $D$ (here p_nom=100 as $D=100$ MW). When the load shedding generator produces 1 unit, it means that 1 unit of demand is not met. Because it has a marginal cost of $V$, the optimizer will only activate it when all available generating options are more expensive than $V$, or when there is not enough physical capacity to meet demand.
If supply becomes insufficient, the model dispatches load shedding at cost $V$, and the market price increases up to exactly $V$. Thus, VOLL acts as an upper bound on market prices and ensures feasibility by allowing controlled load shedding instead of forcing the model to meet demand at arbitrarily high cost.
n.add(
"Generator",
"load-shedding",
bus="electricity",
carrier="load",
marginal_cost=1000,
p_nom=100,
)
n.optimize(solver_name="gurobi")
/tmp/ipykernel_3687/1458849580.py:1: 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(solver_name="gurobi") WARNING:pypsa.consistency:The following generators have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers: Index(['load-shedding'], dtype='object', name='name')
INFO:linopy.model: Solve problem using Gurobi solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.1s
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-9ty6gfqy.lp
Reading time = 0.00 seconds
obj: 2118 rows, 966 columns, 4088 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: 966 primals, 2118 duals Objective: 6.37e+07 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.
('ok', 'optimal')
Now, the peak price is capped at 1000 €/MWh:
get_price_duration(n).plot(ax=ax, label="VOLL", legend=True)
fig
This results in some changes in the cost-optimal capacity mix, in particular in terms of backup capacities.
capacities["VOLL"] = n.statistics.optimal_capacity(round=2)
capacities
| inelastic | VOLL | ||
|---|---|---|---|
| component | carrier | ||
| Generator | solar | 331.40 | 331.27 |
| wind | 278.24 | 257.94 | |
| Link | electrolysis | 9.61 | NaN |
| turbine | 30.40 | NaN | |
| StorageUnit | battery storage | 297.76 | 228.78 |
| Store | hydrogen storage | 912.02 | NaN |
Linear demand curve¶
In reality, electricity demand is at least partially elastic. Consumers would use less electricity if it were more expensive, or more if prices were low.
The utility is the area under the demand curve. For a linear demand curve $p = a - bd$, where $p$ is the price, the utility is quadratic
$U(d) = ad - 0.5 b d^2$.
Maximum consumption occurs when the price is $p=0$, in which case $d_{max} = D = a/b$.
Price elasticity of demand, often expressed as a percentage, is calculated with $\varepsilon = \frac{\%\ \text{change in quantity}}{\%\ \text{change in price}}.$
For a choice of $a=2000$ and $b=20$, the demand curve looks like this:
x = np.linspace(0, 100, 200)
plt.figure(figsize=(8, 4))
plt.plot(x, 2000 - 20 * x)
plt.xlabel("Demand (MW)")
plt.ylabel("Price (€/MWh)")
Text(0, 0.5, 'Price (€/MWh)')
That means, for instance, at a price of 1000 €/MWh, the demand would be only 50 MW. At a price of 400 €/MWh, 80 MW. And so on.
Therefore, the price elasticity in this case would be $\varepsilon=-5\%$. This elasticity is consistent with the empirical findings for the German electricity market reported in Arnold (2023).
Since PyPSA prefers a fixed demand, elasticity is modeled by introducing a variable for unserved demand $g$. So instead of defining a variable demand $d$, we substitute like before
$d = D-g = a/b - g$.
Applying this substitution to the utility function
$U(d) = ad - 0.5 b d^2$
yields
$U(g) = \frac{a^2}{2b} - 0.5 b g^2$.
As the constant term $\frac{a^2}{2b}$ does not affect the optimization, it can be omitted. The remaining term $- 0.5 b g^2$ represents the loss of utility due to unserved demand. When switching from welfare maximization to PyPSA’s cost-minimization convention, this term becomes a positive cost
$b/2$,
which corresponds to a virtual load-shedding generator whose marginal cost (not total cost) increases linearly as $b g$. Each each additional unit of unserved load costs more than the previous one, but the marginal cost grows linearly. This captures the idea that small amounts of unmet load are tolerable, but large amounts become increasingly costly.
n.remove("Generator", "load-shedding")
n.add(
"Generator",
"load-shedding",
bus="electricity",
carrier="load",
marginal_cost_quadratic=20 / 2,
p_nom=100,
)
Due to the quadratic terms in the objective function, this addition transforms the model into a quadratic programming (QP) problem. This quadratic formulation could cause HiGHS to stall or take a long time to iterate. If available, prefer another solver like Gurobi for the following examples about QPs.
n.optimize(solver_name="gurobi")
/tmp/ipykernel_3687/1458849580.py:1: 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(solver_name="gurobi") WARNING:pypsa.consistency:The following generators have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers: Index(['load-shedding'], dtype='object', name='name')
INFO:linopy.model: Solve problem using Gurobi solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.1s
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-yeu4x5ic.lp
Reading time = 0.00 seconds
obj: 2118 rows, 966 columns, 4088 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: 966 primals, 2118 duals Objective: 5.14e+07 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.
('ok', 'optimal')
The price duration curve is considerably smoother with less extreme prices and fewer zero-price hours:
get_price_duration(n).plot(ax=ax, label="linear-elastic", legend=True)
fig
Also, the optimized capacity mix is drastically different. The model cuts down on balancing technologies and rather curtails a lot demand instead.
capacities["linear-elastic"] = n.statistics.optimal_capacity(round=2)
capacities
| inelastic | VOLL | linear-elastic | ||
|---|---|---|---|---|
| component | carrier | |||
| Generator | solar | 331.40 | 331.27 | 204.46 |
| wind | 278.24 | 257.94 | 228.34 | |
| Link | electrolysis | 9.61 | NaN | NaN |
| turbine | 30.40 | NaN | NaN | |
| StorageUnit | battery storage | 297.76 | 228.78 | 111.36 |
| Store | hydrogen storage | 912.02 | NaN | NaN |
The drawback and explanation here is that the linear demand curve becomes unrealistically elastic at higher prices.
Partial demand elasticity¶
It is also possible to mix different demand modelling approaches. For instance, keeping 80% of demand perfectly inelastic, while modelling 20% with a linear demand curve.
This just requires adjusting the capacity and cost terms of the load shedding generator (adjusts the slope of the linear demand curve) This plot shows the 80% inelastic block as a vertical line at the baseline demand and a linear downward segment for the 20% elastic share. The inelastic part is effectively a price-insensitive wall that extends upward (theoretically to infinity), while the elastic portion slopes down as willingness to pay drops with additional curtailed load.
# Parameters of linear demand curve
a, b = 2000, 20
# Set share of elastic demand, here 20%
share_elastic = 0.2
# Get total demand
D = n.loads.loc["demand", "p_set"]
# Set load-shedding parameters according to elasticity
n.generators.at["load-shedding", "p_nom_max"] = share_elastic * D
n.generators.at["load-shedding", "marginal_cost_quadratic"] = b / (2 * share_elastic)
plt.figure(figsize=(8, 4))
plt.plot(
[D * (1 - share_elastic), D], [a, 0], marker="o", label="Linear elastic segment"
)
plt.vlines(
D * (1 - share_elastic),
ymin=a,
ymax=a * 1.2,
linestyles="dashed",
label="Perfectly inelastic segment",
)
plt.xlim(left=0)
plt.ylim(top=a * 1.2)
plt.xlabel("Demand (MW)")
plt.ylabel("Price (€/MWh)")
plt.legend(loc="upper left")
<matplotlib.legend.Legend at 0x7b19101be490>
n.optimize(solver_name="gurobi")
/tmp/ipykernel_3687/1458849580.py:1: 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(solver_name="gurobi") WARNING:pypsa.consistency:The following generators have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers: Index(['load-shedding'], dtype='object', name='name')
INFO:linopy.model: Solve problem using Gurobi solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.1s
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-mqfftgy0.lp
Reading time = 0.00 seconds
obj: 2118 rows, 966 columns, 4088 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: 966 primals, 2118 duals Objective: 6.27e+07 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.
('ok', 'optimal')
get_price_duration(n).plot(ax=ax, label="linear-elastic 20%", legend=True)
fig
Piecewise-linear demand curve¶
It is also possible to model a set of piecewise linear demand curves, e.g. to approximate a log-log demand curve ($\ln p = a - b \ln d$) as shown in Brown, Neumann, Riepin (2025), Section 3.2 and Appendix A.
The following choice of segments reflects an elasticity of $-5\%$ at a price of 100 €/MWh and a demand around 100 MW.
plt.figure(figsize=(8, 4))
x = [0, 95, 100, 110]
y = [8000, 400, 200, 0]
plt.plot(x, y, marker="o", label="Piecewise linear approximation of log-log")
plt.xlim(left=0)
plt.ylim(a * -0.05, a * 1.2)
plt.xlabel("Demand (MW)")
plt.ylabel("Price (€/MWh)")
plt.legend(loc="upper left")
<matplotlib.legend.Legend at 0x7b18fab70e10>
Each piecewise linear segment is modelled as a generator with its own marginal cost and quadratic marginal cost, so dispatch picks the cheapest segments first and traces the stepped demand curve.
The marginal_cost_quadratic $c_i^{(2)}$ is $ \frac{1}{2}$ times the slope $b$ (see examples above). For example for the middle segment
$c_\text{middle}^{(2)} = \frac{1}{2} b = \frac{1}{2}\frac{400-200}{100-95} = 20$.
The marginal_cost $c_i$ is the lower bound of each segment. For example for the middle segment
$c_\text{middle} = c_\text{right} + 2 * c_\text{right}^{(2)} * \overline{P}_\text{right} = 0 + 2 \cdot 10 \cdot 10=200$,
where $\overline{P}_\text{i}$ is the installed capacity.
n.remove("Generator", "load-shedding")
# Add load-shedding generators to model segments from right to left (cheapest first)
p_nom = [10, 5, 95]
# Quadratic marginal costs
mc2 = 0.5 * np.array([20, 40, 80])
# Marginal costs (lower bound of each segment)
mc_right = 0
mc_middle = mc_right + 2 * mc2[0] * p_nom[0]
mc_left = mc_middle + 2 * mc2[1] * p_nom[1]
n.add(
"Generator",
name=["load-shedding-right", "load-shedding-middle", "load-shedding-left"],
bus="electricity",
carrier="load",
p_nom=p_nom,
marginal_cost_quadratic=mc2,
marginal_cost=[mc_right, mc_middle, mc_left],
)
n.optimize(solver_name="gurobi")
/tmp/ipykernel_3687/1458849580.py:1: 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(solver_name="gurobi") WARNING:pypsa.consistency:The following generators have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers: Index(['load-shedding-right', 'load-shedding-middle', 'load-shedding-left'], dtype='object', name='name')
INFO:linopy.model: Solve problem using Gurobi solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.1s
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-26bvrn9i.lp
Reading time = 0.00 seconds
obj: 2502 rows, 1158 columns, 4664 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: 1158 primals, 2502 duals Objective: 5.62e+07 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.
('ok', 'optimal')
get_price_duration(n).plot(ax=ax, label="log-log approximation", legend=True)
fig
