It’s a bit lengthy but hopefully it is valuable. It’s a runnable Dantzig-Wolfe/Lagrangian decomposition algorithm (can converges to the Lagrangian dual bound). The relevant line is commented with a # 🍅🍅🍅
Block Decomposition pglib_opf_case3012wp_k
import LinearAlgebra.norm as norm
import LinearAlgebra.⋅ as ⋅
import SparseArrays, PGLib; RawDict = PGLib.pglib("pglib_opf_case3012wp_k.m");
import Random; Random.seed!(1);
macro load_data()
esc(quote
T, H = 3, 24;
daynumvec = [91, 183, 91] # length(~) == T && sum(~) == 365
# Network's technical data
(B, N) = length.((RawDict["branch"], RawDict["bus"])); # (3572, 3012)
BrateA = [RawDict["branch"]["$b"]["rate_a"] for b = 1:B]; # we'll extract a subvector as Fpmaxy0
N2RawN = collect(1:3013); popat!(N2RawN, 3009); # length(N2RawN) == N != RawN
bft, A = let # load a rectified bft
bft = [findfirst(x -> x == RawDict["branch"]["$r"][c == 1 ? "f_bus" : "t_bus"],
N2RawN) for r in 1:B, c in 1:2]; # extrema(bft) == (1, N) && size(bft) == (B, 2)
# Don't execute to save time assert_is_connected(bft)
rectify_dir!(bft)
A = bft_2_A(bft)
bft, A
end;
D, D2N, Dref = let
D = length(RawDict["load"]); # 2271
D2N = [findfirst(x -> x == RawDict["load"]["$d"]["load_bus"], N2RawN) for d in 1:D];
@assert length(Set(D2N)) == length(D2N) == D
@assert D2N == sort(D2N)
Dref = let
Dref = [RawDict["load"]["$d"]["pd"] for d in 1:D]; Dref[Dref .< 0.0001] .= 0.0001;
Dref += rand(-1e-5:1.7e-12:1e-5, D) # disturb!
Dref[Dref .< 0.0001] *= 10;
Dref
end;
D, D2N, Dref
end;
G, G2N, Gpmaxy0, InvG = let
G_gw = length(RawDict["gen"]) - 6 # we drop 6 gen's due to zero capacity
Gpmaxy0_gw = [RawDict["gen"]["$g"]["pmax"] for g in 1:G_gw] # capacity for each unaggregated gen
@assert minimum(Gpmaxy0_gw) == 0.0011
G_generatorwise2N = [findfirst(x -> x == RawDict["gen"]["$g"]["gen_bus"], N2RawN) for g in 1:G_gw] # length = 496
G2N = sort(collect(Set(G_generatorwise2N))) # length = 341
G = length(G2N) # G == G_buswise
Gpmaxy0 = [sum(Gpmaxy0_gw[findall(x -> x == n, G_generatorwise2N)]) for n in G2N]
Gpmaxy0 += rand(-1e-4:1.7e-11:1e-4, G)
InvG = rand(0.5:1.7e-12:0.67, G) .* Gpmaxy0 # newly_added_capacity := InvG[g] * invg[g]::Int
G, G2N, Gpmaxy0, InvG
end;
# Network's simplifications
DGN2N, BareN2N = D2N ∩ G2N, setdiff(1:N, D2N ∪ G2N);
DonlyN2N, GonlyN2N = setdiff(D2N, DGN2N), setdiff(G2N, DGN2N); # @assert +(length(BareN2N), length(DonlyN2N), length(GonlyN2N), length(DGN2N)) == N
n8, b8 = zeros(Int, N), ones(Int, B); # 📗 n8 = let
# [0]: a normal bus
# [-2]: 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
# [-4]: invalid bus
# end; b8 = let
# [1]: a normal line
# [-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
# [-4]: invalid branch
# end;
# make minor simplification to the network
begin # This code removes all bare nodes whose degree = 1, by set!(-4)
bv, nv, mv = Int[], Int[], Int[]
for n in BareN2N # depth = 0
adj_n_vec = adjnode(n)
length(adj_n_vec) == 0 && error("No here")
if length(adj_n_vec) == 1
m = adj_n_vec[1] # the other node
b = get_b(n, m)
push!(bv, b)
push!(nv, n)
m ∈ BareN2N && push!(mv, m) # make provision for the next deeper layer
end
end; n8[nv] .= -4; b8[bv] .= -4; # mark as invalid
asubsetofBareN = copy(mv)
invalid_node_vec = findall(x -> x == -4, n8)
bv, nv, mv = Int[], Int[], Int[]
for n in asubsetofBareN # depth = 1
adj_n_vec = setdiff(adjnode(n), invalid_node_vec)
length(adj_n_vec) == 0 && error("No here")
if length(adj_n_vec) == 1
m = adj_n_vec[1] # the other node
b = get_b(n, m)
push!(bv, b)
push!(nv, n)
m ∈ BareN2N && push!(mv, m)
end
end; n8[nv] .= -4; b8[bv] .= -4; # mark as invalid
end
begin # This code marks all Gen-only-degree-1 nodes, by set!(-3)
bv, nv = Int[], Int[]
for n in GonlyN2N
adj_n_vec = filter(x -> n8[x] != -4, adjnode(n))
length(adj_n_vec) == 0 && error("No here")
if length(adj_n_vec) == 1
m = adj_n_vec[1] # the other node
b = get_b(n, m)
push!(nv, n)
push!(bv, b)
end
end; b8[bv] .= -3; n8[nv] .= -3;
end
begin # This code marks all Demand-only-degree-1 nodes, by set!(-2)
bv, nv = Int[], Int[]
for n in DonlyN2N
adj_n_vec = filter(x -> n8[x] != -4, adjnode(n))
length(adj_n_vec) == 0 && error("No here")
if length(adj_n_vec) == 1
m = adj_n_vec[1] # the other node
b = get_b(n, m)
push!(nv, n)
push!(bv, b)
end
end; b8[bv] .= -2; n8[nv] .= -2;
end
F2B = findall(x -> x == 1, b8); F = length(F2B); # only need F decisions representing pf
Fpmaxy0 = BrateA[F2B]; # only need Capacity of those (having pf decision) branches
InvF = rand(0.75:1.7e-12:1.25, F) .* Fpmaxy0
KCLnodes = findall(x -> x == 0, n8); # only need KCLconstrs at these nodes
# Network's economical data
Gphcost, Ginvcost, Finvcost = let
Gphcost = rand(0.00017:1.7e-12:0.00223, G)
Ginvcost = rand(117.3:1.3e-10:1115.7, G) .* Gphcost
M = maximum(Ginvcost)
Finvcost = rand(0.7M:1.3e-11:1.7M, F)
Gphcost, Ginvcost, Finvcost
end
Y, iY = 30, 5;
iY2Y = collect(1:(Y ÷ iY):Y); # years that we make investment
Y2iY = repeat(1:iY; inner = Y ÷ iY); # call `iY = Y2iY[y]`
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; function bpf_genonly(b, pgvec)
n, m = bft[b, :] # bpf dir is defined as n -> m
n2g = node -> findfirst(x -> x == node, G2N)
is_source = node -> n8[node] == -3
if is_source(n)
is_source(m) && error("No here ioyhbgopiueryt54y")
return pgvec[n2g(n)]
elseif is_source(m)
return -pgvec[n2g(m)]
else
error("No here -oiysdg0wieuy62udfh")
end
end; function bpf_Donly(b, Dvec)
n, m = bft[b, :] # bpf dir is defined as n -> m
n2d = node -> findfirst(x -> x == node, D2N)
is_load = node -> n8[node] == -2
if is_load(m)
is_load(n) && error("No here fopi7yhw9e8r769")
return Dvec[n2d(m)]
elseif is_load(n)
return -Dvec[n2d(n)]
else
error("No here lkipc709tyerywr76672w46resy")
end
end; function bpf(b, Dvec, pgvec, fvec) # dispatch function of power flow on branch b
b8[b] == -4 && return false
b8[b] == -3 && return bpf_genonly(b, pgvec)
b8[b] == -2 && return bpf_Donly(b, Dvec)
return fvec[findfirst(x -> x == b, F2B)] # [Fallback] the corresponding JuMP's decision
end;
@load_data();
import JuMP, Gurobi; GRB_ENV = Gurobi.Env(); function JuMP_add_constr_KCL(pD, pg, pf) # Demand data vec, pg and pf decision vec
m = JuMP.owner_model(pf[1])
for n ∈ KCLnodes
p_demand = (n ∈ D2N ? pD[findfirst(x -> x == n, D2N)] : false)
p_gen = (n ∈ G2N ? pg[findfirst(x -> x == n, G2N)] : false)
bv, sgv = SparseArrays.findnz(view(A, :, n))
d = b -> bpf(b, pD, pg, pf) # dispatch function
JuMP.@constraint(m, sgv ⋅ d.(bv) + p_gen == p_demand) # at node n
end
end;
Random.seed!(D);
Dref = rand(0.0006878765231343202:1.7345234e-12:0.013646062201588, D);
Gpmaxy0, Fpmaxy0, InvF, iF, iF2F = let
model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(model)
JuMP.@variable(model, false ≤ pg[g = 1:G]);
JuMP.@variable(model, pf[f = 1:F]);
JuMP.@variable(model, false ≤ pfb[f = 1:F]);
JuMP.@constraint(model, pf .<= pfb);
JuMP.@constraint(model, pf .>= -pfb);
D⅁ = (1.0 + 0.0977(T-1)) * (1 + 0.01Y)Dref;
JuMP_add_constr_KCL(D⅁, pg, pf)
JuMP.@objective(model, Min, sum(pfb));
JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false) # Devise data by Gurobi
Gpmaxy0 = 1.4 * JuMP.value.(pg);
Fpmaxy0 = JuMP.value.(pfb);
# 💡 devise a deficient system at initial state
iF = 5; # suppose we select `iF` lines to invest
ind_vec = sortperm(Fpmaxy0);
iF2F = ind_vec[end-(iF-1):end] # ⚠️ Note that this is orderless
drop_ratio = 0.1 # This is all-time feasible!
drop_ratio = 0.2 # This is already INfeasible at (y = 20, t = 3, h = 22)
drop_cap = drop_ratio * Fpmaxy0[iF2F];
Random.seed!(3iF);
InvF = rand(0.57:1.7e-12:0.77, iF) .* drop_cap; # we are not investing only once, therefore define a <1 ratio
Fpmaxy0[iF2F] *= (1 - drop_ratio);
Gpmaxy0, Fpmaxy0, InvF, iF, iF2F
end;
Darray = Array{Float64}(undef, D, H, T, Y); Random.seed!(2D); for y = 1:Y
Dref_y = (1.0 + 0.01Y * (y-1)/(Y-1))Dref
for t = 1:T
le, ue = 0.63 + 0.017(t-1), 1.0 + 0.0977(t-1)
for h = 1:H
Darray[:, h, t, y] = rand(le:1.73e-12:ue, D) .* Dref_y
end
end
end # generate Load data
# Master and SubBlock Assets
function drawmax!(udv, θ) # undrawn_bitvector, JuMP.value.(θ)
still_has, max_i = false, 0
if any(udv)
still_has, max_v = true, -Inf
for (i, undrawn) in enumerate(udv)
if undrawn && max_v < θ[i]
max_v, max_i = θ[i], i
end
end
udv[max_i] = false # write! in place
end
return still_has, max_i # `max_i` is valid only if `still_has` == true
end; function fˈs_pen_4(ȷ, π::Array{Float64}, fpmax) # [idle] for block ȷ
I, Y = size(π)
m = JuMP.owner_model(fpmax[1])
if ȷ == 1
return JuMP.@expression(m, sum(fpmax[ı, 2, ȷ]π[ı, ȷ] for ı = 1:I))
elseif ȷ == Y+1
return JuMP.@expression(m, - sum(fpmax[ı, 1, ȷ]π[ı, ȷ-1] for ı = 1:I))
end
return JuMP.@expression(m, sum(fpmax[ı, 2, ȷ]π[ı, ȷ] - fpmax[ı, 1, ȷ]π[ı, ȷ-1] for ı = 1:I))
end; macro _pen_code()
esc(quote
I, Y = size(π)
if ȷ == 1
return JuMP.@expression(m, sum(fpmax_ȷ[ı, 2]π[ı, ȷ] for ı = 1:I))
elseif ȷ == Y+1
return JuMP.@expression(m, - sum(fpmax_ȷ[ı, 1]π[ı, ȷ-1] for ı = 1:I))
end
return JuMP.@expression(m, sum(fpmax_ȷ[ı, 2]π[ı, ȷ] - fpmax_ȷ[ı, 1]π[ı, ȷ-1] for ı = 1:I))
end)
end; function fˈs_pen_4_(ȷ, π, fpmax_ȷ::Array{JuMP.VariableRef}) # for block ȷ
m = JuMP.owner_model(fpmax_ȷ[1])
@_pen_code()
end; function fˈs_pen_4_(ȷ, π::Array{JuMP.VariableRef}, fpmax_ȷ) # for block ȷ
m = JuMP.owner_model(π[1])
@_pen_code()
end;
# master
COT = 1e-8
an_UB = 6.7357; out = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # 🟣
JuMP.@variable(out, π[ı = 1:iF, ȷ = 1:iY-1]); # Care θ[ȷ]---π Relation
JuMP.@variable(out, θ[ȷ = 1:iY]);
JuMP.@expression(out, out_obj_tbMAX, sum(θ)); JuMP.@constraint(out, out_obj_tbMAX ≤ an_UB);
JuMP.@objective(out, Max, out_obj_tbMAX);
# subblocks
function Min_ȷ_fixed_part()
Min_ȷ = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # ✅ abs_y = iY2Y[ȷ]-1 + y
JuMP.@variable(Min_ȷ, invf_ȷ[ı = 1:iF], Bin);
JuMP.@expression(Min_ȷ, prim_obj_tbMIN_ȷ, sum(invf_ȷ));
JuMP.@variable(Min_ȷ, fpmax_ȷ[ı = 1:iF, b = 1:2]);
JuMP.set_lower_bound.(view(fpmax_ȷ, :, 2), Fpmaxy0[iF2F]); # lower bounded by the zero budget case
JuMP.@constraint(Min_ȷ, view(fpmax_ȷ,:,2) - view(fpmax_ȷ,:,1) == InvF .* invf_ȷ);
JuMP.@variable(Min_ȷ, false ≤ pg_ȷ[g = 1:G, h = 1:H, t = 1:T, y = 1:Y÷iY] ≤ Gpmaxy0[g]);
JuMP.@variable(Min_ȷ, pf_ȷ[f = 1:F, h = 1:H, t = 1:T, y = 1:Y÷iY]);
for f = 1:F # RateA limit
if f ∈ iF2F
⅁ = fpmax_ȷ[findfirst(x -> x == f, iF2F), 2]
JuMP.@constraint(Min_ȷ, view(pf_ȷ, f, :, :, :) .≥ -⅁)
JuMP.@constraint(Min_ȷ, view(pf_ȷ, f, :, :, :) .≤ ⅁)
else
⅁ = Fpmaxy0[f]
JuMP.set_lower_bound.(view(pf_ȷ, f, :, :, :), -⅁)
JuMP.set_upper_bound.(view(pf_ȷ, f, :, :, :), ⅁)
end
end
return Min_ȷ
end
Min_vec = [Min_ȷ_fixed_part() for _ = 1:iY]; # 🍅🍅🍅
for ȷ = 1:iY # finish building all sub_models
Min_ȷ = Min_vec[ȷ]
fpmax_ȷ = Min_ȷ[:fpmax_ȷ]
pg_ȷ = Min_ȷ[:pg_ȷ]
pf_ȷ = Min_ȷ[:pf_ȷ]
for y = 1:Y÷iY, t = 1:T, h = 1:H # Network KCL
JuMP_add_constr_KCL(
view(Darray, :, h, t, iY2Y[ȷ]-1 + y),
view(pg_ȷ, :, h, t, y),
view(pf_ȷ, :, h, t, y)
)
end
# Set necessary bounds
if ȷ == 1
JuMP.fix.(view(fpmax_ȷ, :, 1), Fpmaxy0[iF2F])
else
JuMP.set_upper_bound.(view(fpmax_ȷ, :, 2), Fpmaxy0[iF2F] + (ȷ)InvF)
end
end
# Main loop
JuMP.set_silent(out)
JuMP.set_silent.(Min_vec)
for ite = 1:typemax(Int)
JuMP.optimize!(out); JuMP.assert_is_solved_and_feasible(out; allow_local = false)
ub = JuMP.objective_bound(out)
udv, Θ, Π = trues(iY), JuMP.value.(θ), JuMP.value.(π)
@info "ite = $ite ▶ ub = $ub, Θ is $Θ; Π is $Π"
outer_saturation = true
while true
still_has, ȷ = drawmax!(udv, Θ)
if still_has == false
@info "outer model saturation, you can check the convergence"
return
end
Min_ȷ = Min_vec[ȷ]; @info "subblock $ȷ is selected"
prim_obj_tbMIN_ȷ = Min_ȷ[:prim_obj_tbMIN_ȷ]
fpmax_ȷ = Min_ȷ[:fpmax_ȷ]
JuMP.@objective(Min_ȷ, Min, prim_obj_tbMIN_ȷ + fˈs_pen_4_(ȷ, Π, fpmax_ȷ))
JuMP.optimize!(Min_ȷ); JuMP.assert_is_solved_and_feasible(Min_ȷ; allow_local = false)
ObjVal_ȷ = JuMP.objective_value(Min_ȷ)
Prim_obj_tbMIN_ȷ, Fpmax_ȷ = JuMP.value(prim_obj_tbMIN_ȷ), JuMP.value.(fpmax_ȷ)
if ObjVal_ȷ + COT < Θ[ȷ]
JuMP.@constraint(out, θ[ȷ] ≤ Prim_obj_tbMIN_ȷ + fˈs_pen_4_(ȷ, π, Fpmax_ȷ))
outer_saturation = false
break # Go back to outer model to update Π
end
end
end
