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,
@@ -170,11 +170,13 @@ 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_evaluations, has_subnetworks, get_subsystem
 export set_subsystem!, add_dual!
 export requires_all_branch_models, supports_branch_filtering, ignores_branch_filtering
+export supports_outages
 export validate_network_model
 export BranchReductionOptimizationTracker
 export get_variable_dict, get_constraint_dict, get_constraint_map_by_type
@@ -531,6 +533,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,30 +79,56 @@ 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))
-            row_axis = axes(existing)[1]
-            add_dual_container!(
-                container,
-                constraint_type,
-                D,
-                row_axis,
-                time_steps;
-                meta = meta,
-            )
-        end
+        return
+    end
+    for meta in metas
+        key = ConstraintKey(constraint_type, D, meta)
+        existing = get_constraint(container, key)
+        _assign_dual_from_existing!(container, key, existing, D, time_steps)
     end
     return
 end
 
+# 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.
+function _assign_dual_from_existing!(
+    container::OptimizationContainer,
+    key::ConstraintKey,
+    existing::SparseAxisArray,
+    ::Type{D},
+    time_steps,
+) where {D}
+    dual_container =
+        SparseAxisArray(Dict(k => zero(Float64) for k in keys(existing.data)))
+    _assign_container!(container.duals, key, dual_container)
+    return
+end
+
+# 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.
+function _assign_dual_from_existing!(
+    container::OptimizationContainer,
+    key::ConstraintKey,
+    existing::DenseAxisArray,
+    ::Type{D},
+    time_steps,
+) where {D}
+    row_axis = axes(existing)[1]
+    add_dual_container!(
+        container,
+        get_entry_type(key),
+        D,
+        row_axis,
+        time_steps;
+        meta = key.meta,
+    )
+    return
+end
+
 function _existing_constraint_metas(
     container::OptimizationContainer,
     ::Type{T},

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,11 @@ mutable struct DeviceModel{
     time_series_names::Dict{Type{<:ParameterType}, String}
     attributes::Dict{String, Any}
     subsystem::Union{Nothing, String}
+    # Keyed by UUID to match PNM's `get_registered_contingencies(::VirtualMODF) ::
+    # Dict{UUID, ContingencySpec}` so the consolidation step in network_model.jl can
+    # set-diff directly. UUIDs are also stable across (de)serialization in a way that
+    # live component references aren't.
+    outages::Dict{Base.UUID, Dict{DataType, Set{String}}}
     device_cache::Vector{D}
     function DeviceModel(
         ::Type{D},
@@ -64,6 +79,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 +89,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 +98,39 @@ 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 !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
+
+"""
+    supports_outages(::Type{<:AbstractDeviceFormulation}) -> Bool
+
+Trait declaring whether a device formulation consumes `DeviceModel.outages`.
+Defaults to `false`. POM specializes this to `true` for security-constrained branch
+formulations (`AbstractSecurityConstrainedStaticBranch`).
+"""
+supports_outages(::Type{<:AbstractDeviceFormulation}) = false
+
 get_component_type(
     ::DeviceModel{D, B},
 ) where {D <: IS.InfrastructureSystemsComponent, B <: AbstractDeviceFormulation} = D
@@ -101,6 +147,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

@@ -23,8 +23,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
@@ -41,9 +44,9 @@ 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
-  a ConflictingInputsError if they are inconsistent with `reduce_radial_branches`
-  or `reduce_degree_two_branches`.
+- PTDF and MODF 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`.
 
 # Examples (concrete types like PTDFPowerModel, CopperPlatePowerModel are defined in PowerSimulations)
 # ptdf = PNM.VirtualPTDF(system)
@@ -57,7 +60,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}
@@ -74,7 +77,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}}(),
@@ -86,7 +89,7 @@ mutable struct NetworkModel{T <: AbstractPowerModel}
         new{T}(
             use_slacks,
             PTDF_matrix,
-            LODF_matrix,
+            MODF_matrix,
             subnetworks,
             Dict{IS.InfrastructureSystemsComponent, Int}(),
             duals,
@@ -104,7 +107,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
@@ -130,6 +133,88 @@ 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 `supports_outages`). POM's
+`AbstractSecurityConstrainedStaticBranch` specialization makes that trait
+return `true`; non-SC formulations default to `false`.
+
+`BranchModelContainer` (`Dict{Symbol, DeviceModelForBranches}`) is defined at
+the top of this file and exported from IOM.
+"""
+function _template_has_outage_aware_branch(branch_models::BranchModelContainer)
+    for v in values(branch_models)
+        if 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)
+        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