ADD: Transformer SOC relaxations

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## staged
 
+- Added SOC transformer relaxations
 - 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)

src/core/solution.jl

@@ -67,6 +67,26 @@ function _sol_data_model_w!(solution::Dict{String,<:Any})
                     end
                     delete!(bus, "w")
                 end
+                if haskey(bus, "Wr")
+                    w = LinearAlgebra.diag(bus["Wr"])
+                    if any(w .< 0) # e.g., as allowed by constraint violation settings
+                        bus["vm"] = zeros(length(w))
+                        bus["vm"][w .>= 0.0] .= sqrt.(w[w .>= 0.0])
+                    else
+                        bus["vm"] = sqrt.(w)
+                        if length(w) == 3
+                            t = [-1 1 0; -1 0 1; 0 -1 1]
+                            va = LinearAlgebra.pinv(t)*[atan(bus["Wi"][2,1], bus["Wr"][2,1]);
+                                                        atan(bus["Wi"][3,1], bus["Wr"][3,1]);
+                                                        atan(bus["Wi"][3,2], bus["Wr"][3,2])]
+                            bus["va"] = [va[findmin(abs.(va .- 0))[2]],
+                                         va[findmin(abs.(va .+ 2*pi/3))[2]],
+                                         va[findmin(abs.(va .- 2*pi/3))[2]]] # TODO: better way to get angles in order
+                        end
+                    end
+                    delete!(bus, "Wr")
+                    delete!(bus, "Wi")
+                end
             end
         end
     end

src/form/bf_mx.jl

@@ -141,9 +141,12 @@ end
 function variable_mc_transformer_power(pm::AbstractUBFModels; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
     transformer_arcs = Vector{Tuple{Int,Int,Int}}(ref(pm, nw, :arcs_transformer))
     connections = Dict{Tuple{Int,Int,Int},Vector{Int}}((l,i,j) => connections for (bus,entry) in ref(pm, nw, :bus_arcs_conns_transformer) for ((l,i,j), connections) in entry)
+    tf_del_ids = [id for (id, transformer) in ref(pm, nw, :transformer) if transformer["configuration"]==DELTA] # ids for delta configurations
 
     if bounded
         bound = Dict{eltype(transformer_arcs), Matrix{Real}}()
+        bound_del = Dict{eltype(tf_del_ids), Matrix{Real}}()
+        cmax_del = Dict{eltype(tf_del_ids), Vector{Real}}()
         for (tr, transformer) in ref(pm, nw, :transformer)
             bus_fr = ref(pm, nw, :bus, transformer["f_bus"])
             bus_to = ref(pm, nw, :bus, transformer["t_bus"])
@@ -162,6 +165,11 @@ function variable_mc_transformer_power(pm::AbstractUBFModels; nw::Int=nw_id_defa
                 bound[tuple_fr][idx,idx] = smax_fr[idx]
                 bound[tuple_to][idx,idx] = smax_to[idx]
             end
+
+            if transformer["configuration"]==DELTA
+                bound_del[tr] = bound[tuple_fr]
+                cmax_del[tr] = cmax_fr
+            end
         end
         # create matrix variables
         (Pt,Qt) = variable_mx_complex(pm.model, transformer_arcs, connections, connections; symm_bound=bound, name=("Pt", "Qt"), prefix="$nw")
@@ -180,6 +188,24 @@ function variable_mc_transformer_power(pm::AbstractUBFModels; nw::Int=nw_id_defa
     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)
     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)
+
+    # create auxilary matrix variables for transformers with delta configuration
+    if !isempty(tf_del_ids)
+        conn_del = Dict{Int,Vector{Int}}(id => transformer["f_connections"] for (id,transformer) in ref(pm, nw, :transformer))
+        (Xtr,Xti) = variable_mx_complex(pm.model, tf_del_ids, conn_del, conn_del; symm_bound=bound_del, name="Xt", prefix="$nw")
+        (CCtr, CCti) = variable_mx_hermitian(pm.model, tf_del_ids, conn_del; sqrt_upper_bound=cmax_del, name="CCt", prefix="$nw")
+        # save references
+        var(pm, nw)[:Xtr] = Xtr
+        var(pm, nw)[:Xti] = Xti
+        var(pm, nw)[:CCtr] = CCtr
+        var(pm, nw)[:CCti] = CCti
+
+        report && _IM.sol_component_value(pm, pmd_it_sym, nw, :transformer, :Xtr, tf_del_ids, Xtr)
+        report && _IM.sol_component_value(pm, pmd_it_sym, nw, :transformer, :Xti, tf_del_ids, Xti)
+        report && _IM.sol_component_value(pm, pmd_it_sym, nw, :transformer, :CCtr, tf_del_ids, CCtr)
+        report && _IM.sol_component_value(pm, pmd_it_sym, nw, :transformer, :CCti, tf_del_ids, CCti)
+    end
+
 end
 
 

src/form/bf_mx_soc.jl

@@ -90,3 +90,130 @@ function constraint_mc_voltage_psd(pm::SOCConicUBFModel, nw::Int, i::Int)
 
     relaxation_psd_to_soc_conic(pm.model, Wr, Wi)
 end
+
+
+@doc raw"""
+    constraint_mc_transformer_power_yy(pm::SOCUBFModels, nw::Int, trans_id::Int, f_bus::Int, t_bus::Int, f_idx::Tuple{Int,Int,Int}, t_idx::Tuple{Int,Int,Int}, f_connections::Vector{Int}, t_connections::Vector{Int}, pol::Int, tm_set::Vector{<:Real}, tm_fixed::Vector{Bool}, tm_scale::Real)
+
+Constraints to model a two-winding, wye-wye connected transformer.
+
+```math
+\begin{align}
+    & {W}_{fr} = {T}_{m}{T}_{m}^{H} {W}_{to}  \\
+    & {s}_{fr} + {s}_{to} = 0
+\end{align}
+```
+"""
+function constraint_mc_transformer_power_yy(pm::SOCUBFModels, nw::Int, trans_id::Int, f_bus::Int, t_bus::Int, f_idx::Tuple{Int,Int,Int}, t_idx::Tuple{Int,Int,Int}, f_connections::Vector{Int}, t_connections::Vector{Int}, pol::Int, tm_set::Vector{<:Real}, tm_fixed::Vector{Bool}, tm_scale::Real)
+    transformer = ref(pm, nw, :transformer, trans_id)
+    tm = [tm_fixed[idx] ? tm_set[idx] : var(pm, nw, :tap, trans_id)[idx] for (idx,(fc,tc)) in enumerate(zip(f_connections,t_connections))]
+
+    Wr_fr = var(pm, nw, :Wr, f_bus)
+    Wr_to = var(pm, nw, :Wr, t_bus)
+    Wi_fr = var(pm, nw, :Wi, f_bus)
+    Wi_to = var(pm, nw, :Wi, t_bus)
+
+    p_fr = [var(pm, nw, :pt, f_idx)[c] for c in f_connections]
+    p_to = [var(pm, nw, :pt, t_idx)[c] for c in t_connections]
+    q_fr = [var(pm, nw, :qt, f_idx)[c] for c in f_connections]
+    q_to = [var(pm, nw, :qt, t_idx)[c] for c in t_connections]
+
+    for (f_idx,fc) in enumerate(f_connections)
+        if haskey(transformer,"controls")  # regcontrol settings
+            w_fr = var(pm, nw, :w)[f_bus]
+            w_to = var(pm, nw, :w)[t_bus]
+            v_ref = transformer["controls"]["vreg"][f_idx]
+            δ = transformer["controls"]["band"][f_idx]
+            r = transformer["controls"]["r"][f_idx]
+            x = transformer["controls"]["x"][f_idx]
+        end
+
+        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])
+            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])
+
+                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)
+
+                # with regcontrol
+                if haskey(transformer,"controls")
+                    # voltage squared ignoring losses: w_drop = (2⋅r⋅p+2⋅x⋅q)
+                    w_drop = JuMP.@expression(pm.model, 2*r*p_to[idx] + 2*x*q_to[idx])
+
+                    # (v_ref-δ)^2 ≤ w_fr-w_drop ≤ (v_ref+δ)^2
+                    # w_fr/1.1^2 ≤ w_to ≤ w_fr/0.9^2
+                    JuMP.@constraint(pm.model, w_fr[fc] - w_drop ≥ (v_ref - δ)^2)
+                    JuMP.@constraint(pm.model, w_fr[fc] - w_drop ≤ (v_ref + δ)^2)
+                    JuMP.@constraint(pm.model, w_fr[fc]/1.1^2 ≤ w_to[tc])
+                    JuMP.@constraint(pm.model, w_fr[fc]/0.9^2 ≥ w_to[tc])
+                end
+            end
+        end
+    end
+
+    JuMP.@constraint(pm.model, p_fr + p_to .== 0)
+    JuMP.@constraint(pm.model, q_fr + q_to .== 0)
+end
+
+@doc raw"""
+    constraint_mc_transformer_power_dy(pm::SOCUBFModels, nw::Int, trans_id::Int, f_bus::Int, t_bus::Int, f_idx::Tuple{Int,Int,Int}, t_idx::Tuple{Int,Int,Int}, f_connections::Vector{Int}, t_connections::Vector{Int}, pol::Int, tm_set::Vector{<:Real}, tm_fixed::Vector{Bool}, tm_scale::Real)
+
+Constraints to model a two-winding, delta-wye connected transformer.
+
+```math
+\begin{align}
+    &{W}_{fr}^{ij}-{W}_{fr}^{ik}-{W}_{fr}^{lj}+{W}_{fr}^{lk} = t_m^2{W}_{to}^{ij} ~\forall i,j \in \{1,2,3\}~ \text{and}~ k,l \in \{2,3,1\}   \\
+    &{S}_{fr} = X_tT_t \\
+    &{S}_{fr}^\Delta = T_tX_t \\
+    & {s}_{fr}^\Delta + {s}_{to} = 0\\
+    & {M}_{\Delta} =
+    \begin{bmatrix}
+    {W}_{fr} & {X}_{t} \\
+     {X}_{t}^{\text{H}} &  {L}_{\Delta}
+    \end{bmatrix} \succeq 0
+\end{align}
+```
+"""
+function constraint_mc_transformer_power_dy(pm::SOCUBFModels, nw::Int, trans_id::Int, f_bus::Int, t_bus::Int, f_idx::Tuple{Int,Int,Int}, t_idx::Tuple{Int,Int,Int}, f_connections::Vector{Int}, t_connections::Vector{Int}, pol::Int, tm_set::Vector{<:Real}, tm_fixed::Vector{Bool}, tm_scale::Real)
+    nph = length(tm_set)
+    @assert length(f_connections) == length(t_connections) && nph == 3 "only phases == 3 dy transformers are currently supported"
+    next = Dict(c=>f_connections[idx%nph+1] for (idx,c) in enumerate(f_connections))
+
+    tm = [tm_fixed[idx] ? tm_set[idx] : var(pm, nw, :tap, trans_id)[idx] for (idx,(fc,tc)) in enumerate(zip(f_connections,t_connections))]
+
+    Wr_fr = var(pm, nw, :Wr, f_bus)
+    Wr_to = var(pm, nw, :Wr, t_bus)
+    Wi_fr = var(pm, nw, :Wi, f_bus)
+    Wi_to = var(pm, nw, :Wi, t_bus)
+    Xtr = var(pm, nw, :Xtr, trans_id)
+    Xti = var(pm, nw, :Xti, trans_id)
+    CCtr = var(pm, nw, :CCtr, trans_id)
+    CCti = var(pm, nw, :CCti, trans_id)
+
+    P_fr = var(pm, nw, :Pt, f_idx)
+    Q_fr = var(pm, nw, :Qt, f_idx)
+    p_to = [var(pm, nw, :pt, t_idx)[c] for c in t_connections]
+    q_to = [var(pm, nw, :qt, t_idx)[c] for c in t_connections]
+
+    for (f_idx,fc) in enumerate(f_connections)
+        for (t_idx,tc) in enumerate(t_connections)
+            JuMP.@constraint(pm.model, Wr_fr[fc,tc]-Wr_fr[fc,next[tc]]-Wr_fr[next[fc],tc]+Wr_fr[next[fc],next[tc]] == (pol*tm_scale)^2*tm[f_idx]*tm[t_idx]*Wr_to[fc,tc])
+            JuMP.@constraint(pm.model, Wi_fr[fc,tc]-Wi_fr[fc,next[tc]]-Wi_fr[next[fc],tc]+Wi_fr[next[fc],next[tc]] == (pol*tm_scale)^2*tm[f_idx]*tm[t_idx]*Wi_to[fc,tc])
+        end
+    end
+
+    Tt = [1 -1 0; 0 1 -1; -1 0 1]  # TODO
+    constraint_SWL_psd(pm.model, Xtr, Xti, Wr_fr, Wi_fr, CCtr, CCti) # link W, CCt and Xt
+    # define powers as affine transformations of X
+    JuMP.@constraint(pm.model, P_fr .== Xtr*Tt)
+    JuMP.@constraint(pm.model, Q_fr .== Xti*Tt)
+    JuMP.@constraint(pm.model, LinearAlgebra.diag(Tt*Xtr) + p_to .== 0)
+    JuMP.@constraint(pm.model, LinearAlgebra.diag(Tt*Xti) + q_to .== 0)
+end

test/transformer.jl

@@ -346,4 +346,21 @@
             @test all(isapprox.(result["solution"]["bus"]["tn_1"]["vm"], [1.045, 1.05]; atol=5E-3))
         end
     end
+
+    @testset "transformer SOC relaxations" begin
+        @testset "2w_yy" begin
+            result = solve_mc_opf(ut_trans_2w_yy_bank, SOCNLPUBFPowerModel, ipopt_solver; solution_processors=[sol_data_model!], make_si=false)
+            @test norm(result["solution"]["bus"]["3"]["vm"]-[0.82944, 0.86067, 0.72315], Inf) <= 2E-2
+        end
+
+        @testset "2w_dy_lead" begin
+            result = solve_mc_opf(ut_trans_2w_dy_lead, SOCConicUBFPowerModel, scs_solver; solution_processors=[sol_data_model!], make_si=false)
+            @test norm(result["solution"]["bus"]["3"]["vm"]-[0.87391, 0.86054, 0.85485], Inf) <= 4.2E-2
+        end
+
+        @testset "3w_dyy_1" begin
+            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
+    end
 end