diff --git a/CHANGELOG.md b/CHANGELOG.md index 440f50d0f..2b1b84756 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## staged +- Fixed indexing issue for single-phase delta load models in linear formulations (LinDist3Flow, FOTP, FOTR, FBS) - Added ZIP load model - Updated documentation in `make_multiconductor!` to better indicate its unsupported nature - Added automatic detection of multinetwork data to `instantiate_mc_model` diff --git a/src/data_model/utils.jl b/src/data_model/utils.jl index dae4ba140..45f234004 100644 --- a/src/data_model/utils.jl +++ b/src/data_model/utils.jl @@ -535,7 +535,7 @@ function _pad_properties_delta!(object::Dict{String,<:Any}, properties::Vector{S merge!(tmp, Dict((k[2], k[1])=>sign*v for (k,v) in tmp)) get_val(x,y) = haskey(tmp, (x,y)) ? tmp[(x,y)] : 0.0 - object[property] = [get_val(phases[1], phases[2]), get_val(phases[2], phases[3]), get_val(phases[3], phases[1])] + object[property] = val_length==1 ? [get_val(connections[1], connections[2]), 0.0,0.0] : [get_val(phases[1], phases[2]), get_val(phases[2], phases[3]), get_val(phases[3], phases[1])] end end end diff --git a/src/form/bf_fbs.jl b/src/form/bf_fbs.jl index 5d9d99ff8..86c332037 100644 --- a/src/form/bf_fbs.jl +++ b/src/form/bf_fbs.jl @@ -590,23 +590,33 @@ function constraint_mc_load_power(pm::FBSUBFPowerModel, load_id::Int; nw::Int=nw # zero-order approximation elseif load["configuration"]==DELTA - vr0 = [var(pm, nw, :vr0, bus_id)[findfirst(isequal(c), terminals)] for c in connections] - vi0 = [var(pm, nw, :vi0, bus_id)[findfirst(isequal(c), terminals)] for c in connections] + vr0 = var(pm, nw, :vr0, bus_id) + vi0 = var(pm, nw, :vi0, bus_id) - nph = length(connections) - prev = Dict(i=>(i+nph-2)%nph+1 for i in 1:nph) - next = Dict(i=>i%nph+1 for i in 1:nph) + nph = length(a) - vrd0 = [vr0[i]-vr0[next[i]] for i in 1:nph] - vid0 = [vi0[i]-vi0[next[i]] for i in 1:nph] + prev = Dict(c=>connections[(idx+nph-2)%nph+1] for (idx,c) in enumerate(connections)) + next = Dict(c=>connections[idx%nph+1] for (idx,c) in enumerate(connections)) - crd0 = [a[i]*vrd0[i]*(vrd0[i]^2+vid0[i]^2)^(alpha[i]/2-1)+b[i]*vid0[i]*(vrd0[i]^2+vid0[i]^2)^(beta[i]/2 -1) for i in 1:nph] - cid0 = [a[i]*vid0[i]*(vrd0[i]^2+vid0[i]^2)^(alpha[i]/2-1)-b[i]*vrd0[i]*(vrd0[i]^2+vid0[i]^2)^(beta[i]/2 -1) for i in 1:nph] - crd0_bus = [crd0[i]-crd0[prev[i]] for i in 1:nph] - cid0_bus = [cid0[i]-cid0[prev[i]] for i in 1:nph] + vrd0 = [vr0[idx]-vr0[next[idx]] for (idx, c) in enumerate(connections)] + vid0 = [vi0[idx]-vi0[next[idx]] for (idx, c) in enumerate(connections)] + + crd0 = Array{Any,1}(undef, nph) + cid0 = Array{Any,1}(undef, nph) + for (idx, c) in enumerate(connections) + crd0[c] = a[idx]*vrd0[c]*(vrd0[c]^2+vid0[c]^2)^(alpha[idx]/2-1)+b[idx]*vid0[c]*(vrd0[c]^2+vid0[c]^2)^(beta[idx]/2 -1) + cid0[c] = a[idx]*vid0[c]*(vrd0[c]^2+vid0[c]^2)^(alpha[idx]/2-1)-b[idx]*vrd0[c]*(vrd0[c]^2+vid0[c]^2)^(beta[idx]/2 -1) + end + + crd0_bus = [crd0[idx]-crd0[prev[idx]] for (idx, c) in enumerate(connections)] + cid0_bus = [cid0[idx]-cid0[prev[idx]] for (idx, c) in enumerate(connections)] + + pd_bus = [ vr0[c]*crd0_bus[c]+vi0[c]*cid0_bus[c] for (idx,c) in enumerate(connections)] + qd_bus = [-vr0[c]*cid0_bus[c]+vi0[c]*crd0_bus[c] for (idx,c) in enumerate(connections)] + + pd_bus = JuMP.Containers.DenseAxisArray(pd_bus, connections) + qd_bus = JuMP.Containers.DenseAxisArray(qd_bus, connections) - pd_bus = [ vr0[i]*crd0_bus[i]+vi0[i]*cid0_bus[i] for i in 1:nph] - qd_bus = [-vr0[i]*cid0_bus[i]+vi0[i]*crd0_bus[i] for i in 1:nph] var(pm, nw, :pd_bus)[load_id] = pd_bus var(pm, nw, :qd_bus)[load_id] = qd_bus @@ -614,9 +624,12 @@ function constraint_mc_load_power(pm::FBSUBFPowerModel, load_id::Int; nw::Int=nw sol(pm, nw, :load, load_id)[:pd_bus] = pd_bus sol(pm, nw, :load, load_id)[:qd_bus] = qd_bus - pd = JuMP.@expression(pm.model, [i in 1:nph], a[i]*(vrd0[i]^2+vid0[i]^2)^(alpha[i]/2) ) - qd = JuMP.@expression(pm.model, [i in 1:nph], b[i]*(vrd0[i]^2+vid0[i]^2)^(beta[i]/2) ) - + pd = [] + qd = [] + for (idx,c) in enumerate(connections) + push!(pd, JuMP.@expression(pm.model, a[idx]*(vrd0[c]^2+vid0[c]^2)^(alpha[idx]/2) )) + push!(qd, JuMP.@expression(pm.model, b[idx]*(vrd0[c]^2+vid0[c]^2)^(beta[idx]/2) )) + end sol(pm, nw, :load, load_id)[:pd] = pd sol(pm, nw, :load, load_id)[:qd] = qd end diff --git a/src/form/bf_mx_lin.jl b/src/form/bf_mx_lin.jl index 342a4c6c9..9680fe2ab 100644 --- a/src/form/bf_mx_lin.jl +++ b/src/form/bf_mx_lin.jl @@ -426,7 +426,8 @@ function constraint_mc_load_power(pm::LPUBFDiagModel, load_id::Int; nw::Int=nw_i JuMP.@constraint(pm.model, Xdi[:,idx] .== 0) end end - + pd_bus = JuMP.Containers.DenseAxisArray(pd_bus, connections) + qd_bus = JuMP.Containers.DenseAxisArray(qd_bus, 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 @@ -437,23 +438,23 @@ function constraint_mc_load_power(pm::LPUBFDiagModel, load_id::Int; nw::Int=nw_i JuMP.@constraint(pm.model, qd[idx]==qd0[idx]) end elseif load["model"]==IMPEDANCE - w = var(pm, nw, :w)[bus_id][[c for c in connections]] + w = var(pm, nw, :w)[bus_id] for (idx,c) in enumerate(connections) - JuMP.@constraint(pm.model, pd[idx]==3*a[idx]*w[idx]) - JuMP.@constraint(pm.model, qd[idx]==3*b[idx]*w[idx]) + JuMP.@constraint(pm.model, pd[idx]==3*a[idx]*w[c]) + JuMP.@constraint(pm.model, qd[idx]==3*b[idx]*w[c]) end else - w = var(pm, nw, :w)[bus_id][[c for c in connections]] + w = var(pm, nw, :w)[bus_id] for (idx,c) in enumerate(connections) - JuMP.@constraint(pm.model, pd[idx]==sqrt(3)/2*a[idx]*(w[idx]+1)) - JuMP.@constraint(pm.model, qd[idx]==sqrt(3)/2*b[idx]*(w[idx]+1)) + JuMP.@constraint(pm.model, pd[idx]==sqrt(3)/2*a[idx]*(w[c]+1)) + JuMP.@constraint(pm.model, qd[idx]==sqrt(3)/2*b[idx]*(w[c]+1)) end end ## reporting; for delta these are not available as saved variables! if report - sol(pm, nw, :load, load_id)[:pd] = pd - sol(pm, nw, :load, load_id)[:qd] = qd + sol(pm, nw, :load, load_id)[:pd] = JuMP.Containers.DenseAxisArray(pd, connections) + sol(pm, nw, :load, load_id)[:qd] = JuMP.Containers.DenseAxisArray(qd, connections) sol(pm, nw, :load, load_id)[:pd_bus] = pd_bus sol(pm, nw, :load, load_id)[:qd_bus] = qd_bus end diff --git a/src/form/fotp.jl b/src/form/fotp.jl index 444b78964..29f0facdf 100644 --- a/src/form/fotp.jl +++ b/src/form/fotp.jl @@ -490,25 +490,34 @@ function constraint_mc_load_power(pm::FOTPUPowerModel, load_id::Int; nw::Int=nw_ # zero-order approximation elseif load["configuration"]==DELTA - vm0 = [var(pm, nw, :vm0, bus_id)[findfirst(isequal(c), bus["terminals"])] for c in connections] - va0 = [var(pm, nw, :va0, bus_id)[findfirst(isequal(c), bus["terminals"])] for c in connections] + vm0 = var(pm, nw, :vm0, bus_id) + va0 = var(pm, nw, :va0, bus_id) vr0 = vm0.*cos.(va0) vi0 = vm0.*sin.(va0) + nph = length(a) + + prev = Dict(c=>connections[(idx+nph-2)%nph+1] for (idx,c) in enumerate(connections)) + next = Dict(c=>connections[idx%nph+1] for (idx,c) in enumerate(connections)) + + vrd0 = [vr0[idx]-vr0[next[idx]] for (idx, c) in enumerate(connections)] + vid0 = [vi0[idx]-vi0[next[idx]] for (idx, c) in enumerate(connections)] + + crd0 = Array{Any,1}(undef, nph) + cid0 = Array{Any,1}(undef, nph) + for (idx, c) in enumerate(connections) + crd0[c] = a[idx]*vrd0[c]*(vrd0[c]^2+vid0[c]^2)^(alpha[idx]/2-1)+b[idx]*vid0[c]*(vrd0[c]^2+vid0[c]^2)^(beta[idx]/2 -1) + cid0[c] = a[idx]*vid0[c]*(vrd0[c]^2+vid0[c]^2)^(alpha[idx]/2-1)-b[idx]*vrd0[c]*(vrd0[c]^2+vid0[c]^2)^(beta[idx]/2 -1) + end - nph = length(connections) - prev = Dict(i=>(i+nph-2)%nph+1 for i in 1:nph) - next = Dict(i=>i%nph+1 for i in 1:nph) + crd0_bus = [crd0[idx]-crd0[prev[idx]] for (idx, c) in enumerate(connections)] + cid0_bus = [cid0[idx]-cid0[prev[idx]] for (idx, c) in enumerate(connections)] - vrd0 = [vr0[i]-vr0[next[i]] for i in 1:nph] - vid0 = [vi0[i]-vi0[next[i]] for i in 1:nph] + pd_bus = [ vr0[c]*crd0_bus[c]+vi0[c]*cid0_bus[c] for (idx, c) in enumerate(connections)] + qd_bus = [-vr0[c]*cid0_bus[c]+vi0[c]*crd0_bus[c] for (idx, c) in enumerate(connections)] - crd0 = [a[i]*vrd0[i]*(vrd0[i]^2+vid0[i]^2)^(alpha[i]/2-1)+b[i]*vid0[i]*(vrd0[i]^2+vid0[i]^2)^(beta[i]/2 -1) for i in 1:nph] - cid0 = [a[i]*vid0[i]*(vrd0[i]^2+vid0[i]^2)^(alpha[i]/2-1)-b[i]*vrd0[i]*(vrd0[i]^2+vid0[i]^2)^(beta[i]/2 -1) for i in 1:nph] - crd0_bus = [crd0[i]-crd0[prev[i]] for i in 1:nph] - cid0_bus = [cid0[i]-cid0[prev[i]] for i in 1:nph] + pd_bus = JuMP.Containers.DenseAxisArray(pd_bus, connections) + qd_bus = JuMP.Containers.DenseAxisArray(qd_bus, connections) - pd_bus = [ vr0[i]*crd0_bus[i]+vi0[i]*cid0_bus[i] for i in 1:nph] - qd_bus = [-vr0[i]*cid0_bus[i]+vi0[i]*crd0_bus[i] for i in 1:nph] var(pm, nw, :pd_bus)[load_id] = pd_bus var(pm, nw, :qd_bus)[load_id] = qd_bus @@ -516,9 +525,12 @@ function constraint_mc_load_power(pm::FOTPUPowerModel, load_id::Int; nw::Int=nw_ sol(pm, nw, :load, load_id)[:pd_bus] = pd_bus sol(pm, nw, :load, load_id)[:qd_bus] = qd_bus - pd = JuMP.@expression(pm.model, [i in 1:nph], a[i]*(vrd0[i]^2+vid0[i]^2)^(alpha[i]/2) ) - qd = JuMP.@expression(pm.model, [i in 1:nph], b[i]*(vrd0[i]^2+vid0[i]^2)^(beta[i]/2) ) - + pd = [] + qd = [] + for (idx,c) in enumerate(connections) + push!(pd, JuMP.@expression(pm.model, a[idx]*(vrd0[c]^2+vid0[c]^2)^(alpha[idx]/2) )) + push!(qd, JuMP.@expression(pm.model, b[idx]*(vrd0[c]^2+vid0[c]^2)^(beta[idx]/2) )) + end sol(pm, nw, :load, load_id)[:pd] = pd sol(pm, nw, :load, load_id)[:qd] = qd end diff --git a/src/form/fotr.jl b/src/form/fotr.jl index 63ae88942..0eaae0f4e 100644 --- a/src/form/fotr.jl +++ b/src/form/fotr.jl @@ -518,23 +518,33 @@ function constraint_mc_load_power(pm::FOTRUPowerModel, load_id::Int; nw::Int=nw_ # zero-order approximation elseif load["configuration"]==DELTA - vr0 = [var(pm, nw, :vr0, bus_id)[findfirst(isequal(c), terminals)] for c in connections] - vi0 = [var(pm, nw, :vi0, bus_id)[findfirst(isequal(c), terminals)] for c in connections] + vr0 = var(pm, nw, :vr0, bus_id) + vi0 = var(pm, nw, :vi0, bus_id) - nph = length(connections) - prev = Dict(i=>(i+nph-2)%nph+1 for i in 1:nph) - next = Dict(i=>i%nph+1 for i in 1:nph) + nph = length(a) - vrd0 = [vr0[i]-vr0[next[i]] for i in 1:nph] - vid0 = [vi0[i]-vi0[next[i]] for i in 1:nph] + prev = Dict(c=>connections[(idx+nph-2)%nph+1] for (idx,c) in enumerate(connections)) + next = Dict(c=>connections[idx%nph+1] for (idx,c) in enumerate(connections)) - crd0 = [a[i]*vrd0[i]*(vrd0[i]^2+vid0[i]^2)^(alpha[i]/2-1)+b[i]*vid0[i]*(vrd0[i]^2+vid0[i]^2)^(beta[i]/2 -1) for i in 1:nph] - cid0 = [a[i]*vid0[i]*(vrd0[i]^2+vid0[i]^2)^(alpha[i]/2-1)-b[i]*vrd0[i]*(vrd0[i]^2+vid0[i]^2)^(beta[i]/2 -1) for i in 1:nph] - crd0_bus = [crd0[i]-crd0[prev[i]] for i in 1:nph] - cid0_bus = [cid0[i]-cid0[prev[i]] for i in 1:nph] + vrd0 = [vr0[idx]-vr0[next[idx]] for (idx, c) in enumerate(connections)] + vid0 = [vi0[idx]-vi0[next[idx]] for (idx, c) in enumerate(connections)] + + crd0 = Array{Any,1}(undef, nph) + cid0 = Array{Any,1}(undef, nph) + for (idx, c) in enumerate(connections) + crd0[c] = a[idx]*vrd0[c]*(vrd0[c]^2+vid0[c]^2)^(alpha[idx]/2-1)+b[idx]*vid0[c]*(vrd0[c]^2+vid0[c]^2)^(beta[idx]/2 -1) + cid0[c] = a[idx]*vid0[c]*(vrd0[c]^2+vid0[c]^2)^(alpha[idx]/2-1)-b[idx]*vrd0[c]*(vrd0[c]^2+vid0[c]^2)^(beta[idx]/2 -1) + end + + crd0_bus = [crd0[idx]-crd0[prev[idx]] for (idx, c) in enumerate(connections)] + cid0_bus = [cid0[idx]-cid0[prev[idx]] for (idx, c) in enumerate(connections)] + + pd_bus = [ vr0[c]*crd0_bus[c]+vi0[c]*cid0_bus[c] for (idx,c) in enumerate(connections)] + qd_bus = [-vr0[c]*cid0_bus[c]+vi0[c]*crd0_bus[c] for (idx,c) in enumerate(connections)] + + pd_bus = JuMP.Containers.DenseAxisArray(pd_bus, connections) + qd_bus = JuMP.Containers.DenseAxisArray(qd_bus, connections) - pd_bus = [ vr0[i]*crd0_bus[i]+vi0[i]*cid0_bus[i] for i in 1:nph] - qd_bus = [-vr0[i]*cid0_bus[i]+vi0[i]*crd0_bus[i] for i in 1:nph] var(pm, nw, :pd_bus)[load_id] = pd_bus var(pm, nw, :qd_bus)[load_id] = qd_bus @@ -542,9 +552,12 @@ function constraint_mc_load_power(pm::FOTRUPowerModel, load_id::Int; nw::Int=nw_ sol(pm, nw, :load, load_id)[:pd_bus] = pd_bus sol(pm, nw, :load, load_id)[:qd_bus] = qd_bus - pd = JuMP.@expression(pm.model, [i in 1:nph], a[i]*(vrd0[i]^2+vid0[i]^2)^(alpha[i]/2) ) - qd = JuMP.@expression(pm.model, [i in 1:nph], b[i]*(vrd0[i]^2+vid0[i]^2)^(beta[i]/2) ) - + pd = [] + qd = [] + for (idx,c) in enumerate(connections) + push!(pd, JuMP.@expression(pm.model, a[idx]*(vrd0[c]^2+vid0[c]^2)^(alpha[idx]/2) )) + push!(qd, JuMP.@expression(pm.model, b[idx]*(vrd0[c]^2+vid0[c]^2)^(beta[idx]/2) )) + end sol(pm, nw, :load, load_id)[:pd] = pd sol(pm, nw, :load, load_id)[:qd] = qd end diff --git a/test/loadmodels.jl b/test/loadmodels.jl index 6c8e38542..55962eb9a 100644 --- a/test/loadmodels.jl +++ b/test/loadmodels.jl @@ -7,8 +7,8 @@ vbase = case3_lm_1230["settings"]["vbases_default"]["sourcebus"] @test isapprox(vm(result, "loadbus") ./ vbase, [0.999993, 0.999992, 0.999993]; atol=1E-5) # single-phase delta loads - @test isapprox(pd(result, "d1ph23"), [0.0, 286.6, 113.4], atol=1E-1) - @test isapprox(qd(result, "d1ph23"), [0.0, 34.5, 265.5], atol=1E-1) + @test isapprox(pd(result, "d1ph23"), [286.6, 113.4, 0.0], atol=1E-1) + @test isapprox(qd(result, "d1ph23"), [34.5, 265.5, 0.0], atol=1E-1) # single-phase wye loads @test isapprox(pd(result, "y1ph2"), [400.0], atol=1E-1) @test isapprox(qd(result, "y1ph2"), [300.0], atol=1E-1)