Alternative forward models

This tutorial was generated using Literate.jl. Download the source as a .jl file. Download the source as a .ipynb file.

This example demonstrates how to train convex and non-convex models.

This example uses the following packages:

using SDDP
import Ipopt
import PowerModels
import Test

Formulation

For our model, we build a simple optimal power flow model with a single hydro-electric generator.

The formulation of our optimal power flow problem depends on model_type, which must be one of the PowerModels formulations.

(To run locally, download pglib_opf_case5_pjm.m and update filename appropriately.)

function build_model(model_type)
    filename = joinpath(@__DIR__, "pglib_opf_case5_pjm.m")
    data = PowerModels.parse_file(filename)
    return SDDP.PolicyGraph(
        SDDP.UnicyclicGraph(0.95);
        sense = :Min,
        lower_bound = 0.0,
        optimizer = Ipopt.Optimizer,
    ) do sp, t
        power_model = PowerModels.instantiate_model(
            data,
            model_type,
            PowerModels.build_opf;
            jump_model = sp,
        )
        # Now add hydro power models. Assume that generator 5 is hydro, and the
        # rest are thermal.
        pg = power_model.var[:it][:pm][:nw][0][:pg][5]
        sp[:pg] = pg
        @variable(sp, x >= 0, SDDP.State, initial_value = 10.0)
        @variable(sp, deficit >= 0)
        @constraint(sp, balance, x.out == x.in - pg + deficit)
        @stageobjective(sp, objective_function(sp) + 1e6 * deficit)
        SDDP.parameterize(sp, [0, 2, 5]) do ω
            return SDDP.set_normalized_rhs(balance, ω)
        end
        return
    end
end

build_model (generic function with 1 method)

Training a convex model

We can build and train a convex approximation of the optimal power flow problem.

The problem with the convex model is that it does not accurately simulate the true dynamics of the problem. Therefore, it under-estimates the true cost of operation.

convex = build_model(PowerModels.DCPPowerModel)
SDDP.train(convex; iteration_limit = 10)

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-24
-------------------------------------------------------------------
problem
  nodes           : 1
  state variables : 1
  scenarios       : Inf
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [20, 20]
  AffExpr in MOI.EqualTo{Float64}         : [13, 13]
  AffExpr in MOI.Interval{Float64}        : [6, 6]
  VariableRef in MOI.GreaterThan{Float64} : [14, 14]
  VariableRef in MOI.LessThan{Float64}    : [11, 11]
numerical stability report
  matrix range     [1e+00, 2e+02]
  objective range  [1e+00, 1e+06]
  bounds range     [4e-01, 6e+00]
  rhs range        [5e-01, 5e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
         1   3.520513e+06  4.875069e+04  2.521582e-01        67   1
         6   7.742569e+05  2.311921e+05  1.901419e+00       354   1
         9   1.168861e+06  3.456228e+05  3.162240e+00       559   1
        10   3.474907e+05  3.470271e+05  3.572340e+00       626   1
-------------------------------------------------------------------
status         : iteration_limit
total time (s) : 3.572340e+00
total solves   : 626
best bound     :  3.470271e+05
simulation ci  :  7.944603e+05 ± 6.503414e+05
numeric issues : 0
-------------------------------------------------------------------

To more accurately simulate the dynamics of the problem, a common approach is to write the cuts representing the policy to a file, and then read them into a non-convex model:

SDDP.write_cuts_to_file(convex, "convex.cuts.json")
non_convex = build_model(PowerModels.ACPPowerModel)
SDDP.read_cuts_from_file(non_convex, "convex.cuts.json")

Now we can simulate non_convex to evaluate the policy.

result = SDDP.simulate(non_convex, 1)

1-element Vector{Vector{Dict{Symbol, Any}}}:
 [Dict(:bellman_term => 326259.1596948804, :noise_term => 0, :node_index => 1, :stage_objective => 27420.553408105046, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 323779.11129137216, :noise_term => 2, :node_index => 1, :stage_objective => 25428.67302085477, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 323779.11129053537, :noise_term => 2, :node_index => 1, :stage_objective => 23580.69745569545, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 326615.8628059642, :noise_term => 0, :node_index => 1, :stage_objective => 27420.55340810496, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 323779.1112914929, :noise_term => 2, :node_index => 1, :stage_objective => 25694.76533243011, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 326615.8628072464, :noise_term => 0, :node_index => 1, :stage_objective => 27420.55340810496, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 323779.1087386605, :noise_term => 5, :node_index => 1, :stage_objective => 18341.123411162327, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 323102.6234109749, :noise_term => 5, :node_index => 1, :stage_objective => 17628.650465247185, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 325228.08685101615, :noise_term => 0, :node_index => 1, :stage_objective => 27420.55340810501, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 328064.8383913751, :noise_term => 0, :node_index => 1, :stage_objective => 27420.553408105476, :objective_state => nothing, :belief => Dict(1 => 1.0))  …  Dict(:bellman_term => 325228.0920350857, :noise_term => 0, :node_index => 1, :stage_objective => 27420.553408105006, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 328064.84357544454, :noise_term => 0, :node_index => 1, :stage_objective => 27420.55340810546, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 323779.10921065696, :noise_term => 5, :node_index => 1, :stage_objective => 19055.17274611591, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 323102.6238829171, :noise_term => 5, :node_index => 1, :stage_objective => 17628.65046524719, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 325228.0878191785, :noise_term => 0, :node_index => 1, :stage_objective => 27420.55340810501, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 323779.10802892444, :noise_term => 5, :node_index => 1, :stage_objective => 17658.713182527557, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 323779.11129053537, :noise_term => 2, :node_index => 1, :stage_objective => 23580.692512244987, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 326615.8628059642, :noise_term => 0, :node_index => 1, :stage_objective => 27420.553408104963, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 323779.1112914929, :noise_term => 2, :node_index => 1, :stage_objective => 25694.76533243009, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 326615.86280724645, :noise_term => 0, :node_index => 1, :stage_objective => 27420.553408104963, :objective_state => nothing, :belief => Dict(1 => 1.0))]

A problem with reading and writing the cuts to file is that the cuts have been generated from trial points of the convex model. Therefore, the policy may be arbitrarily bad at points visited by the non-convex model.

Training a non-convex model

We can also build and train a non-convex formulation of the optimal power flow problem.

The problem with the non-convex model is that because it is non-convex, SDDP.jl may find a sub-optimal policy. Therefore, it may over-estimate the true cost of operation.

non_convex = build_model(PowerModels.ACPPowerModel)
SDDP.train(non_convex; iteration_limit = 10)
result = SDDP.simulate(non_convex, 1)

1-element Vector{Vector{Dict{Symbol, Any}}}:
 [Dict(:bellman_term => 390621.75143196946, :noise_term => 0, :node_index => 1, :stage_objective => 17578.223553921172, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 390154.39565383026, :noise_term => 5, :node_index => 1, :stage_objective => 17578.22355391896, :objective_state => nothing, :belief => Dict(1 => 1.0)), Dict(:bellman_term => 392171.4726263264, :noise_term => 0, :node_index => 1, :stage_objective => 25425.091955298987, :objective_state => nothing, :belief => Dict(1 => 1.0))]

Combining convex and non-convex models

To summarize, training with the convex model constructs cuts at points that may never be visited by the non-convex model, and training with the non-convex model may construct arbitrarily poor cuts because a key assumption of SDDP is convexity.

As a compromise, we can train a policy using a combination of the convex and non-convex models; we'll use the non-convex model to generate trial points on the forward pass, and we'll use the convex model to build cuts on the backward pass.

convex = build_model(PowerModels.DCPPowerModel)

A policy graph with 1 nodes.
 Node indices: 1

non_convex = build_model(PowerModels.ACPPowerModel)

A policy graph with 1 nodes.
 Node indices: 1

To do so, we train convex using the SDDP.AlternativeForwardPass forward pass, which simulates the model using non_convex, and we use SDDP.AlternativePostIterationCallback as a post-iteration callback, which copies cuts from the convex model back into the non_convex model.

SDDP.train(
    convex;
    forward_pass = SDDP.AlternativeForwardPass(non_convex),
    post_iteration_callback = SDDP.AlternativePostIterationCallback(non_convex),
    iteration_limit = 10,
)

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-24
-------------------------------------------------------------------
problem
  nodes           : 1
  state variables : 1
  scenarios       : Inf
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [20, 20]
  AffExpr in MOI.EqualTo{Float64}         : [13, 13]
  AffExpr in MOI.Interval{Float64}        : [6, 6]
  VariableRef in MOI.GreaterThan{Float64} : [14, 14]
  VariableRef in MOI.LessThan{Float64}    : [11, 11]
numerical stability report
  matrix range     [1e+00, 2e+02]
  objective range  [1e+00, 1e+06]
  bounds range     [4e-01, 6e+00]
  rhs range        [5e-01, 5e+00]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
         1   5.265564e+04  6.484820e+04  7.032299e-02        12   1
         3   1.163140e+06  1.381730e+05  1.080794e+00       123   1
         4   1.451693e+06  2.748161e+05  2.912240e+00       315   1
         5   1.357166e+06  3.265345e+05  4.612341e+00       468   1
         9   6.049396e+05  3.565627e+05  5.723814e+00       585   1
        10   6.820817e+05  3.642197e+05  6.726104e+00       678   1
-------------------------------------------------------------------
status         : iteration_limit
total time (s) : 6.726104e+00
total solves   : 678
best bound     :  3.642197e+05
simulation ci  :  5.596734e+05 ± 3.599739e+05
numeric issues : 0
-------------------------------------------------------------------

In practice, if we were to simulate non_convex now, we should obtain a better policy than either of the two previous approaches.