I cannot figure out why the JuMP.primal_status(model) is UNKNOWN_RESULT_STATUS here
julia> model # a model that is solved before and modified
A JuMP Model
├ solver: Gurobi
├ objective_sense: MIN_SENSE
│ └ objective_function_type: JuMP.AffExpr
├ num_variables: 284616
├ num_constraints: 431774
│ ├ ...
...
julia> JuMP.optimize!(model)
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
Non-default parameters:
TimeLimit 120
Optimize a model with 154408 rows, 284616 columns and 882830 nonzeros
Model fingerprint: 0x154f3acc
Variable types: 174000 continuous, 110616 integer (110616 binary)
Coefficient statistics:
Matrix range [9e-01, 2e+01]
Objective range [1e+01, 3e+01]
Bounds range [1e+00, 6e+01]
RHS range [1e+00, 1e+05]
Processing MIP start from previous solve: 0 nodes explored in subMIP, total elapsed time 5s
MIP start from previous solve produced solution with objective 3.60589e+06 (5.04s)
MIP start from previous solve produced solution with objective 3.60589e+06 (5.09s)
Loaded MIP start from previous solve with objective 3.60589e+06
Processed MIP start in 5.10 seconds (4.05 work units)
Presolve removed 1905 rows and 5225 columns
Presolve time: 1.55s
Presolved: 152503 rows, 279391 columns, 863091 nonzeros
Variable types: 172314 continuous, 107077 integer (107077 binary)
Deterministic concurrent LP optimizer: primal simplex, dual simplex, and barrier
Showing barrier log only...
Root barrier log...
Ordering time: 0.08s
Barrier statistics:
AA' NZ : 4.064e+05
Factor NZ : 2.694e+06 (roughly 150 MB of memory)
Factor Ops : 1.249e+08 (less than 1 second per iteration)
Threads : 2
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 1.86652623e+10 -1.75398155e+08 2.11e+08 0.00e+00 2.07e+06 8s
......
30 3.53257740e+06 3.53257740e+06 9.80e-12 1.42e-13 1.88e-11 11s
Barrier solved model in 30 iterations and 11.04 seconds (9.03 work units)
Optimal objective 3.53257740e+06
Root crossover log...
11696 DPushes remaining with DInf 0.0000000e+00 11s
0 DPushes remaining with DInf 0.0000000e+00 11s
18701 PPushes remaining with PInf 0.0000000e+00 11s
0 PPushes remaining with PInf 0.0000000e+00 12s
Push phase complete: Pinf 0.0000000e+00, Dinf 9.6187894e-11 12s
Root simplex log...
Iteration Objective Primal Inf. Dual Inf. Time
29843 3.5325774e+06 0.000000e+00 0.000000e+00 12s
Concurrent spin time: 0.61s (can be avoided by choosing Method=3)
Solved with barrier
29843 3.5325774e+06 0.000000e+00 0.000000e+00 13s
Root relaxation: objective 3.532577e+06, 29843 iterations, 5.63 seconds (3.26 work units)
Total elapsed time = 15.49s (DegenMoves)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 3532577.40 0 5404 3605888.90 3532577.40 2.03% - 16s
H 0 0 3548146.6821 3532577.40 0.44% - 17s
0 0 3533464.00 0 5928 3548146.68 3533464.00 0.41% - 22s
H 0 0 3547803.3639 3533464.00 0.40% - 24s
H 0 0 3547721.3895 3533464.00 0.40% - 25s
0 0 3533938.63 0 6123 3547721.39 3533938.63 0.39% - 27s
0 0 3534118.30 0 6225 3547721.39 3534118.30 0.38% - 29s
0 0 3534348.21 0 6312 3547721.39 3534348.21 0.38% - 31s
0 0 3534363.63 0 6393 3547721.39 3534363.63 0.38% - 32s
0 0 3534364.65 0 6419 3547721.39 3534364.65 0.38% - 34s
0 0 3534364.88 0 6433 3547721.39 3534364.88 0.38% - 35s
0 0 3535004.49 0 6191 3547721.39 3535004.49 0.36% - 41s
H 0 0 3544867.0580 3535004.55 0.28% - 42s
0 0 3535167.41 0 6125 3544867.06 3535167.41 0.27% - 45s
0 0 3535261.79 0 6180 3544867.06 3535261.79 0.27% - 47s
0 0 3535316.12 0 6357 3544867.06 3535316.12 0.27% - 48s
0 0 3535322.42 0 6424 3544867.06 3535322.42 0.27% - 49s
0 0 3535323.29 0 6461 3544867.06 3535323.29 0.27% - 50s
0 0 3535490.20 0 6092 3544867.06 3535490.20 0.26% - 60s
H 0 0 3543860.8838 3535490.22 0.24% - 62s
H 0 0 3543833.9358 3535490.22 0.24% - 64s
0 0 3535532.34 0 6307 3543833.94 3535532.34 0.23% - 65s
0 0 3535559.39 0 6304 3543833.94 3535559.39 0.23% - 66s
0 0 3535565.03 0 6378 3543833.94 3535565.03 0.23% - 67s
0 0 3535566.11 0 6439 3543833.94 3535566.11 0.23% - 67s
0 0 3535670.38 0 6126 3543833.94 3535670.38 0.23% - 79s
H 0 0 3542615.0971 3535670.38 0.20% - 81s
H 0 0 3542344.1418 3535670.38 0.19% - 84s
0 0 3535702.72 0 6146 3542344.14 3535702.72 0.19% - 85s
0 0 3535718.84 0 6228 3542344.14 3535718.84 0.19% - 86s
0 0 3535722.77 0 6310 3542344.14 3535722.77 0.19% - 86s
0 0 3535723.49 0 6345 3542344.14 3535723.49 0.19% - 87s
0 0 3535749.91 0 6304 3542344.14 3535749.91 0.19% - 92s
0 0 3535761.89 0 6384 3542344.14 3535761.89 0.19% - 95s
0 0 3535763.94 0 6468 3542344.14 3535763.94 0.19% - 96s
0 0 3535764.43 0 6505 3542344.14 3535764.43 0.19% - 96s
0 0 3535785.44 0 6344 3542344.14 3535785.44 0.19% - 101s
H 0 0 3542148.8408 3535785.45 0.18% - 109s
0 0 3535790.28 0 6352 3542148.84 3535790.28 0.18% - 110s
0 0 3535791.24 0 6435 3542148.84 3535791.24 0.18% - 111s
0 0 3535804.07 0 6325 3542148.84 3535804.07 0.18% - 116s
Cutting planes:
Learned: 1
Gomory: 108
Cover: 247
Implied bound: 129
Clique: 25
MIR: 9707
StrongCG: 9
Flow cover: 6221
GUB cover: 8
Zero half: 7
Mod-K: 1
RLT: 130
Relax-and-lift: 588
BQP: 1
Explored 1 nodes (135082 simplex iterations) in 120.13 seconds (115.52 work units)
Thread count was 8 (of 8 available processors)
Solution count 10: 3.54215e+06 3.54234e+06 3.54262e+06 ... 3.60589e+06
Time limit reached
Warning: max constraint violation (1.9111e-06) exceeds tolerance
Warning: max bound violation (1.2831e-06) exceeds tolerance
Best objective 3.542148840800e+06, best bound 3.535804097091e+06, gap 0.1791%
User-callback calls 8145, time in user-callback 0.03 sec
julia> JuMP.termination_status(model), JuMP.primal_status(model)
(MathOptInterface.TIME_LIMIT, MathOptInterface.UNKNOWN_RESULT_STATUS)
As you can see, although TimeLimit is reached, 10 solutions are already available.
To introduce the model above before optimize!, run this source
Summary
import JuMP, Gurobi
import JuMP.value as ı
import LinearAlgebra.⋅ as ⋅
import LinearAlgebra.norm as norm
import Random
function get_C_and_O()
C = [
# Case 1: Flat midday, strong evening peak
[10, 9, 9, 9, 10, 12, 15, 18, 20, 18, 16, 15, 14, 15, 16, 18, 22, 28, 32, 30, 26, 20, 15, 12],
# Case 2: Two peaks (morning + evening), midday dip
[12, 11, 11, 12, 14, 18, 24, 26, 22, 18, 15, 14, 13, 14, 18, 24, 30, 34, 32, 28, 22, 18, 15, 13],
# Case 3: Midday solar effect (cheapest at noon, peaks morning & evening)
[16, 15, 14, 14, 15, 18, 24, 30, 28, 22, 18, 12, 10, 12, 16, 22, 28, 34, 36, 32, 28, 24, 20, 18],
# Case 4: Steady climb during day, single high plateau evening
[8, 8, 8, 9, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 36, 34, 30, 26, 20, 14],
# Case 5: Inverted (very low off-peak overnight, high midday, gentle evening)
[5, 5, 5, 6, 8, 12, 18, 24, 28, 32, 36, 38, 36, 34, 30, 28, 26, 24, 22, 20, 18, 14, 10, 8]
]
O = [
# Case 1: Typical hot day (peak ~38°C around 15:00)
[28,28,27,27,28,29,31,33,35,36,37,38,38,38,37,36,35,34,32,31,30,29,29,28],
# Case 2: Extremely hot day (peak ~42°C, late afternoon peak)
[29,28,28,28,29,31,33,35,37,39,40,41,42,42,41,40,38,36,34,33,32,31,30,29],
# Case 3: Milder hot day (peak ~35°C, smooth curve)
[27,27,27,27,28,29,30,32,33,34,35,35,35,35,34,33,32,31,30,29,28,28,28,27],
# Case 4: Heatwave night (doesn’t cool much at night, peak 44°C)
[32,32,31,31,32,34,36,38,40,42,43,44,44,44,43,42,41,40,38,37,36,35,34,33],
# Case 5: Cool morning, sharp rise, peak ~39°C
[27,27,27,27,28,29,30,33,35,37,38,39,39,39,38,37,36,34,32,31,30,29,28,27]
]
return C[rand(1:5)], O[rand(1:5)]
end;
function assert_is_approxly_solved(model)
local (t, p) = (JuMP.termination_status(model), JuMP.primal_status(model))
t ∈ [JuMP.TIME_LIMIT, JuMP.OPTIMAL, JuMP.SOLUTION_LIMIT] || error("terminate with $t")
p == JuMP.FEASIBLE_POINT || error("terminate with primal_status = $p")
return nothing
end;
function get_pair_and_self_rng(J)
d = J÷4
# d = J÷2
local (Rng1, Rng2) = (1:d, d+1:J)
end;
function prev(t, d, T) return (n = t-d; n<1 ? T+n : n) end;
function get_P_AC(O, OH, CND, Q_I, COP) return ceil(Int, ((maximum(O) - OH)CND + maximum(Q_I)) / COP) end;
function gen_ac_data()
CND = .5rand(1:7) # thermal conductance of the building wall
INR = rand(6:20) # thermal inertia of the building
COP = rand(2:.5:4) # how many Watts cooling power for 1-Watt e-power
Q_I = rand(3:9, T) # heat generated by indoor human activity, server etc.
Q_BUS = rand(25:35) # RateA of the cooling water/air pipeline
OH = rand(24:29) # the hottest indoor degree Celsius
OΔ = rand(4:9) # the coldest indoor degree Celsius = OH - OΔ
P_AC = get_P_AC(O, OH, CND, Q_I, COP)
return CND, INR, COP, Q_I, Q_BUS, OH, OΔ, P_AC
end;
function get_E_ES(rng)
M = rand(rng) # Max E
i = rand(0:M) # initial SOC _is_ this value
e = rand(0:min(M, 21)) # ending SOC should ≥ this value
local E_ES = (i = i, M = M, e = e) # assume lb = 0 w.l.o.g.
end;
function get_self_P_BUS(D, U, P_EV, E_EV, O, CND, INR, COP, Q_I, OH, OΔ, P_AC)
local model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)) # a temp household
JuMP.set_silent(model)
bU, pU = add_U_module!(model, U)
bEV, pEV = add_self_EV_module!(model, P_EV, E_EV)
o, q, pAC = add_AC_module!(model, O, CND, INR, COP, Q_I, 0, OH, OΔ, P_AC) # Q_BUS = 0
pBus = JuMP.@variable(model)
JuMP.@constraint(model, pBus .≥ D + pU + pEV + pAC) # No G | ES
JuMP.@objective(model, Min, pBus)
JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false)
val = ı(pBus)
val > 0 || error("The self household has P_BUS = $val")
local P_BUS = ceil(Int, val)
end;
function add_ES_module!(model, P_ES, E_ES)
p_ES = (
c = JuMP.@variable(model, [1:T], lower_bound = 0, upper_bound = P_ES.c),
d = JuMP.@variable(model, [1:T], lower_bound = 0, upper_bound = P_ES.d)
); e_ES = JuMP.@variable(model, [t=1:T], lower_bound = t<T ? 0 : E_ES.e, upper_bound = E_ES.M)
JuMP.@constraint(model, [t=1:T], .95p_ES.c[t] - p_ES.d[t]/.95 == e_ES[t]-(t>1 ? e_ES[t-1] : E_ES.i))
return p_ES, e_ES
end;
function add_AC_module!(model, O, CND, INR, COP, Q_I, Q_BUS, OH, OΔ, P_AC)
pAC = JuMP.@variable(model, [1:T], lower_bound = 0, upper_bound = P_AC)
o = JuMP.@variable(model, [1:T], lower_bound = OH-OΔ, upper_bound = OH)
q = JuMP.@variable(model, [1:T], lower_bound = 0, upper_bound = Q_BUS)
JuMP.@constraint(model, [t=1:T], (O[t]-o[t])CND + Q_I[t] -q[t] -pAC[t]COP == (o[t<T ? t+1 : 1]-o[t])INR)
return o, q, pAC
end;
function add_U_module!(model, U)
bU = JuMP.@variable(model, [t = 1:T, i = eachindex(U)], Bin) # Matrix{JuMP.VariableRef}
JuMP.@constraint(model, sum(bU; dims = 1) .≥ true)
# ✅ dependent variable is modeled as @expression
pU = JuMP.@expression(model, [t=1:T], sum(sum(bU[prev(t,φ-1,T), i]P for (φ,P) = enumerate(v)) for (i,v) = enumerate(U))) # Vector{JuMP.AffExpr}
return bU, pU
end;
function add_self_EV_module!(model, P_EV, E_EV) # self means a household in a block with cardinality 1
bEV, pEV = JuMP.@variable(model, [1:T], Bin), JuMP.@variable(model, [1:T])
JuMP.@constraint(model, (P_EV.m)bEV ≤ pEV)
JuMP.@constraint(model, pEV ≤ (P_EV.M)bEV)
JuMP.@constraint(model, sum(pEV) ≥ E_EV)
return bEV, pEV
end;
function add_EV_1_module!(model, P_EV_1, E_EV_1)
bLent, pLent = JuMP.@variable(model, [1:T], Bin), JuMP.@variable(model, [1:T])
bEV_1, pEV_1 = JuMP.@variable(model, [1:T], Bin), JuMP.@variable(model, [1:T])
JuMP.@constraint(model, bEV_1 ≤ bLent)
JuMP.@constraint(model, (1 .- bLent)P_EV_1.m ≤ pLent)
JuMP.@constraint(model, pLent ≤ (1 .- bLent)P_EV_1.M)
JuMP.@constraint(model, (P_EV_1.m)bEV_1 ≤ pEV_1)
JuMP.@constraint(model, pEV_1 ≤ (P_EV_1.M)bEV_1)
JuMP.@constraint(model, sum(pEV_1) ≥ E_EV_1)
return bEV_1, pEV_1, bLent, pLent
end;
function add_EV_2_module!(model, P_EV_2, E_EV_2, bLent, pLent)
bEV_2, pEV_2 = JuMP.@variable(model, [1:T], Bin), JuMP.@variable(model, [1:T])
JuMP.@constraint(model, bEV_2 ≤ bLent)
JuMP.@constraint(model, (P_EV_2.m)bEV_2 ≤ pEV_2)
JuMP.@constraint(model, pEV_2 ≤ (P_EV_2.M)bEV_2)
JuMP.@constraint(model, sum(pEV_2 + pLent) ≥ E_EV_2)
return bEV_2, pEV_2
end;
function add_self_circuit_breaker_module!(model, P_BUS, D, pU, pEV, pAC)
pBus = JuMP.@variable(model, [1:T], lower_bound = 0, upper_bound = P_BUS)
JuMP.@constraint(model, pBus ≥ D + pU + pEV + pAC) # No G | ES
return pBus
end;
function add_circuit_breaker_pair_module!(model, P_BUS_1, P_BUS_2,
p_ES, G, pLent, pEV_1, pU_1, pAC_1, D_1,
pEV_2, pU_2, pAC_2, D_2)
pBus_1 = JuMP.@variable(model, [1:T], lower_bound = -P_BUS_1, upper_bound = P_BUS_1)
pBus_2 = JuMP.@variable(model, [1:T], lower_bound = 0, upper_bound = P_BUS_2)
JuMP.@constraint(model, pBus_1 == p_ES.c -p_ES.d -G + pLent + pEV_1 + pU_1 + pAC_1 + D_1)
JuMP.@constraint(model, pBus_2 ≥ pEV_2 + pU_2 + pAC_2 + D_2)
return pBus_1, pBus_2
end;
function get_a_paired_block(O)::NamedTuple
model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)) # for a block who has a lender and a borrower house
JuMP.set_silent(model)
# 6 lines
G = rand(0:17, T)
D_1 = rand(0:5, T)
P_ES, E_ES = (c = rand(1:6), d = rand(1:6)), get_E_ES(19:55)
U_1 = [rand(1:4, rand(2:5)) for _ = 1:rand(1:4)] # each entry is a cycle vector of an uninterruptible load
P_EV_1, E_EV_1 = (m = rand((1., 1.5)), M = rand(3:7)), rand(10:39)
CND_1, INR_1, COP_1, Q_I_1, Q_BUS_1, OH_1, OΔ_1, P_AC_1 = gen_ac_data()
# lender house
p_ES, e_ES = add_ES_module!(model, P_ES, E_ES) # only for the lender house
bU_1, pU_1 = add_U_module!(model, U_1)
bEV_1, pEV_1, bLent, pLent = add_EV_1_module!(model, P_EV_1, E_EV_1)
o_1, q_1, pAC_1 = add_AC_module!(model, O, CND_1, INR_1, COP_1, Q_I_1, 0, OH_1, OΔ_1, P_AC_1) # Q_BUS = 0
# 4 lines
D_2 = rand(0:5, T) # borrower house
U_2 = [rand(1:4, rand(2:5)) for _ = 1:rand(1:4)] # each entry is a cycle vector of an uninterruptible load
P_EV_2, E_EV_2 = (m = rand((1., 1.5)), M = rand(3:7)), rand(10:39)
CND_2, INR_2, COP_2, Q_I_2, Q_BUS_2, OH_2, OΔ_2, P_AC_2 = gen_ac_data()
# borrower house
bU_2, pU_2 = add_U_module!(model, U_2)
bEV_2, pEV_2 = add_EV_2_module!(model, P_EV_2, E_EV_2, bLent, pLent)
o_2, q_2, pAC_2 = add_AC_module!(model, O, CND_2, INR_2, COP_2, Q_I_2, 0, OH_2, OΔ_2, P_AC_2) # Q_BUS = 0
# determine the circuit breaker limit
P_BUS_2 = let pBus_2 = JuMP.@variable(model)
local s = JuMP.@constraint(model, pBus_2 .≥ pEV_2 + pU_2 + pAC_2 + D_2)
JuMP.@objective(model, Min, pBus_2)
JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false)
local val = ı(pBus_2)
val > 0 || error("pBus_2 has value $val")
val = ceil(Int, val)
JuMP.delete(model, s)
JuMP.delete(model, pBus_2)
val
end
pBus_2 = JuMP.@variable(model, [1:T], lower_bound = 0, upper_bound = P_BUS_2)
JuMP.@constraint(model, pBus_2 ≥ pEV_2 + pU_2 + pAC_2 + D_2)
P_BUS_1 = let pBus_1 = JuMP.@variable(model)
local m = JuMP.@constraint(model, -pBus_1 .≤ p_ES.c -p_ES.d -G + pLent + pEV_1 + pU_1 + pAC_1 + D_1)
local p = JuMP.@constraint(model, pBus_1 .≥ p_ES.c -p_ES.d -G + pLent + pEV_1 + pU_1 + pAC_1 + D_1)
JuMP.@objective(model, Min, pBus_1)
JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false)
local val = ı(pBus_1)
val ≥ -1e-5 || error("pBus_1 has value $val")
val = max(1, ceil(Int, val))
JuMP.delete.(model, (m, p))
JuMP.delete(model, pBus_1)
val
end
local d = (
P_BUS_1 = P_BUS_1,
P_BUS_2 = P_BUS_2,
G = G,
P_ES = P_ES,
E_ES = E_ES,
D_1 = D_1,
D_2 = D_2,
U_1 = U_1,
U_2 = U_2,
P_EV_1 = P_EV_1,
P_EV_2 = P_EV_2,
E_EV_1 = E_EV_1,
E_EV_2 = E_EV_2,
CND_1 = CND_1,
CND_2 = CND_2,
INR_1 = INR_1,
INR_2 = INR_2,
COP_1 = COP_1,
COP_2 = COP_2,
Q_I_1 = Q_I_1,
Q_I_2 = Q_I_2,
Q_BUS_1 = Q_BUS_1,
Q_BUS_2 = Q_BUS_2,
OH_1 = OH_1,
OH_2 = OH_2,
OΔ_1 = OΔ_1,
OΔ_2 = OΔ_2,
P_AC_1 = P_AC_1,
P_AC_2 = P_AC_2
)
end;
function get_a_self_block(O)::NamedTuple
D = rand(0:5, T) # base demand
U = [rand(1:4, rand(2:5)) for _ = 1:rand(1:4)] # each entry is a cycle vector of an uninterruptible load
P_EV, E_EV= (m = rand((1., 1.5)), M = rand(3:7)), rand(10:39)
CND, INR, COP, Q_I, Q_BUS, OH, OΔ, P_AC = gen_ac_data()
P_BUS = get_self_P_BUS(D, U, P_EV, E_EV, O, CND, INR, COP, Q_I, OH, OΔ, P_AC)
local d = (
P_BUS = P_BUS,
D = D ,
U = U ,
P_EV = P_EV ,
E_EV = E_EV ,
CND = CND ,
INR = INR ,
COP = COP ,
Q_I = Q_I ,
Q_BUS = Q_BUS,
OH = OH ,
OΔ = OΔ ,
P_AC = P_AC
)
end;
function add_a_paired_block!(model, d::NamedTuple)::NamedTuple
# lender house
p_ES, e_ES = add_ES_module!(model, d.P_ES, d.E_ES)
bU_1, pU_1 = add_U_module!(model, d.U_1)
bEV_1, pEV_1, bLent, pLent = add_EV_1_module!(model, d.P_EV_1, d.E_EV_1)
o_1, q_1, pAC_1 = add_AC_module!(model, O, d.CND_1, d.INR_1, d.COP_1, d.Q_I_1, 0, d.OH_1, d.OΔ_1, d.P_AC_1) # Q_BUS = 0
# borrower house
bU_2, pU_2 = add_U_module!(model, d.U_2)
bEV_2, pEV_2 = add_EV_2_module!(model, d.P_EV_2, d.E_EV_2, bLent, pLent)
o_2, q_2, pAC_2 = add_AC_module!(model, O, d.CND_2, d.INR_2, d.COP_2, d.Q_I_2, 0, d.OH_2, d.OΔ_2, d.P_AC_2) # Q_BUS = 0
# circuit breaker pair
pBus_1, pBus_2 = add_circuit_breaker_pair_module!(model, d.P_BUS_1, d.P_BUS_2,
p_ES, d.G, pLent, pEV_1, pU_1, pAC_1, d.D_1,
pEV_2, pU_2, pAC_2, d.D_2)
local x = (
pBus_1 = pBus_1,
pBus_2 = pBus_2,
q_1 = q_1,
q_2 = q_2,
pLent = pLent
)
end;
function add_a_self_block!(model, d::NamedTuple)::NamedTuple
bU, pU = add_U_module!(model, d.U)
bEV, pEV = add_self_EV_module!(model, d.P_EV, d.E_EV)
o, q, pAC = add_AC_module!(model, O, d.CND, d.INR, d.COP, d.Q_I, 0, d.OH, d.OΔ, d.P_AC) # Q_BUS = 0
pBus = add_self_circuit_breaker_module!(model, d.P_BUS, d.D, pU, pEV, pAC)
local x = (
pBus = pBus,
q = q
)
end;
function pre_fill!(model, D, X)
for j = Rng1
local d = get_a_paired_block(O)
push!(D, d)
local x = add_a_paired_block!(model, d)
push!(X, x)
print("\rj = $j")
end
for j = Rng2
local d = get_a_self_block(O)
push!(D, d)
local x = add_a_self_block!(model, d)
push!(X, x)
print("\rj = $j")
end
JuMP.@expression(model, pA, sum(+(X[j].pBus_1, X[j].pBus_2) for j = Rng1) + sum(X[j].pBus for j = Rng2)) # a 24-Vector
JuMP.@expression(model, qA, sum(+(X[j].q_1, X[j].q_2) for j = Rng1) + sum(X[j].q for j = Rng2)) # a 24-Vector
println()
end;
function set_qˈs_ub!(D, X)
for j = Rng1
JuMP.set_upper_bound.(X[j].q_1, D[j].Q_BUS_1)
JuMP.set_upper_bound.(X[j].q_2, D[j].Q_BUS_2)
end
for j = Rng2
JuMP.set_upper_bound.(X[j].q, D[j].Q_BUS)
end
end;
function get_P_A(model, X)
local pA_UB = JuMP.@variable(model)
local c = JuMP.@constraint(model, model[:pA] .≤ pA_UB)
JuMP.@objective(model, Min, pA_UB)
JuMP.set_attribute(model, "SolutionLimit", 15)
JuMP.set_attribute(model, "TimeLimit", 120)
@info "doing the 1st optimize!"
JuMP.optimize!(model); assert_is_approxly_solved(model)
local val = ı(pA_UB)
val > 0 || error("pA_UB has value $val")
JuMP.delete(model, c)
JuMP.delete(model, pA_UB)
local P_A = ceil(Int, val)
end;
function get_Q_A(model) return map(x -> rand(0.25:0.05:0.95)x, ı.(model[:qA])) end;
function get_E_Q(model) return round(Int, rand(0.25:0.05:0.85)sum(ı.(model[:qA]))) end;
function get_global_restrictions()
pre_fill!(model, D, X)
P_A::Int = get_P_A(model, X) # 1️⃣ opt without obj
JuMP.@constraint(model, [t = 1:T], model[:pA][t] ≤ P_A); # add pA limit
JuMP.@objective(model, Min, C ⋅ model[:pA]);
JuMP.set_attribute(model, "SolutionLimit", 1) # since the following optimization might be hard
JuMP.set_attribute(model, "TimeLimit", (20)60)
JuMP.set_attribute(model, "MIPFocus", 1) # to find INT-feasible solution ASAP
@info "doing the 2nd optimize!"
JuMP.optimize!(model); assert_is_approxly_solved(model) # 2️⃣ opt with q = 0 (thus costly)
JuMP.set_attribute(model, "MIPFocus", 0)
JuMP.set_attribute(model, "SolutionLimit", typemax(Int32))
JuMP.set_attribute(model, "TimeLimit", 120)
set_qˈs_ub!(D, X) # start sending cooling power
@info "doing the 3rd optimize!"
JuMP.optimize!(model); assert_is_approxly_solved(model) # 3️⃣ opt with unlimited qA
all(map(x -> x>0, ı.(model[:pA]))) || error("there is a `t` when `pA` < 0");
Q_A::Vector{Float64} = get_Q_A(model)
JuMP.@constraint(model, [t = 1:T], model[:qA][t] ≤ Q_A[t]);
@info "doing the 4th optimize!"
JuMP.optimize!(model); assert_is_approxly_solved(model) # 4️⃣ opt with limited qA
E_Q::Int = get_E_Q(model)
JuMP.@constraint(model, sum(model[:qA]) ≤ E_Q);
error("man interrupt")
@info "doing the 5th optimize!"
JuMP.optimize!(model); assert_is_approxly_solved(model) # 5️⃣ opt with limited Eq
println("Global EV lending $(sum([ı.(X[j].pLent) for j = Rng1]))") # see if there is lending behavior at each time step
return P_A, Q_A, E_Q
end;
const GRB_ENV = Gurobi.Env();
const T = 24;
# Input
Random.seed!(324463511960);
const J = 1000;
# Output
const (C, O) = get_C_and_O(); # price and Celsius vector
const model, D, X = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)), NamedTuple[], NamedTuple[];
const (Rng1, Rng2) = get_pair_and_self_rng(J);
const P_A, Q_A, E_Q = get_global_restrictions()
In less than 20 minutes, you can see ERROR: man interrupt (me introduced manual ERROR to facilitate debugging). Then you type model in REPL to recover the aforesaid phenomenon.