Thermal Generation Dispatch Sweep Example
A classic application is to allocate power generation to meet a demand $d$ at minimum cost, subject to capacity limits. We consider a simplified economic dispatch problem:
\[\begin{align*} \min_{g_i \ge 0,\;\phi \ge 0} & \quad \sum_{i=1}^n c_i\,g_i + \sum_{i=1}^n c_{2,i}\,g_i^2 + c_\phi \,\phi \\ \text{s.t.} & \quad \sum_{i=1}^n g_i + \phi \;=\; d, \quad (:\lambda)\\ & \quad 0 \;\le\; g_i \;\le\; G_i, \quad i=1,\dots,n, \end{align*}\]
where $g_i$ is the power generated by plant $i$, each with linear unit cost $c_i$, quadratic cost component $c_{2,i}$ and capacity $G_i$, and $\phi$ is unmet demand with penalty $c_\phi$. We treat $d$ (the system demand) as a parameter. Differentiating through this QP quantifies how optimal generation shifts as demand changes.
When $d$ changes by a small amount, DiffOpt.forward\_differentiate! solves a linear system capturing the KKT conditions, revealing how $g_1^*$, $g_2^*$, and $\phi^*$ shift with respect to $d$. These sensitivities are critical for power systems operators to understand how different plants ramp up or down as demand fluctuates.
Define and solve the Thermal Dispatch Problem for a range of demands
First, import the libraries.
using Test
using JuMP
import DiffOpt
import HiGHS
using Plots
using Plots.MeasuresFixed data
c = [20.0, 30.0] # linear cost ($/MWh)
c2 = [0.2, 0.1] # quadratic cost ($/MWh²)
G = [150.0, 80.0] # capacities (MW)
cφ = 10_000.0 # penalty for unmet demand ($/MWh)
d_range = 11.0:4.0:300.0 # demand sweep (MW) – slightly coarser for clarity
N = length(d_range)73Storage for results
g1, g2, _φ = zeros(N), zeros(N), zeros(N)
dg1_dd, dg2_dd = zeros(N), zeros(N)
dφ_dd, dJ_dd = zeros(N), zeros(N)([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])Sweep & Differentiation
for (k, d_val) in enumerate(d_range)
println("Demand: ", k, " of ", N, " (", d_val, ")")
model = DiffOpt.quadratic_diff_model(HiGHS.Optimizer)
set_silent(model)
@variable(model, d in Parameter(d_val)) # parameter
@variables(model, begin # decisions
0 <= g[i=1:2] <= G[i]
φ >= 0
end)
@objective(
model,
Min, # cost
sum(c[i] * g[i] + c2[i] * g[i]^2 for i in 1:2) + cφ * φ
)
@constraint(model, con, sum(g) + φ == d) # balance
optimize!(model)
# ---------- store primal results ----------
g1[k] = value(g[1])
g2[k] = value(g[2])
_φ[k] = value(φ)
# ---------- forward sensitivities ----------
DiffOpt.empty_input_sensitivities!(model)
DiffOpt.set_forward_parameter(model, d, 1.0)
DiffOpt.forward_differentiate!(model)
dg1_dd[k] = DiffOpt.get_forward_variable(model, g[1])
dg2_dd[k] = DiffOpt.get_forward_variable(model, g[2])
dφ_dd[k] = DiffOpt.get_forward_variable(model, φ)
# marginal cost λ
dJ_dd[k] = dual.(con)
endDemand: 1 of 73 (11.0)
Demand: 2 of 73 (15.0)
Demand: 3 of 73 (19.0)
Demand: 4 of 73 (23.0)
Demand: 5 of 73 (27.0)
Demand: 6 of 73 (31.0)
Demand: 7 of 73 (35.0)
Demand: 8 of 73 (39.0)
Demand: 9 of 73 (43.0)
Demand: 10 of 73 (47.0)
Demand: 11 of 73 (51.0)
Demand: 12 of 73 (55.0)
Demand: 13 of 73 (59.0)
Demand: 14 of 73 (63.0)
Demand: 15 of 73 (67.0)
Demand: 16 of 73 (71.0)
Demand: 17 of 73 (75.0)
Demand: 18 of 73 (79.0)
Demand: 19 of 73 (83.0)
Demand: 20 of 73 (87.0)
Demand: 21 of 73 (91.0)
Demand: 22 of 73 (95.0)
Demand: 23 of 73 (99.0)
Demand: 24 of 73 (103.0)
Demand: 25 of 73 (107.0)
Demand: 26 of 73 (111.0)
Demand: 27 of 73 (115.0)
Demand: 28 of 73 (119.0)
Demand: 29 of 73 (123.0)
Demand: 30 of 73 (127.0)
Demand: 31 of 73 (131.0)
Demand: 32 of 73 (135.0)
Demand: 33 of 73 (139.0)
Demand: 34 of 73 (143.0)
Demand: 35 of 73 (147.0)
Demand: 36 of 73 (151.0)
Demand: 37 of 73 (155.0)
Demand: 38 of 73 (159.0)
Demand: 39 of 73 (163.0)
Demand: 40 of 73 (167.0)
Demand: 41 of 73 (171.0)
Demand: 42 of 73 (175.0)
Demand: 43 of 73 (179.0)
Demand: 44 of 73 (183.0)
Demand: 45 of 73 (187.0)
Demand: 46 of 73 (191.0)
Demand: 47 of 73 (195.0)
Demand: 48 of 73 (199.0)
Demand: 49 of 73 (203.0)
Demand: 50 of 73 (207.0)
Demand: 51 of 73 (211.0)
Demand: 52 of 73 (215.0)
Demand: 53 of 73 (219.0)
Demand: 54 of 73 (223.0)
Demand: 55 of 73 (227.0)
Demand: 56 of 73 (231.0)
Demand: 57 of 73 (235.0)
Demand: 58 of 73 (239.0)
Demand: 59 of 73 (243.0)
Demand: 60 of 73 (247.0)
Demand: 61 of 73 (251.0)
Demand: 62 of 73 (255.0)
Demand: 63 of 73 (259.0)
Demand: 64 of 73 (263.0)
Demand: 65 of 73 (267.0)
Demand: 66 of 73 (271.0)
Demand: 67 of 73 (275.0)
Demand: 68 of 73 (279.0)
Demand: 69 of 73 (283.0)
Demand: 70 of 73 (287.0)
Demand: 71 of 73 (291.0)
Demand: 72 of 73 (295.0)
Demand: 73 of 73 (299.0)Results with Plot graphs
default(;
size = (1150, 350),
legendfontsize = 8,
guidefontsize = 9,
tickfontsize = 7,
)Stacked-area dispatch plot
plt_dispatch = plot(
d_range,
g1; # lower layer
fillrange = 0, # fill down to zero
label = "g₁",
xlabel = "Demand (MWh)",
ylabel = "g (MWh)",
left_margin = 5Measures.mm,
bottom_margin = 5Measures.mm,
title = "Optimal dispatch",
);
plot!(
plt_dispatch,
d_range,
g1 .+ g2; # upper layer
fillrange = g1, # fill down to g₁
label = "g₂",
);
plot!(
plt_dispatch,
d_range,
_φ; # unmet demand as a line
lw = 2,
c = :red,
label = "φ (unserved)",
);
capacity_limit = d_range[findfirst(g2 .== 80.0)]
vline!(
plt_dispatch,
[capacity_limit, sum(G)];
l = :dash,
c = :black,
label = "capacity limit",
);(b) marginal cost λ(d)
plt_mc = plot(
d_range,
dJ_dd;
lw = 2,
c = :black,
xlabel = "Demand (MWh)",
ylabel = "λ (\$/MWh)",
title = "Electricity Price",
yscale = :log10,
left_margin = 5Measures.mm,
bottom_margin = 5Measures.mm,
legend = false,
);
vline!(plt_mc, [capacity_limit, sum(G)]; l = :dash, c = :gray, label = "");(c) generation sensitivities
plt_sens = plot(
d_range,
dg1_dd;
lw = 2,
xlabel = "Demand (MWh)",
ylabel = "∂g/∂d (MWh per MWh)",
title = "Generation Sensitivities",
label = "∂g₁/∂d",
);
plot!(plt_sens, d_range, dg2_dd; lw = 2, label = "∂g₂/∂d");
hline!(plt_sens, [0, 1]; l = :dot, c = :gray, label = "");
vline!(plt_sens, [capacity_limit, sum(G)]; l = :dash, c = :gray, label = "");combine all
plot_all = plot(plt_dispatch, plt_mc, plt_sens; layout = (1, 3))
Optimal dispatch (left), marginal electricity price $\lambda(d)=\partial J/\partial d$ (center, logarithmic scale), and forward sensitivities $\partial g_i/\partial d$ (right) for demands ranging from $\approx 0$ to MWh. The vertical dashed lines mark when the individual plant capacities are reached–-$g_2$ reaches it's capacity at MWh and $g_1$ at MWh.
This page was generated using Literate.jl.