How do you know it can? In that case it will have both more columns and rows.
The code is as follows. There are 2 if-end blocks in the end. The 1st if-end block is using @expression only. The 2nd if-end block is defining additional @variable and @constraint.
Code
import SparseArrays, Random, PGLib; RawDict = PGLib.pglib("pglib_opf_case118_ieee.m"); Random.seed!(1);
function noise() return 1e-6rand(-1:1e-9:1) end;
function noise(N) return 1e-6rand(-1:1e-9:1, N) end;
function perturb!(v) return for i in eachindex(v) v[i] += noise() end end;
function fill_Dmat!(Dmat)
for c in 2:size(Dmat, 2)
a = rand(0.95:7e-9:1.17, D) + rand(-0.017:7e9:0.017, D)
Dmat[:, c] = view(Dmat, :, c-1) .* a
end
end; function write_to!(Dmat, Darray, t)
for c in 1:size(Dmat, 2)
Darray[:, t, c] = view(Dmat, :, c)
end
end; 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]
perturb!(pmax)
perturb!(cost)
G, G2N, pmax, 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]
perturb!(Dref)
D, D2N, Dref
end;
Y, T, Dura, Darray = let
Y, T = 15, 3
Dura = [6, 3, 1]
Darray = Array{Float64, 3}(undef, D, T, Y)
Dmat = Matrix{Float64}(undef, D, Y) # temporary
Random.seed!(2)
Dmiddle = rand(1.1:7e-13:1.3, D) .* Dref
Dlarge = rand(1.4:7e-13:1.7, D) .* Dref
Dmat[:, 1] = Dref; fill_Dmat!(Dmat); write_to!(Dmat, Darray, 1)
Dmat[:, 1] = Dmiddle; fill_Dmat!(Dmat); write_to!(Dmat, Darray, 2)
Dmat[:, 1] = Dlarge; fill_Dmat!(Dmat); write_to!(Dmat, Darray, 3)
Y, T, Dura, Darray
end;
S, Aprob, Avail, 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);
S = 4; Aprob = 1/S * ones(S)
Random.seed!(1)
Avail = sin.(rand(0.37:7e-13:1.3, S, Y))
S, Aprob, Avail, g_hyd, Ginvcost
end;
(B, N) = length.((RawDict["branch"], RawDict["bus"])) # (186, 118)
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[findfirst(x -> x == b, F2B)::Int] # [Fallback] the corresponding JuMP's decision
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
Fra[Fra .> 5] /= 2 # my modification
perturb!(Fra)
F, F2B, Fra
end;
Finvcost = 80abs2.([RawDict["branch"]["$b"]["br_r"] + RawDict["branch"]["$b"]["br_x"]im for b in F2B]);
perturb!(Finvcost);
INVG = 5ones(G) + 1e5noise(G)
INVF = 10ones(F) + 1e5noise(F)
Finvcost *= 5
Ginvcost *= 5
if false # deterministic formulation with NA constr
det = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV));
# 1st-stage variable
JuMP.@variable(det, invf[f = 1:F, s = 1:S], Bin);
JuMP.@variable(det, invg[g = 1:G, y = 1:Y, s = 1:S], Bin);
# Capacity cap
JuMP.@expression(det, fmax[f = 1:F, s = 1:S], Fra[f] + INVF[f]invf[f, s]);
JuMP.@expression(det, pgmax[g = 1:G, y = 1:Y, s = 1:S], invg[g, y, s]INVG[g]);
y = 1; for g = 1:G, s = 1:S
JuMP.add_to_expression!(pgmax[g, y, s], Gcapinit[g])
end; for y = 2:Y, g = 1:G, s = 1:S
JuMP.add_to_expression!(pgmax[g, y, s], pgmax[g, y-1, s])
end;
# operational
JuMP.@variable(det, pf[f = 1:F, t = 1:T, y = 1:Y, s = 1:S]);
JuMP.@variable(det, pg[g = 1:G, t = 1:T, y = 1:Y, s = 1:S] ≥ 0);
JuMP.@constraint(det, [g = 1:G, y = 1:Y, s = 1:S], view(pg, g, :, y, s) .≤ pgmax[g, y, s]);
g = g_hyd; JuMP.@constraint(det, [y = 1:Y, s = 1:S], # 💧
view(pg, g, :, y, s)'Dura ≤ Avail[s, y]pgmax[g, y, s]sum(Dura));
for f = 1:F, s = 1:S
v, B = view(pf, f, :, :, s), fmax[f, s]
JuMP.@constraint(det, -B .≤ v); JuMP.@constraint(det, v .≤ B)
end
# invest cost is w.p.1 due to NA, generation cost is scene-wise
JuMP.@expression(det, invf_cost, Finvcost'view(invf, :, S))
JuMP.@expression(det, invg_cost[y = 1:Y], Ginvcost'view(invg, :, y, S))
JuMP.@expression(det, gen_cost[s = 1:S, y = 1:Y], Gptcost'view(pg, :, :, y, s)Dura)
JuMP.@objective(det, Min, (invf_cost + sum(invg_cost)) + (Aprob'sum(gen_cost; dims = 2))[1])
for s = 1:S-1 # NA
JuMP.@constraint(det, view(invf, :, s) == view(invf, :, S))
JuMP.@constraint(det, view(invg, :, :, s) == view(invg, :, :, S))
end
# network power balance
for t = 1:T, y = 1:Y, s = 1:S, nState = 0:3, n ∈ findall(ø -> ø == nState, n8)
p_gen = (nState < 2 ? false : pg[findfirst(x -> x == n, G2N)::Int, t, y, s])
p_demand = (nState == 1 || nState == 3 ? Darray[findfirst(x -> x == n, D2N)::Int, t, y] : false)
bv, sgv = SparseArrays.findnz(view(A, :, n))
f = b -> bpf(b, view(Darray, :, t, y), view(pg, :, t, y, s), view(pf, :, t, y, s))
JuMP.@constraint(det, sgv'f.(bv) + p_gen == p_demand)
end
JuMP.optimize!(det); JuMP.assert_is_solved_and_feasible(det; allow_local = false);
JuMP.objective_value(det) # ✅ 9177.799530047192
end
if true # deterministic formulation [alternative]
det = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV));
# 1st-stage variable
JuMP.@variable(det, invf[f = 1:F, s = 1:S], Bin);
JuMP.@variable(det, invg[g = 1:G, y = 1:Y, s = 1:S], Bin);
JuMP.@variable(det, Gcapinit[g] ≤ pgmax[g = 1:G, y = 1:Y, s = 1:S] ≤ Gcapinit[g] + INVG[g]y);
y = 1; JuMP.@constraint(det, [g = 1:G, s = 1:S], pgmax[g, y, s] == Gcapinit[g] + INVG[g]invg[g, y, s]);
JuMP.@constraint(det, [g = 1:G, y = 2:Y, s = 1:S], pgmax[g, y, s] == pgmax[g, y-1, s] + INVG[g]invg[g, y, s]);
JuMP.@expression(det, fmax[f = 1:F, s = 1:S], Fra[f] + INVF[f]invf[f, s]);
# operational
JuMP.@variable(det, pf[f = 1:F, t = 1:T, y = 1:Y, s = 1:S]); for f = 1:F, s = 1:S
v, B = view(pf, f, :, :, s), fmax[f, s]
JuMP.@constraint(det, -B .≤ v)
JuMP.@constraint(det, v .≤ B)
end
JuMP.@variable(det, pg[g = 1:G, t = 1:T, y = 1:Y, s = 1:S] ≥ 0);
JuMP.@constraint(det, [g = 1:G, y = 1:Y, s = 1:S], view(pg, g, :, y, s) .≤ pgmax[g, y, s]);
g = g_hyd; JuMP.@constraint(det, [y = 1:Y, s = 1:S], # 💧
view(pg, g, :, y, s)'Dura ≤ Avail[s, y]pgmax[g, y, s]sum(Dura));
# invest cost is w.p.1 due to NA, generation cost is scene-wise
JuMP.@expression(det, invf_cost, Finvcost'view(invf, :, S));
JuMP.@expression(det, invg_cost[y = 1:Y], Ginvcost'view(invg, :, y, S));
JuMP.@expression(det, gen_cost[s = 1:S, y = 1:Y], Gptcost'view(pg, :, :, y, s)Dura);
JuMP.@objective(det, Min, (invf_cost + sum(invg_cost)) + (Aprob'sum(gen_cost; dims = 2))[1]);
for s = 1:S-1 # NA
JuMP.@constraint(det, view(invf, :, s) == view(invf, :, S))
JuMP.@constraint(det, view(invg, :, :, s) == view(invg, :, :, S))
end
# real-time network power balance
for t = 1:T, y = 1:Y, s = 1:S, nState = 0:3, n ∈ findall(ø -> ø == nState, n8)
p_gen = (nState < 2 ? false : pg[findfirst(x -> x == n, G2N)::Int, t, y, s])
p_demand = (nState == 1 || nState == 3 ? Darray[findfirst(x -> x == n, D2N)::Int, t, y] : false)
bv, sgv = SparseArrays.findnz(view(A, :, n))
f = b -> bpf(b, view(Darray, :, t, y), view(pg, :, t, y, s), view(pf, :, t, y, s))
JuMP.@constraint(det, sgv'f.(bv) + p_gen == p_demand)
end
JuMP.optimize!(det); JuMP.assert_is_solved_and_feasible(det; allow_local = false);
JuMP.objective_value(det) # ✅ 9177.799530047192
end
