Piecewise costs and constraints¶
In this example, we demonstrate how to represent non-linear system relationships as piecewise-linear curves in PyPSA. This allows us to describe such relationships as:
- Economies of scale (investment costs per unit capacity reduce as an asset's capacity increases)
- Marginal cost offer curves (marginal cost per unit dispatched energy changes an asset's dispatch increases)
- Part-load efficiencies (efficiency varies when an asset operates below its rated capacity)
Here, we build a model from the single-node capacity expansion example and update the electrolysis component with various piecewise curves. Note that the goal here is not necessarily to showcase a realistic scenario, but rather to demonstrate how it is possible to model piecewise costs and constraints in PyPSA.
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
n = pypsa.examples.model_energy()
# Add a hydrogen load to ensure electrolysis *must* be dispatched in the solution
n.add("Load", "hydrogen-demand", bus="hydrogen", p_set=n.loads_t.p_set["demand"] / 100)
n.components["Link"].static.loc["electrolysis", "marginal_cost"] = 10
Reference case¶
Our reference case will be a model in which the rated values for our CCGT are used.
def update_electrolysis(n: pypsa.Network, **new_attrs) -> pypsa.Network:
n_new = n.copy()
n_new.add(
"Link",
"electrolysis",
**new_attrs,
**n.links.loc["electrolysis"].drop(new_attrs.keys()),
overwrite=True,
)
n_new.sanitize()
return n_new
n_ref = update_electrolysis(n)
n_ref.optimize()
Creating piecewise relationships¶
Let's assume the following relationships:
The marginal investment cost of our electrolysis facility decreases until a tipping point where additional capacity faces a diseconomy. This could be caused by land purchase costs increasing for additional capacity beyond a certain level. We will describe this with a sigmoid-like curve:
$$ \frac{C_{inv}}{C_{inv}^{max}}\ =\frac{\left(\ (\frac{P_{nom}}{P_{nom}^{max}}-1/2)\sqrt{\left(1-k/4\right)}\right)}{\sqrt{\left(1-k(\frac{P_{nom}}{P_{nom}^{max}}-1/2)^{2}\right)}}+1/2 $$
For this example, we will use a scaling factor $k$ of 2.4
The electrolysis component in PyPSA represents a collection of real facilities, each with different marginal costs which will be applied in order of merit. We will represent this marginal cost offer curve as:
$$ C_m = C_0 + \alpha \cdot \frac{P}{P^{max}} + \beta \cdot (\frac{P}{P^{max}})^3 $$
For this example, let's assume $C_0$ = 5, $\alpha$ = 2 €/MWh, $\beta$ = 3€/MW³h
The efficiency of electrolysis degrades at part-load. Let's assume it follows a square root curve bounded by $\eta$ = 70% at full load and $\eta$ = 40% at 0% load:
$$ \eta = (\eta^{max} - \eta^{min}) \cdot \sqrt{\frac{P}{P^{max}}} + \eta^{min} $$
We can approximate these non-linear relationships with a series of straight lines, or segments, which meet at breakpoints.
Here, we will heuristically select points along each curve to derive these segments.
def plot_piecewise_constraints(
x_max,
n_breakpoints,
y_ref,
y_func,
title,
xlabel,
ylabel,
x_min=0,
cumulative: bool = False,
):
x_continuous = np.linspace(x_min, x_max, 500)
x_breakpoints = np.linspace(x_min, x_max, n_breakpoints)
y_breakpoints = y_func(x_breakpoints)
if cumulative:
x_vals = np.repeat(x_breakpoints, 2)[:-1]
y_vals = np.roll(np.repeat(y_breakpoints, 2), -1)[:-1]
else:
x_vals = x_breakpoints
y_vals = y_breakpoints
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(x_continuous, y_func(x_continuous), lw=2, label="Non-linear curve")
ax.plot(x_vals, y_vals, "--o", lw=2, label="Piecewise linear")
ax.plot([x_min, x_max], y_ref, "-.o", lw=2, label="Reference linear")
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
ax.legend()
plt.tight_layout()
plt.show()
Investment cost – economy of scale followed by diseconomy of scale¶
P_max_inv = 5000 # MW
C_max_inv = 190000 # €/MW
k = 3
def get_capital_cost_curve(P: np.array):
p = P / P_max_inv - 0.5
return C_max_inv * (p * np.sqrt(1 - k / 4) / np.sqrt(1 - k * p**2) + 0.5)
plot_piecewise_constraints(
P_max_inv,
4,
(C_max_inv, C_max_inv),
get_capital_cost_curve,
"Capital cost piecewise curve approximation",
"Capacity (MW)",
"Marginal capital cost (€/MW)",
cumulative=True,
)
We integrate the capital cost curve when defining the piecewise constraint in PyPSA. This represents the need to pay the costs for each marginal capacity increase at the capital cost of that increment. This means that, compared to a non-linear curve, the piecewise constraint deviates at high dispatch level as the area under each curve differs. This deviation would reduce with more breakpoints.
def get_capex_curve(P: np.array):
capital_cost = get_capital_cost_curve(P)
return np.cumsum(
(capital_cost * P) - (capital_cost * pd.Series(P).shift(1).fillna(0).values)
)
plot_piecewise_constraints(
P_max_inv,
4,
(0, C_max_inv * P_max_inv),
get_capex_curve,
"CAPEX curve",
"Capacity (MW)",
"Total capital cost (€)",
)
Marginal cost offer curve¶
C0_mc = 5 # €/h
alpha_mc = 2 # €/MWh
beta_mc = 3 # €/MW²h
P_max = 5000
def get_marginal_cost_curve(P_frac):
return C0_mc + alpha_mc * P_frac + beta_mc * P_frac**3
ref_mc = C0_mc + alpha_mc + beta_mc
plot_piecewise_constraints(
1,
4,
(ref_mc, ref_mc),
get_marginal_cost_curve,
"Marginal cost offer curve",
"Load rate (Dispatch / Rated dispatch)",
"Cost (€/MWh)",
cumulative=True,
)
As with capital cost, we consider the integral of the marginal cost curve when defining the piecewise constraint in PyPSA. That is, each segment of the piecewise curve will define different marginal cost rates for increasing dispatch. This means that, compared to a non-linear curve, the piecewise constraint deviates at high dispatch level as the area under each curve is slightly different.
def get_opex_curve(P_out, p_max=P_max):
P_frac = P_out / p_max
marginal_cost = get_marginal_cost_curve(P_frac)
return np.cumsum(
(marginal_cost * P_out)
- (marginal_cost * pd.Series(P_out).shift(1).fillna(0).values)
)
plot_piecewise_constraints(
5000,
4,
(0, ref_mc * P_max),
get_opex_curve,
"OPEX curve",
"Dispatch (MWh)",
"Operating cost (€)",
)
Part-load efficiency¶
eta_rated = 0.7 # efficiency at full load
eta_min = 0.4 # efficiency at zero load
P_frac_min = 0 # minimum load fraction (η is finite at zero load)
P_max = 1000
def get_part_load_eff_curve(P_frac):
return (eta_rated - eta_min) * np.sqrt(P_frac) + eta_min
plot_piecewise_constraints(
1,
6,
(eta_rated, eta_rated),
get_part_load_eff_curve,
"Part-load efficiency",
"Load rate (Dispatch / Rated dispatch)",
"Efficiency",
x_min=P_frac_min,
cumulative=False,
)
We do not consider part-load as a marginal constraint, therefore we do not take the integral of the piecewise curve when applying the constraint:
def get_load_in_curve(P_out):
eff = get_part_load_eff_curve(P_out / P_max)
return P_out / eff
plot_piecewise_constraints(
P_max,
6,
(0, P_max / eta_rated),
get_load_in_curve,
"Energy consumption curve",
"Dispatch (MWh)",
"Energy consumption (MWh)",
x_min=P_frac_min,
)
Piecewise attributes for extendable components¶
Of the above curves, only (dis)economies of scale can be applied to the case where our electrolysis component is extendable.
The other two relate to the load rate of a component (p/p_nom) which is a non-linear relationship of two decision variables when the component is extendable.
Let's see what happens when we add this piecewise constraint:
p_nom_breakpoints = np.array([0, 1000, 2000, 3000, 4000])
capital_costs = get_capital_cost_curve(p_nom_breakpoints)
n_piecewise_capex = update_electrolysis(
n,
capital_cost=dict(zip(p_nom_breakpoints, capital_costs)),
)
n_piecewise_capex.optimize(reformulate_sos="auto")
We can see the results stored in a new component piecewise data table:
display(n_piecewise_capex.c.links.piecewise.capital_cost)
We can see that the capital expenditure in electrolysis capacity sits on our piecewise curve:
capital_cost_opt = n_piecewise_capex.links.loc[
"electrolysis", "capital_cost_piecewise_opt"
]
p_nom_opt = n_piecewise_capex.links.loc["electrolysis", "p_nom_opt"]
capex_opt = capital_cost_opt * p_nom_opt
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(p_nom_opt, capex_opt, label="Optimal capital cost")
ax.plot(
p_nom_breakpoints,
get_capex_curve(p_nom_breakpoints),
"--o",
color="red",
label="Piecewise linear capital cost",
)
ax.legend()
ax.set_xlabel("Nominal capacity (MW)")
ax.set_ylabel("CAPEX (€)")
ax.set_title("Electrolysis capital costs (CAPEX)")
Piecewise attributes for non-extendable components¶
Let's fix the capacity and apply our marginal cost piecewise constraint.
Unlike capital cost, which is defined as capital expenditure at each p_nom level, marginal cost is defined per unit dispatch, namely p_pu.
This means, our breakpoints exist in the range zero to one.
p_pu_breakpoints = np.array([0, 0.25, 0.5, 0.75, 1.0])
marginal_cost = get_marginal_cost_curve(p_pu_breakpoints)
n_piecewise_opex = update_electrolysis(
n,
marginal_cost=dict(zip(p_pu_breakpoints, marginal_cost)),
p_nom_extendable=False,
p_nom=n_ref.links.loc["electrolysis", "p_nom_opt"],
)
n_piecewise_opex.optimize()
As with capital cost in our previous example, we can see our marginal cost defined in the piecewise data table.
This time, however, the values are given per unit, i.e., €/MWh.
In the capital cost case, they were given as total cost at each point along the curve (i.e., €).
display(n_piecewise_opex.c.links.piecewise.marginal_cost)
We can see that the marginal costs of electrolysis follow our piecewise curve:
marginal_cost_opt = n_piecewise_opex.links_t.marginal_cost_piecewise_opt["electrolysis"]
p = n_piecewise_opex.links_t.p["electrolysis"]
opex_opt = marginal_cost_opt * p
p_nom_opt = n_piecewise_opex.links.loc["electrolysis", "p_nom"]
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(p, opex_opt, label="Optimal OPEX")
ax.plot(
p_pu_breakpoints * p_nom_opt,
get_opex_curve(p_pu_breakpoints * p_nom_opt, p_nom_opt),
"--o",
color="red",
label="Piecewise linear OPEX",
)
ax.legend()
ax.set_xlabel("Dispatch (MWh)")
ax.set_ylabel("OPEX (€)")
ax.set_title("Per-snapshot electrolysis operating costs (OPEX)")
Now let's fix the capacity and apply our part-load efficiency piecewise constraint.
As with marginal cost, it requires the capacity to be fixed and for the breakpoints to be defined per unit, i.e., p_pu.
Part-load efficiency curves often make the problem much harder to solve, especially compared to piecewise costs. Therefore, we will only solve the first 200 snapshots of our problem for this example.
p_pu_breakpoints = np.array([0, 0.25, 0.5, 0.75, 1.0])
efficiencies = get_part_load_eff_curve(p_pu_breakpoints)
n_piecewise_efficiency = update_electrolysis(
n,
efficiency=dict(zip(p_pu_breakpoints, efficiencies)),
p_nom_extendable=False,
p_nom=n_ref.links.loc["electrolysis", "p_nom_opt"],
)
n_piecewise_efficiency.optimize(reformulate_sos=True, snapshots=n.snapshots[:200])
We can see that the electrolysis efficiency follow our piecewise curve:
p1_opt = -1 * n_piecewise_efficiency.links_t.p1["electrolysis"]
p = n_piecewise_efficiency.links_t.p["electrolysis"]
p_nom_opt = n_piecewise_efficiency.links.loc["electrolysis", "p_nom"]
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(p, p1_opt, label="Optimal operation")
ax.plot(
p_pu_breakpoints * p_nom_opt,
efficiencies * p_nom_opt * p_pu_breakpoints,
"--o",
color="red",
label="Piecewise linear operation",
)
ax.legend()
ax.set_xlabel("Electricity input (MWh)")
ax.set_ylabel("Hydrogen output (MWh)")
ax.set_title("Per-snapshot electrolysis operating efficiency")