From ce5187d6c67cb3fe21a327a68e53d7dc2cee2714 Mon Sep 17 00:00:00 2001 From: Kshitij I Girigoudar <42916010+keatsig@users.noreply.github.com> Date: Fri, 23 Dec 2022 10:09:28 -0700 Subject: [PATCH] ADD: Transformer SOC relaxations (#412) * Added transformer SOC relaxations * Updated unit tests * Update bf_mx_soc.jl * Update transformer.jl Co-authored-by: David M Fobes --- CHANGELOG.md | 1 + src/core/solution.jl | 20 +++++++ src/form/bf_mx.jl | 26 +++++++++ src/form/bf_mx_soc.jl | 127 ++++++++++++++++++++++++++++++++++++++++++ test/transformer.jl | 17 ++++++ 5 files changed, 191 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 725ffcfa6..eaa0e9117 100644 --- a/CHANGELOG.md +++ b/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) diff --git a/src/core/solution.jl b/src/core/solution.jl index fa724e708..2c85a05f6 100644 --- a/src/core/solution.jl +++ b/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 diff --git a/src/form/bf_mx.jl b/src/form/bf_mx.jl index ae8465f7b..a69886d5a 100644 --- a/src/form/bf_mx.jl +++ b/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 diff --git a/src/form/bf_mx_soc.jl b/src/form/bf_mx_soc.jl index aaedbc05d..5f0cf13cf 100644 --- a/src/form/bf_mx_soc.jl +++ b/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 diff --git a/test/transformer.jl b/test/transformer.jl index bf881d1d5..ea94938d8 100644 --- a/test/transformer.jl +++ b/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