FIX: indexing of single-phase delta loads for linear formulations (#402)

* Fixed indexing issue

Update CHANGELOG.md

* Update bf_mx_lin.jl

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`

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

src/form/bf_fbs.jl

@@ -590,33 +590,46 @@ 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
 
         if report
             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

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

src/form/fotp.jl

@@ -490,35 +490,47 @@ 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
 
         if report
             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

src/form/fotr.jl

@@ -518,33 +518,46 @@ 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
 
         if report
             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

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)