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:
- Solve stage 1
- Solve stage 2
- Solve stage 3
- Add cut to stage 2
- Go to step (3) utill convergence
- Add cut to stage 1
- Go to step (1) until convergence
Or you could do
- Solve stage 1
- Solve stage 2
- Solve stage 3
- Add cut to stage 2
- Add cut to stage 1
- 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?
- sure. I’ll edit.
- 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?