Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

FIX: Bugs in center tap transformer model #407

Merged
merged 11 commits into from
Dec 20, 2022
Merged

FIX: Bugs in center tap transformer model #407

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
b803312
Modified terminals of virtual transformer objects
keatsig Oct 26, 2022
9ed98c3
Corrected center-tap transformer polarity
keatsig Nov 3, 2022
cc79864
Merge branch 'lanl-ansi:main' into virtual-transformer-objects
keatsig Nov 3, 2022
aa33f84
Update utils.jl
keatsig Nov 3, 2022
71540e4
Merge branch 'lanl-ansi:main' into virtual-transformer-objects
keatsig Nov 14, 2022
a829357
Updated initial point for triplex connections
keatsig Nov 15, 2022
59c5643
Updated center tap secondary connections model for ACP formulation
keatsig Nov 28, 2022
683644d
Updated procedure to identify triplex nodes
keatsig Dec 8, 2022
9f4669e
Updated load/generator models connected to triplex nodes for ACR and …
keatsig Dec 13, 2022
26fc80d
Updated load/generator models connected to triplex nodes for linear f…
keatsig Dec 13, 2022
1681b12
Updated documentation
keatsig Dec 15, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

## staged

- Fixed bugs in center-tap transformer modeling
- Add wye-connected CapControl for IVR and FOT (polar) formulations
- Fixed indexing issue for single-phase delta load models in linear formulations (LinDist3Flow, FOTP, FOTR, FBS)
- Added ZIP load model
Expand Down
6 changes: 5 additions & 1 deletion docs/src/manual/load-model.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,11 @@ Firstly, define the three-phase delta transformation matrix

$$M^\Delta_3 = \begin{bmatrix}\;\;\;1 & -1 & \;\;0\\ \;\;\;0 & \;\;\;1 & -1\\ -1 & \;\;\;0 & \;\;\;1\end{bmatrix},$$

which can be extended to more phases in a straight-forward manner. Now,
which can be extended to more phases in a straight-forward manner. For loads connected between split-phase terminals of triplex nodes (usually located on the secondary side of center-tapped transformers), we define a single-phase delta transformation matrix

$$M^\Delta_1 = \begin{bmatrix} 1 & -1 \end{bmatrix}.$$

Now,

$$U^d = M^\Delta U^\text{bus},\;\;\; I^\text{bus} = \left(M^\Delta\right)^T I^d.$$

Expand Down
6 changes: 4 additions & 2 deletions src/core/data.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -448,16 +448,18 @@ The returned bounds are for the pairs 1->2, 2->3, 3->1
function _calc_bus_vm_ll_bounds(bus::Dict; vdmin_eps::Real=0.1)::Tuple
vmax = bus["vmax"]
vmin = bus["vmin"]
is_triplex = length(bus["terminals"]) < 3
if haskey(bus, "vm_ll_max")
vdmax = bus["vm_ll_max"]
else
# implied valid upper bound
vdmax = _mat_mult_rm_nan([1 1 0; 0 1 1; 1 0 1], vmax)
Td = is_triplex ? [1 1] : [1 1 0; 0 1 1; 1 0 1]
vdmax = _mat_mult_rm_nan(Td, vmax)
end
if haskey(bus, "vm_ll_min")
vdmin = bus["vm_ll_min"]
else
vdmin = fill(0.0, length(vmin))
vdmin = is_triplex ? [0.0] : fill(0.0, length(vmin))
end

return (vdmin, vdmax)
Expand Down
29 changes: 23 additions & 6 deletions src/data_model/eng2math.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -481,7 +481,7 @@ function _map_eng2math_transformer!(data_math::Dict{String,<:Any}, data_eng::Dic
z_sc = Dict([(key, im*x_sc[i]) for (i,key) in enumerate([(i,j) for i in 1:nrw for j in i+1:nrw])])

dims = length(eng_obj["tm_set"][1])
transformer_t_bus_w = _build_loss_model!(data_math, name, to_map, r_s, z_sc, y_sh; nphases=dims, status=Int(eng_obj["status"] == ENABLED))
transformer_t_bus_w = _build_loss_model!(data_math, name, to_map, r_s, z_sc, y_sh,eng_obj["connections"][1]; nphases=dims, status=Int(eng_obj["status"] == ENABLED))

for w in 1:nrw
# 2-WINDING TRANSFORMER
Expand All @@ -494,7 +494,7 @@ function _map_eng2math_transformer!(data_math::Dict{String,<:Any}, data_eng::Dic
"t_bus" => transformer_t_bus_w[w],
"tm_nom" => tm_nom,
"f_connections" => eng_obj["connections"][w],
"t_connections" => get(data_math, "is_kron_reduced", false) ? collect(1:dims) : collect(1:dims+1),
"t_connections" => get(data_math, "is_kron_reduced", false) ? eng_obj["connections"][1] : collect(1:dims+1),
"configuration" => eng_obj["configuration"][w],
"polarity" => eng_obj["polarity"][w],
"tm_set" => eng_obj["tm_set"][w],
Expand Down Expand Up @@ -525,6 +525,23 @@ function _map_eng2math_transformer!(data_math::Dict{String,<:Any}, data_eng::Dic
)
data_math["transformer"]["$(transformer_2wa_obj["index"])"]["controls"] = reg_obj
end
if w==3 && eng_obj["polarity"][w]==-1 # identify center-tapped transformer and mark all secondary-side nodes as triplex by adding va_start
default_va = [0, -120, 120][eng_obj["connections"][1][1]]
data_math["bus"]["$(transformer_2wa_obj["f_bus"])"]["va_start"] = haskey(data_eng["bus"][eng_obj["bus"][w]],"va_start") ? data_eng["bus"][eng_obj["bus"][w]]["va_start"] : [default_va, (default_va+180)]
idx = 0
bus_ids = []
t_bus = haskey(data_eng, "line") ? [data["t_bus"] for (_,data) in data_eng["line"] if data["f_bus"] == eng_obj["bus"][w]] : []
while length(t_bus)>0 || idx<length(bus_ids)
for bus_idx in t_bus
bus_id = data_math["bus_lookup"]["$bus_idx"]
push!(bus_ids, bus_id)
default_va = [0, -120, 120][eng_obj["connections"][1][1]]
data_math["bus"]["$bus_id"]["va_start"] = haskey(data_eng["bus"]["$bus_idx"],"va_start") ? data_eng["bus"]["$bus_idx"]["va_start"] : [default_va, (default_va+180)]
end
idx += 1
t_bus = [data["t_bus"] for (_,data) in data_eng["line"] if data["f_bus"] == data_math["bus"]["$(bus_ids[idx])"]["name"]]
end
end

push!(to_map, "transformer.$(transformer_2wa_obj["index"])")

Expand Down Expand Up @@ -699,8 +716,8 @@ function _map_eng2math_load!(data_math::Dict{String,<:Any}, data_eng::Dict{Strin
"dispatchable" => eng_obj["dispatchable"] == NO ? 0 : 1,
"index" => length(data_math["load"])+1,
"pd" => eng_obj["pd_nom"]*eng_obj["zipv"][idx],
)
)

data_math["load"]["$(math_obj["index"])"] = math_obj

push!(to_map, "load.$(math_obj["index"])")
Expand All @@ -721,7 +738,7 @@ function _map_eng2math_load!(data_math::Dict{String,<:Any}, data_eng::Dict{Strin
math_obj["qd"] = eng_obj["qd_nom"]

math_obj["vnom_kv"] = eng_obj["vm_nom"]

data_math["load"]["$(math_obj["index"])"] = math_obj

push!(data_math["map"], Dict{String,Any}(
Expand Down Expand Up @@ -810,7 +827,7 @@ function _map_eng2math_solar!(data_math::Dict{String,<:Any}, data_eng::Dict{Stri
end
end

N = _infer_int_dim_unit(eng_obj, false)
N = eng_obj["configuration"]==DELTA && length(eng_obj["connections"])==1 ? 1 : _infer_int_dim_unit(eng_obj, false) # if solar is delta-connected to triplex node, N can be equal to 1
for (fr_k, to_k, def) in [("pg_lb", "pmin", -Inf), ("pg_ub", "pmax", Inf), ("qg_lb", "qmin", -Inf), ("qg_ub", "qmax", Inf)]
math_obj[to_k] = haskey(eng_obj, fr_k) ? eng_obj[fr_k] : fill(def, N)
end
Expand Down
30 changes: 21 additions & 9 deletions src/data_model/transformations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -702,29 +702,41 @@ function _apply_phase_projection_delta!(data_eng::Dict{String,<:Any})
if haskey(data_eng, "load")
for (_,eng_obj) in data_eng["load"]
if eng_obj["configuration"] == DELTA
_pad_properties_delta!(eng_obj, ["pd_nom", "qd_nom"], eng_obj["connections"], all_conductors)
_pad_connections!(eng_obj, "connections", all_conductors)
bus_terminals[eng_obj["bus"]] = haskey(bus_terminals, eng_obj["bus"]) ? _pad_connections!(bus_terminals, eng_obj["bus"], eng_obj["connections"]) : eng_obj["connections"]
if eng_obj["connections"]==[1, 2] && eng_obj["vm_nom"]==0.24 # check if load is connected between split-phase terminals of triplex node (nominal line-line voltage=240V), TODO: better generalization
bus_terminals[eng_obj["bus"]] = eng_obj["connections"] = [1]
else
_pad_properties_delta!(eng_obj, ["pd_nom", "qd_nom"], eng_obj["connections"], all_conductors)
_pad_connections!(eng_obj, "connections", all_conductors)
bus_terminals[eng_obj["bus"]] = haskey(bus_terminals, eng_obj["bus"]) ? _pad_connections!(bus_terminals, eng_obj["bus"], eng_obj["connections"]) : eng_obj["connections"]
end
end
end
end

if haskey(data_eng, "generator")
for (_,eng_obj) in data_eng["generator"]
if eng_obj["configuration"]==DELTA
_pad_properties_delta!(eng_obj, ["pg", "qg", "vg", "pg_lb", "pg_ub", "qg_lb", "qg_ub"], eng_obj["connections"], all_conductors)
_pad_connections!(eng_obj, "connections", all_conductors)
bus_terminals[eng_obj["bus"]] = haskey(bus_terminals, eng_obj["bus"]) ? _pad_connections!(bus_terminals, eng_obj["bus"], eng_obj["connections"]) : eng_obj["connections"]
if eng_obj["connections"]==[1, 2] && eng_obj["vg"][1]==0.24 # check if generator is connected between split-phase terminals of triplex node (nominal line-line voltage=240V), TODO: better generalization
bus_terminals[eng_obj["bus"]] = eng_obj["connections"] = [1]
else
_pad_properties_delta!(eng_obj, ["pg", "qg", "vg", "pg_lb", "pg_ub", "qg_lb", "qg_ub"], eng_obj["connections"], all_conductors)
_pad_connections!(eng_obj, "connections", all_conductors)
bus_terminals[eng_obj["bus"]] = haskey(bus_terminals, eng_obj["bus"]) ? _pad_connections!(bus_terminals, eng_obj["bus"], eng_obj["connections"]) : eng_obj["connections"]
end
end
end
end

if haskey(data_eng, "solar")
for (_,eng_obj) in data_eng["solar"]
if eng_obj["configuration"]==DELTA
_pad_properties_delta!(eng_obj, ["pg", "qg", "vg", "pg_lb", "pg_ub", "qg_lb", "qg_ub"], eng_obj["connections"], all_conductors)
_pad_connections!(eng_obj, "connections", all_conductors)
bus_terminals[eng_obj["bus"]] = haskey(bus_terminals, eng_obj["bus"]) ? _pad_connections!(bus_terminals, eng_obj["bus"], eng_obj["connections"]) : eng_obj["connections"]
if eng_obj["connections"]==[1, 2] && eng_obj["vg"][1]==0.24 # check if solar is connected between split-phase terminals of triplex node (nominal line-line voltage=240V), TODO: better generalization
bus_terminals[eng_obj["bus"]] = eng_obj["connections"] = [1]
else
_pad_properties_delta!(eng_obj, ["pg", "qg", "vg", "pg_lb", "pg_ub", "qg_lb", "qg_ub"], eng_obj["connections"], all_conductors)
_pad_connections!(eng_obj, "connections", all_conductors)
bus_terminals[eng_obj["bus"]] = haskey(bus_terminals, eng_obj["bus"]) ? _pad_connections!(bus_terminals, eng_obj["bus"], eng_obj["connections"]) : eng_obj["connections"]
end
end
end
end
Expand Down
9 changes: 5 additions & 4 deletions src/data_model/utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -260,7 +260,8 @@ function _build_loss_model!(
to_map::Vector{String},
r_s::Vector{Float64},
zsc::Dict{Tuple{Int,Int},Complex{Float64}},
ysh::Complex{Float64};
ysh::Complex{Float64},
connections::Vector{Int};
nphases::Int=3,
status::Int=1,
)::Vector{Int}
Expand Down Expand Up @@ -342,7 +343,7 @@ function _build_loss_model!(
"vmax" => fill(Inf, nphases),
"vm_pair_lb" => Tuple{Any,Any,Real}[],
"vm_pair_ub" => Tuple{Any,Any,Real}[],
"terminals" => collect(1:nphases),
"terminals" => connections[collect(1:nphases)],
"grounded" => fill(false, nphases),
"base_kv" => 1.0,
"bus_type" => status == 0 ? 4 : 1,
Expand Down Expand Up @@ -395,8 +396,8 @@ function _build_loss_model!(
"br_status"=>status,
"f_bus"=>bus_ids[i],
"t_bus"=>bus_ids[j],
"f_connections"=>collect(1:nphases),
"t_connections"=>collect(1:nphases),
"f_connections"=>data_math["bus"]["$(bus_ids[i])"]["terminals"][collect(1:nphases)],
"t_connections"=>data_math["bus"]["$(bus_ids[j])"]["terminals"][collect(1:nphases)],
"br_r" => LinearAlgebra.diagm(0=>fill(real(z[l]), nphases)),
"br_x" => LinearAlgebra.diagm(0=>fill(imag(z[l]), nphases)),
"g_fr" => LinearAlgebra.diagm(0=>fill(g_fr, nphases)),
Expand Down
47 changes: 31 additions & 16 deletions src/form/acp.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1130,9 +1130,11 @@ function constraint_mc_load_power_delta(pm::AbstractUnbalancedACPModel, nw::Int,
va = var(pm, nw, :va, bus_id)

nph = length(a)
is_triplex = nph < 3
conn_bus = is_triplex ? ref(pm, nw, :bus, bus_id)["terminals"] : connections

prev = Dict(c=>connections[(idx+nph-2)%nph+1] for (idx,c) in enumerate(connections))
next = Dict(c=>connections[idx%nph+1] for (idx,c) in enumerate(connections))
prev = Dict(c=>conn_bus[(idx+nph-2)%nph+1] for (idx,c) in enumerate(conn_bus))
next = is_triplex ? conn_bus[2] : Dict(c=>conn_bus[idx%nph+1] for (idx,c) in enumerate(conn_bus))

vrd = Dict()
vid = Dict()
Expand All @@ -1156,20 +1158,25 @@ function constraint_mc_load_power_delta(pm::AbstractUnbalancedACPModel, nw::Int,

crd_bus = Dict()
cid_bus = Dict()
for (idx, c) in enumerate(connections)
crd_bus[c] = JuMP.@NLexpression(pm.model, crd[c]-crd[prev[c]])
cid_bus[c] = JuMP.@NLexpression(pm.model, cid[c]-cid[prev[c]])
for (idx, c) in enumerate(conn_bus)
if is_triplex
crd_bus[c] = JuMP.@NLexpression(pm.model, (-1.0)^(c-1)*crd[1])
cid_bus[c] = JuMP.@NLexpression(pm.model, (-1.0)^(c-1)*cid[1])
else
crd_bus[c] = JuMP.@NLexpression(pm.model, crd[c]-crd[prev[c]])
cid_bus[c] = JuMP.@NLexpression(pm.model, cid[c]-cid[prev[c]])
end
end

pd_bus = Vector{JuMP.NonlinearExpression}([])
qd_bus = Vector{JuMP.NonlinearExpression}([])
for (idx,c) in enumerate(connections)
for (idx,c) in enumerate(conn_bus)
push!(pd_bus, JuMP.@NLexpression(pm.model, vm[c]*cos(va[c])*crd_bus[c]+vm[c]*sin(va[c])*cid_bus[c]))
push!(qd_bus, JuMP.@NLexpression(pm.model, -vm[c]*cos(va[c])*cid_bus[c]+vm[c]*sin(va[c])*crd_bus[c]))
end

pd_bus = JuMP.Containers.DenseAxisArray(pd_bus, connections)
qd_bus = JuMP.Containers.DenseAxisArray(qd_bus, connections)
pd_bus = JuMP.Containers.DenseAxisArray(pd_bus, conn_bus)
qd_bus = JuMP.Containers.DenseAxisArray(qd_bus, conn_bus)

var(pm, nw, :pd_bus)[id] = pd_bus
var(pm, nw, :qd_bus)[id] = qd_bus
Expand Down Expand Up @@ -1198,8 +1205,11 @@ function constraint_mc_generator_power_delta(pm::AbstractUnbalancedACPModel, nw:
qg = var(pm, nw, :qg, id)

nph = length(connections)
prev = Dict(c=>connections[(idx+nph-2)%nph+1] for (idx,c) in enumerate(connections))
next = Dict(c=>connections[idx%nph+1] for (idx,c) in enumerate(connections))
is_triplex = nph < 3
conn_bus = is_triplex ? ref(pm, nw, :bus, bus_id)["terminals"] : connections

prev = Dict(c=>conn_bus[(idx+nph-2)%nph+1] for (idx,c) in enumerate(conn_bus))
next = is_triplex ? conn_bus[2] : Dict(c=>conn_bus[idx%nph+1] for (idx,c) in enumerate(conn_bus))

vrg = Dict()
vig = Dict()
Expand All @@ -1217,19 +1227,24 @@ function constraint_mc_generator_power_delta(pm::AbstractUnbalancedACPModel, nw:

crg_bus = Dict()
cig_bus = Dict()
for c in connections
crg_bus[c] = JuMP.@NLexpression(pm.model, crg[c]-crg[prev[c]])
cig_bus[c] = JuMP.@NLexpression(pm.model, cig[c]-cig[prev[c]])
for c in conn_bus
if is_triplex
crg_bus[c] = JuMP.@NLexpression(pm.model, (-1.0)^(c-1)*crg[1])
cig_bus[c] = JuMP.@NLexpression(pm.model, (-1.0)^(c-1)*cig[1])
else
crg_bus[c] = JuMP.@NLexpression(pm.model, crg[c]-crg[prev[c]])
cig_bus[c] = JuMP.@NLexpression(pm.model, cig[c]-cig[prev[c]])
end
end

pg_bus = Vector{JuMP.NonlinearExpression}([])
qg_bus = Vector{JuMP.NonlinearExpression}([])
for c in connections
for c in conn_bus
push!(pg_bus, JuMP.@NLexpression(pm.model, vm[c]*cos(va[c])*crg_bus[c]+vm[c]*sin(va[c])*cig_bus[c]))
push!(qg_bus, JuMP.@NLexpression(pm.model, -vm[c]*cos(va[c])*cig_bus[c]+vm[c]*sin(va[c])*crg_bus[c]))
end
pd_bus = JuMP.Containers.DenseAxisArray(pg_bus, connections)
qd_bus = JuMP.Containers.DenseAxisArray(qg_bus, connections)
pg_bus = JuMP.Containers.DenseAxisArray(pg_bus, conn_bus)
qg_bus = JuMP.Containers.DenseAxisArray(qg_bus, conn_bus)

var(pm, nw, :pg_bus)[id] = pg_bus
var(pm, nw, :qg_bus)[id] = qg_bus
Expand Down
45 changes: 29 additions & 16 deletions src/form/acr.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -821,10 +821,11 @@ function constraint_mc_load_power_delta(pm::AbstractUnbalancedACRModel, nw::Int,
vi = var(pm, nw, :vi, bus_id)

nph = length(a)
@assert nph == 3 "only phases == 3 delta loads are currently supported"
is_triplex = nph < 3
conn_bus = is_triplex ? ref(pm, nw, :bus, bus_id)["terminals"] : connections

prev = Dict(c=>connections[(idx+nph-2)%nph+1] for (idx,c) in enumerate(connections))
next = Dict(c=>connections[idx%nph+1] for (idx,c) in enumerate(connections))
next = is_triplex ? conn_bus[2] : Dict(c=>conn_bus[idx%nph+1] for (idx,c) in enumerate(conn_bus))

vrd = Dict()
vid = Dict()
Expand All @@ -842,20 +843,25 @@ function constraint_mc_load_power_delta(pm::AbstractUnbalancedACRModel, nw::Int,

crd_bus = Dict()
cid_bus = Dict()
for (idx, c) in enumerate(connections)
crd_bus[c] = JuMP.@NLexpression(pm.model, crd[c]-crd[prev[c]])
cid_bus[c] = JuMP.@NLexpression(pm.model, cid[c]-cid[prev[c]])
for (idx, c) in enumerate(conn_bus)
if is_triplex
crd_bus[c] = JuMP.@NLexpression(pm.model, (-1.0)^(c-1)*crd[1])
cid_bus[c] = JuMP.@NLexpression(pm.model, (-1.0)^(c-1)*cid[1])
else
crd_bus[c] = JuMP.@NLexpression(pm.model, crd[c]-crd[prev[c]])
cid_bus[c] = JuMP.@NLexpression(pm.model, cid[c]-cid[prev[c]])
end
end

pd_bus = Vector{JuMP.NonlinearExpression}([])
qd_bus = Vector{JuMP.NonlinearExpression}([])
for (idx,c) in enumerate(connections)
for (idx,c) in enumerate(conn_bus)
push!(pd_bus, JuMP.@NLexpression(pm.model, vr[c]*crd_bus[c]+vi[c]*cid_bus[c]))
push!(qd_bus, JuMP.@NLexpression(pm.model, -vr[c]*cid_bus[c]+vi[c]*crd_bus[c]))
end

pd_bus = JuMP.Containers.DenseAxisArray(pd_bus, connections)
qd_bus = JuMP.Containers.DenseAxisArray(qd_bus, connections)
pd_bus = JuMP.Containers.DenseAxisArray(pd_bus, conn_bus)
qd_bus = JuMP.Containers.DenseAxisArray(qd_bus, conn_bus)

var(pm, nw, :pd_bus)[id] = pd_bus
var(pm, nw, :qd_bus)[id] = qd_bus
Expand Down Expand Up @@ -893,10 +899,11 @@ function constraint_mc_generator_power_delta(pm::AbstractUnbalancedACRModel, nw:
qg = var(pm, nw, :qg, id)

nph = length(pmin)
@assert nph == 3 "only phases == 3 delta generators are currently supported"
is_triplex = nph < 3
conn_bus = is_triplex ? ref(pm, nw, :bus, bus_id)["terminals"] : connections

prev = Dict(c=>connections[(idx+nph-2)%nph+1] for (idx,c) in enumerate(connections))
next = Dict(c=>connections[idx%nph+1] for (idx,c) in enumerate(connections))
next = is_triplex ? conn_bus[2] : Dict(c=>conn_bus[idx%nph+1] for (idx,c) in enumerate(conn_bus))

vrg = Dict()
vig = Dict()
Expand All @@ -914,19 +921,25 @@ function constraint_mc_generator_power_delta(pm::AbstractUnbalancedACRModel, nw:

crg_bus = Dict()
cig_bus = Dict()
for c in connections
crg_bus[c] = JuMP.@NLexpression(pm.model, crg[c]-crg[prev[c]])
cig_bus[c] = JuMP.@NLexpression(pm.model, cig[c]-cig[prev[c]])
for c in conn_bus
if is_triplex
crg_bus[c] = JuMP.@NLexpression(pm.model, (-1.0)^(c-1)*crg[1])
cig_bus[c] = JuMP.@NLexpression(pm.model, (-1.0)^(c-1)*cig[1])
else
crg_bus[c] = JuMP.@NLexpression(pm.model, crg[c]-crg[prev[c]])
cig_bus[c] = JuMP.@NLexpression(pm.model, cig[c]-cig[prev[c]])
end
end

pg_bus = Vector{JuMP.NonlinearExpression}([])
qg_bus = Vector{JuMP.NonlinearExpression}([])
for (idx,c) in enumerate(connections)
for (idx,c) in enumerate(conn_bus)
push!(pg_bus, JuMP.@NLexpression(pm.model, vr[c]*crg_bus[c]+vi[c]*cig_bus[c]))
push!(qg_bus, JuMP.@NLexpression(pm.model, -vr[c]*cig_bus[c]+vi[c]*crg_bus[c]))
end
pd_bus = JuMP.Containers.DenseAxisArray(pg_bus, connections)
qd_bus = JuMP.Containers.DenseAxisArray(qg_bus, connections)

pg_bus = JuMP.Containers.DenseAxisArray(pg_bus, conn_bus)
qg_bus = JuMP.Containers.DenseAxisArray(qg_bus, conn_bus)

var(pm, nw, :pg_bus)[id] = pg_bus
var(pm, nw, :qg_bus)[id] = qg_bus
Expand Down
Loading