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.

1 Like

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

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?

cost_1 = [10858, 10686, 11845]
g = [57, 88, 31]
capacity = [63, 61, 123]
cost_2 = [[0 7 7]
          [7 0 5]            
          [9 8 0]
]

model = Model(HiGHS.Optimizer)
set_silent(model)

# Master Problem-I 
@variable(model, x[1:O, 1:O], Bin)
@variable(model, t >=-100)
@constraint(model, [i in 1:O], sum(x[i, j] for j in 1:O) == 1)
@constraint(model, [i in 1:O, j in 1:O], x[i, j] <= x[j,j])

RMP_cost = sum(cost_1[j] * x[j, j] for j in 1:O) 

@objective(model, Min, RMP_cost + sum(t))

subproblem = Model(HiGHS.Optimizer)
set_attribute(subproblem, "presolve", "off")
set_silent(subproblem)
@variable(subproblem, x_copy[i in 1:O, j in 1:O])
@variable(subproblem, y[1:O, 1:O] >=0)

@constraint(subproblem, con1[i in 1:O, j in 1:O], y[i, j] == x_copy[i, j] * g[i]) 
@constraint(subproblem, con2[j in 1:O], sum(y[i, j] for i in 1:O) <= capacity[j])

SP_cost = sum(cost_2[i, j] * y[i, j] for i in 1:O for j in 1:O)

@objective(subproblem, Min, SP_cost)

function solve_subproblem(model, x)
    fix.(model[:x_copy], x)
    optimize!(model)
    if is_solved_and_feasible(model; dual = true)
        return ( is_feasible = true,
        obj = objective_value(model),
        y = value.(model[:y]),
        π = reduced_cost.(model[:x_copy]),
        d1 = dual.(model[:con1]),
        d2 = dual.(model[:con2])
        )
    end
    return (
            is_feasible = false, df1 = dual.(con1), df2 = dual.(con2)
        )
        
end

function print_iteration(k, args...)
    f(x) = Printf.@sprintf("%12.4e", x)
    println(lpad(k, 9), " ", join(f.(args), " "))
    return
end

println("Iteration  Lower Bound  Upper Bound          Gap")
for k in 1:10
    optimize!(model)
    t_k = value.(t)
    x_k = value.(x)
    lower_bound = t_k
    ret = solve_subproblem(subproblem, x_k)
    if ret.is_feasible
        upper_bound = ret.obj
        gap = abs(upper_bound - lower_bound) / abs(upper_bound)
        print_iteration(k, lower_bound, upper_bound, gap)
        if gap < ABSOLUTE_OPTIMALITY_GAP
            println("Terminating with the optimal solution")
            break
        end
        cut = @constraint(model, θ >= sum(ret.d2[j] * capacity[j] for j in 1:O) + 
            sum(sum(ret.d1[i, j] * x[i, j] for j in 1:O) * g[i] for i in 1:O)
        )
        @info "Adding the optimality cut $(cut) at iter $k"
    else
        fcut = @constraint(model, 0 >= sum(ret.df2[j] * capacity[j] for j in 1:O) + 
            sum(sum(ret.df1[i, j] * x[i, j] for j in 1:O) * g[i] for i in 1:O)
        )
        @info "Adding the feasibility cut $(fcut) at iter $k"
    end
end

Could you please give it a look (two-stage decomposition)? why same feasible cuts are generated in all iterations? It works fine when iterative method is used [sub problem in function] but not when I use In-place iterative method. I can’t find out the error that causes this problem.

Thank you.

Iteration  Lower Bound  Upper Bound          Gap
[ Info: Adding the feasibility cut -57 x[1,2] - 88 x[2,2] >= -61 at iter 1
[ Info: Adding the feasibility cut -57 x[1,2] - 88 x[2,2] >= -61 at iter 2
[ Info: Adding the feasibility cut -57 x[1,2] - 88 x[2,2] >= -61 at iter 3
[ Info: Adding the feasibility cut -57 x[1,2] - 88 x[2,2] >= -61 at iter 4
[ Info: Adding the feasibility cut -57 x[1,2] - 88 x[2,2] >= -61 at iter 5
[ Info: Adding the feasibility cut -57 x[1,2] - 88 x[2,2] >= -61 at iter 6
[ Info: Adding the feasibility cut -57 x[1,2] - 88 x[2,2] >= -61 at iter 7
[ Info: Adding the feasibility cut -57 x[1,2] - 88 x[2,2] >= -61 at iter 8
[ Info: Adding the feasibility cut -57 x[1,2] - 88 x[2,2] >= -61 at iter 9
[ Info: Adding the feasibility cut -57 x[1,2] - 88 x[2,2] >= -61 at iter 10

Try import Pkg; Pkg.pkg"add HiGHS_jll@1.8.0" then restart Julia.

There is a bug in the v1.9.0 release of HiGHS that will be fixed in next release.

1 Like

Thank you, I switched to Gurobi and it’s solved.

1 Like

Following up again: I’m about to release support for HiGHS v1.10, which will fix this bug: Add support for HiGHS@1.10.0 by odow · Pull Request #270 · jump-dev/HiGHS.jl · GitHub

Seems Gurobi supports this now.