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:
- 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
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?
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.