Merge pull request #1033 from NREL-Sienna/jd/radial_branch_reduction

Jd/radial branch reduction

Project.toml

@@ -45,13 +45,13 @@ JuMP = "1"
 LinearAlgebra = "1"
 Logging = "1"
 MathOptInterface = "1"
-PowerModels = "~0.20"
-PowerNetworkMatrices = "^0.9"
+PowerModels = "^0.20"
+PowerNetworkMatrices = "^0.10"
 PowerSystems = "^3"
 PrettyTables = "2"
 ProgressMeter = "^1.5"
 SHA = "0.7"
 Serialization = "1"
-TimeSeries = "~0.23"
+TimeSeries = "~0.23, ~0.24"
 TimerOutputs = "~0.5"
 julia = "^1.6"

src/core/definitions.jl

@@ -53,9 +53,6 @@ const KiB = 1024
 const MiB = KiB * KiB
 const GiB = MiB * KiB
 
-const UNSUPPORTED_POWERMODELS =
-    [PM.SOCBFPowerModel, PM.SOCBFConicPowerModel, PM.IVRPowerModel]
-
 const PSI_NAME_DELIMITER = "__"
 
 const M_VALUE = 1e6

src/core/network_model.jl

@@ -20,7 +20,7 @@ Establishes the model for a particular device specified by type.
   - `use_slacks::Bool`: Adds slacks to the network modelings
   - `PTDF::PTDF`: PTDF Array calculated using PowerNetworkMatrices
   - `duals::Vector{DataType}`: Constraint types to calculate the duals
-
+  - `reduce_radial_branches::Bool`: Skips modeling radial branches in the system to reduce problem size
 # Example
 
 ptdf_array = PTDF(system)
@@ -32,21 +32,34 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel}
     subnetworks::Dict{Int, Set{Int}}
     bus_area_map::Dict{PSY.ACBus, Int}
     duals::Vector{DataType}
+    radial_network_reduction::PNM.RadialNetworkReduction
+    reduce_radial_branches::Bool
 
     function NetworkModel(
         ::Type{T};
         use_slacks = false,
         PTDF_matrix = nothing,
+        reduce_radial_branches = false,
         subnetworks = Dict{Int, Set{Int}}(),
         duals = Vector{DataType}(),
     ) where {T <: PM.AbstractPowerModel}
         _check_pm_formulation(T)
-        new{T}(use_slacks, PTDF_matrix, subnetworks, Dict{PSY.ACBus, Int}(), duals)
+        new{T}(
+            use_slacks,
+            PTDF_matrix,
+            subnetworks,
+            Dict{PSY.ACBus, Int}(),
+            duals,
+            PNM.RadialNetworkReduction(),
+            reduce_radial_branches,
+        )
     end
 end
 
 get_use_slacks(m::NetworkModel) = m.use_slacks
 get_PTDF_matrix(m::NetworkModel) = m.PTDF_matrix
+get_reduce_radial_branches(m::NetworkModel) = m.reduce_radial_branches
+get_radial_network_reduction(m::NetworkModel) = m.radial_network_reduction
 get_duals(m::NetworkModel) = m.duals
 get_network_formulation(::NetworkModel{T}) where {T} = T
 get_reference_buses(m::NetworkModel{T}) where {T <: PM.AbstractPowerModel} =
@@ -62,13 +75,26 @@ function add_dual!(model::NetworkModel, dual)
     return
 end
 
+function check_radial_branch_reduction_compatibility(
+    ::Type{T},
+) where {T <: PM.AbstractPowerModel}
+    if T ∈ INCOMPATIBLE_WITH_RADIAL_BRANCHES_POWERMODELS
+        error("Network Model $T is not compatible with radial branch reduction")
+    end
+    return
+end
+
 function instantiate_network_model(
     model::NetworkModel{T},
     sys::PSY.System,
 ) where {T <: PM.AbstractPowerModel}
     if isempty(model.subnetworks)
         model.subnetworks = PNM.find_subnetworks(sys)
     end
+    if model.reduce_radial_branches
+        check_radial_branch_reduction_compatibility(T)
+        model.radial_network_reduction = PNM.RadialNetworkReduction(sys)
+    end
     return
 end
 
@@ -90,7 +116,12 @@ end
 function instantiate_network_model(model::NetworkModel{StandardPTDFModel}, sys::PSY.System)
     if get_PTDF_matrix(model) === nothing
         @info "PTDF Matrix not provided. Calculating using PowerNetworkMatrices.PTDF"
-        model.PTDF_matrix = PNM.PTDF(sys)
+        model.PTDF_matrix =
+            PNM.PTDF(sys; reduce_radial_branches = model.reduce_radial_branches)
+    end
+    if model.reduce_radial_branches
+        @assert !isempty(model.PTDF_matrix.radial_network_reduction)
+        model.radial_network_reduction = model.PTDF_matrix.radial_network_reduction
     end
     get_PTDF_matrix(model).subnetworks
     model.subnetworks = deepcopy(get_PTDF_matrix(model).subnetworks)
@@ -107,19 +138,27 @@ function _assign_subnetworks_to_buses(
 ) where {T <: Union{CopperPlatePowerModel, StandardPTDFModel}}
     subnetworks = model.subnetworks
     temp_bus_map = Dict{Int, Int}()
-    for bus in get_available_components(PSY.ACBus, sys)
+    radial_network_reduction = PSI.get_radial_network_reduction(model)
+    for bus in PSI.get_available_components(PSY.ACBus, sys)
         bus_no = PSY.get_number(bus)
+        mapped_bus_no = PNM.get_mapped_bus_number(radial_network_reduction, bus)
         if haskey(temp_bus_map, bus_no)
             model.bus_area_map[bus] = temp_bus_map[bus_no]
+            continue
         else
+            bus_mapped = false
             for (subnet, bus_set) in subnetworks
-                if bus_no ∈ bus_set
+                if mapped_bus_no ∈ bus_set
                     temp_bus_map[bus_no] = subnet
                     model.bus_area_map[bus] = subnet
+                    bus_mapped = true
                     break
                 end
             end
         end
+        if !bus_mapped
+            error("Bus $(PSY.summary(bus)) not mapped to any reference bus")
+        end
     end
     return
 end

src/core/optimization_container.jl

@@ -431,9 +431,14 @@ function _make_system_expressions!(
     subnetworks::Dict{Int, Set{Int}},
     dc_bus_numbers::Vector{Int},
     ::Type{<:PM.AbstractPowerModel},
+    bus_reduction_map::Dict{Int64, Set{Int64}},
 )
     time_steps = get_time_steps(container)
-    ac_bus_numbers = collect(Iterators.flatten(values(subnetworks)))
+    if isempty(bus_reduction_map)
+        ac_bus_numbers = collect(Iterators.flatten(values(subnetworks)))
+    else
+        ac_bus_numbers = collect(keys(bus_reduction_map))
+    end
     container.expressions = Dict(
         ExpressionKey(ActivePowerBalance, PSY.ACBus) =>
             _make_container_array(ac_bus_numbers, time_steps),
@@ -450,9 +455,14 @@ function _make_system_expressions!(
     subnetworks::Dict{Int, Set{Int}},
     dc_bus_numbers::Vector{Int},
     ::Type{<:PM.AbstractActivePowerModel},
+    bus_reduction_map::Dict{Int64, Set{Int64}},
 )
     time_steps = get_time_steps(container)
-    ac_bus_numbers = collect(Iterators.flatten(values(subnetworks)))
+    if isempty(bus_reduction_map)
+        ac_bus_numbers = collect(Iterators.flatten(values(subnetworks)))
+    else
+        ac_bus_numbers = collect(keys(bus_reduction_map))
+    end
     container.expressions = Dict(
         ExpressionKey(ActivePowerBalance, PSY.ACBus) =>
             _make_container_array(ac_bus_numbers, time_steps),
@@ -467,6 +477,7 @@ function _make_system_expressions!(
     subnetworks::Dict{Int, Set{Int}},
     ::Vector{Int},
     ::Type{CopperPlatePowerModel},
+    bus_reduction_map::Dict{Int64, Set{Int64}},
 )
     time_steps = get_time_steps(container)
     subnetworks_ref_buses = collect(keys(subnetworks))
@@ -482,17 +493,24 @@ function _make_system_expressions!(
     subnetworks::Dict{Int, Set{Int}},
     dc_bus_numbers::Vector{Int},
     ::Type{T},
+    bus_reduction_map::Dict{Int64, Set{Int64}},
 ) where {T <: Union{PTDFPowerModel, StandardPTDFModel}}
     time_steps = get_time_steps(container)
-    ac_bus_numbers = sort!(collect(Iterators.flatten(values(subnetworks))))
+    if isempty(bus_reduction_map)
+        ac_bus_numbers = collect(Iterators.flatten(values(subnetworks)))
+    else
+        ac_bus_numbers = collect(keys(bus_reduction_map))
+    end
     subnetworks = collect(keys(subnetworks))
     container.expressions = Dict(
         ExpressionKey(ActivePowerBalance, PSY.System) =>
             _make_container_array(subnetworks, time_steps),
         ExpressionKey(ActivePowerBalance, PSY.DCBus) =>
             _make_container_array(dc_bus_numbers, time_steps),
         ExpressionKey(ActivePowerBalance, PSY.ACBus) =>
-            _make_container_array(ac_bus_numbers, time_steps),
+        # Bus numbers are sorted to guarantee consistency in the order between the
+        # containers
+            _make_container_array(sort!(ac_bus_numbers), time_steps),
     )
     return
 end
@@ -502,9 +520,10 @@ function initialize_system_expressions!(
     ::Type{T},
     subnetworks::Dict{Int, Set{Int}},
     system::PSY.System,
+    bus_reduction_map::Dict{Int64, Set{Int64}},
 ) where {T <: PM.AbstractPowerModel}
     dc_bus_numbers = [PSY.get_number(b) for b in PSY.get_components(PSY.DCBus, system)]
-    _make_system_expressions!(container, subnetworks, dc_bus_numbers, T)
+    _make_system_expressions!(container, subnetworks, dc_bus_numbers, T, bus_reduction_map)
     return
 end
 
@@ -520,7 +539,7 @@ function build_impl!(
         transmission,
         transmission_model.subnetworks,
         sys,
-    )
+        transmission_model.radial_network_reduction.bus_reduction_map)
 
     # Order is required
     for device_model in values(template.devices)

src/devices_models/device_constructors/branch_constructor.jl

@@ -106,7 +106,7 @@ function construct_device!(
             container,
             devices,
             device_model,
-            get_network_formulation(network_model),
+            network_model,
         )
     end
     return
@@ -273,7 +273,7 @@ function construct_device!(
         container,
         devices,
         device_model,
-        get_network_formulation(network_model),
+        network_model,
     )
     add_constraint_dual!(container, sys, device_model)
     return
@@ -331,7 +331,7 @@ function construct_device!(
 ) where {T <: PSY.ACBranch}
     devices =
         get_available_components(T, sys, get_attribute(model, "filter_function"))
-    branch_rate_bounds!(container, devices, model, get_network_formulation(network_model))
+    branch_rate_bounds!(container, devices, model, network_model)
 
     add_constraints!(container, RateLimitConstraintFromTo, devices, model, network_model)
     add_constraints!(container, RateLimitConstraintToFrom, devices, model, network_model)
@@ -356,7 +356,7 @@ function construct_device!(
 ) where {T <: PSY.ACBranch}
     devices =
         get_available_components(T, sys, get_attribute(model, "filter_function"))
-    branch_rate_bounds!(container, devices, model, get_network_formulation(network_model))
+    branch_rate_bounds!(container, devices, model, network_model)
     add_constraint_dual!(container, sys, model)
     return
 end