Hi guys!
I use Gurobi in Julia in JUMP. The model has linear constraints and second-order cone constraints. However, the “termination_status” is：locally_solved, sometimes is numerous error. Why and how to solve this problem?
Thanks very much!
Gurobi performs relatively inadequate on quadratic models, even if they are convex. If your problem is large scale, a more reliable way out is building an MILP model instead.
Is there any source for this generic statement?
Hi @Z_M,
This generally means that your problem has numerical issues which makes it difficult to solve. Can you share a reproducible example, or a log of the solve?
Gurobi performs relatively inadequate on quadratic models
I don’t agree with this statement, but this thread is not an appropriate place to discuss in further detail. Let’s keep to the original question.
1 Like
From my experience.
Now in retrospect:
I find a case where it seems that large scale MILP may run out of memory during branch-and-bound:
OUT_OF_MEMORY 10001 Available memory was exhausted (quoted from Gurobi’s doc).
Gurobi Optimizer version 12.0.3 build v12.0.3rc0 (linux64 - "Debian GNU/Linux 12 (bookworm)")
CPU model: AMD EPYC 7763 64-Core Processor, instruction set [SSE2|AVX|AVX2]
Thread count: 128 physical cores, 256 logical processors, using up to 32 threads
Non-default parameters:
MemLimit 250
Optimize a model with 3393948 rows, 5933377 columns and 17269092 nonzeros
Model fingerprint: 0x4a259666
Variable types: 3549601 continuous, 2383776 integer (2383776 binary)
Coefficient statistics:
Matrix range [9e-01, 2e+01]
Objective range [1e+01, 3e+01]
Bounds range [1e+00, 3e+05]
RHS range [1e+00, 2e+02]
Presolve removed 12966 rows and 55944 columns (presolve time = 5s)...
Presolve removed 40630 rows and 83175 columns (presolve time = 10s)...
Presolve removed 41595 rows and 87934 columns (presolve time = 15s)...
Presolve removed 44899 rows and 91752 columns (presolve time = 20s)...
Presolve removed 45739 rows and 92437 columns (presolve time = 25s)...
Presolve removed 45992 rows and 92537 columns (presolve time = 30s)...
Presolve removed 45992 rows and 92537 columns (presolve time = 36s)...
Presolve removed 45992 rows and 92544 columns (presolve time = 40s)...
Presolve removed 45992 rows and 92732 columns (presolve time = 46s)...
Presolve removed 45992 rows and 93608 columns (presolve time = 50s)...
Presolve removed 48505 rows and 95421 columns (presolve time = 55s)...
Presolve removed 48621 rows and 95633 columns (presolve time = 61s)...
Presolve removed 48621 rows and 95633 columns
Presolve time: 62.18s
Presolved: 3345327 rows, 5837744 columns, 16964820 nonzeros
Variable types: 3502448 continuous, 2335296 integer (2335295 binary)
Deterministic concurrent LP optimizer: primal simplex, dual simplex, and barrier
Showing barrier log only...
Root barrier log...
Ordering time: 1.09s
Barrier statistics:
AA' NZ : 7.315e+06
Factor NZ : 3.996e+07 (roughly 3.0 GB of memory)
Factor Ops : 1.168e+09 (less than 1 second per iteration)
Threads : 29
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 2.72632712e+10 -3.26416868e+09 7.42e+06 2.80e+01 1.27e+05 107s
1 6.27986724e+09 -3.45379732e+09 1.68e+06 1.03e+02 2.93e+04 109s
2 6.92548755e+08 -3.51052539e+09 1.52e+05 4.55e-13 2.97e+03 111s
3 2.11111142e+08 -2.65351771e+09 2.01e+04 4.55e-13 5.66e+02 112s
4 1.93850002e+08 -2.58674553e+09 1.59e+04 4.55e-13 4.93e+02 114s
5 1.40556522e+08 -1.96058994e+09 6.19e+03 4.38e-12 2.73e+02 115s
6 1.17740242e+08 -8.05020290e+08 2.06e+03 2.73e-12 1.02e+02 116s
7 9.17762445e+07 -3.39505703e+08 4.82e+02 1.36e-12 4.30e+01 117s
8 8.12928226e+07 -1.67809061e+08 3.56e+02 7.39e-13 2.47e+01 118s
9 5.96653530e+07 -4.75508089e+07 1.44e+02 6.14e-12 1.05e+01 119s
10 5.04065968e+07 -1.21775821e+07 8.31e+01 3.64e-12 6.08e+00 120s
11 4.42367486e+07 8.70786390e+06 4.24e+01 2.05e-12 3.44e+00 122s
12 4.19905828e+07 1.91516100e+07 2.89e+01 1.36e-12 2.21e+00 123s
13 3.90162016e+07 2.62676489e+07 1.14e+01 1.14e-12 1.23e+00 124s
14 3.81253940e+07 2.96461242e+07 7.28e+00 6.75e-13 8.17e-01 125s
15 3.74602044e+07 3.27208865e+07 4.49e+00 6.42e-13 4.57e-01 126s
16 3.70012648e+07 3.40309311e+07 2.72e+00 1.08e-12 2.86e-01 127s
17 3.66609120e+07 3.50965392e+07 1.49e+00 5.65e-13 1.51e-01 129s
18 3.65111372e+07 3.54821066e+07 9.79e-01 9.09e-13 9.91e-02 130s
19 3.63737373e+07 3.57864427e+07 5.21e-01 7.96e-13 5.65e-02 131s
20 3.63096547e+07 3.59293547e+07 3.15e-01 7.73e-13 3.66e-02 132s
21 3.62738835e+07 3.60385358e+07 2.02e-01 8.30e-13 2.27e-02 134s
22 3.62389780e+07 3.61250377e+07 9.43e-02 8.46e-13 1.10e-02 135s
23 3.62224874e+07 3.61685221e+07 4.49e-02 9.74e-13 5.19e-03 137s
24 3.62151594e+07 3.61816384e+07 2.54e-02 4.99e-13 3.23e-03 138s
25 3.62104881e+07 3.61911860e+07 1.32e-02 8.52e-13 1.86e-03 139s
26 3.62079154e+07 3.61962113e+07 6.65e-03 9.63e-13 1.13e-03 140s
27 3.62064025e+07 3.62014524e+07 2.86e-03 9.45e-13 4.76e-04 141s
28 3.62055122e+07 3.62040836e+07 8.48e-04 7.54e-13 1.37e-04 143s
29 3.62054168e+07 3.62046331e+07 6.50e-04 1.00e-12 7.54e-05 144s
30 3.62051321e+07 3.62049329e+07 8.10e-05 1.45e-12 1.91e-05 146s
31 3.62050972e+07 3.62050592e+07 1.53e-05 1.57e-12 3.65e-06 147s
32 3.62050860e+07 3.62050850e+07 1.82e-07 9.86e-13 9.74e-08 148s
33 3.62050857e+07 3.62050857e+07 2.84e-09 1.31e-12 7.30e-10 150s
Barrier solved model in 33 iterations and 149.63 seconds (121.34 work units)
Optimal objective 3.62050857e+07
Root crossover log...
427815 DPushes remaining with DInf 0.0000000e+00 153s
149 DPushes remaining with DInf 0.0000000e+00 155s
0 DPushes remaining with DInf 0.0000000e+00 155s
690138 PPushes remaining with PInf 0.0000000e+00 156s
98299 PPushes remaining with PInf 0.0000000e+00 160s
38890 PPushes remaining with PInf 0.0000000e+00 165s
11374 PPushes remaining with PInf 0.0000000e+00 170s
0 PPushes remaining with PInf 0.0000000e+00 172s
Push phase complete: Pinf 0.0000000e+00, Dinf 5.4141903e-10 172s
Root simplex log...
Iteration Objective Primal Inf. Dual Inf. Time
1062977 3.6205086e+07 0.000000e+00 0.000000e+00 173s
Concurrent spin time: 9.39s (can be avoided by choosing Method=3)
Solved with barrier
1062977 3.6205086e+07 0.000000e+00 4.391189e+03 181s
1063445 3.6205086e+07 0.000000e+00 0.000000e+00 183s
Extra simplex iterations after uncrush: 468
Root relaxation: objective 3.620509e+07, 1063445 iterations, 102.37 seconds (72.63 work units)
Total elapsed time = 240.06s (DegenMoves)
Total elapsed time = 254.70s (DegenMoves)
Total elapsed time = 267.67s (DegenMoves)
Total elapsed time = 280.55s (DegenMoves)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 3.6205e+07 0 40110 - 3.6205e+07 - - 294s
0 0 3.6215e+07 0 47119 - 3.6215e+07 - - 1602s
0 0 3.6218e+07 0 47750 - 3.6218e+07 - - 1687s
0 0 3.6218e+07 0 47881 - 3.6218e+07 - - 1728s
0 0 3.6218e+07 0 47992 - 3.6218e+07 - - 1768s
0 0 3.6218e+07 0 47976 - 3.6218e+07 - - 1805s
0 0 3.6224e+07 0 43269 - 3.6224e+07 - - 1888s
0 0 3.6225e+07 0 43027 - 3.6225e+07 - - 1967s
0 0 3.6225e+07 0 43234 - 3.6225e+07 - - 2027s
0 0 3.6225e+07 0 43298 - 3.6225e+07 - - 2091s
0 0 3.6225e+07 0 43288 - 3.6225e+07 - - 2143s
0 0 3.6225e+07 0 43329 - 3.6225e+07 - - 2199s
0 0 3.6227e+07 0 40759 - 3.6227e+07 - - 4317s
0 0 3.6227e+07 0 40272 - 3.6227e+07 - - 4461s
0 0 3.6227e+07 0 40482 - 3.6227e+07 - - 4561s
0 0 3.6227e+07 0 40637 - 3.6227e+07 - - 4681s
0 0 3.6228e+07 0 37435 - 3.6228e+07 - - 6865s
0 0 3.6229e+07 0 36916 - 3.6229e+07 - - 7032s
0 0 3.6229e+07 0 37277 - 3.6229e+07 - - 7149s
0 0 3.6229e+07 0 37451 - 3.6229e+07 - - 7269s
0 0 3.6229e+07 0 36313 - 3.6229e+07 - - 8023s
0 0 3.6229e+07 0 36528 - 3.6229e+07 - - 8103s
0 0 3.6229e+07 0 36696 - 3.6229e+07 - - 8119s
0 0 3.6229e+07 0 35760 - 3.6229e+07 - - 8644s
0 0 3.6229e+07 0 35954 - 3.6229e+07 - - 8689s
0 0 3.6229e+07 0 36085 - 3.6229e+07 - - 8701s
0 0 3.6229e+07 0 35438 - 3.6229e+07 - - 9204s
0 0 3.6229e+07 0 35762 - 3.6229e+07 - - 9254s
0 0 3.6229e+07 0 35803 - 3.6229e+07 - - 9266s
0 0 3.6229e+07 0 35419 - 3.6229e+07 - - 9314s
0 0 3.6229e+07 0 35586 - 3.6229e+07 - - 9353s
0 0 3.6229e+07 0 35358 - 3.6229e+07 - - 9388s
0 0 3.6229e+07 0 28502 - 3.6229e+07 - - 9477s
0 2 3.6229e+07 0 27316 - 3.6229e+07 - - 9773s
1 4 3.6229e+07 1 27317 - 3.6229e+07 - 32700 9804s
3 8 3.6229e+07 2 29150 - 3.6229e+07 - 13143 9831s
7 16 3.6229e+07 3 29651 - 3.6229e+07 - 9042 9874s
11 16 3.6229e+07 3 29769 - 3.6229e+07 - 7010 9875s
15 32 3.6229e+07 4 29747 - 3.6229e+07 - 5959 10365s
Cutting planes:
Learned: 36
Gomory: 116
Cover: 10026
Implied bound: 3153
Clique: 1109
MIR: 32167
StrongCG: 515
Flow cover: 35127
GUB cover: 208
Zero half: 957
Network: 2
RLT: 7645
Relax-and-lift: 2605
BQP: 5
PSD: 6
Explored 31 nodes (2022582 simplex iterations) in 10397.66 seconds (6953.59 work units)
Thread count was 32 (of 256 available processors)
Solution count 0
Solve interrupted (error code 10001)
Best objective -, best bound 3.622914834771e+07, gap -
User-callback calls 718823, time in user-callback 1.31 sec
Note that this is the second time that I run this MILP test case. In a previous test I didn’t set a MemLimit, and the julia REPL is killed by zsh. I was very surprised—since I have 256G RAM (I didn’t monitor it, as it’s slow). So I run this second time with MemLimit = 250 (GB). And it terminates with Error code 10001 without even finding a feasible solution.
code
import JuMP, Gurobi
import JuMP.value as ı
import LinearAlgebra.⋅ as ⋅
import Random
import Dates.now as now
macro get_int_decision(model, expr) return esc(quote
let e = JuMP.@expression($model, $expr), a
a = map(_ -> JuMP.@variable($model, integer = true), e)
JuMP.@constraint($model, a == e)
a
end
end) end;
function get_Bool_value(x)
f = z -> round(Bool, ı(z))
ndims(x) == 0 && return f(x)
map(f, x)
end;
function get_C_and_O()
C = [
[10, 9, 9, 9, 10, 12, 15, 18, 20, 18, 16, 15, 14, 15, 16, 18, 22, 28, 32, 30, 26, 20, 15, 12],
[12, 11, 11, 12, 14, 18, 24, 26, 22, 18, 15, 14, 13, 14, 18, 24, 30, 34, 32, 28, 22, 18, 15, 13],
[16, 15, 14, 14, 15, 18, 24, 30, 28, 22, 18, 12, 10, 12, 16, 22, 28, 34, 36, 32, 28, 24, 20, 18],
[8, 8, 8, 9, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 36, 34, 30, 26, 20, 14],
[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 = [
[28,28,27,27,28,29,31,33,35,36,37,38,38,38,37,36,35,34,32,31,30,29,29,28],
[29,28,28,28,29,31,33,35,37,39,40,41,42,42,41,40,38,36,34,33,32,31,30,29],
[27,27,27,27,28,29,30,32,33,34,35,35,35,35,34,33,32,31,30,29,28,28,28,27],
[32,32,31,31,32,34,36,38,40,42,43,44,44,44,43,42,41,40,38,37,36,35,34,33],
[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 get_honest_model(new = false)
m = new ? JuMP.Model(Gurobi.Optimizer) : JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
JuMP.set_silent(m)
JuMP.set_attribute(m, "MIPGap", 0)
JuMP.set_attribute(m, "MIPGapAbs", 0)
JuMP.set_attribute(m, "Threads", 1) # also for the master LP, since we apply a FULLY async algorithm!
m
end;
function get_pair_and_self_Rng(J)
d = J÷4
1:d, d+1:J # Rng1, Rng2
end;
function prev(t, d, T) return (n = t-d; n<1 ? T+n : n) end;
function pc_P_AC(O, OH, CND, Q_I, COP) return ceil(Int, ((maximum(O) - OH)CND + maximum(Q_I)) / COP) end;
function get_E_ES(Rng)::NamedTuple
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
(; i, M, e)
end;
function add_ES_module!(model, P_ES, E_ES)
pES = ( # eES is dependent variable
c = JuMP.@variable(model, [1:T], lower_bound = 0),
d = JuMP.@variable(model, [1:T], lower_bound = 0)
); bES = JuMP.@variable(model, [1:T], Bin)
JuMP.@constraint(model, pES.c ≤ (bES)P_ES.c)
JuMP.@constraint(model, pES.d ≤ (1 .- bES)P_ES.d)
eES = 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], pES.c[t]*.95 - pES.d[t]/.95 == eES[t]-(t>1 ? eES[t-1] : E_ES.i))
return pES, bES, eES
# NamedTuple(k => ı.(v) for (k, v) = pairs(pES))
end;
function gen_ac_data()::Tuple
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 = pc_P_AC(O, OH, CND, Q_I, COP)
return CND, INR, COP, Q_I, Q_BUS, OH, OΔ, P_AC
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 # pAC is dependent variable
end;
function add_U_module!(model, U)
bU::Matrix{JuMP.VariableRef} = JuMP.@variable(model, [t = 1:T, i = eachindex(U)], Bin)
JuMP.@constraint(model, sum(bU; dims = 1) .≥ true)
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 # bEV can be _inferred from_ pEV
end;
function pc_self_P_BUS(D, U, P_EV, E_EV, O, CND, INR, COP, Q_I, OH, OΔ, P_AC)::Int
model = get_honest_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.termination_status(model) == JuMP.OPTIMAL || error("$(JuMP.termination_status(model))")
val = JuMP.objective_value(model)
val > 0 || error("The self household has P_BUS = $val")
ceil(Int, val) # P_BUS
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 = get_honest_model() # for a block who has a lender and a borrower house
# 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
pES, bES, eES = add_ES_module!(model, P_ES, E_ES)
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
pBus_2 = JuMP.@variable(model, [1:T], lower_bound = 0)
temp_x = JuMP.@variable(model)
temp_c = JuMP.@constraint(model, pBus_2 .== temp_x)
JuMP.@constraint(model, pBus_2 ≥ pEV_2 + pU_2 + pAC_2 + D_2)
JuMP.@objective(model, Min, temp_x)
JuMP.optimize!(model)
JuMP.termination_status(model) == JuMP.OPTIMAL || error("$(JuMP.termination_status(model))")
temp_float64 = ı(temp_x)
temp_float64 > 0 || error("common pBus_2 has value $temp_float64")
P_BUS_2 = ceil(Int, temp_float64)
JuMP.delete(model, temp_c)
JuMP.delete(model, temp_x)
JuMP.set_upper_bound.(pBus_2, P_BUS_2)
temp_x = JuMP.@variable(model) # reuse the local name
JuMP.@constraint(model, -temp_x .≤ pES.c -pES.d -G + pLent + pEV_1 + pU_1 + pAC_1 + D_1)
JuMP.@constraint(model, temp_x .≥ pES.c -pES.d -G + pLent + pEV_1 + pU_1 + pAC_1 + D_1)
JuMP.@objective(model, Min, temp_x)
JuMP.optimize!(model)
JuMP.termination_status(model) == JuMP.OPTIMAL || error("$(JuMP.termination_status(model))")
temp_float64 = ı(temp_x)
temp_float64 > -1e-5 || error("pBus_1 has value $temp_float64")
P_BUS_1 = max(1, ceil(Int, temp_float64))
(;P_BUS_1, P_BUS_2, G, P_ES, E_ES, D_1, D_2, U_1, U_2, P_EV_1, P_EV_2,
E_EV_1, E_EV_2, CND_1, CND_2, INR_1, INR_2, COP_1, COP_2, Q_I_1, Q_I_2,
Q_BUS_1, Q_BUS_2, OH_1, OH_2, OΔ_1, OΔ_2, P_AC_1, 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 = pc_self_P_BUS(D, U, P_EV, E_EV, O, CND, INR, COP, Q_I, OH, OΔ, P_AC)
(;P_BUS, D, U, P_EV, E_EV, CND, INR, COP, Q_I, Q_BUS, OH, OΔ, P_AC)
end;
function add_a_paired_block!(model, d::NamedTuple)::NamedTuple
# lender house
pES, bES, eES = 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,
pES, d.G, pLent, pEV_1, pU_1, pAC_1, d.D_1,
pEV_2, pU_2, pAC_2, d.D_2)
(;pBus_1, pBus_2, bLent, bES, bEV_1, bEV_2, bU_1, bU_2, q_1, q_2)
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)
(;pBus, bEV, bU, q)
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 fill_model_D_X!(monolithic_model, D, X)
for (ȷ, g, a) = zip((Rng1, Rng2), (get_a_paired_block, get_a_self_block), (add_a_paired_block!, add_a_self_block!))
for j = ȷ
d = g(O)
push!(D, d)
push!(X, a(monolithic_model, d))
print("\rj = $j")
end
end
end;
const K = 80;
const LOGIS = 255;
Random.seed!(9012493673963);
const J = (K)LOGIS;
const GRB_ENV = Gurobi.Env();
const T = 24;
const (Rng1, Rng2) = get_pair_and_self_Rng(J);
const (C, O) = get_C_and_O(); # price and Celsius vector
const D, X = NamedTuple[], NamedTuple[];
model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)) # Monolithic Gurobi One Shot
JuMP.set_attribute(model, "MemLimit", 250)
fill_model_D_X!(model, D, X)
e = JuMP.@variable(model);
JuMP.@expression(model, pA, sum(X[j].pBus_1 + X[j].pBus_2 for j = Rng1) + sum(X[j].pBus for j = Rng2));
JuMP.@constraint(model, pA .≤ e);
JuMP.@objective(model, Min, e);
@info "Do the 1st optimize!, you can early break it via ctrl+c"
JuMP.optimize!(model);
if JuMP.primal_status(model) == JuMP.FEASIBLE_POINT
P_AGR_BASE = JuMP.objective_value(model)
end
LIFT = 1.0
P_AGR_BASE = 301000
JuMP.set_upper_bound(e, (LIFT)P_AGR_BASE)
JuMP.@objective(model, Min, C ⋅ pA)
@info "Do the 2nd optimize!, you can early break it via ctrl+c"
JuMP.optimize!(model)
set_qˈs_ub!(D, X)
JuMP.@expression(model, qA, sum(X[j].q_1 + X[j].q_2 for j = Rng1) + sum(X[j].q for j = Rng2))
@info "Do the 3rd optimize!, RUN OUT OF MEMORY ERROR CODE 10001"
JuMP.optimize!(model)
See https://support.gurobi.com/hc/en-us/articles/360013195772-How-do-I-avoid-an-out-of-memory-condition
Your best bet is probably to reduce the number of threads. If you have 32 threads and 256 G memory, then each thread has only 8 G of memory.
If you have more questions on this, I’ll ask that you start a new thread, rather than posting at the bottom of someone else’s question.