diff --git a/CHANGELOG.md b/CHANGELOG.md index a87696708..9249d7ecd 100644 --- 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 index 524b5bef3..0954da408 100644 --- 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,14 +841,16 @@ 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) @@ -856,6 +858,10 @@ function constraint_mc_load_power(pm::AbstractUBFModels, load_id::Int; nw::Int=n 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 index 0c741e217..0e59d9f06 100644 --- 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 index 5f0cf13cf..8413f29d8 100644 --- 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 index ea94938d8..bf28baaee 100644 --- 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