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-25
-------------------------------------------------------------------
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.EqualTo{Float64}     : [1, 1]
  VariableRef in MOI.GreaterThan{Float64} : [5, 5]
  VariableRef in MOI.LessThan{Float64}    : [1, 2]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
         1   1.875000e+04  1.991887e+03  1.255894e-02        18   1
        20   1.875000e+03  8.072917e+03  4.283786e-02       360   1
-------------------------------------------------------------------
status         : iteration_limit
total time (s) : 4.283786e-02
total solves   : 360
best bound     :  8.072917e+03
simulation ci  :  7.062233e+03 ± 3.307524e+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")

This should open a webpage that looks like this one.

Using the mouse, you can highlight individual trajectories by hovering over them. This makes it possible to visualize a single trajectory across multiple dimensions.

If you click on the plot, then trajectories that are close to the mouse pointer are shown darker and those further away are shown lighter.

Publication plots

Instead of the interactive Javascript plots, you can also create some publication ready plots using the SDDP.publication_plot function.

SDDP.publication_plot implements a plot recipe to create ribbon plots of each variable against the stages. The first argument is the vector of simulation dictionaries and the second argument is the dictionary key that you want to plot. Standard Plots.jl keyword arguments such as title and xlabel can be used to modify the look of each plot. By default, the plot displays ribbons of the 0-100, 10-90, and 25-75 percentiles. The dark, solid line in the middle is the median (i.e. 50'th percentile).

import Plots
Plots.plot(
    SDDP.publication_plot(simulations; title = "Outgoing volume") do data
        return data[:volume].out
    end,
    SDDP.publication_plot(simulations; title = "Thermal generation") do data
        return data[:thermal_generation]
    end;
    xlabel = "Stage",
    ylims = (0, 200),
    layout = (1, 2),
)

You can save this plot as a PDF using the Plots.jl function savefig:

Plots.savefig("my_picture.pdf")

Plotting the value function

You can obtain an object representing the value function of a node using SDDP.ValueFunction.

V = SDDP.ValueFunction(model[(1, 1)])

A value function for node (1, 1)

The value function can be evaluated using SDDP.evaluate.

SDDP.evaluate(V; volume = 1)

(23019.270833333332, Dict(:volume => -157.8125))

evaluate returns the height of the value function, and a subgradient with respect to the convex state variables.

You can also plot the value function using SDDP.plot

SDDP.plot(V, volume = 0:200, filename = "value_function.html")

This should open a webpage that looks like this one.

Convergence dashboard

If the text-based logging isn't to your liking, you can open a visualization of the training by passing dashboard = true to SDDP.train.

SDDP.train(model; dashboard = true)

By default, dashboard = false because there is an initial overhead associated with opening and preparing the plot.