From a9d4a95fc5710ccda9812797353ebebc3910ff00 Mon Sep 17 00:00:00 2001 From: Nick Laws Date: Fri, 31 May 2024 20:12:43 -0600 Subject: [PATCH] reasonable voltages for 2 line 1 load multiphase model --- docs/src/methods.md | 1 - src/BranchFlowModel.jl | 1 - src/model_multi_phase_network.jl | 172 +++++++++++----------------- test/data/two_line_multi_phase.yaml | 39 +++++++ test/runtests.jl | 72 +++++++----- 5 files changed, 149 insertions(+), 136 deletions(-) create mode 100644 test/data/two_line_multi_phase.yaml diff --git a/docs/src/methods.md b/docs/src/methods.md index caf5205..13b1797 100644 --- a/docs/src/methods.md +++ b/docs/src/methods.md @@ -17,7 +17,6 @@ add_variables constrain_power_balance constrain_substation_voltage constrain_KVL -constrain_loads constrain_bounds check_rank_one get_bus_values diff --git a/src/BranchFlowModel.jl b/src/BranchFlowModel.jl index 84c6ec4..6d3c5c0 100644 --- a/src/BranchFlowModel.jl +++ b/src/BranchFlowModel.jl @@ -35,7 +35,6 @@ export constrain_power_balance, constrain_substation_voltage, constrain_KVL, - constrain_loads, constrain_bounds, get_edge_values, # recover_voltage_current, # TODO validate this method diff --git a/src/model_multi_phase_network.jl b/src/model_multi_phase_network.jl index bc6c304..2e02386 100644 --- a/src/model_multi_phase_network.jl +++ b/src/model_multi_phase_network.jl @@ -7,14 +7,12 @@ add_variables(m, net) constrain_power_balance(m, net) constrain_substation_voltage(m, net) constrain_KVL(m, net) -constrain_loads(m, net) ``` """ function build_model!(m::JuMP.AbstractModel, net::Network{MultiPhase}) add_variables(m, net) # PSD done in add_variables constrain_power_balance(m, net) constrain_KVL(m, net) - constrain_loads(m, net) end # hack for zero problem in MutableArithmetics @@ -93,7 +91,7 @@ function add_variables(m, net::Network{MultiPhase}) end for t in 1:net.Ntimesteps - m[:w][t] = Dict(net.substation_bus => v0) + m[:w][t] = Dict(net.substation_bus => v0 * cj(v0)) # empty dicts for line values in each time step to fill m[:l][t] = Dict() m[:Sij][t] = Dict() @@ -121,30 +119,6 @@ function add_variables(m, net::Network{MultiPhase}) ) end - # Sj is 3x1 - m[:Sj][t][j] = convert(Matrix{GenericAffExpr{ComplexF64, VariableRef}}, reshape([0.0im; 0im; 0im], 3, 1)) - for phs in phases_into_bus(net, j) # fill in Sj vector with variables for phases going into bus j - real_load = 0.0 - if j in real_load_busses(net) - if phs in 1:3 - real_load = net[j, :kws, phs][t] - end - end - - reactive_load = 0.0 - if j in reactive_load_busses(net) - if phs in 1:3 - reactive_load = net[j, :kvars, phs][t] - end - end - m[:Sj][t][j][phs] = @variable(m, - set = ComplexPlane(), base_name="Sj_" * string(t) *"_"* j *"_"* string(phs), - # lower_bound = p.P_lo_bound + p.Q_lo_bound*im, - # upper_bound = p.P_up_bound + p.Q_up_bound*im, - start = -real_load - reactive_load*im - ) - end - # Hermitian variables m[:l][t][i_j] = convert(Matrix{GenericAffExpr{ComplexF64, VariableRef}}, [0 0im 0im; 0im 0. 0im; 0im 0im 0]) m[:w][t][j] = convert(Matrix{GenericAffExpr{ComplexF64, VariableRef}}, [0 0im 0im; 0im 0. 0im; 0im 0im 0]) @@ -210,10 +184,10 @@ function add_variables(m, net::Network{MultiPhase}) i = net.substation_bus define_vars_downstream(i, t, m, net) - function recursive_variables(i::String, t::Int, m::JuMP.AbstractModel, net::Network) - for j in j_to_k(i, net) - define_vars_downstream(j, t, m, net) - recursive_variables(j, t, m, net) + function recursive_variables(j::String, t::Int, m::JuMP.AbstractModel, net::Network) + for k in j_to_k(j, net) + define_vars_downstream(k, t, m, net) + recursive_variables(k, t, m, net) end end recursive_variables(i, t, m, net) @@ -235,40 +209,66 @@ as the first index. For example `m[:loadbalcons]["busname"]` will give the const from JuMP for all time steps. """ function constrain_power_balance(m, net::Network{MultiPhase}) - Sj = m[:Sj] Sij = m[:Sij] lij = m[:l] m[:loadbalcons] = Dict() - # TODO change Pj and Qj to expressions, make P₀ and Q₀ dv's, which will reduce # of variables - # by (Nnodes - 1)*8760 and number of constraints by 6*(Nnodes - 1)*8760 + for j in busses(net) - if isempty(i_to_j(j, net)) && !isempty(j_to_k(j, net)) # source nodes, injection = flows out - con = @constraint(m, [t in 1:net.Ntimesteps], - Sj[t][j] - sum( diag( Sij[t][(j,k)] ) for k in j_to_k(j, net) ) .== 0 - ) - elseif isempty(i_to_j(j, net)) && isempty(j_to_k(j, net)) # unconnected nodes - @warn "Bus $j has no edges, setting Sj and Qj to zero." - con = @constraint(m, [t in 1:net.Ntimesteps], - diag(Sj[t][j]) .== 0 - ) - elseif !isempty(i_to_j(j, net)) && isempty(j_to_k(j, net)) # leaf nodes / sinks, flows in = draw out - con = @constraint(m, [t in 1:net.Ntimesteps], + + # check for loads + Pj = [zeros(net.Ntimesteps) for _ in 1:3] # first dim is phase, like Pj[phs][t] + Qj = [zeros(net.Ntimesteps) for _ in 1:3] + if j in real_load_busses(net) + for phs in 1:3 # put an Sj method in CommonOPF? to make complex vector of phases + Pj[phs] = -net[j, :kws, phs] * 1e3 / net.Sbase + end + end + if j in reactive_load_busses(net) + for phs in 1:3 + Qj[phs] = -net[j, :kvars, phs] * 1e3 / net.Sbase + end + end + Sj = Pj + im * Qj + Sj = hcat(Sj...) # time X phase + + # source nodes, injection = flows out + if isempty(i_to_j(j, net)) && !isempty(j_to_k(j, net)) + + if j == net.substation_bus # include the slack power variables + m[:loadbalcons][j] = @constraint(m, [t in 1:net.Ntimesteps], + m[:Sj][t][j] + Sj[t, :] - sum( diag( Sij[t][(j,k)] ) for k in j_to_k(j, net) ) .== 0 + ) + else # a source node with known injection + m[:loadbalcons][j] = @constraint(m, [t in 1:net.Ntimesteps], + Sj[t, :] - sum( diag( Sij[t][(j,k)] ) for k in j_to_k(j, net) ) .== 0 + ) + + end + + # unconnected nodes + elseif isempty(i_to_j(j, net)) && isempty(j_to_k(j, net)) + @warn "Bus $j has no edges. Setting loadbalcons to zeros" + m[:loadbalcons][j] = zeros(net.Ntimesteps) + + # leaf nodes / sinks, flows in = draw out + elseif !isempty(i_to_j(j, net)) && isempty(j_to_k(j, net)) + m[:loadbalcons][j] = @constraint(m, [t in 1:net.Ntimesteps], sum( diag( - Sij[t][(i,j)] - zij(i,j,net) * lij[t][(i,j)] + Sij[t][(i,j)] - zij_per_unit(i,j,net) * lij[t][(i,j)] ) for i in i_to_j(j, net) ) - + Sj[t][j] .== 0 + + Sj[t, :] .== 0 ) - else # node with lines in and out - # TODO failing for IEEE13 with i = "632", j = "645" - con = @constraint(m, [t in 1:net.Ntimesteps], + + # node with lines in and out + else + m[:loadbalcons][j] = @constraint(m, [t in 1:net.Ntimesteps], sum( diag( - Sij[t][(i,j)] - zij(i,j,net) * lij[t][(i,j)] + Sij[t][(i,j)] - zij_per_unit(i,j,net) * lij[t][(i,j)] ) for i in i_to_j(j, net) ) - + Sj[t][j] + + Sj[t, :] - sum( diag( Sij[t][(j,k)] ) for k in j_to_k(j, net) ) .== 0 ) end - m[:loadbalcons][j] = con end # TODO add shunts, diag( openDSS-cmatrix * m[:w][t][j] ) ? @@ -282,7 +282,7 @@ end Used in defining the KVL constraints, this method returns `M` for only the indices in `phases`. """ -function matrix_phases_to_vec(M::AbstractMatrix{T}, phases::AbstractSet{Int}) where T +function matrix_phases_to_vec(M::AbstractMatrix{T}, phases::AbstractVector{Int}) where T v = T[] for i in phases, j in phases push!(v, M[i,j]) @@ -305,64 +305,20 @@ function constrain_KVL(m, net::Network{MultiPhase}) m[:kvl] = Dict{String, AbstractArray}() for j in busses(net) # substation_bus in here but has empty i_to_j(j, net) for i in i_to_j(j, net) # for radial network there is only one i in i_to_j - phs = phases_into_bus(net, j) - i_j = (i, j) - z = zij(i,j,net) - # need to slice w[t][i] by phases in i_j + phases = phases_into_bus(net, j) + z = zij_per_unit(i,j,net) + # need to slice w[t][i] by phases in edge (i, j) con = @constraint(m, [t in 1:net.Ntimesteps], - T(w[t][j], phs) .== T(w[t][i], phs) - - T(Sij[t][i_j] * cj(z) + z * cj(Sij[t][i_j]), phs) - + T(z * lij[t][i_j] * cj(z), phs) + T( w[t][j], phases ) .== T( w[t][i], phases ) + - T( + Sij[t][(i, j)] * cj(z) + z * cj(Sij[t][(i, j)]), + phases) + + T( + z * lij[t][(i, j)] * cj(z), + phases) ); m[:kvl][j] = con end end nothing end - - -""" - constrain_loads(m, net::Network{MultiPhase}) - -- set loads to negative of Inputs.Pload and Inputs.Qload, - which are normalized by Sbase when creating Inputs. -- keys of Pload and Qload must match Inputs.busses. Any missing keys have load set to zero. -- Inputs.substation_bus is unconstrained, slack bus - -Each of the power injection constraints are stored in the model under `m[:injectioncons]`. -To acces the constraints you can: -```julia -m[:injectioncons]["busname"][t,phs] -``` -where `t` is the integer time step and `phs` is the integer phase. -""" -function constrain_loads(m, net::Network{BranchFlowModel.MultiPhase}) - Sj = m[:Sj] - m[:injectioncons] = Dict() - for j in setdiff(busses(net), [net.substation_bus]) - - # default to zero injection - real_load = Dict(phs => zeros(net.Ntimesteps) for phs in [1,2,3]) - - if j in real_load_busses(net) - for phs in 1:3 - real_load[phs] = net[j, :kws, phs] * 1e3 - end - end - - # default to zero injection - reactive_load = Dict(phs => zeros(net.Ntimesteps) for phs in [1,2,3]) - - if j in reactive_load_busses(net) - for phs in 1:3 - reactive_load[phs] = net[j, :kvars, phs] * 1e3 - end - end - - con = @constraint(m, [t in 1:net.Ntimesteps, phs in 1:3], - Sj[t][j][phs] == (-real_load[phs][t] - reactive_load[phs][t]im ) / net.Sbase - ) - m[:injectioncons][j] = con - end - nothing -end diff --git a/test/data/two_line_multi_phase.yaml b/test/data/two_line_multi_phase.yaml new file mode 100644 index 0000000..cb56cb9 --- /dev/null +++ b/test/data/two_line_multi_phase.yaml @@ -0,0 +1,39 @@ +Network: + substation_bus: b1 + Sbase: 1e3 + Vbase: 1e3 + Ntimesteps: 1 + +Conductor: + - name: cond1-symmetric + busses: + - b1 + - b2 + r0: 0.766 + x0: 1.944 + r1: 0.301 + x1: 0.627 + length: 10 + phases: [1, 2, 3] + - name: cond2-copy-cond1 + busses: + - b2 + - b3 + template: cond1-symmetric + length: 10 + phases: [1, 2, 3] + +Load: + - bus: b3 + kws1: + - 5.6 + kvars1: + - 1.2 + kws2: + - 5.6 + kvars2: + - 1.2 + kws3: + - 5.6 + kvars3: + - 1.2 diff --git a/test/runtests.jl b/test/runtests.jl index 45e5242..b3c1268 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -208,17 +208,40 @@ end end +@testset "basic two-line multiphase" begin + net = BranchFlowModel.CommonOPF.Network(joinpath("data", "two_line_multi_phase.yaml")) + + m = Model(CSDP.Optimizer) + build_model!(m, net) + + @objective(m, Min, + sum( sum(real.(diag(m[:l][t][i_j]))) for t in 1:net.Ntimesteps, i_j in edges(net)) + ) + + optimize!(m) + + vs = Dict( + k => abs.(sqrt(JuMP.value.(w))) + for (k,w) in m[:w][1] + ) + +end + + @testset "ieee13 unbalanced MultiPhase" begin - # make the dss solution to compare - dssfilepath = "data/ieee13/IEEE13Nodeckt.dss" - OpenDSS.Text.Command("Redirect $dssfilepath") - @test(OpenDSS.Solution.Converged() == true) + # # make the dss solution to compare + # dssfilepath = "data/ieee13/IEEE13Nodeckt.dss" + # OpenDSS.Text.Command("Redirect $dssfilepath") + # @test(OpenDSS.Solution.Converged() == true) - dss_voltages = dss_voltages_pu() + # dss_voltages = dss_voltages_pu() - vbase = 4160/sqrt(3) - net = BranchFlowModel.CommonOPF.dss_to_Network(dssfilepath) + # vbase = 4160/sqrt(3) + # net = BranchFlowModel.CommonOPF.dss_to_Network(dssfilepath) + # net.Sbase = 1_000_000 + # net.Vbase = vbase + # net.Zbase = net.Vbase^2 / net.Sbase # p = Inputs( # joinpath("data", "ieee13", "IEEE13Nodeckt.dss"), @@ -238,30 +261,30 @@ end # p.P_up_bound = 10 # p.Q_up_bound = 10 - m = Model(CSDP.Optimizer) -# set_attribute(m, "printlevel", 0) +# m = Model(CSDP.Optimizer) +# # set_attribute(m, "printlevel", 0) - # m = Model(COSMO.Optimizer) +# # m = Model(COSMO.Optimizer) - # m = Model(SCS.Optimizer) +# # m = Model(SCS.Optimizer) - build_model!(m, net) +# build_model!(m, net) - @objective(m, Min, - sum( sum(real.(diag(m[:l][t][i_j]))) for t in 1:net.Ntimesteps, i_j in edges(net)) - ) +# @objective(m, Min, +# sum( sum(real.(diag(m[:l][t][i_j]))) for t in 1:net.Ntimesteps, i_j in edges(net) ) +# ) - # need bounds to get solution? - # optimize!(m) +# # need bounds to get solution? +# optimize!(m) - # @test termination_status(m) in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL] +# # @test termination_status(m) in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL] - # @test_nowarn(check_rank_one(m,p)) +# # @test_nowarn(check_rank_one(m,p)) - # vs = Dict( - # k => sqrt(JuMP.value.(w)) - # for (k,w) in m[:w][1] - # ) +# vs = Dict( +# k => abs.(sqrt(JuMP.value.(w))) +# for (k,w) in m[:w][1] +# ) # I = get_variable_values(:l, m, p) # TODO method for multiphase results @@ -271,9 +294,6 @@ end # println("$b - $i $(phsv - dss_voltages[b][i])") # end # end - - -# # TODO why are BFM voltages not matching openDSS well? end