FIX: SOC and LinDist3Flow formulations for elements with missing phas…

…es (#428) * Updated SOC and LinDist3Flow formulations * Updated changelog and bf_mx.jl
lanl-ansi · Feb 14, 2023 · c638a9a · c638a9a
1 parent e11aa65
commit c638a9a
Show file tree

Hide file tree

Showing 5 changed files with 141 additions and 16 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,8 @@
 
 ## staged
 
+- Added loads/generator models (240V devices) connected between two secondary terminals of center-tapped transformers for SOC formulation
+- Fixed bug with SOC and LinDist3Flow formulations where diagonal entries of matrix variables were defined with type `Vector{JuMP.VariableRef}` (no information about connections) instead of `JuMP.Containers.DenseAxisArray`, leading to errors when single- or two-phase nodes were present in network
 - Fixed bug in `_calc_bus_vm_ll_bounds` where default min `vdmin_eps` was not being used, leading to invalid `Inf` bounds
 
 ## v0.14.6

diff --git a/src/form/bf_mx.jl b/src/form/bf_mx.jl
@@ -50,7 +50,7 @@ function variable_mc_bus_voltage_prod_hermitian(pm::AbstractUBFModels; nw::Int=n
     var(pm, nw)[:Wr] = Wr
     var(pm, nw)[:Wi] = Wi
     # maintain compatibility
-    var(pm, nw)[:w] = Dict{Int, Any}([(id, LinearAlgebra.diag(Wr[id])) for id in bus_ids])
+    var(pm, nw)[:w] = Dict(id => JuMP.Containers.DenseAxisArray(LinearAlgebra.diag(Wr[id]), terminals[id]) for id in bus_ids)
 
     report && _IM.sol_component_value(pm, pmd_it_sym, nw, :bus, :Wr, ids(pm, nw, :bus), Wr)
     report && _IM.sol_component_value(pm, pmd_it_sym, nw, :bus, :Wi, ids(pm, nw, :bus), Wi)
@@ -129,8 +129,8 @@ function variable_mc_branch_power(pm::AbstractUBFModels; nw::Int=nw_id_default,
     var(pm, nw)[:P] = P
     var(pm, nw)[:Q] = Q
 
-    var(pm, nw)[:p] = Dict([(id,LinearAlgebra.diag(P[id])) for id in branch_arcs])
-    var(pm, nw)[:q] = Dict([(id,LinearAlgebra.diag(Q[id])) for id in branch_arcs])
+    var(pm, nw)[:p] = Dict(id => JuMP.Containers.DenseAxisArray(LinearAlgebra.diag(P[id]), connections[id]) for id in branch_arcs)
+    var(pm, nw)[:q] = Dict(id => JuMP.Containers.DenseAxisArray(LinearAlgebra.diag(Q[id]), connections[id]) for id in branch_arcs)
 
     report && _IM.sol_component_value_edge(pm, pmd_it_sym, nw, :branch, :Pf, :Pt, ref(pm, nw, :arcs_branch_from), ref(pm, nw, :arcs_branch_to), P)
     report && _IM.sol_component_value_edge(pm, pmd_it_sym, nw, :branch, :Qf, :Qt, ref(pm, nw, :arcs_branch_from), ref(pm, nw, :arcs_branch_to), Q)
@@ -191,8 +191,8 @@ function variable_mc_switch_power(pm::AbstractUBFModels; nw::Int=nw_id_default,
     var(pm, nw)[:Psw] = P_aux
     var(pm, nw)[:Qsw] = Q_aux
 
-    var(pm, nw)[:psw] = Dict([(id,LinearAlgebra.diag(P_aux[id])) for id in switch_arcs])
-    var(pm, nw)[:qsw] = Dict([(id,LinearAlgebra.diag(Q_aux[id])) for id in switch_arcs])
+    var(pm, nw)[:psw] = Dict(id => JuMP.Containers.DenseAxisArray(LinearAlgebra.diag(P_aux[id]), connections[id]) for id in switch_arcs)
+    var(pm, nw)[:qsw] = Dict(id => JuMP.Containers.DenseAxisArray(LinearAlgebra.diag(Q_aux[id]), connections[id]) for id in switch_arcs)
 
     report && _IM.sol_component_value_edge(pm, pmd_it_sym, nw, :switch, :Pf, :Pt, ref(pm, nw, :arcs_switch_from), ref(pm, nw, :arcs_switch_to), P_expr)
     report && _IM.sol_component_value_edge(pm, pmd_it_sym, nw, :switch, :Qf, :Qt, ref(pm, nw, :arcs_switch_from), ref(pm, nw, :arcs_switch_to), Q_expr)
@@ -245,8 +245,8 @@ function variable_mc_transformer_power(pm::AbstractUBFModels; nw::Int=nw_id_defa
     var(pm, nw)[:Pt] = Pt
     var(pm, nw)[:Qt] = Qt
 
-    var(pm, nw)[:pt] = pt = Dict([(id,LinearAlgebra.diag(Pt[id])) for id in transformer_arcs])
-    var(pm, nw)[:qt] = qt = Dict([(id,LinearAlgebra.diag(Qt[id])) for id in transformer_arcs])
+    var(pm, nw)[:pt] = pt = Dict(id => JuMP.Containers.DenseAxisArray(LinearAlgebra.diag(Pt[id]), connections[id]) for id in transformer_arcs)
+    var(pm, nw)[:qt] = qt = Dict(id => JuMP.Containers.DenseAxisArray(LinearAlgebra.diag(Qt[id]), connections[id]) for id in transformer_arcs)
 
     report && _IM.sol_component_value_edge(pm, pmd_it_sym, nw, :transformer, :Pf, :Pt, ref(pm, nw, :arcs_transformer_from), ref(pm, nw, :arcs_transformer_to), Pt)
     report && _IM.sol_component_value_edge(pm, pmd_it_sym, nw, :transformer, :Qf, :Qt, ref(pm, nw, :arcs_transformer_from), ref(pm, nw, :arcs_transformer_to), Qt)
@@ -841,21 +841,27 @@ function constraint_mc_load_power(pm::AbstractUBFModels, load_id::Int; nw::Int=n
             sol(pm, nw, :load, load_id)[:qd_bus] = var(pm, nw, :qd_bus)[load_id]
         end
     elseif load["configuration"]==DELTA
+        is_triplex = length(connections)<3
+        conn_bus = is_triplex ? bus["terminals"] : connections
         # link Wy, CCd and X
-        Wr = var(pm, nw, :Wr, bus_id)[[findfirst(isequal(c), terminals) for c in connections],[findfirst(isequal(c), terminals) for c in connections]]
-        Wi = var(pm, nw, :Wi, bus_id)[[findfirst(isequal(c), terminals) for c in connections],[findfirst(isequal(c), terminals) for c in connections]]
+        Wr = var(pm, nw, :Wr, bus_id)[[findfirst(isequal(c), terminals) for c in conn_bus],[findfirst(isequal(c), terminals) for c in conn_bus]]
+        Wi = var(pm, nw, :Wi, bus_id)[[findfirst(isequal(c), terminals) for c in conn_bus],[findfirst(isequal(c), terminals) for c in conn_bus]]
         CCdr = var(pm, nw, :CCdr, load_id)
         CCdi = var(pm, nw, :CCdi, load_id)
         Xdr = var(pm, nw, :Xdr, load_id)
         Xdi = var(pm, nw, :Xdi, load_id)
-        Td = [1 -1 0; 0 1 -1; -1 0 1]  # TODO
+        Td = is_triplex ? [1 -1] : [1 -1 0; 0 1 -1; -1 0 1]  # TODO
         constraint_SWL_psd(pm.model, Xdr, Xdi, Wr, Wi, CCdr, CCdi)
         # define pd/qd and pd_bus/qd_bus as affine transformations of X
         pd_bus = LinearAlgebra.diag(Xdr*Td)
         qd_bus = LinearAlgebra.diag(Xdi*Td)
         pd = LinearAlgebra.diag(Td*Xdr)
         qd = LinearAlgebra.diag(Td*Xdi)
 
+        pd_bus = JuMP.Containers.DenseAxisArray(pd_bus, conn_bus)
+        qd_bus = JuMP.Containers.DenseAxisArray(qd_bus, conn_bus)
+        pd = JuMP.Containers.DenseAxisArray(pd, connections)
+        qd = JuMP.Containers.DenseAxisArray(qd, connections)
         var(pm, nw, :pd_bus)[load_id] = pd_bus
         var(pm, nw, :qd_bus)[load_id] = qd_bus
         var(pm, nw, :pd)[load_id] = pd

diff --git a/src/form/bf_mx_lin.jl b/src/form/bf_mx_lin.jl
@@ -21,6 +21,17 @@ function variable_mc_branch_power(pm::LPUBFDiagModel; nw::Int=nw_id_default, bou
 end
 
 
+"""
+    variable_mc_switch_power(pm::LPUBFDiagModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
+
+Switch power variables.
+"""
+function variable_mc_switch_power(pm::LPUBFDiagModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
+    variable_mc_switch_power_real(pm; nw=nw, bounded=bounded, report=report)
+    variable_mc_switch_power_imaginary(pm; nw=nw, bounded=bounded, report=report)
+end
+
+
 """
     variable_mc_capcontrol(pm::AbstractLPUBFModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true)
 

diff --git a/src/form/bf_mx_soc.jl b/src/form/bf_mx_soc.jl
@@ -1,3 +1,47 @@
+"""
+    variable_mc_generator_power(pm::SOCUBFModels; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
+
+The variable creation for generators in SOC branch flow model.
+Delta generators always need an auxilary power variable (X) and current squared variable (CC) similar to delta loads.
+Wye generators however, don't need any variables.
+"""
+function variable_mc_generator_power(pm::SOCUBFModels; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
+    variable_mc_generator_power_real(pm; nw=nw, bounded=bounded, report=report)
+    variable_mc_generator_power_imaginary(pm; nw=nw, bounded=bounded, report=report)
+    # create auxilary variables for delta generators
+    gen_del_ids = [id for (id, gen) in ref(pm, nw, :gen) if gen["configuration"]==DELTA]
+    variable_mc_generator_power_delta_aux(pm, gen_del_ids; nw=nw)
+    bounded && variable_mc_generator_current(pm, gen_del_ids; nw=nw, bounded=bounded, report=report)
+end
+
+
+"""
+    variable_mc_generator_current(pm::SOCUBFModels, gen_ids::Vector{Int}; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
+
+For the SOC branch-flow formulation, the delta-generator needs an explicit current variable.
+"""
+function variable_mc_generator_current(pm::SOCUBFModels, gen_ids::Vector{Int}; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
+    @assert(bounded)
+
+    connections = Dict{Int,Vector{Int}}(i => gen["connections"] for (i,gen) in ref(pm, nw, :gen))
+    # calculate bounds
+    bound = Dict{eltype(gen_ids), Matrix{Real}}()
+    for (id, gen) in ref(pm, nw, :gen)
+        bus = ref(pm, nw, :bus, gen["gen_bus"])
+        cmax = _calc_gen_current_max(gen, bus)
+        bound[id] = cmax*cmax'
+    end
+    # create matrix variables
+    (CCgr,CCgi) = variable_mx_hermitian(pm.model, gen_ids, connections; symm_bound=bound, name="CCg", prefix="$nw")
+    # save references
+    var(pm, nw)[:CCgr] = CCgr
+    var(pm, nw)[:CCgi] = CCgi
+
+    report && _IM.sol_component_value(pm, pmd_it_sym, nw, :gen, :CCgr, gen_ids, CCgr)
+    report && _IM.sol_component_value(pm, pmd_it_sym, nw, :gen, :CCgi, gen_ids, CCgi)
+end
+
+
 "Defines relationship between branch (series) power flow, branch (series) current and node voltage magnitude"
 function constraint_mc_model_current(pm::SOCUBFModels, nw::Int, i::Int, f_bus::Int, f_idx::Tuple{Int,Int,Int}, g_sh_fr::Matrix{<:Real}, b_sh_fr::Matrix{<:Real})
     branch = ref(pm, nw, :branch, f_idx[1])
@@ -130,17 +174,17 @@ function constraint_mc_transformer_power_yy(pm::SOCUBFModels, nw::Int, trans_id:
 
         for (t_idx,tc) in enumerate(t_connections)
             if tm_fixed[t_idx]
-                JuMP.@constraint(pm.model, Wr_fr[fc,tc] == (pol*tm_scale)^2*tm[f_idx]*tm[t_idx]*Wr_to[fc,tc])
-                JuMP.@constraint(pm.model, Wi_fr[fc,tc] == (pol*tm_scale)^2*tm[f_idx]*tm[t_idx]*Wi_to[fc,tc])
+                JuMP.@constraint(pm.model, Wr_fr[f_idx,t_idx] == (pol*tm_scale)^2*tm[f_idx]*tm[t_idx]*Wr_to[f_idx,t_idx])
+                JuMP.@constraint(pm.model, Wi_fr[f_idx,t_idx] == (pol*tm_scale)^2*tm[f_idx]*tm[t_idx]*Wi_to[f_idx,t_idx])
             else
                 PolyhedralRelaxations.construct_univariate_relaxation!(pm.model, x->x^2, tm[idx], tmsqr[idx], [JuMP.has_lower_bound(tm[idx]) ? JuMP.lower_bound(tm[idx]) : 0.9, JuMP.has_upper_bound(tm[idx]) ? JuMP.upper_bound(tm[idx]) : 1.1], false)
                 tmsqr_Wr = JuMP.@variable(pm.model, base_name="$(nw)_tmsqr_Wr_to_$(trans_id)_$(f_bus)_$(fc)_$(t_bus)_$(tc)")
                 tmsqr_Wi = JuMP.@variable(pm.model, base_name="$(nw)_tmsqr_Wi_to_$(trans_id)_$(f_bus)_$(fc)_$(t_bus)_$(tc)")
-                PolyhedralRelaxations.construct_bilinear_relaxation!(pm.model, tmsqr[idx], Wr_to[fc,tc], tmsqr_Wr, [JuMP.lower_bound(tmsqr[idx]), JuMP.upper_bound(tmsqr[idx])], [JuMP.has_lower_bound(Wr_to[fc,tc]) ? JuMP.lower_bound(Wr_to[fc,tc]) : -(1.1^2), JuMP.has_upper_bound(Wr_to[fc,tc]) ? JuMP.upper_bound(Wr_to[fc,tc]) : 1.1^2])
-                PolyhedralRelaxations.construct_bilinear_relaxation!(pm.model, tmsqr[idx], Wi_to[fc,tc], tmsqr_Wi, [JuMP.lower_bound(tmsqr[idx]), JuMP.upper_bound(tmsqr[idx])], [JuMP.has_lower_bound(Wi_to[fc,tc]) ? JuMP.lower_bound(Wi_to[fc,tc]) : -(1.1^2), JuMP.has_upper_bound(Wi_to[fc,tc]) ? JuMP.upper_bound(Wi_to[fc,tc]) : 1.1^2])
+                PolyhedralRelaxations.construct_bilinear_relaxation!(pm.model, tmsqr[idx], Wr_to[f_idx,t_idx], tmsqr_Wr, [JuMP.lower_bound(tmsqr[idx]), JuMP.upper_bound(tmsqr[idx])], [JuMP.has_lower_bound(Wr_to[f_idx,t_idx]) ? JuMP.lower_bound(Wr_to[f_idx,t_idx]) : -(1.1^2), JuMP.has_upper_bound(Wr_to[f_idx,t_idx]) ? JuMP.upper_bound(Wr_to[f_idx,t_idx]) : 1.1^2])
+                PolyhedralRelaxations.construct_bilinear_relaxation!(pm.model, tmsqr[idx], Wi_to[f_idx,t_idx], tmsqr_Wi, [JuMP.lower_bound(tmsqr[idx]), JuMP.upper_bound(tmsqr[idx])], [JuMP.has_lower_bound(Wi_to[f_idx,t_idx]) ? JuMP.lower_bound(Wi_to[f_idx,t_idx]) : -(1.1^2), JuMP.has_upper_bound(Wi_to[f_idx,t_idx]) ? JuMP.upper_bound(Wi_to[f_idx,t_idx]) : 1.1^2])
 
-                JuMP.@constraint(pm.model, Wr_fr[fc,tc] == (pol*tm_scale)^2*tmsqr_Wr_to)
-                JuMP.@constraint(pm.model, Wi_fr[fc,tc] == (pol*tm_scale)^2*tmsqr_Wi_to)
+                JuMP.@constraint(pm.model, Wr_fr[f_idx,t_idx] == (pol*tm_scale)^2*tmsqr_Wr_to)
+                JuMP.@constraint(pm.model, Wi_fr[f_idx,t_idx] == (pol*tm_scale)^2*tmsqr_Wi_to)
 
                 # with regcontrol
                 if haskey(transformer,"controls")
@@ -217,3 +261,54 @@ function constraint_mc_transformer_power_dy(pm::SOCUBFModels, nw::Int, trans_id:
     JuMP.@constraint(pm.model, LinearAlgebra.diag(Tt*Xtr) + p_to .== 0)
     JuMP.@constraint(pm.model, LinearAlgebra.diag(Tt*Xti) + q_to .== 0)
 end
+
+
+@doc raw"""
+    constraint_mc_generator_power_delta(pm::SOCUBFModels, nw::Int, gen_id::Int, bus_id::Int, connections::Vector{Int}, pmin::Vector{<:Real}, pmax::Vector{<:Real}, qmin::Vector{<:Real}, qmax::Vector{<:Real}; report::Bool=true, bounded::Bool=true)
+
+Adds constraints for delta-connected generators similar to delta-connected loads (using auxilary variable X).
+
+```math
+\begin{align}
+&\text{Three-phase delta transformation matrix: }  T^\Delta = \begin{bmatrix}\;\;\;1 & -1 & \;\;0\\ \;\;\;0 & \;\;\;1 & -1\\ -1 & \;\;\;0 & \;\;\;1\end{bmatrix} \\
+&\text{Single-phase delta transformation matrix (triple nodes): }  T^\Delta = \begin{bmatrix}\;1 & -1 \end{bmatrix} \\
+&\text{Line-neutral generation power: }  S_{bus} = diag(T^\Delta X_g) \\
+&\text{Line-line generation power: }  S^\Delta = diag(X_g T^\Delta)
+\end{align}
+```
+"""
+function constraint_mc_generator_power_delta(pm::SOCUBFModels, nw::Int, gen_id::Int, bus_id::Int, connections::Vector{Int}, pmin::Vector{<:Real}, pmax::Vector{<:Real}, qmin::Vector{<:Real}, qmax::Vector{<:Real}; report::Bool=true, bounded::Bool=true)
+    gen = ref(pm, nw, :gen, gen_id)
+    bus_id = gen["gen_bus"]
+    bus = ref(pm, nw, :bus, bus_id)
+    terminals = bus["terminals"]
+    is_triplex = length(connections)<3
+    conn_bus = is_triplex ? bus["terminals"] : connections
+
+    pg = var(pm, nw, :pg, gen_id)
+    qg = var(pm, nw, :qg, gen_id)
+    Xgr = var(pm, nw, :Xgr, gen_id)
+    Xgi = var(pm, nw, :Xgi, gen_id)
+    CCgr = var(pm, nw, :CCgr, gen_id)
+    CCgi = var(pm, nw, :CCgi, gen_id)
+    Wr = var(pm, nw, :Wr, bus_id)[[findfirst(isequal(c), terminals) for c in conn_bus],[findfirst(isequal(c), terminals) for c in conn_bus]]
+    Wi = var(pm, nw, :Wi, bus_id)[[findfirst(isequal(c), terminals) for c in conn_bus],[findfirst(isequal(c), terminals) for c in conn_bus]]
+
+    Tg = is_triplex ? [1 -1] : [1 -1 0; 0 1 -1; -1 0 1]  # TODO
+    constraint_SWL_psd(pm.model, Xgr, Xgi, Wr, Wi, CCgr, CCgi)
+    # define pg/qg and pg_bus/qg_bus as affine transformations of X
+    JuMP.@constraint(pm.model, pg .== LinearAlgebra.diag(Tg*Xgr))
+    JuMP.@constraint(pm.model, qg .== LinearAlgebra.diag(Tg*Xgi))
+
+    pg_bus = LinearAlgebra.diag(Xgr*Tg)
+    qg_bus = LinearAlgebra.diag(Xgi*Tg)
+    pg_bus = JuMP.Containers.DenseAxisArray(pg_bus, conn_bus)
+    qg_bus = JuMP.Containers.DenseAxisArray(qg_bus, conn_bus)
+    var(pm, nw, :pg_bus)[gen_id] = pg_bus
+    var(pm, nw, :qg_bus)[gen_id] = qg_bus
+
+    if report
+        sol(pm, nw, :gen, gen_id)[:pg_bus] = pg_bus
+        sol(pm, nw, :gen, gen_id)[:qg_bus] = qg_bus
+    end
+end
diff --git a/test/transformer.jl b/test/transformer.jl
@@ -362,5 +362,16 @@
             result = solve_mc_opf(ut_trans_3w_dyy_1, SOCConicUBFPowerModel, scs_solver; solution_processors=[sol_data_model!], make_si=false)
             @test norm(result["solution"]["bus"]["3"]["vm"]-[0.93180, 0.88827, 0.88581], Inf) <= 7.2E-2
         end
+
+        @testset "3w_center_tap" begin
+            apply_voltage_bounds!(trans_3w_center_tap; vm_lb=0.95, vm_ub=1.05)
+            result = solve_mc_opf(trans_3w_center_tap, SOCConicUBFPowerModel, scs_solver; solution_processors=[sol_data_model!], make_si=false)
+
+            sbase = trans_3w_center_tap["settings"]["sbase_default"]
+
+            @test all(isapprox.(result["solution"]["load"]["l2a"]["pd"]*sbase, 12.0; atol=1E-5))
+            @test all(isapprox.(result["solution"]["load"]["l3"]["pd"]*sbase, [10.0, 10.0]; atol=1E-5))
+            @test all(isapprox.(result["solution"]["bus"]["tn_1"]["vm"], [1.045, 1.05]; atol=5E-3))
+        end
     end
 end