julia> JuMP.unset_silent(outer)
Set parameter OutputFlag to value 1
julia> JuMP.optimize!(outer)
Gurobi Optimizer version 12.0.2 build v12.0.2rc0 (win64 - Windows 11.0 (26100.2))
CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 39 rows, 41 columns and 127 nonzeros
Model fingerprint: 0x6abf7598
Model has 38 quadratic objective terms
Coefficient statistics:
Matrix range [1e+00, 3e+02]
Objective range [9e-18, 5e+00]
QObjective range [2e+00, 2e+00]
Bounds range [1e+03, 1e+03]
RHS range [4e+02, 2e+03]
Presolve removed 11 rows and 2 columns
Presolve time: 0.00s
Presolved: 28 rows, 39 columns, 80 nonzeros
Presolved model has 38 quadratic objective terms
Ordering time: 0.00s
Barrier statistics:
Free vars : 39
AA' NZ : 6.900e+01
Factor NZ : 1.060e+02
Factor Ops : 8.340e+02 (less than 1 second per iteration)
Threads : 1
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 1.29618707e+01 1.29610499e+01 4.29e+02 1.02e+01 1.08e+06 0s
1 7.79452936e+03 -9.29391938e+03 4.59e+01 4.59e-01 8.38e+04 0s
2 2.34735483e+03 -3.67336088e+03 3.99e+00 4.00e-02 1.54e+04 0s
3 9.28074635e+01 -8.89553051e+02 5.28e-01 5.42e-03 3.32e+03 0s
4 2.67970847e+00 -2.69217523e+02 6.93e-02 1.81e-04 4.56e+02 0s
5 4.93346415e+00 -2.33632458e+02 6.85e-02 1.78e-04 4.42e+02 0s
6 2.61337406e+00 -2.70695461e+02 3.15e-01 1.50e-04 4.63e+02 0s
7 1.09725171e+01 -1.34247786e+02 5.05e-01 1.37e-04 3.67e+02 0s
8 6.15979734e+00 -7.90459574e+01 2.02e-01 6.79e-05 1.87e+02 0s
9 2.03422156e-02 -8.54163247e+01 1.06e-01 4.22e-05 1.57e+02 0s
10 6.28528487e-01 -3.31827999e+01 2.99e-02 4.33e-11 3.35e+01 0s
11 1.33921978e-03 -2.70800877e+01 4.09e-13 2.12e-12 1.23e+00 0s
12 2.55952604e-04 -3.06308511e+01 1.99e-13 2.79e-12 1.10e-02 0s
13 3.33498799e-04 -3.42217263e+01 3.09e-13 2.70e-12 1.19e-04 0s
14 5.57483325e-04 -3.51843270e+01 9.52e-14 1.82e-12 4.34e-07 0s
15 2.39292335e-04 -3.74834526e+01 2.37e-13 1.82e-12 2.28e-09 0s
16 2.43578332e-04 -3.89336401e+01 1.87e-13 2.41e-12 1.43e-10 0s
17 2.11498031e-03 -5.43399078e+01 4.19e-02 3.42e-12 7.46e-11 0s
18 2.96317192e-04 -5.53941734e+01 1.40e-13 5.61e-12 2.39e-12 0s
19 4.12094691e-04 -5.67513066e+01 5.10e-14 4.84e-12 3.94e-13 0s
20 6.75682559e-04 -5.83702637e+01 1.42e-13 1.61e-12 4.12e-16 0s
21 8.90546712e-04 -5.95038795e+01 2.26e-13 3.58e-12 1.14e-16 0s
22 9.87057970e-04 -6.26382723e+01 1.59e-13 3.55e-12 2.54e-19 0s
Barrier solved model in 22 iterations and 0.01 seconds (0.00 work units)
Optimal objective 9.87057970e-04
User-callback calls 113, time in user-callback 0.00 sec
julia> JuMP.solution_summary(outer)
solution_summary(; result = 1, verbose = false)
├ solver_name : Gurobi
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ ├ raw_status : Model was solved to optimality (subject to tolerances), and an optimal solution is available.
│ └ objective_bound : -2.70801e+01
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : 9.87058e-04
│ └ dual_objective_value : -4.97810e+01
└ Work counters
├ solve_time (sec) : 7.99990e-03
├ simplex_iterations : 0
├ barrier_iterations : 22
└ node_count : 0
julia> vio = JuMP.MOI.get(JuMP.backend(outer), Gurobi.ModelAttribute("MaxVio"))
43.97206084166428
julia> Gurobi.GRBwrite(JuMP.unsafe_backend(outer), "outer.mps")
0
This “OPTIMAL” solution is very weird—noticing the MaxVio, and the termination gap. Is this normal? (@odow I remember that there was also a gap displayed here)
The code to reproduce this model outer is here. (A lightweight PGlib.jl data package is required.)
code
import SparseArrays, PGLib; RawDict = PGLib.pglib("pglib_opf_case118_ieee.m");
import Random
function rectify_dir!(bft)
for r in 1:size(bft, 1)
if bft[r, 1] > bft[r, 2]
bft[r, 1], bft[r, 2] = bft[r, 2], bft[r, 1]
elseif bft[r, 1] == bft[r, 2]
error("self-branch")
end
end
end; function assert_is_connected(bft)
B = size(bft, 1)
pb = Vector(1:B) # primal bs
sn = [1] # subnet
for ite in 1:B # only du a full iteration as a conservative choice
progress = false
for (i, b) in enumerate(pb) # take out branch b
if bft[b, 1] in sn
on = bft[b, 2] # the other node
on ∉ sn && push!(sn, on)
elseif bft[b, 2] in sn
on = bft[b, 1] # the other node
on ∉ sn && push!(sn, on)
else
continue
end
popat!(pb, i); progress = true; break # fathom branch b
end
if progress == false
error("ite = $ite. All the rest branches are not connected to the subnet being investigated. Check the current subnet")
end
1:maximum(bft) ⊆ sn && return # The graph is proved to be connected
end
error("here shouldn't be reached")
end; function bft_2_A(bft)
B, N = size(bft, 1), maximum(bft)
return SparseArrays.sparse([Vector(1:B); Vector(1:B)], vec(bft), [-ones(Int, B); ones(Int, B)], B, N)
end; function adjnode(n)
nv = Int[]
for ci in findall(x -> x == n, bft)
b, c = ci.I # row(= branch), col
on = bft[b, 3 - c]
on ∉ nv && push!(nv, on)
end
sort(nv)
end; function get_b(n, m) # use this after `rectify_dir!(bft)` had been done
n > m && ((n, m) = (m, n))
for r in 1:size(bft, 1)
(bft[r, 1] == n && bft[r, 2] == m) && return r
end
error("no such branch exist")
end;
G, G2N, Gcapinit, Gptcost = let
G = length(RawDict["gen"])
G2N = [RawDict["gen"]["$g"]["gen_bus"] for g in 1:G]
cost = [RawDict["gen"]["$g"]["cost"] for g in 1:G]
pmax = [RawDict["gen"]["$g"]["pmax"] for g in 1:G]
gi = findall(x -> x > 1e-5, pmax)
G2N, cost, pmax, G = G2N[gi], cost[gi], pmax[gi], length(gi)
cost = [.001i[1] for i in cost]
f = x -> round(x; digits = 6)
G, G2N, pmax, f.(cost)
end;
D, D2N, Dref = let
D = length(RawDict["load"])
D2N = [RawDict["load"]["$d"]["load_bus"] for d in 1:D]
Dref = [RawDict["load"]["$d"]["pd"] for d in 1:D]
D, D2N, Dref
end;
T, Dura, Dmat = let
Random.seed!(1); Dlarge = rand(1.1:7e-13:1.7, D) .* Dref
Random.seed!(6); Dmiddle = rand(.7:7e-13:1.7, D) .* Dref
Dmat = [Dref Dmiddle Dlarge]
Dura = [ 6, 3, 1]
T = size(Dmat, 2)
T, Dura, Dmat
end;
S, prob_vec, avai_vec, g_hyd, Ginvcost = let
g_hyd = 17; Gptcost[g_hyd] = 0 # 💧 pg[g = 17] is a hydro generator
Ginvcost = 2.0Gptcost;
Ginvcost[g_hyd] = maximum(Gptcost);
avai_vec = [.3, .79] # availability factor for 💧 generator
prob_vec = [.5, .5] # The associated probability vector
S = length(prob_vec)
S, prob_vec, avai_vec, g_hyd, Ginvcost
end;
# feasible cost = 959.7425433749277
# SP solution cost = 904.3816383661987
# NA-relaxed cost = 884.6469843243801
B, N = length.([RawDict["branch"], RawDict["bus"]]) # (186, 118), fixed for this network
bft = let # load a rectified bft
bft = [RawDict["branch"]["$r"][c == 1 ? "f_bus" : "t_bus"] for r in 1:B, c in 1:2]; # raw bft
assert_is_connected(bft);
rectify_dir!(bft);
bft
end; A = bft_2_A(bft);
import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
function bpf(b, Dvec, pgvec, pfvec) # dispatch function of power flow on branch b
b == 9 && return -pgvec[findfirst(x -> x == 10, G2N)::Int] # b8[9] = -3
b == 7 && return bpf(9, Dvec, pgvec, pfvec) # b8[7] = -1
b == 94 && return bpf(93, Dvec, pgvec, pfvec) # b8[94] = -1
b == 127 && return -bpf(126, Dvec, pgvec, pfvec) # b8[127] = -1
b == 134 && return -pgvec[findfirst(x -> x == 87, G2N)::Int] # b8[134] = -3
b == 176 && return -pgvec[findfirst(x -> x == 111, G2N)::Int] # b8[176] = -3
b == 184 && return Dvec[findfirst(x -> x == 117, D2N)::Int] # b8[184] = -2
dup_bvec = [67, 76, 99, 86, 124, 139, 142]
b ∈ dup_bvec && return false # b8[dup_bvec] .= 0; b8[dup_bvec .- 1] .= 2
return pfvec[(f = findfirst(x -> x == b, F2B);)] # [Fallback] the corresponding JuMP.@variable
end
n8 = let # elements being
# [1]: a bus with load only, and whose degree is larger than 2
# [2]: a bus with gen only, and whose degree is larger than 2
# [3]: a bus with both power injection and load
# [-1]: a bus with load only, and whose degree is 1, we won't add KCL constr here
# [-3]: a bus with gen only, and whose degree is 1, we won't add KCL constr here
# [-2]: a bare bus whose degree is 2, we won't add KCL constr here
# [0]: a bare joint bus
n8 = zeros(Int8, N); n8[D2N] .+= 1; n8[G2N] .+= 2
n8[117] = -1
n8[[10, 87, 111]] .= -3
n8[[9, 63, 81]] .= -2
n8
end; b8 = let # Note this procedure is specific to case118, therefore executing once is sufficient
# elements being
# [2]: this line represents a 2-parallel branch, works as a normal line
# [1]: a normal line
# [0]: this line was banned (pf ≡ 0) due to parallel redundancy, we won't allocate pf variable here
# [-1]: power flow of this line hinges on one of its adjacent line, we won't allocate pf variable here
# [-2]: power flow of this line hinges on the corresponding pure load, we won't allocate pf variable here
# [-3]: power flow of this line hinges on the corresponding pure gen, we won't allocate pf variable here
b8 = ones(Int8, size(bft, 1))
b8[9] = -3
b8[7] = -1
b8[94] = -1
b8[127] = -1
b8[134] = -3
b8[176] = -3
b8[184] = -2
dup_bvec = [67, 76, 99, 86, 124, 139, 142]
b8[dup_bvec] .= 0; b8[dup_bvec .- 1] .= 2
b8
end;
F, F2B, fra = let
F2B = findall(s -> s ∈ 1:2, b8); F = length(F2B);
fra = [RawDict["branch"]["$b"]["rate_a"] for b in F2B]; # f_rate_a is a used part of b_rate_a, under Dref, if fra *= .59, then INFEASIBLE
F, F2B, fra
end;
min_s = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(min_s);
JuMP.@variable(min_s, 0 ≤ invg_s[1:G] ≤ 300); # ⚠️ set an artificial ub here
JuMP.@variable(min_s, 0 ≤ pg_s[1:G, 1:T]);
JuMP.@expression(min_s, power_max, invg_s + Gcapinit)
JuMP.@constraint(min_s, [g = 1:G, t = 1:T], pg_s[g, t] ≤ power_max[g]);
JuMP.@variable(min_s, -fra[f] ≤ pf_s[f = 1:F, t = 1:T] ≤ fra[f]);
for t in 1:T, nState in 0:3, n in findall(ø -> ø == nState, n8) # real-time everynode KCL
p_gen = (nState < 2 ? false : pg_s[findfirst(x -> x == n, G2N)::Int, t])
p_demand = (nState == 1 || nState == 3 ? Dmat[findfirst(x -> x == n, D2N)::Int, t] : false)
bv, sgv = SparseArrays.findnz(view(A, :, n))
f = b -> bpf(b, view(Dmat, :, t), view(pg_s, :, t), view(pf_s, :, t))
JuMP.@constraint(min_s, sgv'f.(bv) + p_gen == p_demand)
end;
JuMP.@expression(min_s, hyd_energy, Dura'view(pg_s, g_hyd, :))
JuMP.@expression(min_s, hyd_energy_max, power_max[g_hyd]sum(Dura))
JuMP.@expression(min_s, cost_vec_s, [Ginvcost'invg_s, (Gptcost'pg_s)Dura])
constr_s = JuMP.@constraint(min_s, 0 == 1); # initial dumb
function solve_min_s(min_s, constr_s, β_s, avai_vec_s, prob_vec_s)
JuMP.delete(min_s, constr_s)
constr_s = JuMP.@constraint(min_s, min_s[:hyd_energy] ≤ min_s[:hyd_energy_max]avai_vec_s)
invg_s, cost_vec_s = min_s[:invg_s], min_s[:cost_vec_s]
bias_s = JuMP.@expression(min_s, sum(cost_vec_s)prob_vec_s)
JuMP.@objective(min_s, Min, β_s'invg_s + bias_s) # ✅ θs-β cut
JuMP.optimize!(min_s); JuMP.assert_is_solved_and_feasible(min_s; allow_local = false)
slope_s, bias_s = JuMP.value.(invg_s), JuMP.value(bias_s)
return constr_s, JuMP.objective_value(min_s), slope_s, bias_s
end
outer = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(outer); # Linear Programming
JuMP.@variable(outer, outer_β[1:G, 1:S]); JuMP.@constraint(outer, sum(outer_β; dims = 2) .== 0);
JuMP.@variable(outer, outer_O ≤ 1e3); # ⚠️ this ub needs to be valid
JuMP.@variable(outer, outer_θ[1:S]);
JuMP.@constraint(outer, outer_O == sum(outer_θ))
good_objc = JuMP.@constraint(outer, 0 ≤ 1) # initial dumb
function solve_outer_1(outer, good_objc)
JuMP.delete(outer, good_objc)
JuMP.@objective(outer, Max, outer[:outer_O])
JuMP.optimize!(outer); JuMP.assert_is_solved_and_feasible(outer; allow_local = false)
return dual_upper_bound = JuMP.objective_bound(outer)
end;
function solve_outer_2(outer, β, ub, lb)
good_objc = JuMP.@constraint(outer, .7ub + .3lb ≤ outer[:outer_O])
outer_β = outer[:outer_β]
Δ = JuMP.@expression(outer, outer_β - β)
JuMP.@objective(outer, Min, sum(x -> (x)x, Δ))
JuMP.optimize!(outer); JuMP.assert_is_solved_and_feasible(outer; allow_local = false)
vio = JuMP.MOI.get(JuMP.backend(outer), Gurobi.ModelAttribute("MaxVio"))
vio > 9e-7 && @error "outer_QP: vio = $vio"
return good_objc, JuMP.value.(outer_β)
end;
β = zeros(G, S); # start from zero
lbv = similar(prob_vec);
for ite in 1:8
global good_objc
for s in 1:S
global constr_s
constr_s, lbv[s], slope_s, bias_s = solve_min_s(min_s, constr_s, view(β, :, s), avai_vec[s], prob_vec[s])
JuMP.@constraint(outer, outer_θ[s] ≤ bias_s + slope_s'view(outer_β, :, s))
end; lb = sum(lbv)
ub = solve_outer_1(outer, good_objc)
@info "ite = $ite, (P) lb = $lb | $ub = ub (D)"
good_objc, β = solve_outer_2(outer, β, ub, lb)
end
global good_objc # ite == 9
for s in 1:S
global constr_s
constr_s, lbv[s], slope_s, bias_s = solve_min_s(min_s, constr_s, view(β, :, s), avai_vec[s], prob_vec[s])
JuMP.@constraint(outer, outer_θ[s] ≤ bias_s + slope_s'view(outer_β, :, s))
end; lb = sum(lbv)
ub = solve_outer_1(outer, good_objc)
@info "(P) lb = $lb | $ub = ub (D)"
good_objc = JuMP.@constraint(outer, .7ub + .3lb ≤ outer[:outer_O])
outer_β = outer[:outer_β]
Δ = JuMP.@expression(outer, outer_β - β)
JuMP.@objective(outer, Min, sum(x -> (x)x, Δ))
JuMP.unset_silent(outer)
JuMP.optimize!(outer)
JuMP.solution_summary(outer)
vio = JuMP.MOI.get(JuMP.backend(outer), Gurobi.ModelAttribute("MaxVio"))
Gurobi.GRBwrite(JuMP.unsafe_backend(outer), "outer.mps")