Release v0.21.2 (#919)

* Improve performance of calc_connected_components (#914) --------- Co-authored-by: mtanneau3 <[email protected]> * Add in-place conversion to basic data format (#915) * Add in-place conversion to basic data format * Add unit test for in-place make_basic_network --------- Co-authored-by: mtanneau3 <[email protected]> * No deepcopy when renumbering components (#916) * No deepcopy when renumbering components * In-place _renumber_components and unit tests * Fix typo --------- Co-authored-by: mtanneau3 <[email protected]> Co-authored-by: Carleton Coffrin <[email protected]> * add notes to changelog and prep for release * Update psse parser for 3winding transformers (#917) * fix psse parser 3 winding tranformers * update readme * add note to changelog and readme * fix quote counting check in psse parser --------- Co-authored-by: Mathieu Tanneau <[email protected]> Co-authored-by: mtanneau3 <[email protected]> Co-authored-by: Rahmat Heidari <[email protected]>
lanl-ansi · Jul 4, 2024 · 30e3d39 · 30e3d39
1 parent 094e1b9
commit 30e3d39
Show file tree

Hide file tree

Showing 14 changed files with 327 additions and 40 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -4,6 +4,12 @@ PowerModels.jl Change Log
 ### Staged
 - nothing
 
+### v0.21.2
+- In place building of basic network data (#915,#916)
+- Performance improvements to `calc_connected_components` (#914)
+- Fix three winding transformer parsing in psse data (#917)
+- Fix quote counting check in psse parser (#920)
+
 ### v0.21.1
 - Fix bug in `calc_theta_delta_bounds` (#907)
 

diff --git a/Project.toml b/Project.toml
@@ -2,7 +2,7 @@ name = "PowerModels"
 uuid = "c36e90e8-916a-50a6-bd94-075b64ef4655"
 authors = ["Carleton Coffrin"]
 repo = "https://github.com/lanl-ansi/PowerModels.jl"
-version = "0.21.1"
+version = "0.21.2"
 
 [deps]
 InfrastructureModels = "2030c09a-7f63-5d83-885d-db604e0e9cc0"

diff --git a/README.md b/README.md
@@ -61,11 +61,12 @@ The primary developer is Carleton Coffrin (@ccoffrin) with support from the foll
 - David Fobes (@pseudocubic) LANL, PSS(R)E v33 data support
 - Rory Finnegan (@rofinn) Invenia, Memento Logging
 - Frederik Geth (@frederikgeth) CSIRO, storage modeling advise, Branch Flow and current-voltage formulation
+- Rahmat Heidari (@hei06j) CSIRO, improved PSS(R)E data support
 - Jonas Kersulis (@kersulis) University of Michigan, Sparse SDP formulation
 - Miles Lubin (@mlubin) MIT, Julia/JuMP advise
 - Yeesian Ng (@yeesian) MIT, Documenter.jl setup
 - Kaarthik Sundar (@kaarthiksundar) LANL, OBBT utility
-- Mathieu Tanneau (@mtanneau) Georgia Tech, PTDF matrix computation
+- Mathieu Tanneau (@mtanneau) Georgia Tech, PTDF matrix computation, performance and memory improvements
 - Byron Tasseff (@tasseff) LANL, multi-infrastructure updates
 
 

diff --git a/src/core/data.jl b/src/core/data.jl
@@ -2526,16 +2526,43 @@ function calc_connected_components(data::Dict{String,<:Any}; edges=["branch", "d
         end
     end
 
-    component_lookup = Dict(i => Set{Int}([i]) for i in active_bus_ids)
-    touched = Set{Int}()
-
-    for i in active_bus_ids
-        if !(i in touched)
-            _cc_dfs(i, neighbors, component_lookup, touched)
+    sorted_bus_ids = sort(collect(active_bus_ids))
+    i0 = sorted_bus_ids[1] - 1  # this is to track un-visited buses
+    # The two dictionaries below are used to track the connected components
+    # `component_lookup` maps each bus to the ID of the connected component it belongs to, which is the smallest bus ID in said components
+    # `components` maps the connected component ID to a `Set{Int}` of all bus IDs in that component
+    component_lookup = Dict(i => i0 for i in active_bus_ids)
+    components = Dict{Int,Set{Int}}()
+
+    # ⚠️ it is important to iterate over _sorted_ bus IDs to ensure that components are labeled correctly  
+    for i in sorted_bus_ids
+        if component_lookup[i] != i0
+            continue  # bus already flagged; skip
+        end
+
+        # We have a new connected component!
+        component_lookup[i] = i
+        components[i] = Set([i])
+        V = [i]  # list of all bus IDs in this connected component
+        while length(V) > 0
+            j = pop!(V)
+            for k in neighbors[j]
+                if component_lookup[k] == i0
+                    # Bus `k` is connected to a bus `j` that's connected to `i`
+                    push!(V, k)
+                    component_lookup[k] = i
+                    push!(components[i], k)
+                else
+                    # Bus was already visited
+                    # Check that it's in the right connected component, and move on
+                    @assert component_lookup[k] == i "Unexpected error with bus $k while computing connected components; please report this as a bug."
+                    continue
+                end
+            end
         end
     end
 
-    ccs = (Set(values(component_lookup)))
+    ccs = (Set(values(components)))
 
     return ccs
 end

diff --git a/src/core/data_basic.jl b/src/core/data_basic.jl
@@ -17,17 +17,38 @@ users requiring any of the features listed above for their analysis should use
 the non-basic PowerModels routines.
 """
 function make_basic_network(data::Dict{String,<:Any})
+    # These initial checks are redundant with the checks in make_basic_network!
+    # We keep them here so they run _before_ we create a deepcopy of the data
+    # The checks are fast so this redundancy has little performance impact
     if _IM.ismultiinfrastructure(data)
         Memento.error(_LOGGER, "make_basic_network does not support multiinfrastructure data")
     end
-
+    
     if _IM.ismultinetwork(data)
         Memento.error(_LOGGER, "make_basic_network does not support multinetwork data")
     end
 
-    # make a copy of data so that modifications do not change the input data
     data = deepcopy(data)
 
+    make_basic_network!(data)
+
+    return data
+end
+
+"""
+given a powermodels data dict, modifies it in-place to conform to basic network model requirements.
+
+See [`make_basic_network`](@ref) for more information.
+"""
+function make_basic_network!(data::Dict{String,<:Any})
+    if _IM.ismultiinfrastructure(data)
+        Memento.error(_LOGGER, "make_basic_network does not support multiinfrastructure data")
+    end
+
+    if _IM.ismultinetwork(data)
+        Memento.error(_LOGGER, "make_basic_network does not support multinetwork data")
+    end
+
     # TODO transform PWL costs into linear costs
     for (i,gen) in data["gen"]
         if get(gen, "cost_model", 2) != 2
@@ -110,7 +131,7 @@ function make_basic_network(data::Dict{String,<:Any})
     # re-number non-bus component ids
     for comp_key in keys(pm_component_status)
         if comp_key != "bus"
-            data[comp_key] = _renumber_components(data[comp_key])
+            data[comp_key] = _renumber_components!(data[comp_key])
         end
     end
 
@@ -148,18 +169,18 @@ end
 given a component dict returns a new dict where components have been renumbered
 from 1-to-n ordered by the increasing values of the orginal component id.
 """
-function _renumber_components(comp_dict::Dict{String,<:Any})
-    renumbered_dict = Dict{String,Any}()
-
-    comp_ordered = sort([comp for (i,comp) in comp_dict], by=(x) -> x["index"])
+function _renumber_components!(comp_dict::Dict{String,<:Any})
+    comp_ordered = sort([(comp["index"], comp) for (i, comp) in comp_dict], by=(x) -> x[1])
 
-    for (i,comp) in enumerate(comp_ordered)
-        comp = deepcopy(comp)
-        comp["index"] = i
-        renumbered_dict["$i"] = comp
+    # Delete existing keys
+    empty!(comp_dict)
+    # Update component indices and re-build the dict keys
+    for (i_new, (i_old, comp)) in enumerate(comp_ordered)
+        comp["index"] = i_new
+        comp_dict["$(i_new)"] = comp
     end
 
-    return renumbered_dict
+    return comp_dict
 end
 
 

diff --git a/src/io/psse.jl b/src/io/psse.jl
@@ -434,8 +434,13 @@ function _psse2pm_transformer!(pm_data::Dict, pti_data::Dict, import_all::Bool)
                     else
                         br_r, br_x = transformer["R1-2"], transformer["X1-2"]
                     end
-                    br_r *= (transformer["NOMV1"]^2 / _get_bus_value(transformer["I"], "base_kv", pm_data)^2) * (pm_data["baseMVA"] / transformer["SBASE1-2"])
-                    br_x *= (transformer["NOMV1"]^2 / _get_bus_value(transformer["I"], "base_kv", pm_data)^2) * (pm_data["baseMVA"] / transformer["SBASE1-2"])
+                    if isapprox(transformer["NOMV1"], 0.0)
+                        br_r *= (pm_data["baseMVA"] / transformer["SBASE1-2"])
+                        br_x *= (pm_data["baseMVA"] / transformer["SBASE1-2"])
+                    else
+                        br_r *= (transformer["NOMV1"]^2 / _get_bus_value(transformer["I"], "base_kv", pm_data)^2) * (pm_data["baseMVA"] / transformer["SBASE1-2"])
+                        br_x *= (transformer["NOMV1"]^2 / _get_bus_value(transformer["I"], "base_kv", pm_data)^2) * (pm_data["baseMVA"] / transformer["SBASE1-2"])
+                    end
                 end
 
                 # Zeq scaling for tap2 (see eq (4.21b) in PROGRAM APPLICATION GUIDE 1 in PSSE installation folder)
@@ -448,8 +453,13 @@ function _psse2pm_transformer!(pm_data::Dict, pti_data::Dict, import_all::Bool)
                         br_r *= (transformer["WINDV2"]/_get_bus_value(transformer["J"], "base_kv", pm_data))^2
                         br_x *= (transformer["WINDV2"]/_get_bus_value(transformer["J"], "base_kv", pm_data))^2
                     else  # "for off-nominal turns ratio in pu of nominal winding voltage, NOMV1, NOMV2 and NOMV3."
-                        br_r *= (transformer["WINDV2"]*(transformer["NOMV2"]/_get_bus_value(transformer["J"], "base_kv", pm_data)))^2
-                        br_x *= (transformer["WINDV2"]*(transformer["NOMV2"]/_get_bus_value(transformer["J"], "base_kv", pm_data)))^2
+                        if isapprox(transformer["NOMV2"], 0.0)
+                            br_r *= transformer["WINDV2"]^2
+                            br_x *= transformer["WINDV2"]^2
+                        else
+                            br_r *= (transformer["WINDV2"]*(transformer["NOMV2"]/_get_bus_value(transformer["J"], "base_kv", pm_data)))^2
+                            br_x *= (transformer["WINDV2"]*(transformer["NOMV2"]/_get_bus_value(transformer["J"], "base_kv", pm_data)))^2
+                        end
                     end
                 end
 
@@ -474,23 +484,31 @@ function _psse2pm_transformer!(pm_data::Dict, pti_data::Dict, import_all::Bool)
                 if sub_data["rate_c"] == 0.0
                     delete!(sub_data, "rate_c")
                 end
-
+                
                 if import_all
                     sub_data["windv1"] = transformer["WINDV1"]
                     sub_data["windv2"] = transformer["WINDV2"]
                     sub_data["nomv1"] = transformer["NOMV1"]
                     sub_data["nomv2"] = transformer["NOMV2"]
                 end
 
+                # Assumes CW = 1, namely, for off-nominal turns ratio in pu of windning bus base voltage
                 sub_data["tap"] = pop!(transformer, "WINDV1") / pop!(transformer, "WINDV2")
                 sub_data["shift"] = pop!(transformer, "ANG1")
 
                 # Unit Transformations
-                if transformer["CW"] != 1  # NOT "for off-nominal turns ratio in pu of winding bus base voltage"
+                if transformer["CW"] == 1       # "for off-nominal turns ratio in pu of winding bus base voltage"
+                    # do nothing
+                elseif transformer["CW"] == 2   # "for winding voltage in kv"
                     sub_data["tap"] *= _get_bus_value(transformer["J"], "base_kv", pm_data) / _get_bus_value(transformer["I"], "base_kv", pm_data)
-                    if transformer["CW"] == 3  # "for off-nominal turns ratio in pu of nominal winding voltage, NOMV1, NOMV2 and NOMV3."
-                        sub_data["tap"] *= transformer["NOMV1"] / transformer["NOMV2"]
+                elseif transformer["CW"] == 3   # "for off-nominal turns ratio in pu of nominal winding voltage, NOMV1, NOMV2 and NOMV3."
+                    if !iszero(transformer["NOMV1"]) && !iszero(transformer["NOMV2"])
+                        sub_data["tap"] *= (transformer["NOMV1"] / transformer["NOMV2"]) * (_get_bus_value(transformer["J"], "base_kv", pm_data) / _get_bus_value(transformer["I"], "base_kv", pm_data))
+                    else
+                        # do nothing
                     end
+                else
+                    error(_LOGGER, "psse data parsing error, unsupported value for `CW` on transformer")
                 end
 
                 if import_all
@@ -544,13 +562,29 @@ function _psse2pm_transformer!(pm_data::Dict, pti_data::Dict, import_all::Bool)
 
                 # Unit Transformations
                 if transformer["CZ"] != 1  # NOT "for resistance and reactance in pu on system MVA base and winding voltage base"
-                    br_r12 *= (transformer["NOMV1"] / _get_bus_value(bus_id1, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE1-2"])
-                    br_r23 *= (transformer["NOMV2"] / _get_bus_value(bus_id2, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE2-3"])
-                    br_r31 *= (transformer["NOMV3"] / _get_bus_value(bus_id3, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE3-1"])
+                    if isapprox(transformer["NOMV1"], 0.0)
+                        br_r12 *= (pm_data["baseMVA"] / transformer["SBASE1-2"])
+                        br_x12 *= (pm_data["baseMVA"] / transformer["SBASE1-2"])
+                    else
+                        br_r12 *= (transformer["NOMV1"] / _get_bus_value(bus_id1, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE1-2"])
+                        br_x12 *= (transformer["NOMV1"] / _get_bus_value(bus_id1, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE1-2"])
+                    end
 
-                    br_x12 *= (transformer["NOMV1"] / _get_bus_value(bus_id1, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE1-2"])
-                    br_x23 *= (transformer["NOMV2"] / _get_bus_value(bus_id2, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE2-3"])
-                    br_x31 *= (transformer["NOMV3"] / _get_bus_value(bus_id3, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE3-1"])
+                    if isapprox(transformer["NOMV2"], 0.0)
+                        br_r23 *= (pm_data["baseMVA"] / transformer["SBASE2-3"])
+                        br_x23 *= (pm_data["baseMVA"] / transformer["SBASE2-3"])
+                    else
+                        br_r23 *= (transformer["NOMV2"] / _get_bus_value(bus_id2, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE2-3"])
+                        br_x23 *= (transformer["NOMV2"] / _get_bus_value(bus_id2, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE2-3"])
+                    end
+
+                    if isapprox(transformer["NOMV3"], 0.0)
+                        br_r31 *= (pm_data["baseMVA"] / transformer["SBASE3-1"])
+                        br_x31 *= (pm_data["baseMVA"] / transformer["SBASE3-1"])
+                    else
+                        br_r31 *= (transformer["NOMV3"] / _get_bus_value(bus_id3, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE3-1"])
+                        br_x31 *= (transformer["NOMV3"] / _get_bus_value(bus_id3, "base_kv", pm_data))^2 * (pm_data["baseMVA"] / transformer["SBASE3-1"])
+                    end                    
                 end
 
                 # See "Power System Stability and Control", ISBN: 0-07-035958-X, Eq. 6.72
@@ -600,11 +634,18 @@ function _psse2pm_transformer!(pm_data::Dict, pti_data::Dict, import_all::Bool)
                     sub_data["shift"] = pop!(transformer, "ANG$m")
 
                     # Unit Transformations
-                    if transformer["CW"] != 1  # NOT "for off-nominal turns ratio in pu of winding bus base voltage"
+                    if transformer["CW"] == 1       # "for off-nominal turns ratio in pu of winding bus base voltage"
+                        # do nothing
+                    elseif transformer["CW"] == 2   # "for winding voltage in kv"
                         sub_data["tap"] /= _get_bus_value(bus_id, "base_kv", pm_data)
-                        if transformer["CW"] == 3  # "for off-nominal turns ratio in pu of nominal winding voltage, NOMV1, NOMV2 and NOMV3."
-                            sub_data["tap"] *= transformer["NOMV$m"]
+                    elseif transformer["CW"] == 3   # "for off-nominal turns ratio in pu of nominal winding voltage, NOMV1, NOMV2 and NOMV3."
+                        if !iszero(transformer["NOMV$m"])
+                            sub_data["tap"] *= transformer["NOMV$m"] / _get_bus_value(bus_id, "base_kv", pm_data)
+                        else
+                            # do nothing
                         end
+                    else
+                        error(_LOGGER, "psse data parsing error, unsupported value for `CW` on transformer")
                     end
 
                     if import_all

diff --git a/src/io/pti.jl b/src/io/pti.jl
@@ -607,8 +607,8 @@ Comments, typically indicated at the end of a line with a `'/'` character,
 are also extracted separately, and `Array{Array{String}, String}` is returned.
 """
 function _get_line_elements(line::AbstractString)
-    if count(i->(i=="'"), line) % 2 == 1
-        throw(Memento.error(_LOGGER, "There are an uneven number of single-quotes in \"{line}\", the line cannot be parsed."))
+    if count(i->(i=='\''), line) % 2 == 1
+        throw(Memento.error(_LOGGER, "There are an uneven number of single-quotes in \"$line\", the line cannot be parsed."))
     end
 
     line_comment = split(line, _comment_split, limit=2)

diff --git a/test/data-basic.jl b/test/data-basic.jl
@@ -35,6 +35,57 @@
         @test isapprox(result["objective"], 16551.7; atol=1e0)
     end
 
+    @testset "Re-number components in-place" begin
+        data = PowerModels.parse_file("../test/data/matpower/case7_tplgy.m")
+
+        # Scramble indices before renumbering
+        v = [data["gen"]["$i"] for i in 1:3]
+        data["gen"]["1"]["index"] = 2
+        data["gen"]["2"]["index"] = 3
+        data["gen"]["3"]["index"] = 1
+
+        d = PowerModels._renumber_components!(data["gen"])
+        @test length(data["gen"]) == 3
+        # Check index consistency
+        @test data["gen"]["1"]["index"] == 1
+        @test data["gen"]["2"]["index"] == 2
+        @test data["gen"]["3"]["index"] == 3
+        # Check that modifications were done in-place
+        @test d === data["gen"]
+        @test data["gen"]["1"] === v[3]
+        @test data["gen"]["2"] === v[1]
+        @test data["gen"]["3"] === v[2]
+    end
+
+    @testset "In-place make_basic_network!" begin
+        data = PowerModels.parse_file("../test/data/matpower/case7_tplgy.m")
+        data_ = deepcopy(data)
+
+        # Test #1: `make_basic_network` is non-destructive
+        data_basic = PowerModels.make_basic_network(data)
+        @test data == data_        # input data not modified
+        @test data_basic !== data  # we created a new dict
+
+        # Test #2: `make_basic_network!` is destructive and works in-place
+        data_basic = PowerModels.make_basic_network!(data)
+        @test data_basic === data  # in-place modification
+
+        # Test #3: double-check the output of `make_basic_network!`
+        #   This may be redundant with the above tests, unless `make_basic_network` does not
+        #   call the in-place version under the hood
+        @test length(data_basic["bus"]) == 4
+        @test length(data_basic["load"]) == 2
+        @test length(data_basic["shunt"]) == 2
+        @test length(data_basic["gen"]) == 1
+        @test length(data_basic["branch"]) == 2
+        @test length(data_basic["dcline"]) == 0
+        @test length(data_basic["switch"]) == 0
+        @test length(data_basic["storage"]) == 0
+
+        result = solve_opf(data_basic, ACPPowerModel, nlp_solver)
+        @test isapprox(result["objective"], 1036.52; atol=1e0)
+    end
+
 end