FIX: Bugs in center tap transformer model

CHANGELOG.md

@@ -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

docs/src/manual/load-model.md

@@ -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.$$
 

src/core/data.jl

@@ -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)

src/data_model/eng2math.jl

@@ -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
@@ -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],
@@ -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"])")
 
@@ -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"])")
@@ -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}(
@@ -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

src/data_model/transformations.jl

@@ -702,9 +702,13 @@ 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
@@ -722,9 +726,13 @@ function _apply_phase_projection_delta!(data_eng::Dict{String,<:Any})
     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

src/data_model/utils.jl

@@ -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}
@@ -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,
@@ -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)),

src/form/acp.jl

@@ -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()
@@ -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
@@ -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()
@@ -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

src/form/acr.jl

@@ -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()
@@ -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
@@ -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()
@@ -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