[Do not review yet] IOM-side port of PSI#1579 (SCUC branch formulation using MODF)

Project.toml

@@ -39,7 +39,7 @@ JuMP = "^1.28"
 LinearAlgebra = "1"
 Logging = "1"
 MathOptInterface = "1"
-PowerNetworkMatrices = "^0.19"
+PowerNetworkMatrices = "^0.20"
 PrettyTables = "3.1"
 Random = "^1.10"
 Serialization = "1"

src/InfrastructureOptimizationModels.jl

@@ -16,7 +16,7 @@ import LinearAlgebra
 import JSON3
 import InfrastructureSystems
 import PowerNetworkMatrices
-import PowerNetworkMatrices: PTDF, VirtualPTDF, LODF, VirtualLODF
+import PowerNetworkMatrices: PTDF, VirtualPTDF, LODF, VirtualLODF, VirtualMODF
 import InfrastructureSystems: @assert_op, TableFormat, list_recorder_events, get_name
 import InfrastructureSystems:
     get_value_curve, get_power_units, get_function_data, get_proportional_term,
@@ -168,7 +168,8 @@ export InitialCondition
 
 # Network Relevant Exports
 export NetworkModel
-export get_PTDF_matrix, get_LODF_matrix, get_reduce_radial_branches
+export get_PTDF_matrix, get_MODF_matrix, get_reduce_radial_branches
+export get_outages
 export get_duals, get_reference_buses, get_subnetworks, get_bus_area_map
 export get_power_flow_evaluation, has_subnetworks, get_subsystem
 export set_subsystem!, add_dual!
@@ -507,6 +508,7 @@ export PTDF
 export VirtualPTDF
 export LODF
 export VirtualLODF
+export VirtualMODF
 export get_name
 export get_model_base_power
 export get_optimizer_stats

src/common_models/add_constraint_dual.jl

@@ -79,16 +79,27 @@ function assign_dual_variable!(
     if isempty(metas)
         device_names = IS.get_name.(devices)
         add_dual_container!(container, constraint_type, D, device_names, time_steps)
-    else
-        # Reuse the existing constraint container's row axis so the dual axis
-        # matches the constraint exactly. Network reductions (radial /
-        # degree-two) drop branches that pass the device-model filter, so the
-        # constraint axis is a strict subset of IS.get_name.(devices). Sizing
-        # the dual from the device list would leave the dual broadcast in
-        # process_duals incompatible with the constraint matrix.
-        for meta in metas
-            existing =
-                get_constraint(container, ConstraintKey(constraint_type, D, meta))
+        return
+    end
+    for meta in metas
+        key = ConstraintKey(constraint_type, D, meta)
+        existing = get_constraint(container, key)
+        if existing isa SparseAxisArray
+            # Sparse constraints (e.g. post-contingency flow-rate constraints
+            # keyed by (outage_id, name, t)) have no `axes`. Mirror the
+            # constraint's exact sparse keys into a Float64 dual container so
+            # the dual matches the constraint storage one-to-one.
+            dual_container =
+                SparseAxisArray(Dict(k => zero(Float64) for k in keys(existing.data)))
+            _assign_container!(container.duals, key, dual_container)
+        else
+            # Reuse the existing constraint container's row axis so the dual
+            # axis matches the constraint exactly. Network reductions (radial /
+            # degree-two) drop branches that pass the device-model filter, so
+            # the constraint axis is a strict subset of IS.get_name.(devices).
+            # Sizing the dual from the device list would leave the dual
+            # broadcast in process_duals incompatible with the constraint
+            # matrix.
             row_axis = axes(existing)[1]
             add_dual_container!(
                 container,

src/core/device_model.jl

@@ -24,6 +24,7 @@ end
         duals::Vector{DataType},
         services::Vector{ServiceModel}
         attributes::Dict{String, Any}
+        outages::AbstractVector{<:IS.InfrastructureSystemsComponent}
     )
 
 Establishes the model for a particular device specified by type. Uses the keyword argument
@@ -38,6 +39,15 @@ feedforward to enable passing values between operation model at simulation time
   - `duals::Vector{DataType} = Vector{DataType}()`: use to pass constraint type to calculate the duals. The DataType needs to be a valid ConstraintType
   - `time_series_names::Dict{Type{<:TimeSeriesParameter}, String} = get_default_time_series_names(D, B)` : use to specify time series names associated to the device`
   - `attributes::Dict{String, Any} = get_default_attributes(D, B)` : use to specify attributes to the device
+  - `outages::AbstractVector{<:IS.InfrastructureSystemsComponent} = IS.InfrastructureSystemsComponent[]` :
+    N-1 contingencies to model when the formulation is security-constrained. The
+    constructor stores the `IS.get_uuid(outage)` of each entry as a key in the model's
+    `outages::Dict{UUID, Dict{DataType, Set{String}}}` field with empty inner maps;
+    template validation in downstream packages fills the inner maps with the per-type
+    set of monitored component names that each outage carries. Power-specific
+    validation (e.g. checking that entries are `PSY.Outage` subtypes) lives in
+    `PowerOperationsModels`. If `B` is not security-constrained, a non-empty value is
+    dropped with a warning.
 
 # Example
 ```julia
@@ -55,6 +65,7 @@ mutable struct DeviceModel{
     time_series_names::Dict{Type{<:ParameterType}, String}
     attributes::Dict{String, Any}
     subsystem::Union{Nothing, String}
+    outages::Dict{Base.UUID, Dict{DataType, Set{String}}}
     device_cache::Vector{D}
     function DeviceModel(
         ::Type{D},
@@ -64,6 +75,8 @@ mutable struct DeviceModel{
         duals = Vector{DataType}(),
         time_series_names = get_default_time_series_names(D, B),
         attributes = Dict{String, Any}(),
+        outages::AbstractVector{<:IS.InfrastructureSystemsComponent} =
+        IS.InfrastructureSystemsComponent[],
     ) where {D <: IS.InfrastructureSystemsComponent, B <: AbstractDeviceFormulation}
         attributes_ = get_default_attributes(D, B)
         for (k, v) in attributes
@@ -72,6 +85,7 @@ mutable struct DeviceModel{
 
         _check_device_formulation(D)
         _check_device_formulation(B)
+        outages_field = _add_device_model_outages(D, B, outages)
         new{D, B}(
             feedforwards,
             use_slacks,
@@ -80,11 +94,35 @@ mutable struct DeviceModel{
             time_series_names,
             attributes_,
             nothing,
+            outages_field,
             Vector{D}(),
         )
     end
 end
 
+function _add_device_model_outages(
+    ::Type{D},
+    ::Type{B},
+    outages::AbstractVector{<:IS.InfrastructureSystemsComponent},
+) where {D <: IS.InfrastructureSystemsComponent, B <: AbstractDeviceFormulation}
+    field = Dict{Base.UUID, Dict{DataType, Set{String}}}()
+    isempty(outages) && return field
+    if !_formulation_supports_outages(B)
+        @warn "DeviceModel{$D, $B}: 'outages' kwarg ignored — formulation does \
+               not support N-1 contingencies."
+        return field
+    end
+    for outage in outages
+        field[IS.get_uuid(outage)] = Dict{DataType, Set{String}}()
+    end
+    return field
+end
+
+# Multi-dispatch flag for formulations that consume `DeviceModel.outages`.
+# Default: false. Security-constrained branch formulations live in
+# PowerOperationsModels and specialize this trait to `true` there.
+_formulation_supports_outages(::Type{<:AbstractDeviceFormulation}) = false
+
 get_component_type(
     ::DeviceModel{D, B},
 ) where {D <: IS.InfrastructureSystemsComponent, B <: AbstractDeviceFormulation} = D
@@ -101,6 +139,7 @@ get_attributes(m::DeviceModel) = m.attributes
 get_attribute(::Nothing, ::String) = nothing
 get_attribute(m::DeviceModel, key::String) = get(m.attributes, key, nothing)
 get_subsystem(m::DeviceModel) = m.subsystem
+get_outages(m::DeviceModel) = m.outages
 get_device_cache(m::DeviceModel) = m.device_cache
 
 set_subsystem!(m::DeviceModel, id::String) = m.subsystem = id

src/core/dual_processing.jl

@@ -1,3 +1,18 @@
+# DenseAxisArray duals broadcast over the backing array. Post-contingency
+# duals are SparseAxisArray (Dict-backed), where `.data .= …` is undefined, so
+# copy per key instead.
+function _copy_dual_values!(dual::DenseAxisArray, constraint::DenseAxisArray)
+    dual.data .= jump_value.(constraint).data
+    return
+end
+
+function _copy_dual_values!(dual::SparseAxisArray, constraint::SparseAxisArray)
+    for (k, cref) in constraint.data
+        dual.data[k] = jump_value(cref)
+    end
+    return
+end
+
 function process_duals(container::OptimizationContainer, lp_optimizer)
     var_container = get_variables(container)
     for (k, v) in var_container
@@ -68,7 +83,7 @@ function process_duals(container::OptimizationContainer, lp_optimizer)
     if JuMP.has_duals(jump_model)
         for (key, dual) in get_duals(container)
             constraint = get_constraint(container, key)
-            dual.data .= jump_value.(constraint).data
+            _copy_dual_values!(dual, constraint)
         end
     end
 

src/core/network_model.jl

@@ -28,8 +28,11 @@ Establishes the NetworkModel for a given AC network formulation type.
     Adds slack buses to the network modeling.
 - `PTDF_matrix::Union{PNM.PowerNetworkMatrix, Nothing}` = nothing
     PTDF/VirtualPTDF matrix produced by PowerNetworkMatrices (optional).
-- `LODF_matrix::Union{PNM.PowerNetworkMatrix, Nothing}` = nothing
-    LODF/VirtualLODF matrix produced by PowerNetworkMatrices (optional).
+- `MODF_matrix::Union{PNM.VirtualMODF, Nothing}` = nothing
+    VirtualMODF matrix for security-constrained models (N-k contingencies).
+    If `nothing` and the template includes a security-constrained branch
+    formulation, the matrix is constructed from the system during
+    `instantiate_network_model!` (same pattern as PTDF).
 - `reduce_radial_branches::Bool` = false
     Enable radial branch reduction when building network matrices.
 - `reduce_degree_two_branches::Bool` = false
@@ -45,7 +48,7 @@ Establishes the NetworkModel for a given AC network formulation type.
 # Notes
 - `modeled_branch_types` and `reduced_branch_tracker` are internal fields managed by the model.
 - `subsystem` can be set after construction via `set_subsystem!(model, id)`.
-- PTDF/LODF inputs are validated against the requested reduction flags and may raise
+- PTDF inputs are validated against the requested reduction flags and may raise
   a ConflictingInputsError if they are inconsistent with `reduce_radial_branches`
   or `reduce_degree_two_branches`.
 
@@ -59,7 +62,7 @@ Establishes the NetworkModel for a given AC network formulation type.
 mutable struct NetworkModel{T <: AbstractPowerModel}
     use_slacks::Bool
     PTDF_matrix::Union{Nothing, PNM.PowerNetworkMatrix}
-    LODF_matrix::Union{Nothing, PNM.PowerNetworkMatrix}
+    MODF_matrix::Union{Nothing, PNM.VirtualMODF}
     subnetworks::Dict{Int, Set{Int}}
     bus_area_map::Dict{IS.InfrastructureSystemsComponent, Int}
     duals::Vector{DataType}
@@ -76,7 +79,7 @@ mutable struct NetworkModel{T <: AbstractPowerModel}
         ::Type{T};
         use_slacks = false,
         PTDF_matrix = nothing,
-        LODF_matrix = nothing,
+        MODF_matrix = nothing,
         reduce_radial_branches = false,
         reduce_degree_two_branches = false,
         subnetworks = Dict{Int, Set{Int}}(),
@@ -91,7 +94,7 @@ mutable struct NetworkModel{T <: AbstractPowerModel}
         new{T}(
             use_slacks,
             PTDF_matrix,
-            LODF_matrix,
+            MODF_matrix,
             subnetworks,
             Dict{IS.InfrastructureSystemsComponent, Int}(),
             duals,
@@ -109,7 +112,7 @@ end
 
 get_use_slacks(m::NetworkModel) = m.use_slacks
 get_PTDF_matrix(m::NetworkModel) = m.PTDF_matrix
-get_LODF_matrix(m::NetworkModel) = m.LODF_matrix
+get_MODF_matrix(m::NetworkModel) = m.MODF_matrix
 get_reduce_radial_branches(m::NetworkModel) = m.reduce_radial_branches
 get_network_reduction(m::NetworkModel) = m.network_reduction
 get_duals(m::NetworkModel) = m.duals
@@ -135,6 +138,85 @@ function add_dual!(model::NetworkModel, dual)
     return
 end
 
+function _build_network_reductions(
+    model::NetworkModel,
+    irreducible_buses::Vector{Int64},
+)
+    reductions = PNM.NetworkReduction[]
+    if model.reduce_radial_branches
+        push!(reductions, PNM.RadialReduction(; irreducible_buses = irreducible_buses))
+    end
+    if model.reduce_degree_two_branches
+        push!(
+            reductions,
+            PNM.DegreeTwoReduction(; irreducible_buses = irreducible_buses),
+        )
+    end
+    return reductions
+end
+
+# Verify a user-provided MODF Matrix was built with the same network reduction
+# as the active reduction (derived from the PTDF Matrix). Equality of the bus
+# reduction map is the decisive check: it fixes the reduced bus/arc numbering
+# the post-contingency builder uses to index `modf_matrix[arc, outage_spec]`.
+function _validate_provided_modf_reduction!(
+    modf::PNM.VirtualMODF,
+    network_reduction::PNM.NetworkReductionData,
+)
+    if PNM.get_bus_reduction_map(modf.network_reduction_data) !=
+       PNM.get_bus_reduction_map(network_reduction)
+        throw(
+            IS.ConflictingInputsError(
+                "The provided MODF Matrix was built with a different network \
+                reduction than the active reduction derived from the PTDF \
+                Matrix. Rebuild the MODF with a consistent network reduction, \
+                or omit it so it is recalculated automatically.",
+            ),
+        )
+    end
+    return
+end
+
+"""
+True if any branch DeviceModel in `branch_models` uses a formulation that
+consumes `DeviceModel.outages` (per `_formulation_supports_outages`). POM's
+`AbstractSecurityConstrainedStaticBranch` specialization makes that trait
+return `true`; non-SC formulations default to `false`.
+"""
+function _template_has_outage_aware_branch(branch_models::BranchModelContainer)
+    for v in values(branch_models)
+        if _formulation_supports_outages(get_formulation(v))
+            return true
+        end
+    end
+    return false
+end
+
+"""
+Drop outages from each outage-aware-branch `DeviceModel` whose UUID isn't
+registered on `modf_matrix`; without this they'd `KeyError` downstream in
+post-contingency expression construction. PNM's `_register_outages!` silently
+skips outages it can't convert to a `NetworkModification`, so the IOM-side
+view of `m.outages` can be a strict superset of what's actually usable.
+"""
+function _consolidate_device_model_outages_with_modf!(
+    branch_models::BranchModelContainer,
+    modf_matrix::PNM.VirtualMODF,
+)
+    registered = PNM.get_registered_contingencies(modf_matrix)
+    for m in values(branch_models)
+        _formulation_supports_outages(get_formulation(m)) || continue
+        for uuid in setdiff(keys(m.outages), keys(registered))
+            @warn "Outage $(uuid) (DeviceModel{$(get_component_type(m)), \
+                   $(get_formulation(m))}) is not registered on the MODF \
+                   matrix and will not contribute any post-contingency \
+                   constraints." _group = LOG_GROUP_MODELS_VALIDATION
+            delete!(m.outages, uuid)
+        end
+    end
+    return
+end
+
 # Default implementations for network model compatibility checks
 # These can be extended in PowerOperationsModels for specific network formulations
 requires_all_branch_models(::Type{<:AbstractPowerModel}) = true

src/core/network_reductions.jl

@@ -37,14 +37,16 @@ Base.isempty(
 ) =
     isempty(reduction_tracker.variable_dict) &&
     isempty(reduction_tracker.parameter_dict) &&
-    isempty(reduction_tracker.constraint_dict)
+    isempty(reduction_tracker.constraint_dict) &&
+    isempty(reduction_tracker.constraint_map_by_type)
 
 Base.empty!(
     reduction_tracker::BranchReductionOptimizationTracker,
 ) = begin
     empty!(reduction_tracker.variable_dict)
     empty!(reduction_tracker.parameter_dict)
     empty!(reduction_tracker.constraint_dict)
+    empty!(reduction_tracker.constraint_map_by_type)
 end
 
 function BranchReductionOptimizationTracker()

src/core/optimization_container.jl

@@ -1292,9 +1292,19 @@ function _calculate_dual_variable_value!(
     T <: ConstraintType,
     D <: Union{IS.InfrastructureSystemsComponent, IS.InfrastructureSystemsContainer},
 }
-    constraint_duals = jump_value.(get_constraint(container, key))
+    constraint_container = get_constraint(container, key)
     dual_variable_container = get_duals(container)[key]
 
+    if constraint_container isa SparseAxisArray
+        # SparseAxisArray (Dict-backed, e.g. post-contingency constraints keyed
+        # by (outage_id, name, t)) has no `axes`, so `Iterators.product` is
+        # undefined. Copy per stored key instead; the dual container was
+        # mirrored from the same keys in `assign_dual_variable!`.
+        _copy_dual_values!(dual_variable_container, constraint_container)
+        return
+    end
+
+    constraint_duals = jump_value.(constraint_container)
     # Needs to loop since the container ordering might not match in the DenseAxisArray
     for index in Iterators.product(axes(constraint_duals)...)
         dual_variable_container[index...] = constraint_duals[index...]