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

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...]

src/operation/decision_model_store.jl

@@ -125,6 +125,28 @@ function write_output!(
     return
 end
 
+# Sparse expressions (e.g., post-contingency flows keyed by
+# `(outage_id, branch_name, t)`) are pre-allocated as 2D dense storage with
+# the non-time tuple flattened into encoded `"a__b"` columns by
+# `get_column_names_from_axis_array(::SparseAxisArray)`. `to_matrix` returns
+# `(n_time, n_cols)`; transpose to match the `(cols, time)` layout the dense
+# storage expects.
+function write_output!(
+    store::DecisionModelStore,
+    name::Symbol,
+    key::OptimizationContainerKey,
+    index::DecisionModelIndexType,
+    update_timestamp::Dates.DateTime,
+    array::SparseAxisArray{T},
+) where {T}
+    columns = get_column_names_from_axis_array(array)[1]
+    matrix = to_matrix(array)
+    container = getfield(store, get_store_container_type(key))
+    container[key][index] =
+        DenseAxisArray(permutedims(matrix), columns, 1:size(matrix, 1))
+    return
+end
+
 function read_outputs(
     store::DecisionModelStore,
     key::OptimizationContainerKey;