Plotting tools
This tutorial was generated using Literate.jl. Download the source as a .jl
file. Download the source as a .ipynb
file.
In our previous tutorials, we formulated, solved, and simulated multistage stochastic optimization problems. However, we haven't really investigated what the solution looks like. Luckily, SDDP.jl
includes a number of plotting tools to help us do that. In this tutorial, we explain the tools and make some pretty pictures.
Preliminaries
The next two plot types help visualize the policy. Thus, we first need to create a policy and simulate some trajectories. So, let's take the model from Markovian policy graphs, train it for 20 iterations, and then simulate 100 Monte Carlo realizations of the policy.
using SDDP, HiGHS
Ω = [
(inflow = 0.0, fuel_multiplier = 1.5),
(inflow = 50.0, fuel_multiplier = 1.0),
(inflow = 100.0, fuel_multiplier = 0.75),
]
model = SDDP.MarkovianPolicyGraph(;
transition_matrices = Array{Float64,2}[
[1.0]',
[0.75 0.25],
[0.75 0.25; 0.25 0.75],
],
sense = :Min,
lower_bound = 0.0,
optimizer = HiGHS.Optimizer,
) do subproblem, node
t, markov_state = node
@variable(subproblem, 0 <= volume <= 200, SDDP.State, initial_value = 200)
@variables(subproblem, begin
thermal_generation >= 0
hydro_generation >= 0
hydro_spill >= 0
inflow
end)
@constraints(
subproblem,
begin
volume.out == volume.in + inflow - hydro_generation - hydro_spill
thermal_generation + hydro_generation == 150.0
end
)
probability =
markov_state == 1 ? [1 / 6, 1 / 3, 1 / 2] : [1 / 2, 1 / 3, 1 / 6]
fuel_cost = [50.0, 100.0, 150.0]
SDDP.parameterize(subproblem, Ω, probability) do ω
JuMP.fix(inflow, ω.inflow)
@stageobjective(
subproblem,
ω.fuel_multiplier * fuel_cost[t] * thermal_generation
)
end
end
SDDP.train(model; iteration_limit = 20, run_numerical_stability_report = false)
simulations = SDDP.simulate(
model,
100,
[:volume, :thermal_generation, :hydro_generation, :hydro_spill],
)
println("Completed $(length(simulations)) simulations.")
-------------------------------------------------------------------
SDDP.jl (c) Oscar Dowson and contributors, 2017-24
-------------------------------------------------------------------
problem
nodes : 5
state variables : 1
scenarios : 1.08000e+02
existing cuts : false
options
solver : serial mode
risk measure : SDDP.Expectation()
sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
VariableRef : [7, 7]
AffExpr in MOI.EqualTo{Float64} : [2, 2]
VariableRef in MOI.GreaterThan{Float64} : [5, 5]
VariableRef in MOI.LessThan{Float64} : [1, 2]
-------------------------------------------------------------------
iteration simulation bound time (s) solves pid
-------------------------------------------------------------------
1 1.562500e+04 1.991887e+03 1.465201e-02 18 1
20 3.375000e+04 8.072917e+03 4.963994e-02 360 1
-------------------------------------------------------------------
status : iteration_limit
total time (s) : 4.963994e-02
total solves : 360
best bound : 8.072917e+03
simulation ci : 6.710737e+03 ± 3.916781e+03
numeric issues : 0
-------------------------------------------------------------------
Completed 100 simulations.
Great! Now we have some data in simulations
to visualize.
Spaghetti plots
The first plotting utility we discuss is a spaghetti plot (you'll understand the name when you see the graph).
To create a spaghetti plot, begin by creating a new SDDP.SpaghettiPlot
instance as follows:
plt = SDDP.SpaghettiPlot(simulations)
A spaghetti plot with 100 scenarios and 3 stages.
We can add plots to plt
using the SDDP.add_spaghetti
function.
SDDP.add_spaghetti(plt; title = "Reservoir volume") do data
return data[:volume].out
end
In addition to returning values from the simulation, you can compute things:
SDDP.add_spaghetti(plt; title = "Fuel cost", ymin = 0, ymax = 250) do data
if data[:thermal_generation] > 0
return data[:stage_objective] / data[:thermal_generation]
else # No thermal generation, so return 0.0.
return 0.0
end
end
Note that there are many keyword arguments in addition to title
. For example, we fixed the minimum and maximum values of the y-axis using ymin
and ymax
. See the SDDP.add_spaghetti
documentation for all the arguments.
Having built the plot, we now need to display it using SDDP.plot
.
SDDP.plot(plt, "spaghetti_plot.html")