Nested Benders with MIP/LP solvers

I want to solve a nested optimization problem using Benders decomposition, where the first stage relies on the second stage, and the second stage relies on the third stage. The challenge is that the second stage is updated with Benders cuts from the third stage, which changes its dual formulation in each iteration. As a result, the cuts passed to the master problem are also updated. Could someone provide a small, reproducible example of such a case? I am looking for a solution method based on the JuMP Benders decomposition structure.

I have a package SDDP.jl which implements nested Benders.

If you don’t want to use SDDP.jl directly, here’s a tutorial that implements nested Benders:

There’s also this tutorial for regular two-stage Benders: Benders decomposition · JuMP

A nice realization is that there is nothing special about three stages or two stages. Stages 1 and 2 form a two-stage Benders, and stages 2 and 3 form a two-stage Benders. All that matters is how you iterate on the graph.

You could do:

  1. Solve stage 1
  2. Solve stage 2
  3. Solve stage 3
  4. Add cut to stage 2
  5. Go to step (3) utill convergence
  6. Add cut to stage 1
  7. Go to step (1) until convergence

Or you could do

  1. Solve stage 1
  2. Solve stage 2
  3. Solve stage 3
  4. Add cut to stage 2
  5. Add cut to stage 1
  6. Go to step (1) until convergence
2 Likes

Hi @odow

Currently, I am trying to figure it by myself without any external package, but thanks for [SDDP.jl]. I understand your steps, but my concern is with step 5 (add cut to stage 1). In case of two stage Benders, this is straightforward, since the number of constraints are fixed. But if I understand correctly, in three stage, additional constraints are added to stage 2, so I am not sure if following the [Benders decomposition · JuMP], all three stages can be synced or is it something like manually adding cuts from one Benders to another one.

Stripping the code from Introductory theory · SDDP.jl, removing the stuff related to stochasticity, and removing the arbitrary graph in favor of a linear graph gives:

import HiGHS
import JuMP

struct Node
    subproblem::JuMP.Model
    states::Dict{Symbol,Tuple{JuMP.VariableRef,JuMP.VariableRef}}
    cost_to_go::JuMP.VariableRef
end

struct PolicyGraph
    nodes::Vector{Node}
end

function PolicyGraph(
    subproblem_builder::Function;
    stages::Int,
    lower_bound::Float64,
    optimizer,
)
    nodes = Node[]
    for t in 1:stages
        model = JuMP.Model(optimizer)
        JuMP.set_silent(model)
        states = subproblem_builder(model, t)
        JuMP.@variable(model, cost_to_go >= lower_bound)
        obj = JuMP.objective_function(model)
        JuMP.@objective(model, Min, obj + cost_to_go)
        if t == stages
            JuMP.fix(cost_to_go, 0.0; force = true)
        end
        push!(nodes, Node(model, states, cost_to_go))
    end
    return PolicyGraph(nodes)
end

function forward_pass(model::PolicyGraph, io::IO = stdout)
    println(io, "| Forward Pass")
    incoming_state =
        Dict(k => JuMP.fix_value(v[1]) for (k, v) in model.nodes[1].states)
    trajectory, simulation_cost = Dict{Symbol,Float64}[], 0.0
    for t in 1:length(model.nodes)
        node = model.nodes[t]
        println(io, "| | Visiting node $(t)")
        println(io, "| | | x = ", incoming_state)
        for (k, v) in incoming_state
            JuMP.fix(node.states[k][1], v; force = true)
        end
        JuMP.optimize!(node.subproblem)
        JuMP.assert_is_solved_and_feasible(node.subproblem)
        outgoing_state = Dict(k => JuMP.value(v[2]) for (k, v) in node.states)
        println(io, "| | | x′ = ", outgoing_state)
        stage_cost =
            JuMP.objective_value(node.subproblem) - JuMP.value(node.cost_to_go)
        simulation_cost += stage_cost
        println(io, "| | | C(x, u) = ", stage_cost)
        incoming_state = outgoing_state
        push!(trajectory, outgoing_state)
    end
    return trajectory, simulation_cost
end

function backward_pass(
    model::PolicyGraph,
    trajectory::Vector{Dict{Symbol,Float64}},
    io::IO = stdout,
)
    println(io, "| Backward pass")
    for t in length(trajectory)-1:-1:1
        outgoing_states = trajectory[t]
        node = model.nodes[t]
        println(io, "| | Visiting node $t")
        cut_expression = zero(JuMP.AffExpr)
        next_node = model.nodes[t+1]
        for (k, v) in outgoing_states
            JuMP.fix(next_node.states[k][1], v; force = true)
        end
        JuMP.optimize!(next_node.subproblem)
        V = JuMP.objective_value(next_node.subproblem)
        dVdx = Dict(k => JuMP.reduced_cost(v[1]) for (k, v) in next_node.states)
        println(io, "| | | | V, dVdx′ = ", V, ", ", dVdx)
        cut_expression += JuMP.@expression(
            node.subproblem,
            V + sum(dVdx[k] * (x[2] - outgoing_states[k]) for (k, x) in node.states)
        )
        c = JuMP.@constraint(node.subproblem, node.cost_to_go >= cut_expression)
        println(io, "| | | Adding cut : ", c)
    end
    return nothing
end

function lower_bound(model::PolicyGraph)
    JuMP.optimize!(model.nodes[1].subproblem)
    return JuMP.objective_value(model.nodes[1].subproblem)
end

upper_bound(model::PolicyGraph) = last(forward_pass(model, devnull))

function train(
    model::PolicyGraph;
    iteration_limit::Int,
    io::IO = stdout,
)
    for i in 1:iteration_limit
        println(io, "Starting iteration $(i)")
        outgoing_states, _ = forward_pass(model, io)
        backward_pass(model, outgoing_states, io)
        println(io, "| Finished iteration")
        println(io, "| | lower_bound = ", lower_bound(model))
        println(io, "| | upper_bound = ", upper_bound(model))
    end
    return
end

function evaluate_policy(
    model::PolicyGraph;
    node::Int,
    incoming_state::Dict{Symbol,Float64},
)
    for (k, v) in incoming_state
        JuMP.fix(model.nodes[node].states[k][1], v; force = true)
    end
    JuMP.optimize!(model.nodes[node].subproblem)
    return Dict(
        k => JuMP.value.(v) for
        (k, v) in JuMP.object_dictionary(model.nodes[node].subproblem)
    )
end

model = PolicyGraph(;
    stages = 3,
    lower_bound = 0.0,
    optimizer = HiGHS.Optimizer,
) do sp, t
    JuMP.@variables(sp, begin
        volume_in == 200
        0 <= volume_out <= 200
        thermal_generation >= 0
        hydro_generation >= 0
        hydro_spill >= 0
        inflow == 50
    end)
    JuMP.@constraints(sp, begin
        volume_out == volume_in + inflow - hydro_generation - hydro_spill
        demand_constraint, thermal_generation + hydro_generation == 150.0
    end)
    fuel_cost = [50.0, 100.0, 150.0]
    JuMP.@objective(sp, Min, fuel_cost[t] * thermal_generation)
    return Dict(:volume => (volume_in, volume_out))
end;

train(model; iteration_limit = 3)

evaluate_policy(model; node = 1, incoming_state = Dict(:volume => 150.0))

Here’s the output

julia> model = PolicyGraph(;
           stages = 3,
           lower_bound = 0.0,
           optimizer = HiGHS.Optimizer,
       ) do sp, t
           JuMP.@variables(sp, begin
               volume_in == 200
               0 <= volume_out <= 200
               thermal_generation >= 0
               hydro_generation >= 0
               hydro_spill >= 0
               inflow == 50
           end)
           JuMP.@constraints(sp, begin
               volume_out == volume_in + inflow - hydro_generation - hydro_spill
               demand_constraint, thermal_generation + hydro_generation == 150.0
           end)
           fuel_cost = [50.0, 100.0, 150.0]
           JuMP.@objective(sp, Min, fuel_cost[t] * thermal_generation)
           return Dict(:volume => (volume_in, volume_out))
       end;

julia> train(model; iteration_limit = 3)
Starting iteration 1
| Forward Pass
| | Visiting node 1
| | | x = Dict(:volume => 200.0)
| | | x′ = Dict(:volume => 0.0)
| | | C(x, u) = 0.0
| | Visiting node 2
| | | x = Dict(:volume => 0.0)
| | | x′ = Dict(:volume => 0.0)
| | | C(x, u) = 10000.0
| | Visiting node 3
| | | x = Dict(:volume => 0.0)
| | | x′ = Dict(:volume => 0.0)
| | | C(x, u) = 15000.0
| Backward pass
| | Visiting node 2
| | | | V, dVdx′ = 15000.0, Dict(:volume => -150.0)
| | | Adding cut : 150 volume_out + cost_to_go ≥ 15000
| | Visiting node 1
| | | | V, dVdx′ = 22500.0, Dict(:volume => -150.0)
| | | Adding cut : 150 volume_out + cost_to_go ≥ 22500
| Finished iteration
| | lower_bound = 2500.0
| | upper_bound = 7500.0
Starting iteration 2
| Forward Pass
| | Visiting node 1
| | | x = Dict(:volume => 200.0)
| | | x′ = Dict(:volume => 150.0)
| | | C(x, u) = 2500.0
| | Visiting node 2
| | | x = Dict(:volume => 150.0)
| | | x′ = Dict(:volume => 100.0)
| | | C(x, u) = 5000.0
| | Visiting node 3
| | | x = Dict(:volume => 100.0)
| | | x′ = Dict(:volume => 0.0)
| | | C(x, u) = 0.0
| Backward pass
| | Visiting node 2
| | | | V, dVdx′ = 0.0, Dict(:volume => -150.0)
| | | Adding cut : 150 volume_out + cost_to_go ≥ 15000
| | Visiting node 1
| | | | V, dVdx′ = 5000.0, Dict(:volume => -100.0)
| | | Adding cut : 100 volume_out + cost_to_go ≥ 20000
| Finished iteration
| | lower_bound = 5000.0
| | upper_bound = 5000.000000000003
Starting iteration 3
| Forward Pass
| | Visiting node 1
| | | x = Dict(:volume => 200.0)
| | | x′ = Dict(:volume => 199.99999999999997)
| | | C(x, u) = 5000.0
| | Visiting node 2
| | | x = Dict(:volume => 199.99999999999997)
| | | x′ = Dict(:volume => 100.0)
| | | C(x, u) = 2.8421709430404007e-12
| | Visiting node 3
| | | x = Dict(:volume => 100.0)
| | | x′ = Dict(:volume => 0.0)
| | | C(x, u) = 0.0
| Backward pass
| | Visiting node 2
| | | | V, dVdx′ = 0.0, Dict(:volume => -150.0)
| | | Adding cut : 150 volume_out + cost_to_go ≥ 15000
| | Visiting node 1
| | | | V, dVdx′ = 2.8421709430404007e-12, Dict(:volume => -100.0)
| | | Adding cut : 100 volume_out + cost_to_go ≥ 20000
| Finished iteration
| | lower_bound = 5000.0
| | upper_bound = 5000.0

I tried to see if I could write nested Benders in the fewest possible lines. I’m fairly pleased with this:

using JuMP, HiGHS
function forward_pass(nodes)
    x′s, cost = Vector{Float64}[fix_value.(first.(nodes[1][:states]))], 0.0
    for (t, sp) in enumerate(nodes)
        fix.(first.(sp[:states]), x′s[end]; force = true)
        optimize!(sp)
        assert_is_solved_and_feasible(sp)
        cost += objective_value(sp) - value(sp[:θ])
        push!(x′s, value.(last.(sp[:states])))
    end
    return x′s, cost
end
function backward_pass(nodes, x′s)
    for (sp, x′, sp′) in reverse(collect(zip(nodes, x′s[2:end], nodes[2:end])))
        fix.(first.(sp′[:states]), x′; force = true)
        optimize!(sp′)
        assert_is_solved_and_feasible(sp′; dual = true)
        V, π = objective_value(sp′), reduced_cost.(first.(sp′[:states]))
        @constraint(sp, sp[:θ] >= V + π' * (last.(sp[:states]) .- x′))
    end
    optimize!(nodes[1])
    assert_is_solved_and_feasible(nodes[1])
    return objective_value(nodes[1])
end
function train(nodes::Vector{Model}; atol = 1e-6)
    lb, ub = -Inf, Inf
    while abs(ub - lb) > atol
        x, ub = forward_pass(nodes)
        lb = backward_pass(nodes, x)
        println("Iteration: lower_bound = $lb, upper_bound = $ub")
    end
end
nodes = map(1:3) do t
    sp = Model(HiGHS.Optimizer)
    set_silent(sp)
    @variable(sp, x_in == 200)
    @variable(sp, 0 <= x_out <= 200)
    sp[:states] = [(x_in, x_out)]
    @variable(sp, u[i in [:g, :h]] >= 0)
    @variable(sp, 0 <= θ <= (t == 3 ? 0 : Inf))
    @constraint(sp, x_out <= x_in - u[:h] + 50)
    @constraint(sp, sum(u) == 150)
    @objective(sp, Min, 50 * t * u[:g] + θ)
    return sp
end;
train(nodes)
1 Like

Technically speaking, in your function backward_pass():
1.after optimize!(nodes[1]), there should be an @assert OPTIMAL
2. the returned object should be objective_bound instead of objective_value?

  1. sure. I’ll edit.
  2. no. It’s an LP. The bound is related to the missing cuts that define the value function; not the solution of the first-stage to sub-optimality.

I think the objective_bound is valid. And when it is solved to optimal, the 2 values coincides.
I find a new problem, please see my new topic.Gurobi12 reports obj_value < obj_bound in a Min-Program?