ADD: functions and capabilitites that automatically initialize the bo…

…undary power flow variables. (#15) - Added functions and capabilitites that automatically initialize the boundary power flow variables (i.e., `pbound_fr`, `pbound_to`, `qbound_fr`, and `qbound_to`). - `p/qbound_fr` variables are initialized based on the data parsed from the transmission system data. - `p/qbound_to` variables are initialized by adding all the loads in the respective distribution systems and dividing them by 3 (assuming the initial state of the distribution system is balanced). - Refactored multiple functions in `ref.jl` to avoid unnecessary loops. These unnecessary loops were making the refs processes unnecessary long.
lanl-ansi · Jun 14, 2023 · 0dd3c42 · 0dd3c42
1 parent bdece4d
commit 0dd3c42
Show file tree

Hide file tree

Showing 13 changed files with 151 additions and 85 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -4,6 +4,10 @@
 
 - Fixed implementation of polynomial nl costs above quadratic in `objective.jl`.
 - Bumped PMITD compatibility of `IM`, `PMD` and `PM` to the latest versions (i.e., V0.7.7, v0.14.9, and v0.19.9).
+- Added functions and capabilitites that automatically initialize the boundary power flow variables (i.e., `pbound_fr`, `pbound_to`, `qbound_fr`, and `qbound_to`).
+- `p/qbound_fr` variables are initialized based on the data parsed from the tranmission system data.
+- `p/qbound_to` variables are initialized by adding all the loads in the respective distribution systems and dividing them by 3 (assuming the initial state of the distribution system is balanced).
+- Refactored multiple functions in `ref.jl` to avoid unnecesary loops. These unnecesary loops were making the refs processes unccessary long.
 
 ## v0.7.7
 

diff --git a/src/core/ref.jl b/src/core/ref.jl
@@ -46,15 +46,19 @@ end
 Removes/filters-out the loads at buses (i.e., boundary buses) where distribution systems are going to be integrated/connected.
 """
 function _ref_filter_transmission_integration_loads!(ref::Dict{Symbol,<:Any})
-    # Loops over all T-D pmitd available
-    for (nw_it, nw_ref_it) in ref[:it][:pmitd][:nw]
+    # Loops over all nws
+    for (nw, nw_ref) in ref[:it][:pm][:nw]
         # Filters only the ones that have the "transmission_boundary" key
-        for (i, conn) in filter(x -> "transmission_boundary" in keys(x.second), nw_ref_it)
-            # Filters out only the loads connected to the transmission-boundary bus
-            for (nw, nw_ref) in ref[:it][:pm][:nw]
-                nw_ref[:load] = Dict(x for x in nw_ref[:load] if x.second["load_bus"] != conn["transmission_boundary"] )
-                nw_ref[:bus_loads][conn["transmission_boundary"]] = []
-            end
+        for (i, conn) in filter(x -> "transmission_boundary" in keys(x.second), ref[:it][:pmitd][:nw][nw])
+            # Get init (start) values before deleting the boundary load info.
+            pbound_fr_start = nw_ref[:load][nw_ref[:bus_loads][conn["transmission_boundary"]][1]]["pd"]
+            qbound_fr_start = nw_ref[:load][nw_ref[:bus_loads][conn["transmission_boundary"]][1]]["qd"]
+            conn["pbound_fr_start"] = pbound_fr_start
+            conn["qbound_fr_start"] = qbound_fr_start
+
+            # Remove loads
+            nw_ref[:load] = Dict(x for x in nw_ref[:load] if x.second["load_bus"] != conn["transmission_boundary"] )
+            nw_ref[:bus_loads][conn["transmission_boundary"]] = []
         end
     end
 end
@@ -68,29 +72,27 @@ end
 Removes/filters-out the slack generators at buses/nodes where the transmission system is going to be integrated/connected.
 """
 function _ref_filter_distribution_slack_generators!(ref::Dict{Symbol,<:Any})
-     # Loops over all T-D pmitd available
-    for (nw_it, nw_ref_it) in ref[:it][:pmitd][:nw]
+    # Loops over all nws
+    for (nw, nw_ref) in ref[:it][:pmd][:nw]
         # Filters only the ones that have the "distribution_boundary" key
-        for (i, conn) in filter(x -> "distribution_boundary" in keys(x.second), nw_ref_it)
+        for (i, conn) in filter(x -> "distribution_boundary" in keys(x.second), ref[:it][:pmitd][:nw][nw])
             # Filters out only the gens connected to the distribution-boundary bus (virtual slack bus from opendss)
-            for (nw, nw_ref) in ref[:it][:pmd][:nw]
-                nw_ref[:gen] = Dict(x for x in nw_ref[:gen] if x.second["gen_bus"] != conn["distribution_boundary"] )
-                nw_ref[:bus_conns_gen][conn["distribution_boundary"]] = []
-                nw_ref[:bus_gens][conn["distribution_boundary"]] = []
-                # Unrestrict buspairs connected to reference bus
-                for (j, bus_pair) in nw_ref[:buspairs]
-                    if (bus_pair["vm_fr_min"] == 1.0) # only need to check one
-                        bus_pair["vm_fr_min"] = 0.0
-                        bus_pair["vm_fr_max"] = Inf
-                    end
+            nw_ref[:gen] = Dict(x for x in nw_ref[:gen] if x.second["gen_bus"] != conn["distribution_boundary"] )
+            nw_ref[:bus_conns_gen][conn["distribution_boundary"]] = []
+            nw_ref[:bus_gens][conn["distribution_boundary"]] = []
+            # Unrestrict buspairs connected to reference bus
+            for (j, bus_pair) in nw_ref[:buspairs]
+                if (bus_pair["vm_fr_min"] == 1.0) # only need to check one
+                    bus_pair["vm_fr_min"] = 0.0
+                    bus_pair["vm_fr_max"] = Inf
                 end
-                # Modify v_min and v_max, remove va and vm, and change bus type for reference bus
-                nw_ref[:bus][conn["distribution_boundary"]]["vmin"] = [0.0, 0.0, 0.0]
-                nw_ref[:bus][conn["distribution_boundary"]]["vmax"] = [Inf, Inf, Inf]
-                # nw_ref[:bus][conn["distribution_boundary"]]["va"] = [0.0, -120.0, 120.0]
-                # nw_ref[:bus][conn["distribution_boundary"]]["vm"] = [1.0,1.0,1.0]
-                nw_ref[:bus][conn["distribution_boundary"]]["bus_type"] = 1
             end
+            # Modify v_min and v_max, remove va and vm, and change bus type for reference bus
+            nw_ref[:bus][conn["distribution_boundary"]]["vmin"] = [0.0, 0.0, 0.0]
+            nw_ref[:bus][conn["distribution_boundary"]]["vmax"] = [Inf, Inf, Inf]
+            # nw_ref[:bus][conn["distribution_boundary"]]["va"] = [0.0, -120.0, 120.0]
+            # nw_ref[:bus][conn["distribution_boundary"]]["vm"] = [1.0,1.0,1.0]
+            nw_ref[:bus][conn["distribution_boundary"]]["bus_type"] = 1
         end
     end
 end
@@ -104,43 +106,56 @@ end
 Creates the boundary `refs` that integrate/connect the transmission and distribution system bus(es).
 """
 function _ref_connect_transmission_distribution!(ref::Dict{Symbol,<:Any})
-    # Loops over all T-D pmitd available
-    for (nw_it, nw_ref_it) in ref[:it][:pmitd][:nw]
-
-        for i in 1:length(nw_ref_it)   # loop through all boundary objects
-            boundary_number = BOUNDARY_NUMBER - 1 + i # boundary number index
-
-            for (nw, nw_ref) in ref[:it][:pmd][:nw]
-                # create :boundary structure if does not exists; inserts to dictionary if it already exists
-                if !haskey(nw_ref_it, :boundary)
-                    nw_ref_it[:boundary] = Dict(boundary_number => Dict("f_bus" => 0, "t_bus" => 0, "index" => 0, "name" => "empty", "f_connections" => [1], "t_connections" => [1, 2, 3]))
-                else
-                    nw_ref_it[:boundary][boundary_number] = Dict("f_bus" => 0, "t_bus" => 0, "index" => 0, "name" => "empty", "f_connections" => [1], "t_connections" => [1, 2, 3])
-                end
+    # Loops over all nws
+    for (nw, nw_ref) in ref[:it][:pmd][:nw]
+        # get the specific nw pmitd data
+        nw_ref_it = ref[:it][:pmitd][:nw][nw]
+        # Loops over all boundary objects
+        for i in 1:length(nw_ref_it)
+            # boundary number index
+            boundary_number = BOUNDARY_NUMBER - 1 + i
 
-                # modify default values with actual values coming from linking file information
-                nw_ref_it[:boundary][boundary_number]["f_bus"] = nw_ref_it[Symbol(boundary_number)]["transmission_boundary"]
-                nw_ref_it[:boundary][boundary_number]["t_bus"] = nw_ref_it[Symbol(boundary_number)]["distribution_boundary"]
-                nw_ref_it[:boundary][boundary_number]["index"] = boundary_number
-                nw_ref_it[:boundary][boundary_number]["name"] = "_itd_boundary_$boundary_number"
-
-                # Add bus reference from transmission (pm)
-                # The dictionary represents Dict(original bus_index => boundary # that belongs to)
-                trans_bus = nw_ref_it[Symbol(boundary_number)]["transmission_boundary"]
-                if !haskey(nw_ref_it, :bus_from)
-                    nw_ref_it[:bus_from] = Dict(trans_bus => Dict("boundary" => boundary_number))
-                else
-                    nw_ref_it[:bus_from][trans_bus] = Dict("boundary" => boundary_number)
-                end
+            # create :boundary structure if does not exists; inserts to dictionary if it already exists
+            if !haskey(nw_ref_it, :boundary)
+                nw_ref_it[:boundary] = Dict(boundary_number => Dict("f_bus" => 0, "t_bus" => 0, "index" => 0, "name" => "empty", "f_connections" => [1], "t_connections" => [1, 2, 3], "pbound_fr_start" => 0, "qbound_fr_start" => 0, "pbound_to_start" => 0, "qbound_to_start" => 0))
+            else
+                nw_ref_it[:boundary][boundary_number] = Dict("f_bus" => 0, "t_bus" => 0, "index" => 0, "name" => "empty", "f_connections" => [1], "t_connections" => [1, 2, 3], "pbound_fr_start" => 0, "qbound_fr_start" => 0, "pbound_to_start" => 0, "qbound_to_start" => 0)
+            end
 
-                # Add bus reference from distribution (pmd)
-                # The dictionary represents Dict(original bus_index => boundary # that belongs to)
-                dist_bus = nw_ref_it[Symbol(boundary_number)]["distribution_boundary"]
-                if !haskey(nw_ref_it, :bus_to)
-                    nw_ref_it[:bus_to] = Dict(dist_bus => Dict("boundary" => boundary_number))
-                else
-                    nw_ref_it[:bus_to][dist_bus] = Dict("boundary" => boundary_number)
-                end
+            # modify default values with actual values coming from linking file information
+            nw_ref_it[:boundary][boundary_number]["f_bus"] = nw_ref_it[Symbol(boundary_number)]["transmission_boundary"]
+            nw_ref_it[:boundary][boundary_number]["t_bus"] = nw_ref_it[Symbol(boundary_number)]["distribution_boundary"]
+            nw_ref_it[:boundary][boundary_number]["index"] = boundary_number
+            nw_ref_it[:boundary][boundary_number]["name"] = "_itd_boundary_$boundary_number"
+            nw_ref_it[:boundary][boundary_number]["pbound_fr_start"] = nw_ref_it[Symbol(boundary_number)]["pbound_fr_start"]
+            nw_ref_it[:boundary][boundary_number]["qbound_fr_start"] = nw_ref_it[Symbol(boundary_number)]["qbound_fr_start"]
+
+            # Compute pbound_to and qbound_to start values for specific nw
+            pload_totals = _compute_boundary_active_power_start_values_distribution(nw_ref)
+            qload_totals = _compute_boundary_reactive_power_start_values_distribution(nw_ref)
+            # Get the ckt_name related to the boundary number
+            source_id = nw_ref[:bus][nw_ref_it[:boundary][boundary_number]["t_bus"]]["source_id"]
+            ckt_name = split(source_id, ".")[2]
+            # Assumes balance power initiliazation
+            nw_ref_it[:boundary][boundary_number]["pbound_to_start"] = pload_totals[ckt_name][1]/3
+            nw_ref_it[:boundary][boundary_number]["qbound_to_start"] = qload_totals[ckt_name][1]/3
+
+            # Add bus reference from transmission (pm)
+            # The dictionary represents Dict(original bus_index => boundary # that belongs to)
+            trans_bus = nw_ref_it[Symbol(boundary_number)]["transmission_boundary"]
+            if !haskey(nw_ref_it, :bus_from)
+                nw_ref_it[:bus_from] = Dict(trans_bus => Dict("boundary" => boundary_number))
+            else
+                nw_ref_it[:bus_from][trans_bus] = Dict("boundary" => boundary_number)
+            end
+
+            # Add bus reference from distribution (pmd)
+            # The dictionary represents Dict(original bus_index => boundary # that belongs to)
+            dist_bus = nw_ref_it[Symbol(boundary_number)]["distribution_boundary"]
+            if !haskey(nw_ref_it, :bus_to)
+                nw_ref_it[:bus_to] = Dict(dist_bus => Dict("boundary" => boundary_number))
+            else
+                nw_ref_it[:bus_to][dist_bus] = Dict("boundary" => boundary_number)
             end
 
             # :arcs_boundary_from for boundary
@@ -194,3 +209,63 @@ function _ref_remove_refbus_distribution!(ref::Dict{Symbol,<:Any})
         nw_ref[:ref_buses] = Dict{Int,Any}()
     end
 end
+
+
+"""
+    function _compute_boundary_active_power_start_values_distribution(
+        nw_ref::Dict{Symbol,<:Any}
+    )
+
+Computes the starting values for `pbound_to` variables. Returns dictionary with summation of the activate power loads.
+Returns dictionary with the summation of the active power loads for each dist. system.
+"""
+function _compute_boundary_active_power_start_values_distribution(nw_ref::Dict{Symbol,<:Any})
+
+    # Dicts to store summation of total load in dist. system
+    pload_totals = Dict()
+
+    # loop through all loads to add them up
+    for (_, load_info) in nw_ref[:load]
+        load_name = load_info["name"]
+        ckt_name = split(load_name, ".")
+        pd = load_info["pd"]
+
+        if !haskey(pload_totals, ckt_name[1])
+            pload_totals[ckt_name[1]] = sum(pd)
+        else
+            pload_totals[ckt_name[1]] += sum(pd)
+        end
+    end
+
+    return pload_totals
+end
+
+
+"""
+    function _compute_boundary_reactive_power_start_values_distribution(
+        nw_ref::Dict{Symbol,<:Any}
+    )
+
+Computes the starting values for `qbound_to` and adds them to `ref`.
+Returns dictionary with the summation of the reactive power loads for each dist. system.
+"""
+function _compute_boundary_reactive_power_start_values_distribution(nw_ref::Dict{Symbol,<:Any})
+
+    # Dicts to store summation of total load in dist. system
+    qload_totals = Dict()
+
+    # loop through all loads to add them up
+    for (_, load_info) in nw_ref[:load]
+        load_name = load_info["name"]
+        ckt_name = split(load_name, ".")
+        qd = load_info["qd"]
+
+        if !haskey(qload_totals, ckt_name[1])
+            qload_totals[ckt_name[1]] = sum(qd)
+        else
+            qload_totals[ckt_name[1]] += sum(qd)
+        end
+    end
+
+    return qload_totals
+end
diff --git a/src/core/variable.jl b/src/core/variable.jl
@@ -53,7 +53,7 @@ function variable_boundary_power_real_to(pm::AbstractPowerModelITD; nw::Int=nw_i
     # creating the variables
     pbound_to = Dict{Any,Any}((l,j,i) => JuMP.@variable(pm.model,
             [c in connections[(l,j,i)]], base_name="$(nw)_pbound_to_$((l,j,i))",
-            start = _PMD.comp_start_value(ref(pm, nw, :boundary, l), "pbound_to_start", c, 0.0)
+            start = _PMD.comp_start_value(ref(pm, nw, :boundary, l), "pbound_to_start", c)
         ) for (l,j,i) in ref(pm, nw, :arcs_boundary_to)
     )
 
@@ -93,7 +93,7 @@ function variable_boundary_power_imaginary_to(pm::AbstractPowerModelITD; nw::Int
     # creating the variables
     qbound_to = Dict{Any,Any}((l,j,i) => JuMP.@variable(pm.model,
             [c in connections[(l,j,i)]], base_name="$(nw)_qbound_to_$((l,j,i))",
-            start = _PMD.comp_start_value(ref(pm, nw, :boundary, l), "qbound_to_start", c, 0.0)
+            start = _PMD.comp_start_value(ref(pm, nw, :boundary, l), "qbound_to_start", c)
         ) for (l,j,i) in ref(pm, nw, :arcs_boundary_to)
     )
 

diff --git a/test/data/transmission/pglib_opf_case500_goc.m b/test/data/transmission/pglib_opf_case500_goc.m
@@ -1733,15 +1733,15 @@
 % INFO    : === Translation Options ===
 % INFO    : Phase Angle Bound:           30.0 (deg.)
 % INFO    : Setting Flat Start
-% INFO    : 
+% INFO    :
 % INFO    : === Generator Bounds Update Notes ===
-% INFO    : 
+% INFO    :
 % INFO    : === Base KV Replacement Notes ===
-% INFO    : 
+% INFO    :
 % INFO    : === Transformer Setting Replacement Notes ===
-% INFO    : 
+% INFO    :
 % INFO    : === Line Capacity Monotonicity Notes ===
-% INFO    : 
+% INFO    :
 % INFO    : === Voltage Setpoint Replacement Notes ===
 % INFO    : Bus 1	: V=0.98540926462, theta=-8.36784152229 -> V=1.0, theta=0.0
 % INFO    : Bus 2	: V=0.974020345007, theta=-11.1447179681 -> V=1.0, theta=0.0
@@ -2243,7 +2243,7 @@
 % INFO    : Bus 498	: V=1.00866690287, theta=-2.66209509196 -> V=1.0, theta=0.0
 % INFO    : Bus 499	: V=1.00340669673, theta=-3.59395611678 -> V=1.0, theta=0.0
 % INFO    : Bus 500	: V=1.003371832, theta=-3.59691858681 -> V=1.0, theta=0.0
-% INFO    : 
+% INFO    :
 % INFO    : === Generator Setpoint Replacement Notes ===
 % INFO    : Gen at bus 272	: Pg=46.8993571031, Qg=18.0091132505 -> Pg=34.463, Qg=7.856
 % INFO    : Gen at bus 272	: Vg=1.0343582 -> Vg=1.0
@@ -2640,5 +2640,5 @@
 % INFO    : Gen at bus 499	: Vg=1.039689 -> Vg=1.0
 % INFO    : Gen at bus 499	: Pg=1.73259069124, Qg=0.665908015041 -> Pg=1.2695, Qg=0.2905
 % INFO    : Gen at bus 499	: Vg=1.039689 -> Vg=1.0
-% INFO    : 
+% INFO    :
 % INFO    : === Writing Matpower Case File Notes ===
diff --git a/test/io.jl b/test/io.jl
@@ -22,7 +22,6 @@
         @test haskey(pmitd_data["it"], _PMD.pmd_it_name)
     end
 
-
     @testset "parse_link_file (invalid extension)" begin
         path = joinpath(dirname(dist_path), "case3_balanced.raw")
         @test_throws ErrorException parse_link_file(path)

diff --git a/test/opfitd.jl b/test/opfitd.jl
@@ -18,7 +18,6 @@
         @test result["termination_status"] == LOCALLY_SOLVED
     end
 
-
     @testset "solve_model (with network inputs): Balanced case5-case3 ACR-ACR with polynomial nl terms" begin
         pm_file = joinpath(dirname(trans_path), "case5_withload.m")
         pmd_file = joinpath(dirname(dist_path), "case3_balanced.dss")
@@ -128,7 +127,6 @@
         @test isapprox(result["objective"], 18005.97151; atol = 1e-4)
     end
 
-
     @testset "solve_model (with network inputs): Balanced case5-case13 ACR-ACR Without Dist. Generator ACR-ACR" begin
         pm_file = joinpath(dirname(trans_path), "case5_withload_ieee13.m")
         pmd_file = joinpath(dirname(dist_path), "caseIEEE13_balanced_withoutgen.dss")
@@ -294,7 +292,6 @@
         @test isapprox(result["objective"], 18005.97151; atol = 1e-4)
     end
 
-
     @testset "solve_model (with network inputs): Balanced case5-case3 IVR-IVR with polynomial nl terms" begin
         pm_file = joinpath(dirname(trans_path), "case5_withload.m")
         pmd_file = joinpath(dirname(dist_path), "case3_balanced.dss")
@@ -320,7 +317,7 @@
         @test result["termination_status"] == LOCALLY_SOLVED
     end
 
-     @testset "solve_model (with network inputs): Balanced case5-case3 IVR-IVR with piecewise linear terms" begin
+    @testset "solve_model (with network inputs): Balanced case5-case3 IVR-IVR with piecewise linear terms" begin
         pm_file = joinpath(dirname(trans_path), "case5_pwlc_withload.m")
         pmd_file = joinpath(dirname(dist_path), "case3_balanced.dss")
         pmitd_file = joinpath(dirname(bound_path), "case5_case3_bal.json")

diff --git a/test/opfitd_hybrids.jl b/test/opfitd_hybrids.jl
@@ -42,7 +42,6 @@
         @test result["termination_status"] == LOCALLY_SOLVED
     end
 
-
     @testset "solve_model (with network inputs): Unbalanced case5-case3 Without Dist. Generator BFA-LinDist3FlowPowerModel" begin
         pm_file = joinpath(dirname(trans_path), "case5_withload.m")
         pmd_file = joinpath(dirname(dist_path), "case3_unbalanced_withoutgen.dss")
@@ -63,7 +62,6 @@
         @test result["termination_status"] == OPTIMAL
     end
 
-
     @testset "solve_model (with network inputs): Balanced case5-case3 Without Dist. Generator SDPWRM-SOCConicUBF" begin
         pm_file = joinpath(dirname(trans_path), "case5_withload.m")
         pmd_file = joinpath(dirname(dist_path), "case3_balanced_notransformer_withoutgen.dss")

diff --git a/test/opfitd_mn.jl b/test/opfitd_mn.jl
@@ -32,7 +32,6 @@
         @test isapprox(result["objective"], 58399.9993; atol = 1e-4)
     end
 
-
     @testset "solve_mn_opfitd: Multinetwork case5-case3 x2 Without Dist. Generator ACP-ACP" begin
         pm_file = joinpath(dirname(trans_path), "case5_with2loads.m")
         pmd_file1 = joinpath(dirname(dist_path), "case3_unbalanced_withoutgen_mn.dss")

diff --git a/test/opfitd_ms.jl b/test/opfitd_ms.jl
@@ -101,7 +101,6 @@
         @test all(isapprox.(result["solution"]["it"]["pmd"]["bus"]["3bus_bal_nogen.primary"]["va"][2], -121.0428; atol=1e-3))
     end
 
-
     @testset "solve_model opfitd (with network inputs): Multi-System Balanced case5-case3x2 With PV ACP-ACP" begin
         pm_file = joinpath(dirname(trans_path), "case5_with2loads.m")
         pmd_file1 = joinpath(dirname(dist_path), "case3_balanced_withPV.dss")

diff --git a/test/opfitd_solution.jl b/test/opfitd_solution.jl
@@ -15,7 +15,6 @@
         @test all(isapprox.(result["solution"]["it"]["pmd"]["bus"]["3bus_unbal.primary"]["vm"][1], 0.9352; atol=1e-3))
     end
 
-
     @testset "solve_model (with network inputs): Balanced case5-case3 Without Dist. Generator IVR-IVR" begin
         pm_file = joinpath(dirname(trans_path), "case5_withload.m")
         pmd_file = joinpath(dirname(dist_path), "case3_balanced_withoutgen.dss")