Adaptive robust optimization (ARO) is typically NP-hard.
It’s known that in the RHS uncertainty setting, the subproblem constitutes a disjointed bilinear programming (DBP).
A naive way to solve DBP is employing Gurobi’s MINLP capability.
Nonetheless, the computational performance is unsatisfying.
Relaxation linearization technique (RLT) might be applicable, but it can be questionable when we do NOT have relatively complete recourse (RCR).
The reason is told in this tutorial, in which the example is from doi.org/10.1016/j.orl.2013.05.003.
import JuMP, MosekTools, Ipopt
import Gurobi; const GRB_ENV = Gurobi.Env();
import LinearAlgebra.tr as tr; function sm(p, x) return sum(p .* x) end
# The inappropriateness of applying RLT to solve the DBP in ARO who lacks RCR
function optimise(m) return (JuMP.optimize!(m); (JuMP.termination_status(m), JuMP.primal_status(m))) end
macro set_objective_function(m, f) return esc(:(JuMP.set_objective_function($m, JuMP.@expression($m, $f)))) end
function Model(solver, mode)
solver == "Ipopt" && (m = JuMP.Model(Ipopt.Optimizer))
solver == "Mosek" && (m = JuMP.Model(() -> MosekTools.Optimizer()))
solver == "Gurobi" && (m = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)))
occursin("m", mode) && JuMP.set_objective_sense(m, JuMP.MIN_SENSE)
occursin("M", mode) && JuMP.set_objective_sense(m, JuMP.MAX_SENSE)
occursin("s", mode) && JuMP.set_silent(m)
m
end
macro opt_ass_opt(m)
n = string(m, ": ")
return esc(quote
let (s, p) = optimise($m)
(
s == ifelse(JuMP.solver_name($m) == "Ipopt", JuMP.LOCALLY_SOLVED, JuMP.OPTIMAL)
&& p == JuMP.FEASIBLE_POINT
) || error($n * string(s))
end
end)
end
# 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎
C_x, C_y, C_z, B_d = let # Data
[22 33 24; 33 23 30; 20 25 27], [400, 414, 326], [18, 25, 20], [206, 274, 220]
end
st1 = Model("Mosek", "ms") # The 1st-stage problem who picks a trial `z`
JuMP.@variable(st1, st1_y[1:3], Bin)
JuMP.@variable(st1, st1_z[1:3] >= 0); JuMP.@constraint(st1, st1_z .<= 800 * st1_y)
@set_objective_function(st1, sm(C_y, st1_y) + sm(C_z, st1_z))
function stage2_problem_primal(z, d)
st2 = Model("Mosek", "ms")
JuMP.@variable(st2, x[i = 1:3, j = 1:3] >= 0)
JuMP.@constraint(st2, DI[i = 1:3], z[i] >= sum(x[i, :]))
JuMP.@constraint(st2, DJ[j = 1:3], sum(x[:, j]) >= d[j] )
@set_objective_function(st2, sm(C_x, x))
@opt_ass_opt(st2)
return JuMP.objective_value(st2)
end
function stage2_problem_dual(z, d)
Dst2 = Model("Mosek", "Ms")
JuMP.@variable(Dst2, DI[i = 1:3] >= 0)
JuMP.@variable(Dst2, DJ[j = 1:3] >= 0)
JuMP.@constraint(Dst2, x[i = 1:3, j = 1:3], C_x[i, j] + DI[i] - DJ[j] >= 0)
@set_objective_function(Dst2, sm(DJ, d) - sm(DI, z))
@opt_ass_opt(Dst2)
return JuMP.objective_value(Dst2)
end
function slacked_stage2_problem_primal(z, d)
sl2 = Model("Mosek", "ms")
JuMP.@variable(sl2, s[i = 1:3] >= 0) # nonnegative slack variable
JuMP.@variable(sl2, x[i = 1:3, j = 1:3] >= 0)
JuMP.@constraint(sl2, Ds[i = 1:3], z[i] + s[i] >= sum(x[i, :]))
JuMP.@constraint(sl2, DJ[j = 1:3], sum(x[:, j]) >= d[j] )
@set_objective_function(sl2, sum(s))
@opt_ass_opt(sl2)
return JuMP.objective_value(sl2)
end
function slacked_stage2_problem_dual_objf(z, d, DJ, Ds) return sm(DJ, d) - sm(Ds, z) end
function slacked_stage2_problem_dual(z, d)
Dsl2 = Model("Mosek", "Ms")
JuMP.@variable(Dsl2, 0 <= Ds[i = 1:3] ≤ 1) # "≤ 1" ⇐ [s >= 0]
JuMP.@variable(Dsl2, DJ[j = 1:3] >= 0)
JuMP.@constraint(Dsl2, x[i = 1:3, j = 1:3], Ds[i] - DJ[j] >= 0)
@set_objective_function(Dsl2, slacked_stage2_problem_dual_objf(z, d, DJ, Ds))
@opt_ass_opt(Dsl2)
objval = JuMP.objective_value(Dsl2)
cn = JuMP.value(sm(DJ, d))
pz = JuMP.value.(-Ds)
return objval, cn, pz
end
function slacked_subproblem(z) # cheat via Gurobi
BL = Model("Gurobi", "Ms") # disjointed bilinear program
JuMP.@variable(BL, 0 <= g[1:3] <= 1)
JuMP.@constraint(BL, g[1] + g[2] <= 1.2); JuMP.@constraint(BL, sum(g) <= 1.8);
JuMP.@expression(BL, d, B_d .+ 40 * g);
# 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎
JuMP.@variable(BL, 0 <= Ds[i = 1:3] ≤ 1) # "≤ 1" ⇐ [s >= 0]
JuMP.@variable(BL, DJ[j = 1:3] >= 0)
JuMP.@constraint(BL, x[i = 1:3, j = 1:3], Ds[i] - DJ[j] >= 0)
# 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎
@set_objective_function(BL, slacked_stage2_problem_dual_objf(z, d, DJ, Ds))
@opt_ass_opt(BL)
@info JuMP.value.(d) # the most adverse scene
JuMP.objective_bound(BL) # if this >0, then the arg `z` should then be cut off by a feas. cut
end
function slacked_subproblem(z, g0, Ds0, DJ0) # suffix `0` means start values for the variables
BL = Model("Ipopt", "Ms")
JuMP.@variable(BL, 0 <= g[1:3] <= 1)
JuMP.@constraint(BL, g[1] + g[2] <= 1.2); JuMP.@constraint(BL, sum(g) <= 1.8);
JuMP.@expression(BL, d, B_d .+ 40 * g);
# 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎
JuMP.@variable(BL, 0 <= Ds[i = 1:3] ≤ 1) # "≤ 1" ⇐ [s >= 0]
JuMP.@variable(BL, DJ[j = 1:3] >= 0)
JuMP.@constraint(BL, x[i = 1:3, j = 1:3], Ds[i] - DJ[j] >= 0)
# 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎
@set_objective_function(BL, slacked_stage2_problem_dual_objf(z, d, DJ, Ds))
JuMP.set_start_value.(g, g0)
JuMP.set_start_value.(Ds, Ds0)
JuMP.set_start_value.(DJ, DJ0)
@opt_ass_opt(BL) # local optimal
JuMP.objective_value(BL), JuMP.value.(d) # Ipopt can only provide a primal value
end
function slacked_subproblem(z, ::String) # TODO the dual bound is very loose currently due to partial RLT
CR = Model("Mosek", "Ms")
JuMP.@variable(CR, 0 <= g[1:3] <= 1)
JuMP.@constraint(CR, g[1] + g[2] <= 1.2); JuMP.@constraint(CR, sum(g) <= 1.8);
JuMP.@expression(CR, d, B_d .+ 40 * g);
# 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎
JuMP.@variable(CR, 0 <= Ds[i = 1:3] ≤ 1) # "≤ 1" ⇐ [s >= 0]
JuMP.@variable(CR, DJ[j = 1:3] >= 0)
JuMP.@constraint(CR, x[i = 1:3, j = 1:3], Ds[i] - DJ[j] >= 0)
# 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎
JuMP.@variable(CR, 0 <= gDJ[1:3, 1:3] <= 1)
@set_objective_function(CR, sm(DJ, B_d) + 40 * tr(gDJ) - sm(Ds, z))
@opt_ass_opt(CR)
db = JuMP.objective_bound(CR)
pb, dt = slacked_subproblem(z, JuMP.value.(g), JuMP.value.(Ds), JuMP.value.(DJ)) # via Ipopt
@info "lb = $pb < $db = ub"
return pb, db, dt
end
@opt_ass_opt(st1)
zt = JuMP.value.(st1_z) # generate the first trial `z`
pb, db, dt = slacked_subproblem(zt, "CR") # we see lb = 772, thus this `zt` need to be cut off by a feas. cut
# the following 2 lines add the feas. cut
_, cn, cz = slacked_stage2_problem_dual(zt, dt)
JuMP.@constraint(st1, cn + sm(cz, st1_z) <= 0)
@opt_ass_opt(st1)
zt = JuMP.value.(st1_z) # generate the second trial `z`
pb, db, dt = slacked_subproblem(zt, "CR") # we see lb = 7.653811934460464e-6 < 120.0 = ub
# this result is uninformative in that
# 1. lb is not strictly positive
# 2. ub is not zero
# 3. since the gap is not closed, `dt` may not be optimal
# By comparison
slacked_subproblem(zt) # if we cheat via Gurobi, we get an ub = 0.000238, which in this example can be deemed zero
# This suggest that this trial `zt` is essentially feasible, then we could leave the "slacked" 2nd-stage problem
# TODO write the 3 `subproblem::Function` resembling `slacked_subproblem`
# use the dual bound to update upper bound, and use the primal feas. solution to add opt. cut to update lower bound