An unexpected new finding is that This Warning persists, even if I opt not to use mapreduce. (I only used LinearAlgebra.dot)
The whole code have 310 lines (copy paste to REPL if you want to try). The Warning should appear within 2 minutes.
import JuMP
import Gurobi
import Random
import Distributions: Uniform as UDistri
import LinearAlgebra: dot as ip
import LinearAlgebra: norm as norm
const UPDTOL, argZ_TimeLimit = 0.01, 8.0
GRB_ENV = Gurobi.Env()
Random.seed!(3)
function GurobiDirectModel() return JuMP.direct_model(Gurobi.Optimizer(GRB_ENV)) end
function jvr6(x) return round(JuMP.value(x); digits = 6) end
function jdr6(x) return round( JuMP.dual(x); digits = 6) end
function optimise(m) return (JuMP.optimize!(m); JuMP.termination_status(m)) end
# function my_ip(x::AbstractVector{JuMP.VariableRef}, y)
# JuMP.@expression(JuMP.owner_model(first(x)), sum(xi * yi for (xi, yi) in zip(x, y)))
# end
function opt_ass_opt(m, str)
sta = optimise(m)
sta == JuMP.OPTIMAL || error(str * ": " * string(sta))
end
function opt_ass_optortime(m, str)
sta = optimise(m)
sta == JuMP.OPTIMAL && return
sta != JuMP.TIME_LIMIT && error(str * ": " * string(sta))
JuMP.has_values(m) || error(str * ": no primal solution")
end
function brcs(v) return ones(T) * transpose(v) end
function gen_load_pattern(fLM)
fL = length(fLM)
fLP = rand(UDistri(0.7 * fLM[1], 1.4 * fLM[1]), T, 1)
for i in 2:fL
fLP = [fLP rand(UDistri(0.7 * fLM[i], 1.4 * fLM[i]), T)]
end
return fLP
end
begin # base case: such that the first stage solution is to turn on all the generators
G = 2;
UT = DT = 3;
PI = [0.45, 0.375, 0.5];
PS = [4.5, 4, 5.5];
ZS = [0, 0, 1];
CST = [0.63, 0.60, 0.72];
CG = [0.8, 0.68, 0.72];
CL_BASE = [1.6, 1.3776, 1.4442];
NW = [2, 3];
W = 2;
NL = [4, 5, 6];
L = 3;
LM = [4.0, 3.5, 3.0];
T = 24;
(CL = brcs(CL_BASE); CL[end, :] *= 2);
end
MZ = gen_load_pattern(LM)
MY, yM, NY = let
YABSMAX = [1.6, 2.3]
Pr = UDistri.(0.2 * rand(W), YABSMAX / 1.1)
Ybaseline = [rand(Pr[w]) for t in 1:T, w in 1:W]
YMIN, YMAX = 0.9 * Ybaseline, 1.1 * Ybaseline # used to construct vertices
all0() = Bool.(zero(YMIN))
vertexList = [all0() for _ in eachindex(YMIN)]
for i in eachindex(YMIN)
vertexList[i][i] = true
end
pushfirst!(vertexList, all0())
for i in 1:1+length(YMIN)
push!(vertexList, .!(vertexList[i]))
end
ref = vertexList[1]
yM = [(ref[i,j] ? YMAX[i,j] : YMIN[i,j]) for i in eachindex(eachrow(YMIN)), j in eachindex(eachcol(YMIN))]
yM = reshape(yM, (T, W, 1))
popfirst!(vertexList)
for ref in vertexList
yM = cat(yM, [(ref[i,j] ? YMAX[i,j] : YMIN[i,j]) for i in eachindex(eachrow(YMIN)), j in eachindex(eachcol(YMIN))]; dims = 3)
end # âś… after this line, yM is already gened
NY = size(yM, 3)
PrVec = rand(UDistri(.4, .6), NY)
PrVec = PrVec / sum(PrVec)
MY = sum(yM[:, :, i] * PrVec[i] for i in 1:NY)
@assert all(MY .< YMAX) && all(MY .> YMIN)
function is_strict_in(MY, yM)
NY = size(yM, 3)
ø = GurobiDirectModel()
JuMP.set_silent(ø)
JuMP.@variable(ø, c[1:NY] >= 1e-5)
JuMP.@constraint(ø, sum(c) == 1)
JuMP.@constraint(ø, sum(yM[:, :, i] * c[i] for i in 1:NY) .== MY)
status = optimise(ø)
status == JuMP.OPTIMAL && return true
return false
end
@assert is_strict_in(MY, yM)
MY, yM, size(yM)[3]
end;
Z_list = [MZ]; # must involve MZ to ensure feasibility
function recruit(Z_list, Z_new)
for Z in Random.shuffle(Z_list)
norm(Z .- Z_new, Inf) < 5e-5 && return false
end
push!(Z_list, Z_new); return true
end
trh = GurobiDirectModel(); JuMP.set_objective_sense(trh, JuMP.MIN_SENSE); JuMP.@variable(trh, oLh); JuMP.set_silent(trh);
JuMP.@variable(trh, trh_u[t = 1:T, g = 1:G+1], Bin);
JuMP.@variable(trh, trh_v[t = 1:T, g = 1:G+1], Bin);
JuMP.@variable(trh, trh_x[t = 1:T, g = 1:G+1], Bin);
let # 1st stage feasible region
JuMP.@constraint(trh, trh_x .- [transpose(ZS); trh_x[1:T-1, :]] .== trh_u .- trh_v)
JuMP.@constraint(trh, [g = 1:G+1, t = 1:T-UT+1], sum(trh_x[i, g] for i in t:t+UT-1) >= UT * trh_u[t, g])
JuMP.@constraint(trh, [g = 1:G+1, t = T-UT+1:T], sum(trh_x[i, g] - trh_u[t, g] for i in t:T) >= 0)
JuMP.@constraint(trh, [g = 1:G+1, t = 1:T-DT+1], sum(1 - trh_x[i, g] for i in t:t+DT-1) >= DT * trh_v[t, g])
JuMP.@constraint(trh, [g = 1:G+1, t = T-DT+1:T], sum(1 - trh_x[i, g] - trh_v[t, g] for i in t:T) >= 0)
end;
u, v, x = let
JuMP.set_objective_function(trh, ip(brcs(CST), trh_u)) # to derive a feas (u, v, x)
opt_ass_opt(trh, "trh")
u = jvr6.(trh_u);
v = jvr6.(trh_v);
x = jvr6.(trh_x);
u, v, x
end
fx2Z = GurobiDirectModel(); JuMP.set_objective_sense(fx2Z, JuMP.MAX_SENSE); JuMP.set_silent(fx2Z);
f_c_x, f_c_Y, f_c_Z = let
JuMP.@variable(fx2Z, Dbal[t = 1:T]);
JuMP.@variable(fx2Z, Dvp[t = 1:T, w = 1:W] >= 0);
JuMP.@variable(fx2Z, Dzt[t = 1:T, l = 1:L] >= 0);
JuMP.@variable(fx2Z, Dvr[t = 1:T, g = 1:G+1] >= 0);
JuMP.@variable(fx2Z, Dpi[t = 1:T, g = 1:G+1] >= 0);
JuMP.@variable(fx2Z, Dps[t = 1:T, g = 1:G+1] >= 0);
JuMP.@constraint(fx2Z, [t = 1:T, g = 1:G+1], Dps[t, g] - Dpi[t, g] - Dvr[t, g] + CG[g] == 0);
JuMP.@constraint(fx2Z, [t = 1:T, g = 1:G+1], Dbal[t] - CG[g] + Dvr[t, g] >= 0);
JuMP.@constraint(fx2Z, [t = 1:T, w = 1:W], Dvp[t, w] + Dbal[t] >= 0);
JuMP.@constraint(fx2Z, [t = 1:T, l = 1:L], Dzt[t, l] - CL[t, l] - Dbal[t] >= 0);
JuMP.@expression(fx2Z, f_c_x, brcs(PI) .* Dpi .- brcs(PS) .* Dps); # in x2
JuMP.@expression(fx2Z, f_c_Y, -Dvp); # in x2
JuMP.@expression(fx2Z, f_c_Z, CL .- Dzt)
f_c_x, f_c_Y, f_c_Z
end;
function set_f_objective(x, Y, Z)
JuMP.set_objective_function(fx2Z,
ip(f_c_x, x)
+ ip(f_c_Y, Y)
+ ip(f_c_Z, Z)
)
end
function get_f_value(x, Y, Z)
set_f_objective(x, Y, Z)
opt_ass_opt(fx2Z, "fx2Z")
JuMP.objective_value(fx2Z)
end
function get_f_cut(x, Y, Z)
set_f_objective(x, Y, Z)
opt_ass_opt(fx2Z, "fx2Z")
cn = ip(jvr6.(f_c_Z), Z) + ip(jvr6.(f_c_Y), Y)
px = jvr6.(f_c_x)
return cn, px
end
function add_cut_for_oLh(cn, px)
JuMP.@constraint(trh, oLh >= cn + ip(px, trh_x))
end
function get_cut_oLh_coeff(x, Q, P)
cn, px = 0., zero(x)
for j in 1:Q[1]
for i in 1:P[1]
p = Q[2][j] * P[2][i]
this_cn, this_px = get_f_cut(x, Q[3][:, :, j], P[3][i])
cn += p * this_cn
px += p * this_px
end
end
return cn, px
end
function add_cut_for_oLh(x, Q, P, oLh_check)
cn, px = get_cut_oLh_coeff(x, Q, P)
if oLh_check < cn + ip(px, x) - UPDTOL
add_cut_for_oLh(cn, px)
return true
end
return false
end
argZ = GurobiDirectModel(); JuMP.set_objective_sense(argZ, JuMP.MAX_SENSE); JuMP.set_attribute(argZ, "TimeLimit", argZ_TimeLimit); JuMP.set_silent(argZ);
JuMP.@variable(argZ, argZ_Z[t = 1:T, l = 1:L] == MZ[t, l]); # TBD
argZ_c_x, argZ_c_Y, argZ_c_b2, argZ_cn = let
JuMP.@variable(argZ, Dbal[t = 1:T]);
JuMP.@variable(argZ, Dvp[t = 1:T, w = 1:W] >= 0);
JuMP.@variable(argZ, Dzt[t = 1:T, l = 1:L] >= 0);
JuMP.@variable(argZ, Dvr[t = 1:T, g = 1:G+1] >= 0);
JuMP.@variable(argZ, Dpi[t = 1:T, g = 1:G+1] >= 0);
JuMP.@variable(argZ, Dps[t = 1:T, g = 1:G+1] >= 0);
JuMP.@constraint(argZ, [t = 1:T, g = 1:G+1], Dps[t, g] - Dpi[t, g] - Dvr[t, g] + CG[g] == 0);
JuMP.@constraint(argZ, [t = 1:T, g = 1:G+1], Dbal[t] - CG[g] + Dvr[t, g] >= 0);
JuMP.@constraint(argZ, [t = 1:T, w = 1:W], Dvp[t, w] + Dbal[t] >= 0);
JuMP.@constraint(argZ, [t = 1:T, l = 1:L], Dzt[t, l] - CL[t, l] - Dbal[t] >= 0);
JuMP.@expression(argZ, argZ_c_x, brcs(PI) .* Dpi .- brcs(PS) .* Dps); # in x2
JuMP.@expression(argZ, argZ_c_Y, -Dvp); # in x2
JuMP.@expression(argZ, argZ_c_b2, -argZ_Z); # is beta2
JuMP.@expression(argZ, argZ_cn, ip(CL .- Dzt, argZ_Z)); # the rest
argZ_c_x, argZ_c_Y, argZ_c_b2, argZ_cn
end;
function set_argZ_objective(x, Y, b2) # (x2, b2)
JuMP.set_objective_function(argZ, argZ_cn
+ ip(argZ_c_x, x)
+ ip(argZ_c_Y, Y)
+ ip(argZ_c_b2, b2)
)
end
tr2 = GurobiDirectModel(); JuMP.@variable(tr2, oΛ2); JuMP.set_silent(tr2); # ✅
JuMP.@variable(tr2, tr2_b2[eachindex(eachrow(MZ)), eachindex(eachcol(MZ))]);
JuMP.@variable(tr2, tr2_x[t = 1:T, g = 1:G+1] == x[t, g]);
JuMP.@variable(tr2, tr2_Y[t = 1:T, w = 1:W] == MY[t, w]);
function add_cut_for_oΛ2(cn, px, pY, pb)
JuMP.@constraint(tr2, oΛ2 >= cn + ip(px, tr2_x) + ip(pY, tr2_Y) + ip(pb, tr2_b2))
end
let
set_argZ_objective(x, (local Y = MY), (local b2 = zero(MZ)));
opt_ass_optortime(argZ, "argZ")
cn, px, pY = JuMP.value(argZ_cn), jvr6.(argZ_c_x), jvr6.(argZ_c_Y)
add_cut_for_oΛ2(cn, px, pY, (local pb2 = -JuMP.fix_value.(argZ_Z)))
JuMP.@objective(tr2, Min, ip(MZ, tr2_b2) + oΛ2); # ❄️ stable obj
try
opt_ass_opt(tr2, "tr2")
finally # finish establishing Z's support
JuMP.unfix.(argZ_Z)
(JuMP.set_lower_bound.(argZ_Z, 0.9 * MZ); JuMP.set_upper_bound.(argZ_Z, 1.1 * MZ))
JuMP.@variable(argZ, argZ_a[eachindex(eachrow(MZ)), eachindex(eachcol(MZ))])
JuMP.@constraint(argZ, argZ_a .>= argZ_Z .- MZ)
JuMP.@constraint(argZ, argZ_a .>= MZ .- argZ_Z)
JuMP.@constraint(argZ, [t = 1:T], sum(argZ_a[t, :] ./ (0.1 * MZ[t, :])) <= L/3)
end
end
function psi_via_19(x, i)
JuMP.fix.(tr2_x, x)
JuMP.fix.(tr2_Y, yM[:, :, i])
opt_ass_opt(tr2, "tr2")
JuMP.objective_value(tr2)
end
function get_proper_Q(x)
JuMP.set_objective_function(fr9, ip(fr9_p, [psi_via_19(x, i) for i in 1:NY]));
opt_ass_opt(fr9, "fr9")
prob_Y = jvr6.(fr9_p)
itmp = prob_Y .> 1e-6
proper_Q_p = prob_Y[itmp]
proper_Q_n = length(proper_Q_p)
proper_Q_Y = yM[:, :, itmp]
proper_Q_n, proper_Q_p, proper_Q_Y
end
function get_proper_P(x, Y, Z_list)
fr10 = GurobiDirectModel(); JuMP.set_silent(fr10); JuMP.set_objective_sense(fr10, JuMP.MAX_SENSE); # JuMP.set_attribute(fr10, "Presolve", 0);
JuMP.@variable(fr10, p[eachindex(Z_list)] >= 0); JuMP.@constraint(fr10, sum(p) == 1);
# JuMP.@constraint(fr10, my_ip(p, Z_list) .== MZ) # this style triggers Warning
JuMP.@constraint(fr10, sum(p[i] * Z_list[i] for i in eachindex(p)) .== MZ) # this equivalent style also triggers
JuMP.set_objective_function(fr10, ip(p, [get_f_value(x, Y, Z) for Z in Z_list]));
opt_ass_opt(fr10, "fr10")
prob_Z = jvr6.(p)
itmp = prob_Z .> 1e-6
proper_P_p = prob_Z[itmp]
proper_P_n = length(proper_P_p)
proper_P_Z = Z_list[itmp]
proper_P_n, proper_P_p, proper_P_Z
end
fr9 = GurobiDirectModel(); JuMP.set_objective_sense(fr9, JuMP.MAX_SENSE); JuMP.set_silent(fr9);
JuMP.@variable(fr9, fr9_p[1:NY] >= 0); JuMP.@constraint(fr9, sum(fr9_p) == 1);
# JuMP.@constraint(fr9, my_ip(fr9_p, [yM[:, :, i] for i in 1:size(yM, 3)]) .== MY);
JuMP.@constraint(fr9, sum(fr9_p[i] * yM[:, :, i] for i in 1:NY) .== MY);
proper_Q = get_proper_Q(x) # âś…
JuMP.fix.(tr2_x, x);
i = 1; # for one Y
Y = proper_Q[3][:, :, i];
JuMP.fix.(tr2_Y, Y);
opt_ass_opt(tr2, "tr2");
b2 = jvr6.(tr2_b2); # 🥑
set_argZ_objective(x, Y, b2);
opt_ass_optortime(argZ, "argZ");
Z_new = jvr6.(argZ_Z);
recruit(Z_list, Z_new);
proper_P = get_proper_P(x, Y, Z_list) # âś…
cn, px = get_cut_oLh_coeff(x, proper_Q, proper_P); add_cut_for_oLh(cn, px)
JuMP.set_objective_function(trh, ip(brcs(CST), trh_u) + oLh) # ❄️ stable obj
opt_ass_opt(trh, "trh")
@info "trh's initialization succeed, main loop begins"
updVec = trues(2)
while true
opt_ass_opt(trh, "trh")
x = jvr6.(trh_x); oLh_check = JuMP.value(oLh); lb = JuMP.objective_bound(trh);
JuMP.fix.(tr2_x, x) # 🖌️
proper_Q = get_proper_Q(x) # âś…
updVec .= false
for j in 1:proper_Q[1]
Y = proper_Q[3][:, :, j]
JuMP.fix.(tr2_Y, Y) # 🖌️
opt_ass_opt(tr2, "tr2")
b2 = jvr6.(tr2_b2); oΛ2_check = JuMP.value(oΛ2) # 🥑✂️
set_argZ_objective(x, Y, b2); opt_ass_optortime(argZ, "argZ")
Z_new = jvr6.(argZ_Z)
cn, px, pY, pb = JuMP.value(argZ_cn), jvr6.(argZ_c_x), jvr6.(argZ_c_Y), jvr6.(argZ_c_b2)
if oΛ2_check < cn + ip(px, x) + ip(pY, Y) + ip(pb, b2) - UPDTOL
add_cut_for_oΛ2(cn, px, pY, pb)
updVec[2] = true
end
recruit(Z_list, Z_new)
proper_P = get_proper_P(x, Y, Z_list) # âś…
add_cut_for_oLh(x, proper_Q, proper_P, oLh_check) && (updVec[1] = true)
end
@info "$updVec, lb = $lb"
all(updVec .== false) && break
end
The result should be like
[ Info: Bool[1, 1], lb = -208.90723332481141
┌ Warning: The addition operator has been used on JuMP expressions a large number of times. This warning is safe to ignore but may indicate that model generation is slower than necessary. For performance reasons, you should not add expressions in a loop. Instead of x += y, use add_to_expression!(x,y) to modify x in place. If y is a single variable, you may also use add_to_expression!(x, coef, y) for x += coef*y.
â”” @ JuMP K:\judepot1112\packages\JuMP\i68GU\src\operators.jl:282
[ Info: Bool[1, 1], lb = 0.6
...
