Jd/radial branch reduction

LICENSE

@@ -1,6 +1,6 @@
 BSD 3-Clause License
 
-Copyright (c) 2018, 2023 Alliance for Sustainable Energy, LLC
+Copyright (c) 2018, 2023 Alliance for Sustainable Energy, LLC and The Regents of the University of California
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without

Project.toml

@@ -1,7 +1,7 @@
 name = "PowerSimulations"
 uuid = "e690365d-45e2-57bb-ac84-44ba829e73c4"
 authors = ["Jose Daniel Lara", "Clayton Barrows", "Daniel Thom", "Dheepak Krishnamurthy", "Sourabh Dalvi"]
-version = "0.25.2"
+version = "0.26.0"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -45,13 +45,13 @@ JuMP = "1"
 LinearAlgebra = "1"
 Logging = "1"
 MathOptInterface = "1"
-PowerModels = "~0.19"
-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