start transition to CommonOPF.Network

CHANGELOG.md

@@ -1,5 +1,11 @@
 # BranchFlowModel Changelog
 
+## dev
+
+### transition to CommonOPF.Network (from CommonOPF.Inputs)
+- add model_single_phase_network.jl and test build and solve
+    - this will replace model_single_phase.jl in next release
+
 ## v0.4.3
 - upgrade CommonOPF to v0.3.8
 

src/BranchFlowModel.jl

@@ -12,6 +12,15 @@ const DEFAULT_AMP_LIMIT = 1000.0
 
 export 
     Inputs,
+    # new CommonOPF methods
+    Network,
+    edges,
+    busses,
+    i_to_j,
+    j_to_k,
+    rij, 
+    xij,
+    # end new methods TODO distinguish all exports by module
     singlephase38linesInputs,
     dsstxt_to_sparse_array, 
     dss_files_to_dict,
@@ -29,10 +38,6 @@ export
     get_ijlinelength,
     get_ij_idx,
     # recover_voltage_current,  # TODO validate this method
-    i_to_j, 
-    j_to_k, 
-    rij, 
-    xij,
     zij,
     check_soc_inequalities,
     check_connected_graph,
@@ -72,6 +77,7 @@ include("utils.jl")
 include("results.jl")
 include("checks.jl")
 include("model_single_phase.jl")
+include("model_single_phase_network.jl")
 include("model_multi_phase.jl")
 include("decomposition.jl")
 

src/model_single_phase.jl

@@ -217,6 +217,8 @@ end
 - set net injections Pj/Qj to negative of Inputs.Pload/Qload, which are normalized by Sbase when creating Inputs
 - keys of P/Qload must match Inputs.busses. Any missing keys have load set to zero.
 - Inputs.substation_bus is unconstrained, slack bus
+
+TODO this method is same as LinDistFlow single phase: should it be moved to CommonOPF?
 """
 function constrain_loads(m, p::Inputs)
     Pj = m[:Pj]

src/model_single_phase_network.jl

@@ -0,0 +1,254 @@
+"""
+    build_model!(m::JuMP.AbstractModel, net::Network{SinglePhase})
+
+Add variables and constraints to `m` using the values in `net`. Calls the following functions:
+```julia
+add_variables(m, net)
+constrain_power_balance(m, net)
+constrain_substation_voltage(m, net)
+constrain_KVL(m, net)
+constrain_cone(m, net)
+constrain_loads(m, net)
+```
+"""
+function build_model!(m::JuMP.AbstractModel, net::Network{SinglePhase}; relaxed::Bool=true)
+
+    add_variables(m, net)
+    constrain_power_balance(m, net)
+    constrain_substation_voltage(m, net)
+    constrain_KVL(m, net)
+    constrain_loads(m, net)
+    if relaxed
+        constrain_cone(m, net)
+    else
+        constrain_psd(m, net)
+    end
+end
+
+
+function add_variables(m, net::Network)
+    T = 1:net.Ntimesteps
+    # bus injections are expressions, defined in constrain_power_balance
+    # @variables m begin
+    #     p.P_lo_bound <= Pj[p.busses, T] <= p.P_up_bound
+    #     p.Q_lo_bound <= Qj[p.busses, T] <= p.Q_up_bound
+    # end
+    # TODO warn when applying power injection lower bounds and SDP, radial b/c for radial
+    # networks with power injections unbounded below (voltage angle relaxation) the SOCP (more
+    # efficiently) yields the same solution as the SDP and the (non-convex) exact equations.
+
+    # voltage squared
+   @variable(m, vsqrd[busses(net), T] >= (net.v_lolim)^2) #<= p.v_uplim^2 ) 
+
+    # line flows, netower sent from i to j
+    # @variable(m, net.P_lo_bound <= Pij[edges(net), T] <= p.P_up_bound )
+    # @variable(m, net.Q_lo_bound <= Qij[edges(net), T] <= p.Q_up_bound )
+    @variable(m, Pij[edges(net), T])
+    @variable(m, Qij[edges(net), T])
+
+    # current squared (non-negative)
+    @variable(m, lij[edges(net), T] >= 0)
+
+    # @constraint(m, [edge in net.edge_keys, t in T],
+    #     lij[edge, t] <= net.amps_limit[edge]
+    # )
+
+    nothing
+end
+
+
+"""
+    function constrain_power_balance(m, net::Network)
+
+Pij in - losses == sum of line flows out + net injection
+NOTE: using sum over Pij for future expansion to mesh grids
+i -> j -> k
+"""
+function constrain_power_balance(m, net::Network)
+    Pij = m[:Pij]
+    Qij = m[:Qij]
+    lij = m[:lij]
+    Pj = m[:Pj] = Dict()  # we fill these in with expression indexed on bus and time
+    Qj = m[:Qj] = 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
+    # AND? using expressions for net injections means that we do not have to delete and redifine
+    # constraints? i.e. can just define another expression for Pj when it is a function of a DER ?
+    for j in busses(net)
+
+        # source nodes, injection = flows out
+        if isempty(i_to_j(j, net)) && !isempty(j_to_k(j, net))
+            substation_Pload = zeros(net.Ntimesteps)
+            if j in real_load_busses(net)
+                substation_Pload = net[j][:Load][:kws1] / net.Sbase
+            end
+            Pj[j] = @expression(m, [t = 1:net.Ntimesteps],
+                - substation_Pload[t] - sum( Pij[(j,k), t] for k in j_to_k(j, net) )
+            )
+            substation_Qload = zeros(net.Ntimesteps)
+            if j in reactive_load_busses(net)
+                substation_Qload = net[j][:Load][:kvars1] / net.Sbase
+            end
+            Qj[j] = @expression(m, [t = 1:net.Ntimesteps],
+                - substation_Qload[t] - sum( Qij[(j,k), t] for k in j_to_k(j, net) )
+            )
+
+        # unconnected nodes
+        elseif isempty(i_to_j(j, net)) && isempty(j_to_k(j, net))
+            @warn "Bus $j has no edges, setting Pj and Qj to zero."
+            Pj[j] = zeros(net.Ntimesteps)
+            Qj[j] = zeros(net.Ntimesteps)
+
+        # leaf nodes / sinks, flows in = draw out
+        elseif !isempty(i_to_j(j, net)) && isempty(j_to_k(j, net))
+            Pj[j] = @expression(m, [t = 1:net.Ntimesteps],
+                -sum( Pij[(i,j), t] for i in i_to_j(j, net) ) +
+                sum( lij[(i,j), t] * rij(i,j,net) for i in i_to_j(j, net) ) 
+            )
+            Qj[j] = @expression(m, [t = 1:net.Ntimesteps],
+                -sum( Qij[(i,j), t] for i in i_to_j(j, net) ) +
+                sum( lij[(i,j), t] * xij(i,j,net) for i in i_to_j(j, net) )
+                # + p.shunt_susceptance[j] * m[:vsqrd][j,t]
+            )
+
+        # intermediate nodes
+        else
+            Pj[j] = @expression(m, [t = 1:net.Ntimesteps],
+                -sum( Pij[(i,j), t] for i in i_to_j(j, net) )
+                + sum( lij[(i,j), t] * rij(i,j,net) for i in i_to_j(j, net) ) 
+                + sum( Pij[(j,k), t] for k in j_to_k(j, net) )
+            )
+            Qj[j] = @expression(m, [t = 1:net.Ntimesteps],
+                -sum( Qij[(i,j), t] for i in i_to_j(j, net) ) 
+                + sum( lij[(i,j), t] * xij(i,j,net) for i in i_to_j(j, net) )
+                + sum( Qij[(j,k), t] for k in j_to_k(j, net) ) 
+                # + p.shunt_susceptance[j] * m[:vsqrd][j,t]
+            )
+        end
+    end
+    # TODO Farivar and Low have b*v and q*v in these equations for shunts? neglected for now
+
+    nothing
+end
+
+
+function constrain_substation_voltage(m, net::Network)
+    if typeof(net.v0) <: Real
+        @constraint(m, con_substationV[t = 1:net.Ntimesteps],
+            m[:vsqrd][net.substation_bus, t] == net.v0^2
+        )
+    else  # vector of time
+        @constraint(m, con_substationV[t = 1:net.Ntimesteps],
+            m[:vsqrd][net.substation_bus, t] == net.v0[t]^2
+        )
+    end
+    nothing
+end
+
+
+function constrain_KVL(m, net::Network)
+    w = m[:vsqrd]
+    P = m[:Pij]
+    Q = m[:Qij]
+    l = m[:lij]
+    m[:vcons] = Dict()
+    for j in busses(net)
+        for i in i_to_j(j, net)  # for radial network there is only one i in i_to_j
+            # if !( (i,j) in keys(p.regulators) )
+                rᵢⱼ = rij(i,j,net)
+                xᵢⱼ = xij(i,j,net)
+                m[:vcons][j] = @constraint(m, [t = 1:net.Ntimesteps],
+                    w[j,t] == w[i,t]
+                        - 2*(rᵢⱼ * P[(i,j),t] + xᵢⱼ * Q[(i,j),t])
+                        + (rᵢⱼ^2 + xᵢⱼ^2) * l[(i,j), t]
+                )
+            # else
+            #     if has_vreg(p, j)
+            #         m[:vcons][j] = @constraint(m, [t = 1:net.Ntimesteps],
+            #             w[j,t] == p.regulators[(i,j)][:vreg]^2
+            #         )
+            #     else  # default turn_ratio is 1.0
+            #         m[:vcons][j] = @constraint(m, [t = 1:net.Ntimesteps],
+            #             w[j,t] == w[i,t] * p.regulators[(i,j)][:turn_ratio]^2 
+            #         )
+            #     end
+            # end
+        end
+    end
+    nothing
+end
+
+
+function constrain_cone(m, net::Network)
+    w = m[:vsqrd]
+    P = m[:Pij]
+    Q = m[:Qij]
+    l = m[:lij]
+    for j in busses(net)
+        for i in i_to_j(j, net)  # for radial network there is only one i in i_to_j
+            # # equivalent but maybe not handled as well in JuMP ?
+            # @constraint(m, [t = 1:net.Ntimesteps],
+            #     w[i,t] * l[(i,j), t] ≥ P[(i,j),t]^2 + Q[(i,j),t]^2
+            # )
+            @constraint(m, [t = 1:net.Ntimesteps], 
+                [w[i,t]/2, l[(i,j), t], P[(i,j),t], Q[(i,j),t]] in JuMP.RotatedSecondOrderCone()
+            )
+        end
+    end
+end
+
+
+function constrain_psd(m, net::Network)
+    w = m[:vsqrd]
+    P = m[:Pij]
+    Q = m[:Qij]
+    l = m[:lij]
+    for j in busses(net)
+        for i in i_to_j(j, net)  # for radial network there is only one i in i_to_j
+            # need to use PSD?
+            @constraint(m, [t = 1:net.Ntimesteps],
+                w[i,t] * l[(i,j), t] == P[(i,j),t]^2 + Q[(i,j),t]^2
+            )
+        end
+    end
+end
+
+
+"""
+    constrain_loads(m, net::Network)
+
+- set net injections Pj/Qj to negative of Inputs.Pload/Qload, which are normalized by Sbase when creating Inputs
+- keys of P/Qload must match Inputs.busses. Any missing keys have load set to zero.
+- Inputs.substation_bus is unconstrained, slack bus
+
+TODO this method is same as LinDistFlow single phase: should it be moved to CommonOPF?
+"""
+function constrain_loads(m, net::Network)
+    Pj = m[:Pj]
+    Qj = m[:Qj]
+    m[:injectioncons] = Dict()
+    for j in setdiff(busses(net), [net.substation_bus])
+        m[:injectioncons][j] = Dict()
+        if j in real_load_busses(net)
+            con = @constraint(m, [t = 1:net.Ntimesteps],
+                Pj[j][t] == -net[j][:Load][:kws1][t] * 1e3 / net.Sbase
+            )
+        else
+            con = @constraint(m, [t = 1:net.Ntimesteps],
+                Pj[j][t] == 0
+            )
+        end
+        m[:injectioncons][j]["p"] = con
+        if j in reactive_load_busses(net)
+            con = @constraint(m, [t = 1:net.Ntimesteps],
+                Qj[j][t] == -net[j][:Load][:kvars1][t] * 1e3 / net.Sbase
+            )
+        else
+            con = @constraint(m, [t = 1:net.Ntimesteps],
+                Qj[j][t] == 0
+            )
+        end
+        m[:injectioncons][j]["q"] = con
+    end
+    nothing
+end

test/runtests.jl

@@ -9,7 +9,7 @@ using LinearAlgebra
 using COSMO
 using CSDP
 using Ipopt
-using Graphs, MetaGraphs  # maybe export what is needed from MetaGraphs in BranchFlowModel ?
+using Graphs, MetaGraphs  # TODO export what is needed from MetaGraphs in BranchFlowModel
 
 # # hack for local testing
 # using Pkg
@@ -43,6 +43,42 @@ end
 
 @testset "BranchFlowModel.jl" begin
 
+
+@testset "CommonOPF.Network" begin
+
+    # dss("clear")
+    # dss("Redirect data/ieee13_makePosSeq/Master.dss")
+    # dss("Solve")
+    # @test(OpenDSSDirect.Solution.Converged() == true)
+    # dss_voltages = dss_voltages_pu()
+
+
+    fp = joinpath("data", "ieee13", "ieee13_single_phase.yaml")
+    net = Network(fp)
+    # net.v0 = dss_voltages["rg60"][1]
+    m = Model(Ipopt.Optimizer)
+    build_model!(m, net; relaxed=false)
+    @objective(m, Min, 
+        sum( m[:lij][edge,t] for t in 1:net.Ntimesteps, edge in edges(net))
+    )
+    optimize!(m)
+
+    vs = get_variable_values(:vsqrd, m, net)
+
+    # NEXT the opendss model has a lot more in it. rather than catch up to opendss to validate use
+    # https://jso.dev/NLPModelsJuMP.jl/dev/tutorial/ to newton-solve the SDP (and check rank 1?) and
+    # compare voltages between the optimal model and the newton solution
+    # for bus in busses(net)
+    #     if bus in keys(vs)
+    #         println(rpad(bus,5), 
+    #             rpad(round(vs[bus][1] -  dss_voltages[bus][1], digits=4), 8)
+    #         )
+    #     end
+    # end
+
+end
+
+
 @testset "Papavasiliou 2018 with shunts" begin
     #=
     Copied network from paper "Analysis of DLMPs"
@@ -569,7 +605,7 @@ end
     );
     @test check_connected_graph(p) == true
     g = make_graph(p.busses, p.edges)
-    @test g.graph.ne == length(g.vprops) - 1  # a radial network has n_edges = n_vertices - 1
+    @test g.graph.ne == Graphs.nv(g) - 1  # a radial network has n_edges = n_vertices - 1
     @test_warn "The per unit impedance values should be much less than one" check_unique_solution_conditions(p)
     p.regulators[("650", "rg60")][:turn_ratio] = dss_voltages["rg60"][1]