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.9, ^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

@@ -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,24 @@ function add_dual!(model::NetworkModel, dual)
     return
 end
 
+function check_radial_branch_redution_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_redution_compatibility(T)
+        model.radial_network_reduction = PNM.RadialNetworkReduction(sys)
+    end
     return
 end
 
@@ -90,14 +114,19 @@ 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
     get_PTDF_matrix(model).subnetworks
     model.subnetworks = deepcopy(get_PTDF_matrix(model).subnetworks)
     if length(model.subnetworks) > 1
         @debug "System Contains Multiple Subnetworks. Assigning buses to subnetworks."
         _assign_subnetworks_to_buses(model, sys)
     end
+    if model.reduce_radial_branches
+        @assert !isempty(model.PTDF_matrix.radial_branches)
+        model.radial_network_reduction = model.PTDF_matrix.radial_branches
+    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,22 @@ 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),
+            _make_container_array(sort!(ac_bus_numbers), time_steps),
     )
     return
 end
@@ -502,9 +518,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 +537,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

src/devices_models/devices/AC_branches.jl

@@ -104,11 +104,18 @@ function branch_rate_bounds!(
     container::OptimizationContainer,
     devices::IS.FlattenIteratorWrapper{B},
     ::DeviceModel{B, <:AbstractBranchFormulation},
-    ::Type{<:PM.AbstractDCPModel},
+    network_model::NetworkModel{<:PM.AbstractDCPModel},
 ) where {B <: PSY.ACBranch}
     var = get_variable(container, FlowActivePowerVariable(), B)
+
+    radial_network_reduction = get_radial_network_reduction(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_network_reduction)
+
     for d in devices
         name = PSY.get_name(d)
+        if name ∈ radial_branches_names
+            continue
+        end
         for t in get_time_steps(container)
             JuMP.set_upper_bound(var[name, t], PSY.get_rate(d))
             JuMP.set_lower_bound(var[name, t], -1.0 * PSY.get_rate(d))
@@ -121,15 +128,23 @@ function branch_rate_bounds!(
     container::OptimizationContainer,
     devices::IS.FlattenIteratorWrapper{B},
     ::DeviceModel{B, <:AbstractBranchFormulation},
-    ::Type{<:PM.AbstractPowerModel},
+    network_model::NetworkModel{<:PM.AbstractPowerModel},
 ) where {B <: PSY.ACBranch}
     vars = [
         get_variable(container, FlowActivePowerFromToVariable(), B),
         get_variable(container, FlowActivePowerToFromVariable(), B),
     ]
+
+    time_steps = get_time_steps(container)
+    radial_network_reduction = get_radial_network_reduction(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_network_reduction)
+
     for d in devices
         name = PSY.get_name(d)
-        for t in get_time_steps(container), var in vars
+        if name ∈ radial_branches_names
+            continue
+        end
+        for t in time_steps, var in vars
             JuMP.set_upper_bound(var[name, t], PSY.get_rate(d))
             JuMP.set_lower_bound(var[name, t], -1.0 * PSY.get_rate(d))
         end
@@ -168,14 +183,55 @@ function add_constraints!(
     container::OptimizationContainer,
     cons_type::Type{RateLimitConstraint},
     devices::IS.FlattenIteratorWrapper{T},
-    model::DeviceModel{T, U},
-    ::NetworkModel{X},
+    ::DeviceModel{T, U},
+    network_model::NetworkModel{V},
 ) where {
     T <: PSY.ACBranch,
     U <: AbstractBranchFormulation,
-    X <: PM.AbstractActivePowerModel,
+    V <: PM.AbstractActivePowerModel,
 }
-    add_range_constraints!(container, cons_type, FlowActivePowerVariable, devices, model, X)
+    time_steps = get_time_steps(container)
+    radial_network_reduction = get_radial_network_reduction(network_model)
+    if isempty(radial_network_reduction)
+        device_names = [PSY.get_name(d) for d in devices]
+    else
+        device_names = PNM.get_meshed_branches(radial_network_reduction)
+    end
+
+    con_lb =
+        add_constraints_container!(
+            container,
+            cons_type(),
+            T,
+            device_names,
+            time_steps;
+            meta = "lb",
+        )
+    con_ub =
+        add_constraints_container!(
+            container,
+            cons_type(),
+            T,
+            device_names,
+            time_steps;
+            meta = "ub",
+        )
+
+    array = get_variable(container, FlowActivePowerVariable(), T)
+
+    for device in devices
+        ci_name = PSY.get_name(device)
+        if ci_name ∈ PNM.get_radial_branches(radial_network_reduction)
+            continue
+        end
+        limits = get_min_max_limits(device, RateLimitConstraint, U) # depends on constraint type and formulation type
+        for t in time_steps
+            con_ub[ci_name, t] =
+                JuMP.@constraint(container.JuMPmodel, array[ci_name, t] <= limits.max)
+            con_lb[ci_name, t] =
+                JuMP.@constraint(container.JuMPmodel, array[ci_name, t] >= limits.min)
+        end
+    end
     return
 end
 
@@ -215,7 +271,7 @@ function add_constraints!(
     cons_type::Type{RateLimitConstraintFromTo},
     devices::IS.FlattenIteratorWrapper{B},
     ::DeviceModel{B, <:AbstractBranchFormulation},
-    ::NetworkModel{T},
+    network_model::NetworkModel{T},
 ) where {B <: PSY.ACBranch, T <: PM.AbstractPowerModel}
     rating_data = [(PSY.get_name(h), PSY.get_rate(h)) for h in devices]
 
@@ -231,7 +287,13 @@ function add_constraints!(
     )
     constraint = get_constraint(container, cons_type(), B)
 
+    radial_network_reduction = get_radial_network_reduction(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_network_reduction)
+
     for r in rating_data
+        if r[1] ∈ radial_branches_names
+            continue
+        end
         for t in time_steps
             constraint[r[1], t] = JuMP.@constraint(
                 get_jump_model(container),
@@ -250,7 +312,7 @@ function add_constraints!(
     cons_type::Type{RateLimitConstraintToFrom},
     devices::IS.FlattenIteratorWrapper{B},
     ::DeviceModel{B, <:AbstractBranchFormulation},
-    ::NetworkModel{T},
+    network_model::NetworkModel{T},
 ) where {B <: PSY.ACBranch, T <: PM.AbstractPowerModel}
     rating_data = [(PSY.get_name(h), PSY.get_rate(h)) for h in devices]
 
@@ -266,7 +328,13 @@ function add_constraints!(
     )
     constraint = get_constraint(container, cons_type(), B)
 
+    radial_network_reduction = get_radial_network_reduction(network_model)
+    radial_branches_names = PNM.get_radial_branches(radial_network_reduction)
+
     for r in rating_data
+        if r[1] ∈ radial_branches_names
+            continue
+        end
         for t in time_steps
             constraint[r[1], t] = JuMP.@constraint(
                 get_jump_model(container),