import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import SparseArrays, Random; import LinearAlgebra.diag as diag
function sm(p, x) return sum(p .* x) end;
function optimise(m)
JuMP.optimize!(m)
return (
JuMP.termination_status(m),
JuMP.primal_status(m),
JuMP.dual_status(m)
)
end;
macro set_objective_function(m, f) return esc(:(JuMP.set_objective_function($m, JuMP.@expression($m, $f)))) end;
function Model(name) # e.g. "st2_lms" or "st1_bMD" or "supp_lMs"
m = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
JuMP.MOI.set(m, JuMP.MOI.Name(), name); s = name[end-1]
s == 'm' && JuMP.set_objective_sense(m, JuMP.MIN_SENSE)
s == 'M' && JuMP.set_objective_sense(m, JuMP.MAX_SENSE)
last(name) == 's' && JuMP.set_silent(m)
m
end;
function solve_to_normality(m)
t, p, d = optimise(m)
name = JuMP.name(m)
t == JuMP.OPTIMAL || error("$name: $(string(t))")
p == JuMP.FEASIBLE_POINT || error("$name: (primal) $(string(p))")
(name[end-2] == 'l' && d != p) && error("$name: (dual) $(string(d))")
end;
function safe_upper_bound(ObjBnd) return max(0., ObjBnd + 9.99, 1.01 * ObjBnd) end;
function stage2_problem_primal(z, d)
st2 = Model("st2_primal_lms")
JuMP.@variable(st2, x[i = 1:M, j = 1:N] ≥ 0)
JuMP.@variable(st2, ζ[i = 1:M] ≥ 0) # an expensive substitute for `z`
JuMP.@constraint(st2, DI[i = 1:M], z[i] + ζ[i] ≥ sum(x[i, :]))
JuMP.@constraint(st2, DJ[j = 1:N], sum(x[:, j]) ≥ d[j] )
@set_objective_function(st2, sm(R_x, x) + sm(R_ζ, ζ))
solve_to_normality(st2)
return JuMP.objective_value(st2)
end;
function st2_dual_obj_f(z, d, DI, DJ) return sm(DJ, d) - sm(DI, z) end; # 1st-stage decision; scene; 2nd-stage dual decision
function stage2_problem_dual(z, d) # ✅
st2 = Model("st2_dual_lMs")
JuMP.@variable(st2, 0 <= DI[i = 1:M] <= R_ζ[i]) # c3 and c4 (in R2)
JuMP.@variable(st2, 0 <= DJ[j = 1:N]) # c2 (in R2)
JuMP.@constraint(st2, x[i = 1:M, j = 1:N], R_x[i, j] + DI[i] - DJ[j] >= 0) # c1 (in R2)
@set_objective_function(st2, st2_dual_obj_f(z, d, DI, DJ))
solve_to_normality(st2)
return JuMP.objective_value(st2), JuMP.value.(DI), JuMP.value.(DJ)
end;
function get_st1_oz_cut(z, d) # `d` is the [quasi]-most-hazardous scene
_, DI, DJ = stage2_problem_dual(z, d)
return sm(DJ, d), -DI # cn, pz
end;
function g2d(g) return B_d + V_d .* g end;
function g_2_DI_DJ(z, g) return stage2_problem_dual(z, g2d(g)) end;
function DI_DJ_2_g(z, DI, DJ)
oppo = Model("oppo_lMs")
JuMP.@variable(oppo, 0 ≤ g[1:N] ≤ 1) # r2 and r3 (in R1)
JuMP.@constraint(oppo, C_g_poly * g .≤ b_g_poly) # r1 (in R1)
@set_objective_function(oppo, st2_dual_obj_f(z, g2d(g), DI, DJ))
solve_to_normality(oppo)
return JuMP.objective_value(oppo), JuMP.value.(g)
end;
function improve_g_by_BCA(z, g) # Block Coordinate Ascent
_, DI, DJ = g_2_DI_DJ(z, g)
return lb, g = DI_DJ_2_g(z, DI, DJ)
end;
function get_A1_A2() # A1 is the uncertainty's; A2 is the st2_dual's
A1 = SparseArrays.spzeros(R1, 1 + C1); A2 = SparseArrays.spzeros(R2, 1 + C2)
A1[eachindex(b_g_poly), 1] .= b_g_poly; A1[eachindex(b_g_poly), 2:end] .= -C_g_poly # r1
A1[length(b_g_poly) .+ (1:N), 2:end] .= SparseArrays.spdiagm(ones(N)) # r2
A1[(length(b_g_poly) + N) .+ (1:N), 2:end] .= -SparseArrays.spdiagm(ones(N)) # r3
A1[(length(b_g_poly) + N) .+ (1:N), 1] .= 1 # r3
for i in 1:M, j in 1:N # c1
A2[N*(i-1)+j, [1, 1+ i, (1+M)+ j]] .= [R_x[i, j], 1, -1]
end
A2[M*N .+ (1:N), (1+M) .+ (1:N)] .= SparseArrays.spdiagm(ones(N)); # c2
A2[(M*N + N) .+ (1:M), 1 .+ (1:M)] .= SparseArrays.spdiagm(ones(M)); # c3
A2[(M*N + N + M) .+ (1:M), 1 .+ (1:M)] .= -SparseArrays.spdiagm(ones(M)); # c4
A2[(M*N + N + M) .+ (1:M), 1] .= R_ζ # c4
return A1, A2
end;
# 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎
M, N = 69, 100; # it was (9, 10)
# C_g_poly = SparseArrays.sparse([1, 4, 6, 3, 4, 6, 1, 4, 7, 8, 2, 4, 5, 6, 7, 8, 4, 8, 1, 3, 8, 2, 4, 5, 6, 3, 1, 2, 6, 8, 1, 2, 4, 5, 7, 3, 5, 6, 7, 2, 3, 5, 7, 2, 4, 5, 6, 8, 2, 3, 4, 5, 8, 2, 6, 7, 1, 3, 1, 2, 3, 4, 5, 7, 8, 2, 5, 1, 6, 7, 5, 6, 8, 1, 2, 4, 5, 6, 7, 1, 2, 3, 4, 8, 1, 2, 3, 4, 5, 6, 8, 3, 4, 5, 1, 2, 5, 7, 2, 4, 7, 1, 2, 6, 7, 8, 5, 8, 3, 8, 1, 5, 1, 8, 3, 4, 6, 8, 1, 2, 3, 6, 1, 3, 8, 1, 2, 6, 7, 4, 5, 6, 8, 3, 5, 6, 8, 2, 6, 7, 4, 6, 7, 8, 5, 8, 2, 3, 7, 4, 7, 1, 2, 3, 4, 7, 1, 2, 4, 5, 7, 8, 4, 5, 8, 1, 3, 5, 8, 1, 4, 6, 7, 8, 1, 6, 3, 4, 5, 6, 7, 1, 5, 7, 1, 2, 5, 8, 1, 3, 8, 4, 5, 6, 7, 3, 5, 8, 1, 2, 3, 4, 5, 7, 8, 6, 1, 2, 4, 5, 6, 7, 8, 1, 2, 3, 7, 2, 4, 8, 2, 3, 4, 5, 6, 8, 5, 7, 8, 2, 3, 5, 1, 3, 4, 5, 6, 8, 1, 4, 5, 4, 6, 7, 3, 6, 1, 2, 3, 6, 8, 1, 3, 4, 6, 7, 4, 6, 8, 2, 4, 5, 3, 5, 8, 1, 2, 4, 7, 8, 1, 3, 6, 3, 5, 6, 8, 2, 4, 5, 7, 8, 1, 2, 4, 5, 7, 8, 1, 2, 3, 4, 5, 6, 2, 3, 5, 3, 5, 6, 7, 8, 2, 4, 6, 1, 6, 1, 3, 4, 6, 8, 4, 7, 1, 2, 4, 5, 6, 7, 8, 2, 3, 4, 1, 2, 3, 5, 6, 8, 1, 3, 4, 8, 2, 4, 5, 7, 3, 8, 2, 4, 7, 4, 7, 8, 2, 3, 4, 6, 8, 1, 2, 6, 7, 6, 8, 2, 6, 8, 1, 2, 4, 6, 7, 8, 1, 4, 6, 8, 3, 4, 5, 6, 8], [1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 17, 17, 17, 17, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23, 24, 24, 24, 25, 25, 25, 25, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 32, 32, 33, 33, 33, 33, 34, 34, 34, 35, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 37, 38, 38, 38, 39, 39, 39, 39, 40, 40, 41, 41, 41, 42, 42, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 47, 47, 47, 47, 47, 48, 48, 49, 49, 49, 49, 49, 50, 50, 50, 51, 51, 51, 51, 52, 52, 52, 53, 53, 53, 53, 54, 54, 54, 55, 55, 55, 55, 55, 55, 55, 56, 57, 57, 57, 57, 57, 57, 57, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 61, 61, 62, 62, 62, 63, 63, 63, 64, 64, 64, 64, 64, 64, 65, 65, 65, 66, 66, 66, 67, 67, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 71, 71, 71, 72, 72, 72, 73, 73, 73, 73, 73, 74, 74, 74, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 77, 78, 78, 78, 78, 78, 78, 79, 79, 79, 80, 80, 80, 80, 80, 81, 81, 81, 82, 82, 83, 83, 83, 83, 83, 84, 84, 85, 85, 85, 85, 85, 85, 85, 86, 86, 86, 87, 87, 87, 87, 87, 87, 88, 88, 88, 88, 89, 89, 89, 89, 90, 90, 91, 91, 91, 92, 92, 92, 93, 93, 93, 93, 93, 94, 94, 94, 94, 95, 95, 96, 96, 96, 97, 97, 97, 98, 98, 98, 99, 99, 99, 99, 100, 100, 100, 100, 100], Bool[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], 8, 100);
# b_g_poly = [40.5, 35.88, 37.41, 46.8, 38.64, 41.16, 32.4, 48.6];
C_g_poly = ones(1, N)
b_g_poly = N/2 * ones(1)
B_d = [0.11, 0.57, 0.77, 0.35, 0.49, 0.25, 0.25, 0.45, 0.05, 0.77, 0.39, 0.83, 0.47, 0.55, 0.43, 0.47, 0.29, 0.49, 0.47, 0.75, 0.09, 0.59, 0.11, 0.09, 0.11, 0.49, 0.29, 0.61, 0.23, 0.53, 0.15, 0.85, 0.17, 0.23, 0.25, 0.57, 0.65, 0.47, 0.65, 0.11, 0.09, 0.13, 0.57, 0.29, 0.19, 0.39, 0.31, 0.63, 0.61, 0.51, 0.75, 0.41, 0.57, 0.59, 0.47, 0.05, 0.29, 0.49, 0.53, 0.79, 0.21, 0.83, 0.13, 0.31, 0.63, 0.43, 0.05, 0.33, 0.77, 0.43, 0.11, 0.41, 0.19, 0.61, 0.89, 0.43, 0.15, 0.35, 0.43, 0.63, 0.59, 0.89, 0.23, 0.23, 0.33, 0.63, 0.69, 0.35, 0.79, 0.31, 0.69, 0.81, 0.15, 0.19, 0.27, 0.29, 0.17, 0.83, 0.85, 0.45];
V_d = [0.56, 0.26, 0.5, 0.5, 0.71, 0.77, 0.32, 0.53, 0.38, 0.2, 0.2, 0.29, 0.68, 0.59, 0.38, 0.29, 0.38, 0.41, 0.26, 0.35, 0.56, 0.68, 0.8, 0.8, 0.62, 0.44, 0.47, 0.53, 0.68, 0.44, 0.26, 0.35, 0.44, 0.59, 0.59, 0.41, 0.53, 0.8, 0.29, 0.62, 0.29, 0.59, 0.23, 0.56, 0.68, 0.2, 0.44, 0.74, 0.65, 0.62, 0.8, 0.71, 0.35, 0.56, 0.23, 0.44, 0.65, 0.35, 0.38, 0.38, 0.44, 0.47, 0.5, 0.38, 0.74, 0.29, 0.2, 0.29, 0.62, 0.23, 0.41, 0.8, 0.77, 0.2, 0.35, 0.65, 0.32, 0.23, 0.74, 0.38, 0.68, 0.56, 0.59, 0.32, 0.8, 0.74, 0.38, 0.2, 0.5, 0.68, 0.74, 0.8, 0.68, 0.53, 0.44, 0.59, 0.26, 0.68, 0.44, 0.38];
# @assert sum(K_y) >= sum(g2d(ones(100)))
K_y = [2.42, 2.42, 1.78, 2.7, 2.06, 1.78, 1.14, 3.5, 3.1, 2.02, 1.82, 3.42, 2.42, 0.9, 3.06, 1.18, 1.7, 3.46, 1.46, 1.62, 1.26, 1.3, 3.46, 2.66, 1.82, 2.34, 1.54, 1.82, 2.58, 2.9, 2.66, 3.14, 2.5, 1.62, 2.78, 3.3, 1.38, 1.26, 2.74, 1.82, 3.14, 1.38, 2.74, 0.94, 3.06, 3.22, 2.62, 2.42, 3.02, 3.18, 1.58, 2.7, 3.06, 1.7, 3.22, 2.74, 2.54, 3.26, 1.3, 2.38, 3.3, 1.02, 2.42, 2.38, 3.1, 1.3, 2.02, 0.94, 3.46];
Random.seed!(1); R_x = rand(0.3:.00117:8.7, M, N);
C_z = rand(2:0.00117:5, M);
C_y = rand(9:0.017:15, M) .* C_z;
R_ζ = 20 * C_z;
R1 = length(b_g_poly) + N + N; C1 = N; # 1 ⇒ uncertainty_side, 2 ⇒ st2_dual side
R2 = M*N + N + M + M; C2 = M + N;
# 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎
st1 = Model("st1_bms"); JuMP.@variable(st1, st1_o);
JuMP.@variables(st1, begin st1_y[1:M], Bin; st1_z[1:M] ≥ 0 end); JuMP.@constraint(st1, st1_z .≤ K_y .* st1_y);
JuMP.set_objective_coefficient(st1, [st1_y; st1_z], [C_y; C_z])
solve_to_normality(st1); z = JuMP.value.(st1_z); # 🍏 initial trial `z`
# 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎 🥎
best_lpr_ul = zeros(2); best_g = zeros(N); # 📙 global book
A1, A2 = get_A1_A2(); old_feas_val = zeros(1); # if the heuristic solution of BLP reaches a different value, then we deem it new
lpr = Model("lpr_lMD");
JuMP.@variables(lpr, begin lpr_g[1:C1]; lpr_Dv[1:C2]; g_x_Dv[1:C1, 1:C2] end);
JuMP.@expression(lpr, cross_multrix, [1 lpr_Dv'; lpr_g g_x_Dv]); JuMP.@expression(lpr, lpr_DI, lpr_Dv[1:M]);
JuMP.@expression(lpr, lpr_DJ, lpr_Dv[M .+ (1:N)]); JuMP.@expression(lpr, g_x_DJ, g_x_Dv[:, M .+ (1:N)]);
JuMP.@constraint(lpr, [A1 * [1; lpr_g]; A2 * [1; lpr_Dv]] .>= 0); # [uncertainty_side; st2_dual_side] constraints
ite = 0
t0 = time()
for r in eachrow(A1)
ite += 1
println("ite = $ite <= 201, time_elapsed = $(time() - t0)")
JuMP.@constraint(lpr, [c in eachrow(A2)], transpose(r) * cross_multrix * c >= 0)
end
The bottom for-loop is slow.