From 83e0de2be98caa9af6bf86e07faf2756b991027e Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 19 May 2026 15:07:28 -0400 Subject: [PATCH 1/2] IOM-side port of PSI#1579 (SCUC branch formulation using MODF) Bumps PNM compat to ^0.20 for VirtualMODF. Renames the NetworkModel LODF_matrix field to MODF_matrix (PNM.VirtualMODF). Adds an outages field/kwarg/getter on DeviceModel plus a _formulation_supports_outages trait that POM will specialize for AbstractSecurityConstrainedStaticBranch. Adds sparse-axis paths for post-contingency duals across dual_processing, optimization_container, decision_model_store, and add_constraint_dual. Extends Base.isempty/empty! on BranchReductionOptimizationTracker to cover constraint_map_by_type. PSI's power-specific changes (devices_models, network_models, services_models, contingency_model, template_validation, etc.) remain out of scope here and will follow in POM. Co-Authored-By: Claude Opus 4.7 --- Project.toml | 2 +- src/InfrastructureOptimizationModels.jl | 6 ++-- src/common_models/add_constraint_dual.jl | 31 +++++++++++++------ src/core/device_model.jl | 39 ++++++++++++++++++++++++ src/core/dual_processing.jl | 17 ++++++++++- src/core/network_model.jl | 17 ++++++----- src/core/network_reductions.jl | 4 ++- src/core/optimization_container.jl | 12 +++++++- src/operation/decision_model_store.jl | 22 +++++++++++++ 9 files changed, 127 insertions(+), 23 deletions(-) diff --git a/Project.toml b/Project.toml index 4419271f..936f17f5 100644 --- a/Project.toml +++ b/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" diff --git a/src/InfrastructureOptimizationModels.jl b/src/InfrastructureOptimizationModels.jl index 1fa5f860..1a387cb6 100644 --- a/src/InfrastructureOptimizationModels.jl +++ b/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 diff --git a/src/common_models/add_constraint_dual.jl b/src/common_models/add_constraint_dual.jl index fa56fc3e..da237ca0 100644 --- a/src/common_models/add_constraint_dual.jl +++ b/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, diff --git a/src/core/device_model.jl b/src/core/device_model.jl index 45b612ef..e4c172b9 100644 --- a/src/core/device_model.jl +++ b/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 diff --git a/src/core/dual_processing.jl b/src/core/dual_processing.jl index 86c95af3..085af597 100644 --- a/src/core/dual_processing.jl +++ b/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 diff --git a/src/core/network_model.jl b/src/core/network_model.jl index 928fa857..402c8631 100644 --- a/src/core/network_model.jl +++ b/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 diff --git a/src/core/network_reductions.jl b/src/core/network_reductions.jl index 47418b0a..435da942 100644 --- a/src/core/network_reductions.jl +++ b/src/core/network_reductions.jl @@ -37,7 +37,8 @@ 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, @@ -45,6 +46,7 @@ Base.empty!( 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() diff --git a/src/core/optimization_container.jl b/src/core/optimization_container.jl index 36f1979b..687cbf3c 100644 --- a/src/core/optimization_container.jl +++ b/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...] diff --git a/src/operation/decision_model_store.jl b/src/operation/decision_model_store.jl index e9f65c64..ac5e0577 100644 --- a/src/operation/decision_model_store.jl +++ b/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; From 877325a89343e569d9c481da37014d3408ac0720 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 19 May 2026 16:42:13 -0400 Subject: [PATCH 2/2] Round 2: pull PNM-only helpers from PSI#1579 into IOM Adds four MODF helpers that PSI placed in core/network_model.jl but are PSY-free, so they belong in IOM rather than POM: - _build_network_reductions: pure PNM, drop-in. - _validate_provided_modf_reduction!: pure PNM, drop-in. - _template_has_outage_aware_branch: recast from PSI's _template_uses_security_constrained_branch to dispatch on the _formulation_supports_outages trait (already in IOM) rather than the POM-resident AbstractSecurityConstrainedStaticBranch type. - _consolidate_device_model_outages_with_modf!: same recast; name kept. The POM port will call these via IOM.* directly. The remaining new helpers in PSI's network_model.jl (irreducible-bus discovery for monitored components, instantiate_network_model! updates) stay POM-bound because they need PSY.System / PSY.Branch / PSY.has_time_series / etc. Co-Authored-By: Claude Opus 4.7 --- src/core/network_model.jl | 79 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) diff --git a/src/core/network_model.jl b/src/core/network_model.jl index 402c8631..9b0b4d4f 100644 --- a/src/core/network_model.jl +++ b/src/core/network_model.jl @@ -138,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