Skip to content

Commit

Permalink
ADD: Transformer SOC relaxations (#412)
Browse files Browse the repository at this point in the history 
* Added transformer SOC relaxations

* Updated unit tests

* Update bf_mx_soc.jl

* Update transformer.jl

Co-authored-by: David M Fobes <[email protected]>
  • Loading branch information
@keatsig @pseudocubic
keatsig and pseudocubic authored Dec 23, 2022
1 parent 7188b0f commit ce5187d
Show file tree
Hide file tree
Showing 5 changed files with 191 additions and 0 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
20 changes: 20 additions & 0 deletions src/core/solution.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
26 changes: 26 additions & 0 deletions src/form/bf_mx.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"])
Expand All @@ -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")
Expand All @@ -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


Expand Down
127 changes: 127 additions & 0 deletions src/form/bf_mx_soc.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
17 changes: 17 additions & 0 deletions test/transformer.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit ce5187d

Please sign in to comment.