Add creation of Y matrix from Network, multiphase breakout

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.0"
 
 [deps]
 CommonOPF = "dad7ea18-2b21-482e-81c1-e84414cc4b03"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 

src/BusInjectionModel.jl

@@ -4,7 +4,7 @@ module BusInjectionModel
 using CommonOPF
 using JuMP
 using LinearAlgebra
-
+using Graphs
 
 """
     ModelType
@@ -19,10 +19,13 @@ end
 
 export
     build_bim!,
-    FixedPointLinear
+    FixedPointLinear,
+    admittance_builder,
 
 
 include("single_phase_model.jl")
 include("single_phase_fpl_model.jl")
+include("multi_phase_fpl_model.jl")
+include("utilities.jl")
 
 end

src/multi_phase_fpl_model.jl

@@ -0,0 +1,166 @@
+
+
+"""
+    build_bim!(m::JuMP.AbstractModel, net::Network{MultiPhase}, ::Val{FixedPointLinear})
+
+Add variables and constraints to `m` using the values in `net` to make an FixedPointLinear branch flow
+model. Calls the following functions:
+- [`add_bim_variables`](@ref)
+- [`constrain_voltage`](@ref)
+```
+"""
+function build_bim!(m::JuMP.AbstractModel, net::Network{MultiPhase}, ::Val{FixedPointLinear})
+    add_bim_variables(m, net)
+    constrain_voltage(m, net)
+end
+
+
+function build_lpf(m::JuMP.AbstractModel, net::Network{MultiPhase})
+    add_variables(m, p)
+    constrain_voltage(m, p)
+    constrain_substation_power(m, p)
+    constrain_bounds(m, p)
+end
+
+
+function add_bim_variables(m::JuMP.AbstractModel, net::Network{MultiPhase})
+    N = length(busses(net))
+    T = net.Ntimesteps
+
+    @variable(m, P₀[1:T])
+    @variable(m, Q₀[1:T])
+    # TODO Pj Qj variables and change p.P/Qref to P/Qload dicts ?
+
+    @variable(m, v[1:N, 1:T] >=0)  # vlo and vhi applied in constraints s.t. they end up in C matrix
+end
+
+
+function constrain_voltage(m::JuMP.AbstractModel, net::Network{MultiPhase})
+    v = m[:v]
+    N = length(busses(net))
+    Y = admittance_builder(net)
+
+    # Bernstein & Dall'anese 2017 equ. 5b, 9b, 9c
+    for t = 1:net.Ntimesteps
+        # TODO should we update the M matrices with each time step? Depends on if Vref is updated with time
+        invDiagConjVref = inv(Diagonal(vec(conj(p.vref[:, t]))))  # TODO mv to Inputs so only calculated once (and Mwye)
+        Mwye = [p.invYLL * invDiagConjVref  -im * p.invYLL * invDiagConjVref]
+        Kwye = inv(Diagonal(abs.(p.vref[:,t]))) * real( Diagonal(conj(p.vref[:,t])) * Mwye )
+        b = abs.(p.vref[:,t]) - Kwye * vec([p.Pref[:,t] / p.Sbase; p.Qref[:,t]] / p.Sbase)
+
+        KL = Kwye[:, 1:N]
+        KR = Kwye[:, N+1:end]
+
+        @constraint(m, [j = 1:N],
+            v[j,t] == sum( KL[j,n] * p.Pref[n,t] / p.Sbase for n in N) 
+                    + sum( KR[j,n] * p.Qref[n,t] / p.Sbase for n in N) 
+                    + b[j]
+        )
+    end
+    #= 
+    equ (10) distributable voltage estimation, less accurate
+
+    v₀ = [complex(net.v0, 0)]
+    w = -p.invYLL * p.YL0 * v₀
+    W = Diagonal(w)
+    big_one = vec(ones(N))
+
+    abs.(W) * (big_one + real(inv(W) * Mwye) * vec([p.Pref[:,t] / p.Sbase; p.Qref[:,t]] / p.Sbase))
+    =# 
+    p.Nequality_cons += N * net.Ntimesteps
+    nothing
+end
+
+
+
+"""
+    constrain_substation_power(m::JuMP.AbstractModel, params::YPQVarraysNodeByTime)
+
+P₀ & Q₀ definitions ∀ t ∈ T
+"""
+function constrain_substation_power(m::JuMP.AbstractModel, net::Network{MultiPhase})
+    P₀ = m[:P₀]
+    Q₀ = m[:Q₀]
+    v₀ = [complex(net.v0, 0)]
+    # parameters using notation from Bernstein and Dall'anese 2017 equ.s 5c, 13a, 13b
+    w = -p.invYLL * p.YL0 * v₀  # can be replaced with ntwk.vnom?
+    c = Diagonal(v₀) * (conj(p.Y₀₀) * conj(v₀) + conj(p.Y0L) * conj(w))
+    creal = real(c)[1]
+    cimag = imag(c)[1]
+
+    # s₀ = G*x + c
+    for t = 1:net.Ntimesteps
+        invDiagConjVref = inv(Diagonal(vec(conj(p.vref[:, t]))))
+        Mwye = [p.invYLL * invDiagConjVref  -im * p.invYLL * invDiagConjVref]
+        gwye = Diagonal(v₀) * conj(p.Y0L) * conj(Mwye)
+
+        @constraint(m, 
+            P₀[t] == 1/p.Sbase * dot(real(gwye), [p.Pref[:,t]; p.Qref[:,t]]) + creal
+        )
+        @constraint(m, 
+            Q₀[t] == 1/p.Sbase * dot(imag(gwye), [p.Pref[:,t]; p.Qref[:,t]]) + cimag
+        )
+    end
+
+    p.Nequality_cons += 2 * net.Ntimesteps
+    nothing
+end
+
+
+function constrain_bounds(m::JuMP.AbstractModel, net::Network{MultiPhase})
+    N = length(busses(net))
+
+    @constraint(m, con_vhi[n=1:N, t=1:net.Ntimesteps],
+        m[:v][n, t] ≤ net.bounds.v_upper_mag
+    )
+    p.Nlte_cons += N * net.Ntimesteps
+
+    @constraint(m, con_vlo[n=1:N, t=1:net.Ntimesteps],
+        p.vlo - m[:v][n, t] ≤ 0
+    )
+    p.Nlte_cons += N * net.Ntimesteps
+
+    @constraint(m, con_P0hi[t=1:net.Ntimesteps],
+        m[:P₀][t] ≤ p.P0hi
+    )
+    p.Nlte_cons += net.Ntimesteps
+
+    @constraint(m, con_P0lo[t=1:net.Ntimesteps],
+        p.P0lo - m[:P₀][t] ≤ 0
+    )
+    p.Nlte_cons += net.Ntimesteps
+
+    @constraint(m, con_Q0hi[t=1:net.Ntimesteps],
+        m[:Q₀][t] ≤ p.Q0hi
+    )
+    p.Nlte_cons += net.Ntimesteps
+
+    @constraint(m, con_Q0lo[t=1:net.Ntimesteps],
+        p.Q0lo - m[:Q₀][t] ≤ 0
+    )
+    p.Nlte_cons += net.Ntimesteps
+    nothing
+end
+
+
+function check_existence_condition(m, net::Network{MultiPhase})
+    # TODO: still applies? if so update the P,Q vectors
+    N = length(busses(net))
+    upperbound = net.v0^2 / (4 * matrixnormstar(p.Z))
+    Snorm = zeros(Float64, net.Ntimesteps)
+    for t in 1:net.Ntimesteps
+        P = value.(m[:P][:,t])
+        Q = value.(m[:P][:,t])
+        Snorm[t] = LinearAlgebra.norm( sqrt(sum( P[j]^2 + Q[j]^2 for j in 0:N)) )
+    end
+    return all(Snorm[t] < upperbound for t in 1:net.Ntimesteps)
+end
+
+
+function matrixnormstar(A::AbstractArray{<:Complex{Float64}, 2})
+    rownorms = zeros(Float64, size(A, 1))
+    for r = 1:size(A,1)
+        rownorms[r] += LinearAlgebra.norm(abs.(A[r,:]))
+    end
+    return maximum(rownorms) 
+end

src/utilities.jl

@@ -0,0 +1,62 @@
+
+
+function admittance_builder(net::Network{MultiPhase})
+
+    order = bfs_order(net)
+    Y = zeros(ComplexF64, (3*length(order), 3*length(order)))
+
+    for edge in CommonOPF.edges(net)
+
+        z_prim = CommonOPF.zij(edge[1],edge[2],net)
+
+        phases = [1,2,3]
+        missing_phases = []
+        for i in [1,2,3]
+            if z_prim[i,:] == ComplexF64[0.0 + 0.0im for i = 1:3]                
+                filter!(!=(i),phases)
+                append!(missing_phases,i)
+            end
+        end
+
+        if length(phases) == 2
+            z_prim = z_prim[1:end .!= missing_phases[1], 1:end .!= missing_phases[1]]
+        elseif  length(phases) == 1
+            z_prim = z_prim[1:end .== phases[1], 1:end .== phases[1]]
+        end
+
+        y_prim = inv(z_prim)
+
+        position1 = findfirst(x -> x == edge[1],order)
+        position2 = findfirst(x -> x == edge[2],order)
+
+        for (row,phs1) in enumerate(phases)
+            for (col,phs2) in enumerate(phases)
+                Y[3*(position1-1)+phs1, 3*(position2-1)+phs2] += y_prim[row,col]
+                Y[3*(position2-1)+phs2, 3*(position1-1)+phs1] += y_prim[row,col]
+            end
+        end
+    end
+
+    non_zero_rows = any(Y .!= 0.0 + 0.0im, dims=2)
+    non_zero_indices = findall(non_zero_rows) .|> x -> x[1]
+    Y = Y[non_zero_indices, non_zero_indices]
+
+    return Y
+
+end
+
+function bfs_order(net::Network{MultiPhase})
+    b = Graphs.bfs_tree(net.graph, 1)
+    order = Any[1]
+    i = 1
+    while i <= length(order)
+        append!(order, neighbors(b, order[i]))
+        i += 1
+    end
+
+    for (i,num) in enumerate(order)
+        order[i] = net.graph.vertex_labels[num]
+    end
+
+    return order
+end