From 5cd2596a880e2d18c37b87829f78b3ba009aeb6f Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Fri, 29 May 2026 18:24:23 -0600 Subject: [PATCH 01/18] update dependency --- Project.toml | 4 +- src/core/network_model.jl | 339 ++++++++++++++++++++------------------ 2 files changed, 179 insertions(+), 164 deletions(-) diff --git a/Project.toml b/Project.toml index 49e401c84..123380dfd 100644 --- a/Project.toml +++ b/Project.toml @@ -49,9 +49,9 @@ JuMP = "^1.28" LinearAlgebra = "1" Logging = "1" MathOptInterface = "1" -PowerFlows = "^0.18" +PowerFlows = "^0.19" PowerModels = "^0.21.5" -PowerNetworkMatrices = "^0.21" +PowerNetworkMatrices = "^0.22" PowerSystems = "^5.10" PrettyTables = "2.4, 3.1" ProgressMeter = "^1.5" diff --git a/src/core/network_model.jl b/src/core/network_model.jl index 55252755a..a9e80ecc7 100644 --- a/src/core/network_model.jl +++ b/src/core/network_model.jl @@ -169,67 +169,59 @@ function _consolidate_device_model_outages_with_modf!( return end +""" +Build fresh `NetworkReduction` specs for one matrix construction. Each reduction +gets an independent `copy` of `irreducible_buses`: PNM mutates that vector in +place while building a matrix (a `DegreeTwoReduction` grows it to the full +retained set), so a spec object must never be shared across two builds or the +second matrix would pin the first's retained set and reduce inconsistently. +Always call this fresh per matrix; never reuse the returned vector. +""" 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)) + push!( + reductions, + PNM.RadialReduction(; irreducible_buses = copy(irreducible_buses)), + ) end if model.reduce_degree_two_branches - push!(reductions, PNM.DegreeTwoReduction(; irreducible_buses = irreducible_buses)) + push!( + reductions, + PNM.DegreeTwoReduction(; irreducible_buses = copy(irreducible_buses)), + ) end return reductions end """ -Buses retained by `nrd` (reduction representatives — i.e. keys of the bus -reduction map). This is the matrix's bus dimension. PNM's invariant is that -irreducible buses are never eliminated, so they remain keys of the bus -reduction map; the `@assert` here makes that invariant load-bearing instead of -silently papered over with a `union`. +Buses retained by `nrd` — the keys of the bus reduction map, which PNM +initializes with every available non-isolated bus and prunes as buses are +eliminated (`Ybus` build). This key set is therefore exactly the matrix's bus +dimension. + +This is *not* the same as `get_irreducible_buses(nrd)`: that field accumulates +(via `union!`) every bus each reduction stage flagged as protected — including +all injection buses added by a `DegreeTwoReduction`. A protected injection bus +that sits on a radial leaf can still be eliminated by a `RadialReduction` stage, +so `irreducible_buses` is not a subset of the retained keys and must not be +treated as one. """ function _retained_buses(nrd::PNM.NetworkReductionData) - retained = Set(keys(PNM.get_bus_reduction_map(nrd))) - @assert issubset(PNM.get_irreducible_buses(nrd), retained) "irreducible buses are not a subset of bus_reduction_map keys; PNM reduction invariant violated" - return retained + return Set(keys(PNM.get_bus_reduction_map(nrd))) end """ -Rebuild PTDF and MODF onto the union of their retained buses when they diverge. -Returns `true` if a rebuild happened; throws if one pass fails to converge them, -since mismatched reductions break the nodal-balance vs. MODF-column dimensions. +True when two network matrices carry the same reduction, compared by full +`bus_reduction_map` equality (keys *and* values). This is the invariant the +post-contingency builder relies on: it resolves monitored arcs from the reduced +bus/arc numbering and indexes `modf[arc, outage_spec]`, so the PTDF (which +dimensions the nodal balance) and the MODF must agree on the entire map, not just +the retained-bus key set. """ -function _reconcile_ptdf_modf_reduction!( - model::NetworkModel{<:AbstractPTDFModel}, - sys::PSY.System, -) - ptdf_nrd = PNM.get_network_reduction_data(model.PTDF_matrix) - modf_nrd = PNM.get_network_reduction_data(get_MODF_matrix(model)) - retained_ptdf = _retained_buses(ptdf_nrd) - retained_modf = _retained_buses(modf_nrd) - retained_ptdf == retained_modf && return false - - @warn "PTDF and MODF reduced to different bus sets \ - (|PTDF retained|=$(length(retained_ptdf)), \ - |MODF retained|=$(length(retained_modf))). Reconciling both onto \ - the cohesive union of retained buses so the nodal-balance and \ - post-contingency dimensions agree." - cohesive = collect(union(retained_ptdf, retained_modf)) - reductions = _build_network_reductions(model, cohesive) - model.PTDF_matrix = - PNM.VirtualPTDF(sys; tol = PTDF_ZERO_TOL, network_reductions = reductions) - model.MODF_matrix = - PNM.VirtualMODF(sys; tol = PTDF_ZERO_TOL, network_reductions = reductions) - - if _retained_buses(PNM.get_network_reduction_data(model.PTDF_matrix)) != - _retained_buses(PNM.get_network_reduction_data(get_MODF_matrix(model))) - throw( - IS.ConflictingInputsError( - "PTDF and MODF reductions remain dimensionally inconsistent \ - after one reconciliation pass; aborting build.", - ), - ) - end - return true +function _reductions_match(a, b) + return PNM.get_bus_reduction_map(PNM.get_network_reduction_data(a)) == + PNM.get_bus_reduction_map(PNM.get_network_reduction_data(b)) end function _get_filters(branch_models::BranchModelContainer) @@ -509,87 +501,73 @@ function instantiate_network_model!( return 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, +function _build_virtual_ptdf( + model::NetworkModel{<:AbstractPTDFModel}, + sys::PSY.System, + irreducible_buses::Vector{Int64}, ) - if PNM.get_bus_reduction_map(PNM.get_network_reduction_data(modf)) != - 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 + return PNM.VirtualPTDF( + sys; + tol = PTDF_ZERO_TOL, + network_reductions = _build_network_reductions(model, irreducible_buses), + ) end -function instantiate_network_model!( +function _build_virtual_modf( model::NetworkModel{<:AbstractPTDFModel}, - branch_models::BranchModelContainer, - number_of_steps::Int, sys::PSY.System, + irreducible_buses::Vector{Int64}, ) - irreducible_buses = - _get_irreducible_buses_due_to_monitored_components(sys, model, branch_models) - if isnothing(get_PTDF_matrix(model)) || !isempty(irreducible_buses) - if !isnothing(get_PTDF_matrix(model)) - @warn "Provided PTDF Matrix is being ignored since irreducible buses were identified from monitored components (TimeSeriesBounds and/or outage-monitored devices). Recalculating PTDF Matrix with PowerNetworkMatrices.VirtualPTDF and the identified irreducible buses." - else - @info "No PTDF Matrix provided. Calculating using PowerNetworkMatrices.VirtualPTDF" - end - - if model.reduce_radial_branches && model.reduce_degree_two_branches - @info "Applying both radial and degree two reductions" - ptdf = PNM.VirtualPTDF( - sys; - tol = PTDF_ZERO_TOL, - network_reductions = PNM.NetworkReduction[ - PNM.RadialReduction(; irreducible_buses = irreducible_buses), - PNM.DegreeTwoReduction(; irreducible_buses = irreducible_buses), - ], - ) - elseif model.reduce_radial_branches - @info "Applying radial reduction" - if !isempty(irreducible_buses) - @warn "Irreducible buses identified from monitored components. The reduction of any radial branch between 2 irreducible buses will be ignored" - end - ptdf = PNM.VirtualPTDF( - sys; - tol = PTDF_ZERO_TOL, - network_reductions = PNM.NetworkReduction[PNM.RadialReduction(; - irreducible_buses = irreducible_buses, - )], - ) - elseif model.reduce_degree_two_branches - @info "Applying degree two reduction" - ptdf = PNM.VirtualPTDF( - sys; - tol = PTDF_ZERO_TOL, - network_reductions = PNM.NetworkReduction[PNM.DegreeTwoReduction(; - irreducible_buses = irreducible_buses, - )], - ) - else - ptdf = PNM.VirtualPTDF(sys; tol = PTDF_ZERO_TOL) - end - model.PTDF_matrix = ptdf - model.network_reduction = deepcopy(PNM.get_network_reduction_data(ptdf)) - else - model.network_reduction = - deepcopy(PNM.get_network_reduction_data(model.PTDF_matrix)) - end + return PNM.VirtualMODF( + sys; + tol = PTDF_ZERO_TOL, + network_reductions = _build_network_reductions(model, irreducible_buses), + ) +end - if !model.reduce_radial_branches && PNM.has_radial_reduction( - PNM.get_reductions(PNM.get_network_reduction_data(model.PTDF_matrix)), +""" +Keep a provided reduced matrix only when it already retains every bus in +`irreducible_buses` (the buses pinned by monitored components — time-series bounds +and/or outage-monitored devices); otherwise rebuild it with `build` so those buses +survive the reduction. Keeping a sufficient matrix avoids discarding valid user +input. Shared by the PTDF and MODF paths; `label` (`"PTDF"` / `"MODF"`) only names +the matrix in the log messages, and `build` is a zero-arg thunk that recomputes it +with the identified irreducible buses pinned. +""" +function _resolve_reduced_matrix( + provided::Union{Nothing, PNM.PowerNetworkMatrix}, + label::String, + irreducible_buses::Vector{Int64}, + build, +) + if isnothing(provided) + @info "No $label Matrix provided. Calculating using PowerNetworkMatrices.Virtual$label" + return build() + end + n_missing = length( + setdiff( + Set(irreducible_buses), + _retained_buses(PNM.get_network_reduction_data(provided)), + ), ) + if n_missing > 0 + @warn "Provided $label Matrix does not retain $n_missing bus(es) pinned by \ + monitored components (time-series bounds and/or outage-monitored \ + devices). Recalculating $label Matrix with \ + PowerNetworkMatrices.Virtual$label and the identified irreducible buses." + return build() + end + return provided +end + +""" +Validate a kept PTDF matrix's embedded reductions against the model's reduce +flags. A rebuilt matrix satisfies these by construction; the checks matter when a +provided matrix is kept. +""" +function _validate_ptdf_reduction_flags(model::NetworkModel{<:AbstractPTDFModel}) + reductions = PNM.get_reductions(PNM.get_network_reduction_data(model.PTDF_matrix)) + if !model.reduce_radial_branches && PNM.has_radial_reduction(reductions) throw( IS.ConflictingInputsError( "The provided PTDF Matrix has reduced radial branches and mismatches the network \ @@ -598,9 +576,7 @@ function instantiate_network_model!( ), ) end - if !model.reduce_degree_two_branches && PNM.has_degree_two_reduction( - PNM.get_reductions(PNM.get_network_reduction_data(model.PTDF_matrix)), - ) + if !model.reduce_degree_two_branches && PNM.has_degree_two_reduction(reductions) throw( IS.ConflictingInputsError( "The provided PTDF Matrix has reduced degree two branches and mismatches the network \ @@ -609,10 +585,7 @@ function instantiate_network_model!( ), ) end - if model.reduce_radial_branches && - PNM.has_ward_reduction( - PNM.get_reductions(PNM.get_network_reduction_data(model.PTDF_matrix)), - ) + if model.reduce_radial_branches && PNM.has_ward_reduction(reductions) throw( IS.ConflictingInputsError( "The provided PTDF Matrix has a ward reduction specified and the keyword argument \\ @@ -621,52 +594,94 @@ function instantiate_network_model!( ), ) end - if model.reduce_radial_branches @assert !isempty(PNM.get_network_reduction_data(model.PTDF_matrix)) end + return +end + +function _set_subnetworks_from_ptdf!( + model::NetworkModel{<:AbstractPTDFModel}, + sys::PSY.System, +) model.subnetworks = _make_subnetworks_from_subnetwork_axes(model.PTDF_matrix) if length(model.subnetworks) > 1 @debug "System Contains Multiple Subnetworks. Assigning buses to subnetworks." _assign_subnetworks_to_buses(model, sys) end + return +end + +""" +Guarantee the PTDF and MODF share an identical `bus_reduction_map` so the +nodal-balance rows and post-contingency MODF columns are dimensionally aligned. +When the two disagree, rebuild them consistently. `VirtualMODF` auto-protects +monitored/outaged buses (`VirtualPTDF` does not), so the MODF's retained set is +the maximal one: rebuild the MODF, then rebuild the PTDF pinning the MODF's +retained buses, which forces the PTDF to reproduce exactly the MODF reduction. +Throw only if the two still diverge, which indicates a PowerNetworkMatrices +reduction defect rather than mismatched inputs. +""" +function _ensure_ptdf_modf_consistent!( + model::NetworkModel{<:AbstractPTDFModel}, + sys::PSY.System, + irreducible_buses::Vector{Int64}, +) + _reductions_match(get_PTDF_matrix(model), get_MODF_matrix(model)) && return + @warn "PTDF and MODF were built with different network reductions. Rebuilding \ + the PTDF to match the MODF's retained buses so the nodal-balance and \ + post-contingency dimensions agree." + model.MODF_matrix = _build_virtual_modf(model, sys, irreducible_buses) + modf_retained = + collect(_retained_buses(PNM.get_network_reduction_data(get_MODF_matrix(model)))) + model.PTDF_matrix = _build_virtual_ptdf(model, sys, modf_retained) + model.network_reduction = deepcopy(PNM.get_network_reduction_data(model.PTDF_matrix)) + _set_subnetworks_from_ptdf!(model, sys) + _reductions_match(get_PTDF_matrix(model), get_MODF_matrix(model)) && return + np = length(_retained_buses(PNM.get_network_reduction_data(get_PTDF_matrix(model)))) + nm = length(_retained_buses(PNM.get_network_reduction_data(get_MODF_matrix(model)))) + throw( + IS.ConflictingInputsError( + "PTDF and MODF resolved to different network reductions even after \ + rebuilding the PTDF to match the MODF's retained buses (PTDF retains \ + $(np) buses, MODF retains $(nm)). The nodal-balance and post-contingency \ + dimensions cannot agree; this indicates a PowerNetworkMatrices reduction defect.", + ), + ) +end + +function instantiate_network_model!( + model::NetworkModel{<:AbstractPTDFModel}, + branch_models::BranchModelContainer, + number_of_steps::Int, + sys::PSY.System, +) + irreducible_buses = + _get_irreducible_buses_due_to_monitored_components(sys, model, branch_models) + # Resolve the PTDF first: it dimensions the nodal balance and supplies the + # active network reduction. Each matrix is (re)built from its own fresh + # reduction spec (see `_build_network_reductions`); reduction objects are + # never shared across the PTDF and MODF builds. + model.PTDF_matrix = _resolve_reduced_matrix( + get_PTDF_matrix(model), + "PTDF", + irreducible_buses, + () -> _build_virtual_ptdf(model, sys, irreducible_buses), + ) + _validate_ptdf_reduction_flags(model) + model.network_reduction = deepcopy(PNM.get_network_reduction_data(model.PTDF_matrix)) + _set_subnetworks_from_ptdf!(model, sys) if _template_uses_security_constrained_branch(branch_models) - if isnothing(get_MODF_matrix(model)) - @info "No MODF Matrix provided. Calculating using PowerNetworkMatrices.VirtualMODF" - model.MODF_matrix = PNM.VirtualMODF( - sys; - tol = PTDF_ZERO_TOL, - network_reductions = _build_network_reductions(model, irreducible_buses), - ) - elseif !isempty(irreducible_buses) - @warn "Provided MODF Matrix is being ignored since irreducible buses were identified from monitored components (TimeSeriesBounds and/or outage-monitored devices). Recalculating MODF Matrix with PowerNetworkMatrices.VirtualMODF and the identified irreducible buses." - model.MODF_matrix = PNM.VirtualMODF( - sys; - tol = PTDF_ZERO_TOL, - network_reductions = _build_network_reductions(model, irreducible_buses), - ) - else - # The provided MODF is kept verbatim. The post-contingency - # expression builder resolves monitored arcs from - # `model.network_reduction` (taken from the PTDF) and then indexes - # `modf_matrix[arc, outage_spec]`; a MODF built with a different - # network reduction or bus numbering would silently mis-key or - # KeyError, so reject the mismatch up front. - _validate_provided_modf_reduction!( - get_MODF_matrix(model), - model.network_reduction, - ) - end - # Reconcile PTDF/MODF reductions before outage consolidation populates - # the branch maps. - if _reconcile_ptdf_modf_reduction!(model, sys) - model.network_reduction = - deepcopy(PNM.get_network_reduction_data(model.PTDF_matrix)) - model.subnetworks = _make_subnetworks_from_subnetwork_axes(model.PTDF_matrix) - if length(model.subnetworks) > 1 - _assign_subnetworks_to_buses(model, sys) - end - end + model.MODF_matrix = _resolve_reduced_matrix( + get_MODF_matrix(model), + "MODF", + irreducible_buses, + () -> _build_virtual_modf(model, sys, irreducible_buses), + ) + # Both matrices must share one reduction so the nodal-balance rows and the + # post-contingency MODF columns line up; reconcile (and fail fast) before + # outage consolidation populates the branch maps. + _ensure_ptdf_modf_consistent!(model, sys, irreducible_buses) _consolidate_device_model_outages_with_modf!(branch_models, get_MODF_matrix(model)) end PNM.populate_branch_maps_by_type!(model.network_reduction, _get_filters(branch_models)) From 34479faa10e84fe2ae3251380a9b2aea374d6e51 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Fri, 29 May 2026 18:24:38 -0600 Subject: [PATCH 02/18] update the reduction management logic --- src/core/network_model.jl | 62 ++++++++++++++++++--------------------- 1 file changed, 28 insertions(+), 34 deletions(-) diff --git a/src/core/network_model.jl b/src/core/network_model.jl index a9e80ecc7..c0f185259 100644 --- a/src/core/network_model.jl +++ b/src/core/network_model.jl @@ -170,26 +170,20 @@ function _consolidate_device_model_outages_with_modf!( end """ -Build fresh `NetworkReduction` specs for one matrix construction. Each reduction -gets an independent `copy` of `irreducible_buses`: PNM mutates that vector in -place while building a matrix (a `DegreeTwoReduction` grows it to the full -retained set), so a spec object must never be shared across two builds or the -second matrix would pin the first's retained set and reduce inconsistently. -Always call this fresh per matrix; never reuse the returned vector. +Build the `NetworkReduction` specs for one matrix construction. As of +PowerNetworkMatrices 0.22 the reduction specs are field-less configuration +objects; the set of buses to protect from elimination is supplied separately via +the `irreducible_buses` keyword on the matrix constructor (`Ybus` / `VirtualPTDF` +/ `VirtualMODF`), which seeds it once into the reduction container and shares it +across all reduction stages. """ -function _build_network_reductions(model::NetworkModel, irreducible_buses::Vector{Int64}) +function _build_network_reductions(model::NetworkModel) reductions = PNM.NetworkReduction[] if model.reduce_radial_branches - push!( - reductions, - PNM.RadialReduction(; irreducible_buses = copy(irreducible_buses)), - ) + push!(reductions, PNM.RadialReduction()) end if model.reduce_degree_two_branches - push!( - reductions, - PNM.DegreeTwoReduction(; irreducible_buses = copy(irreducible_buses)), - ) + push!(reductions, PNM.DegreeTwoReduction()) end return reductions end @@ -426,9 +420,10 @@ function instantiate_network_model!( @info "Applying both radial and degree two reductions" ybus = PNM.Ybus( sys; + irreducible_buses = irreducible_buses, network_reductions = PNM.NetworkReduction[ - PNM.RadialReduction(; irreducible_buses = irreducible_buses), - PNM.DegreeTwoReduction(; irreducible_buses = irreducible_buses), + PNM.RadialReduction(), + PNM.DegreeTwoReduction(), ], ) elseif model.reduce_radial_branches @@ -438,17 +433,15 @@ function instantiate_network_model!( end ybus = PNM.Ybus( sys; - network_reductions = PNM.NetworkReduction[PNM.RadialReduction(; - irreducible_buses = irreducible_buses, - )], + irreducible_buses = irreducible_buses, + network_reductions = PNM.NetworkReduction[PNM.RadialReduction()], ) elseif model.reduce_degree_two_branches @info "Applying degree two reduction" ybus = PNM.Ybus( sys; - network_reductions = PNM.NetworkReduction[PNM.DegreeTwoReduction(; - irreducible_buses = irreducible_buses, - )], + irreducible_buses = irreducible_buses, + network_reductions = PNM.NetworkReduction[PNM.DegreeTwoReduction()], ) else ybus = PNM.Ybus(sys) @@ -509,7 +502,8 @@ function _build_virtual_ptdf( return PNM.VirtualPTDF( sys; tol = PTDF_ZERO_TOL, - network_reductions = _build_network_reductions(model, irreducible_buses), + irreducible_buses = irreducible_buses, + network_reductions = _build_network_reductions(model), ) end @@ -521,7 +515,8 @@ function _build_virtual_modf( return PNM.VirtualMODF( sys; tol = PTDF_ZERO_TOL, - network_reductions = _build_network_reductions(model, irreducible_buses), + irreducible_buses = irreducible_buses, + network_reductions = _build_network_reductions(model), ) end @@ -615,9 +610,9 @@ end """ Guarantee the PTDF and MODF share an identical `bus_reduction_map` so the nodal-balance rows and post-contingency MODF columns are dimensionally aligned. -When the two disagree, rebuild them consistently. `VirtualMODF` auto-protects -monitored/outaged buses (`VirtualPTDF` does not), so the MODF's retained set is -the maximal one: rebuild the MODF, then rebuild the PTDF pinning the MODF's +`VirtualMODF` auto-protects monitored/outaged buses (`VirtualPTDF` does not), so +the MODF's retained set is the maximal one and is taken as the source of truth. +When the two disagree, keep the MODF and rebuild only the PTDF pinning the MODF's retained buses, which forces the PTDF to reproduce exactly the MODF reduction. Throw only if the two still diverge, which indicates a PowerNetworkMatrices reduction defect rather than mismatched inputs. @@ -625,13 +620,11 @@ reduction defect rather than mismatched inputs. function _ensure_ptdf_modf_consistent!( model::NetworkModel{<:AbstractPTDFModel}, sys::PSY.System, - irreducible_buses::Vector{Int64}, ) _reductions_match(get_PTDF_matrix(model), get_MODF_matrix(model)) && return @warn "PTDF and MODF were built with different network reductions. Rebuilding \ the PTDF to match the MODF's retained buses so the nodal-balance and \ post-contingency dimensions agree." - model.MODF_matrix = _build_virtual_modf(model, sys, irreducible_buses) modf_retained = collect(_retained_buses(PNM.get_network_reduction_data(get_MODF_matrix(model)))) model.PTDF_matrix = _build_virtual_ptdf(model, sys, modf_retained) @@ -659,9 +652,10 @@ function instantiate_network_model!( irreducible_buses = _get_irreducible_buses_due_to_monitored_components(sys, model, branch_models) # Resolve the PTDF first: it dimensions the nodal balance and supplies the - # active network reduction. Each matrix is (re)built from its own fresh - # reduction spec (see `_build_network_reductions`); reduction objects are - # never shared across the PTDF and MODF builds. + # active network reduction. The buses to protect from elimination are passed + # as `irreducible_buses` to each matrix constructor (see + # `_build_virtual_ptdf` / `_build_virtual_modf`), which seeds them into the + # reduction container; the reduction specs themselves are field-less config. model.PTDF_matrix = _resolve_reduced_matrix( get_PTDF_matrix(model), "PTDF", @@ -681,7 +675,7 @@ function instantiate_network_model!( # Both matrices must share one reduction so the nodal-balance rows and the # post-contingency MODF columns line up; reconcile (and fail fast) before # outage consolidation populates the branch maps. - _ensure_ptdf_modf_consistent!(model, sys, irreducible_buses) + _ensure_ptdf_modf_consistent!(model, sys) _consolidate_device_model_outages_with_modf!(branch_models, get_MODF_matrix(model)) end PNM.populate_branch_maps_by_type!(model.network_reduction, _get_filters(branch_models)) From e559fe022f91079c8f337476ebc4d996c029a65d Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Sun, 31 May 2026 10:08:26 +0100 Subject: [PATCH 03/18] debug CI/CD failures --- test/includes.jl | 1 + ...ransmission_security_constrained_models.jl | 7 ++ test/test_utils/klu_ci_debug.jl | 72 +++++++++++++++++++ 3 files changed, 80 insertions(+) create mode 100644 test/test_utils/klu_ci_debug.jl diff --git a/test/includes.jl b/test/includes.jl index c2c49d6aa..7f44b5db6 100644 --- a/test/includes.jl +++ b/test/includes.jl @@ -56,6 +56,7 @@ include("test_utils/mbc_system_utils.jl") include("test_utils/mbc_simulation_utils.jl") include("test_utils/events_simulation_utils.jl") include("test_utils/iec_simulation_utils.jl") +include("test_utils/klu_ci_debug.jl") # TEMPORARY CI diagnostic — remove before merge ENV["RUNNING_PSI_TESTS"] = "true" ENV["SIENNA_RANDOM_SEED"] = 1234 # Set a fixed seed for reproducibility in tests diff --git a/test/test_ac_transmission_security_constrained_models.jl b/test/test_ac_transmission_security_constrained_models.jl index a0ea8cbc9..4c2de4350 100644 --- a/test/test_ac_transmission_security_constrained_models.jl +++ b/test/test_ac_transmission_security_constrained_models.jl @@ -4,6 +4,10 @@ # `comment out unfeasible test` (commit 5fe9232bc) was a real modeling issue, # not KLU instability — re-commented and left as a follow-up. +# TEMPORARY CI diagnostic (PR #1619) — trace every KLU solve in this file so the +# last line before a SIGSEGV pinpoints the crashing call. Remove before merge. +PowerNetworkMatrices.KLUWrapper.PNM_KLU_CI_TRACE[] = true + @testset "Security Constrained branch formulation Network DC-PF with VirtualPTDF + auto-MODF" begin # Guards against regressions on the threaded Woodbury code path: combining # VirtualPTDF with MODF contingency solves has shown KLU-solver instability, @@ -1249,3 +1253,6 @@ end nodal = PSI.get_expression(container, PSI.ActivePowerBalance(), PSY.ACBus) @test size(nodal.data, 1) == length(modf_retained) end + +# TEMPORARY CI diagnostic (PR #1619) — disable KLU solve tracing. Remove before merge. +PowerNetworkMatrices.KLUWrapper.PNM_KLU_CI_TRACE[] = false diff --git a/test/test_utils/klu_ci_debug.jl b/test/test_utils/klu_ci_debug.jl new file mode 100644 index 000000000..91db5ec8e --- /dev/null +++ b/test/test_utils/klu_ci_debug.jl @@ -0,0 +1,72 @@ +# TEMPORARY CI DIAGNOSTIC — remove before merge. +# +# PR #1619 fails on Linux/KLU with `ArgumentError: KLU klu_solve failed: invalid +# argument` (KLU_INVALID) followed by SIGSEGV, in the security-constrained MODF +# Woodbury build. The crash does not reproduce on macOS (heap-layout luck), so we +# capture the deterministic Linux signal directly in CI. +# +# CI builds PSI against the *registered* PowerNetworkMatrices 0.22.0, so we cannot +# edit PNM source for CI. Instead we redefine `KLUWrapper.solve!` at runtime to: +# - print the full klu_solve precondition snapshot when the ccall returns FALSE +# (which {Numeric, Symbolic, B} is NULL? what is ldim / nrhs / status?), then +# throw exactly as the original does, and +# - optionally trace every solve's dimensions (gated by PNM_KLU_CI_TRACE[]) so +# the LAST line before a SIGSEGV pinpoints the crashing call. +# +# All output goes to stderr with an explicit flush so it survives a segfault. + +@eval PowerNetworkMatrices.KLUWrapper begin + # Toggled true only around the crashing testset (see + # test_ac_transmission_security_constrained_models.jl) to bound trace volume. + const PNM_KLU_CI_TRACE = Base.Ref(false) + + function solve!( + cache::KLULinSolveCache{Tv, Ti}, + B::StridedVecOrMat{Tv}, + ) where {Tv, Ti} + is_factored(cache) || error("KLULinSolveCache: not factored yet.") + n = _dim(cache) + size(B, 1) == Int(n) || throw(DimensionMismatch( + "size(B, 1) = $(size(B, 1)), cache n = $(Int(n))", + )) + stride(B, 1) == 1 || throw(ArgumentError( + "B must have unit stride in the first dimension.", + )) + nrhs = size(B, 2) + nrhs == 0 && return B + if PNM_KLU_CI_TRACE[] + Base.println( + Base.stderr, + "[KLU-CI-TRACE] pre-solve n(ldim)=$(Int(n)) nrhs=$(Int(nrhs)) ", + "sym_null=$(cache.symbolic == C_NULL) ", + "num_null=$(cache.numeric == C_NULL) ", + "cache_id=$(objectid(cache)) ", + "sym=$(UInt(cache.symbolic)) num=$(UInt(cache.numeric)) ", + "B=$(UInt(pointer(B))) size(B)=$(size(B))", + ) + Base.flush(Base.stderr) + end + ok = _solve_call( + Tv, Ti, cache.symbolic, cache.numeric, n, nrhs, pointer(B), cache.common, + ) + if ok == 0 + Base.println( + Base.stderr, + "[KLU-CI-FAIL] klu_solve returned FALSE: ", + "status=$(Int(cache.common[].status)) ", + "n(ldim)=$(Int(n)) nrhs=$(Int(nrhs)) ", + "sym_null=$(cache.symbolic == C_NULL) ", + "num_null=$(cache.numeric == C_NULL) ", + "B_null=$(pointer(B) == C_NULL) ", + "cache_id=$(objectid(cache)) ", + "sym=$(UInt(cache.symbolic)) num=$(UInt(cache.numeric)) ", + "size(B)=$(size(B)) tid=$(Threads.threadid())", + ) + Base.flush(Base.stderr) + klu_throw(cache.common[], "klu_solve") + end + return B + end +end + +@info "[KLU-CI-DEBUG] KLUWrapper.solve! instrumented for CI diagnostics." From eeb1be9b7acd2178ce821d5ae04c6f6b9d9c3222 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Sun, 31 May 2026 10:28:57 +0100 Subject: [PATCH 04/18] add more instrumentation to debug --- test/test_utils/klu_ci_debug.jl | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/test/test_utils/klu_ci_debug.jl b/test/test_utils/klu_ci_debug.jl index 91db5ec8e..b3ad18d91 100644 --- a/test/test_utils/klu_ci_debug.jl +++ b/test/test_utils/klu_ci_debug.jl @@ -35,9 +35,20 @@ nrhs = size(B, 2) nrhs == 0 && return B if PNM_KLU_CI_TRACE[] + # Read the handles' internal dimension directly to detect a corrupted + # or mismatched factorization before the (possibly fatal) ccall. + # klu_l_symbolic: n is at byte offset 40 (4 doubles + 1 ptr). + # klu_l_numeric: n is at byte offset 0. + sym_n = cache.symbolic == C_NULL ? -1 : + Base.unsafe_load(Base.reinterpret(Base.Ptr{Int64}, cache.symbolic), 6) + num_n = cache.numeric == C_NULL ? -1 : + Base.unsafe_load(Base.reinterpret(Base.Ptr{Int64}, cache.numeric), 1) Base.println( Base.stderr, - "[KLU-CI-TRACE] pre-solve n(ldim)=$(Int(n)) nrhs=$(Int(nrhs)) ", + "[KLU-CI-TRACE] pre-solve Ti=$(Ti) n(ldim)=$(Int(n)) nrhs=$(Int(nrhs)) ", + "Symbolic_n=$(sym_n) Numeric_n=$(num_n) ", + "sizeof_common=$(Base.sizeof(eltype(typeof(cache.common)))) ", + "status=$(Int(cache.common[].status)) ", "sym_null=$(cache.symbolic == C_NULL) ", "num_null=$(cache.numeric == C_NULL) ", "cache_id=$(objectid(cache)) ", From f642037cdf4b62c00f5b243ceb3ca698fb058b39 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Sun, 31 May 2026 11:01:44 +0100 Subject: [PATCH 05/18] add more debug --- test/test_utils/klu_ci_debug.jl | 45 +++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/test/test_utils/klu_ci_debug.jl b/test/test_utils/klu_ci_debug.jl index b3ad18d91..f4685e17f 100644 --- a/test/test_utils/klu_ci_debug.jl +++ b/test/test_utils/klu_ci_debug.jl @@ -81,3 +81,48 @@ end @info "[KLU-CI-DEBUG] KLUWrapper.solve! instrumented for CI diagnostics." + +# --- Parallel AA (Apple Accelerate) solve-path instrumentation --------------- +# The AA backend crashes in the same Woodbury/PTDF solve path on macOS. AA is +# macOS-only, so this block only redefines on Apple. It dumps the opaque factor's +# internal state (status, dimensions, the numeric/symbolic factor pointers, and +# the required-vs-allocated solve-workspace bytes) right before the libSparse +# ccall — so a run under `MallocPreScribble` pinpoints which field is corrupt. +# Shares the same trace toggle as the KLU path. +if Sys.isapple() + @eval PowerNetworkMatrices.AccelerateWrapper begin + function solve!(cache::AAFactorCache, b::StridedVector{Cdouble}) + is_factored(cache) || error("AAFactorCache: not factored yet.") + n = cache.n + length(b) == n || throw(DimensionMismatch( + "length(b) = $(length(b)), cache n = $(n)", + )) + stride(b, 1) == 1 || throw(ArgumentError("b must have unit stride.")) + ws = _ensure_solve_workspace!(cache, 1) + if Main.PowerNetworkMatrices.KLUWrapper.PNM_KLU_CI_TRACE[] + f = cache.numeric + sf = f.symbolicFactorization + req = _solve_workspace_bytes(f, 1) + Base.println( + Base.stderr, + "[AA-CI-TRACE] pre-solve n=$(n) len_b=$(length(b)) ", + "fac_status=$(Int(f.status)) sym_status=$(Int(sf.status)) ", + "sym_rows=$(sf.rowCount) sym_cols=$(sf.columnCount) ", + "numFac=$(UInt(f.numericFactorization)) ", + "symFac=$(UInt(sf.factorization)) ", + "userFactorStorage=$(f.userFactorStorage) ", + "wsStatic=$(f.solveWorkspaceRequiredStatic) ", + "wsPerRHS=$(f.solveWorkspaceRequiredPerRHS) ", + "ws_bytes_req=$(req) ws_have_bytes=$(length(cache.solve_workspace) * 8) ", + "ws_ptr=$(UInt(ws)) nnz=$(cache.nnz) ", + "colStarts_end=$(cache.columnStarts[end]) ", + "n_rowIdx=$(length(cache.rowIndices))", + ) + Base.flush(Base.stderr) + end + GC.@preserve cache _sparse_solve_vector_ws!(cache.numeric, _dense_vector(b), ws) + return b + end + end + @info "[AA-CI-DEBUG] AccelerateWrapper.solve! instrumented for AA diagnostics." +end From 88b280687c9f49a1499e4e0ea262fc649004cc73 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 31 May 2026 10:42:42 +0000 Subject: [PATCH 06/18] fix use-after-free SIGSEGV from deepcopying VirtualPTDF/VirtualMODF in model construction MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit DecisionModel deepcopied the caller's template, which bit-copied the raw libklu factorization Ptr held by a user-provided VirtualPTDF/VirtualMODF. The copy aliased the original's handle while the cleanup finalizer stayed attached only to the original; when GC reclaimed the original template during build!, the finalizer freed the handle out from under the copy. The next klu_solve then read freed memory (Numeric.n == 0) -> KLU_INVALID -> SIGSEGV. Linux crashed deterministically; macOS survived on heap-layout luck. Exposed by PNM 0.22, whose Virtual matrices carry the KLU factorization through a finalizer-guarded Ptr. Add _copy_template_for_build: a deepcopy that shares PTDF_matrix/MODF_matrix by reference (by seeding the deepcopy memo) so the factorization handle is never copied, while still deep-copying the device/branch/service models and network reduction so the model can mutate them independently. Use it in DecisionModel (replacing the unsafe deepcopy) and EmulationModel (which previously mutated the caller's template in place — the same contract violation). https://claude.ai/code/session_01KNFhS4Z85haH3HZFT8VJC4 --- src/operation/decision_model.jl | 2 +- src/operation/emulation_model.jl | 5 +++-- src/operation/problem_template.jl | 33 +++++++++++++++++++++++++++++++ 3 files changed, 37 insertions(+), 3 deletions(-) diff --git a/src/operation/decision_model.jl b/src/operation/decision_model.jl index def56435c..5e1a3f3b9 100644 --- a/src/operation/decision_model.jl +++ b/src/operation/decision_model.jl @@ -70,7 +70,7 @@ function DecisionModel{M}( OptimizationContainer(sys, settings, jump_model, ts_type), ) - template_ = deepcopy(template) + template_ = _copy_template_for_build(template) finalize_template!(template_, sys) model = DecisionModel{M}( name, diff --git a/src/operation/emulation_model.jl b/src/operation/emulation_model.jl index 5d955b341..1cbb1588f 100644 --- a/src/operation/emulation_model.jl +++ b/src/operation/emulation_model.jl @@ -59,13 +59,14 @@ mutable struct EmulationModel{M <: EmulationProblem} <: OperationModel elseif name isa String name = Symbol(name) end - finalize_template!(template, sys) + template_ = _copy_template_for_build(template) + finalize_template!(template_, sys) internal = ISOPT.ModelInternal( OptimizationContainer(sys, settings, jump_model, PSY.SingleTimeSeries), ) new{M}( name, - template, + template_, sys, internal, SimulationInfo(), diff --git a/src/operation/problem_template.jl b/src/operation/problem_template.jl index edae4d5dd..d8fdb5b15 100644 --- a/src/operation/problem_template.jl +++ b/src/operation/problem_template.jl @@ -349,6 +349,39 @@ function _populate_aggregated_service_model!(template::ProblemTemplate, sys::PSY return end +""" + _copy_template_for_build(template::ProblemTemplate) + +Return an independently-mutable copy of `template` for an operation model to own, +finalize, and instantiate during `build!` without ever mutating the caller's +template. + +This is a `deepcopy` with one deliberate exception: the network model's +`PTDF_matrix` and `MODF_matrix` are *shared by reference* instead of copied. Those +`VirtualPTDF` / `VirtualMODF` objects wrap a libklu (and, on Apple, an Accelerate) +factorization through a raw `Ptr` released by a finalizer. `deepcopy` bit-copies +the `Ptr`, creating a second owner of the same handle while the finalizer stays +attached only to the original; when the original is garbage-collected during +`build!`, it frees the handle out from under the copy and the next `klu_solve` +reads freed memory (`Numeric.n == 0`) → SIGSEGV. + +Sharing is safe because the matrices are immutable, lazily-populated caches: +`instantiate_network_model!` only ever *reassigns* `model.PTDF_matrix` / +`model.MODF_matrix` (it never mutates the objects in place), and a reassignment +swaps the field on the copy without disturbing the shared original. +""" +function _copy_template_for_build(template::ProblemTemplate) + network_model = get_network_model(template) + stackdict = IdDict() + for matrix in (get_PTDF_matrix(network_model), get_MODF_matrix(network_model)) + # Seed the deepcopy memo so each live matrix maps to itself: deepcopy then + # reuses the original in place instead of bit-copying (and recursing into) + # its raw factorization handle. + isnothing(matrix) || (stackdict[matrix] = matrix) + end + return Base.deepcopy_internal(template, stackdict)::ProblemTemplate +end + function finalize_template!(template::ProblemTemplate, sys::PSY.System) _populate_aggregated_service_model!(template, sys) _populate_contributing_devices!(template, sys) From c25cbbb21a0d347ceb13c875f4cca185b983cce5 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 31 May 2026 11:21:44 +0000 Subject: [PATCH 07/18] ci: disable fail-fast in the PR test matrix A failure on one platform (currently macOS, which runs Julia x64 under Rosetta on the arm64 runner) was cancelling the Linux and Windows jobs before they could report. Disabling fail-fast gives a clean per-platform verdict so the Linux/KLU result is visible while the macOS failure is investigated separately. https://claude.ai/code/session_01KNFhS4Z85haH3HZFT8VJC4 --- .github/workflows/pr_testing.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/pr_testing.yml b/.github/workflows/pr_testing.yml index bf0324f34..07b12f4c1 100644 --- a/.github/workflows/pr_testing.yml +++ b/.github/workflows/pr_testing.yml @@ -9,6 +9,7 @@ jobs: name: Julia ${{ matrix.julia-version }} - ${{ matrix.os }} runs-on: ${{ matrix.os }} strategy: + fail-fast: false matrix: julia-version: ['1'] julia-arch: [x64] From 8334a726ff597f583a0ef27452bf37e3426ff555 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Sun, 31 May 2026 14:16:40 +0100 Subject: [PATCH 08/18] additional debugging --- test/test_utils/klu_ci_debug.jl | 71 +++++++++++++++++++++++++++++---- 1 file changed, 64 insertions(+), 7 deletions(-) diff --git a/test/test_utils/klu_ci_debug.jl b/test/test_utils/klu_ci_debug.jl index f4685e17f..281f1cd4c 100644 --- a/test/test_utils/klu_ci_debug.jl +++ b/test/test_utils/klu_ci_debug.jl @@ -37,12 +37,23 @@ if PNM_KLU_CI_TRACE[] # Read the handles' internal dimension directly to detect a corrupted # or mismatched factorization before the (possibly fatal) ccall. - # klu_l_symbolic: n is at byte offset 40 (4 doubles + 1 ptr). - # klu_l_numeric: n is at byte offset 0. - sym_n = cache.symbolic == C_NULL ? -1 : - Base.unsafe_load(Base.reinterpret(Base.Ptr{Int64}, cache.symbolic), 6) - num_n = cache.numeric == C_NULL ? -1 : - Base.unsafe_load(Base.reinterpret(Base.Ptr{Int64}, cache.numeric), 1) + # `n` lives at byte offset 0 in {klu_numeric, klu_l_numeric} and byte + # offset 40 (4 doubles + 1 ptr) in {klu_symbolic, klu_l_symbolic}; its + # width follows Ti (Int32 for the `klu_*` family, Int64/SuiteSparse_long + # for `klu_l_*`). Read with the matching width so Int32 caches don't + # report `n`+`nblocks` packed as one Int64. + if Ti === Int64 + sym_n = cache.symbolic == C_NULL ? -1 : + Int(Base.unsafe_load(Base.reinterpret(Base.Ptr{Int64}, cache.symbolic), 6)) + num_n = cache.numeric == C_NULL ? -1 : + Int(Base.unsafe_load(Base.reinterpret(Base.Ptr{Int64}, cache.numeric), 1)) + else + # Int32 family: symbolic n at byte 40 → Int32 index 11; numeric n at byte 0 → index 1. + sym_n = cache.symbolic == C_NULL ? -1 : + Int(Base.unsafe_load(Base.reinterpret(Base.Ptr{Int32}, cache.symbolic), 11)) + num_n = cache.numeric == C_NULL ? -1 : + Int(Base.unsafe_load(Base.reinterpret(Base.Ptr{Int32}, cache.numeric), 1)) + end Base.println( Base.stderr, "[KLU-CI-TRACE] pre-solve Ti=$(Ti) n(ldim)=$(Int(n)) nrhs=$(Int(nrhs)) ", @@ -80,7 +91,53 @@ end end -@info "[KLU-CI-DEBUG] KLUWrapper.solve! instrumented for CI diagnostics." +# Instrument the handle-free path to identify the `Numeric_n=0` mechanism: +# - If a `[KLU-FREE] cache_id=X` (especially `from_finalizer=true`) appears and a +# `[KLU-CI-TRACE] ... cache_id=X` solve runs on the same X → the numeric was +# FREED out from under an in-use cache (use-after-free; GC/lifetime bug). +# - If `Numeric_n` flips 4→0 for some cache_id with NO `[KLU-FREE]` for it in +# between → the struct was OVERWRITTEN/zeroed without a free (heap overrun). +# The backtrace shows whether the free came from the GC finalizer or an explicit +# call (`symbolic_factor!` / `_recover_factorization!`). Body mirrors the original +# `_free_klu_handles!` exactly; logging is gated by the same trace flag. +@eval PowerNetworkMatrices.KLUWrapper begin + function _free_klu_handles!(cache::KLULinSolveCache{Tv, Ti}) where {Tv, Ti} + if PNM_KLU_CI_TRACE[] && (cache.numeric != C_NULL || cache.symbolic != C_NULL) + bt = Base.stacktrace(Base.backtrace()) + from_finalizer = Base.any( + fr -> Base.occursin( + r"finaliz|run_finalizer|jl_gc|_gc_|gc.c", + Base.string(fr.func) * "|" * Base.string(fr.file), + ), + bt, + ) + Base.println( + Base.stderr, + "[KLU-FREE] cache_id=$(objectid(cache)) ", + "numeric=$(UInt(cache.numeric)) symbolic=$(UInt(cache.symbolic)) ", + "from_finalizer=$(from_finalizer) tid=$(Threads.threadid()) ", + "nframes=$(length(bt))", + ) + for fr in bt[1:min(14, length(bt))] + Base.println(Base.stderr, " @ $(fr.func) $(fr.file):$(fr.line)") + end + Base.flush(Base.stderr) + end + if cache.numeric != C_NULL + num_ref = Base.Ref(cache.numeric) + _free_numeric!(Tv, Ti, num_ref, cache.common) + cache.numeric = num_ref[] + end + if cache.symbolic != C_NULL + sym_ref = Base.Ref(cache.symbolic) + _free_symbolic!(Ti, sym_ref, cache.common) + cache.symbolic = sym_ref[] + end + return nothing + end +end + +@info "[KLU-CI-DEBUG] KLUWrapper.solve! + _free_klu_handles! instrumented for CI diagnostics." # --- Parallel AA (Apple Accelerate) solve-path instrumentation --------------- # The AA backend crashes in the same Woodbury/PTDF solve path on macOS. AA is From 50c555673e3936f6896956e781e1962716c5f3c0 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Sun, 31 May 2026 15:18:03 +0100 Subject: [PATCH 09/18] clean up --- src/operation/problem_template.jl | 50 ++--- test/includes.jl | 1 - ...ransmission_security_constrained_models.jl | 7 - test/test_utils/klu_ci_debug.jl | 185 ------------------ 4 files changed, 25 insertions(+), 218 deletions(-) delete mode 100644 test/test_utils/klu_ci_debug.jl diff --git a/src/operation/problem_template.jl b/src/operation/problem_template.jl index d8fdb5b15..3fbac9596 100644 --- a/src/operation/problem_template.jl +++ b/src/operation/problem_template.jl @@ -352,34 +352,34 @@ end """ _copy_template_for_build(template::ProblemTemplate) -Return an independently-mutable copy of `template` for an operation model to own, -finalize, and instantiate during `build!` without ever mutating the caller's -template. - -This is a `deepcopy` with one deliberate exception: the network model's -`PTDF_matrix` and `MODF_matrix` are *shared by reference* instead of copied. Those -`VirtualPTDF` / `VirtualMODF` objects wrap a libklu (and, on Apple, an Accelerate) -factorization through a raw `Ptr` released by a finalizer. `deepcopy` bit-copies -the `Ptr`, creating a second owner of the same handle while the finalizer stays -attached only to the original; when the original is garbage-collected during -`build!`, it frees the handle out from under the copy and the next `klu_solve` -reads freed memory (`Numeric.n == 0`) → SIGSEGV. - -Sharing is safe because the matrices are immutable, lazily-populated caches: -`instantiate_network_model!` only ever *reassigns* `model.PTDF_matrix` / -`model.MODF_matrix` (it never mutates the objects in place), and a reassignment -swaps the field on the copy without disturbing the shared original. +Return an independently-mutable copy of `template` for `build!` to own, sharing the +network model's `PTDF_matrix`/`MODF_matrix` by reference. Those matrices wrap a +libklu/libSparse factorization behind a finalizer-owned raw pointer that is not +safe to `deepcopy` (copying it aliases one handle across two objects). The matrices +are only ever reassigned during build, never mutated in place, so sharing is safe. """ function _copy_template_for_build(template::ProblemTemplate) - network_model = get_network_model(template) - stackdict = IdDict() - for matrix in (get_PTDF_matrix(network_model), get_MODF_matrix(network_model)) - # Seed the deepcopy memo so each live matrix maps to itself: deepcopy then - # reuses the original in place instead of bit-copying (and recursing into) - # its raw factorization handle. - isnothing(matrix) || (stackdict[matrix] = matrix) + src_network = get_network_model(template) + ptdf = get_PTDF_matrix(src_network) + modf = get_MODF_matrix(src_network) + # Detach the matrices so the copy never deepcopies their solver caches; restore + # the caller's template even if the copy throws. + src_network.PTDF_matrix = nothing + src_network.MODF_matrix = nothing + new_network = try + deepcopy(src_network) + finally + src_network.PTDF_matrix = ptdf + src_network.MODF_matrix = modf end - return Base.deepcopy_internal(template, stackdict)::ProblemTemplate + new_network.PTDF_matrix = ptdf + new_network.MODF_matrix = modf + + new_template = ProblemTemplate(new_network) + new_template.devices = deepcopy(get_device_models(template)) + new_template.branches = deepcopy(get_branch_models(template)) + new_template.services = deepcopy(get_service_models(template)) + return new_template end function finalize_template!(template::ProblemTemplate, sys::PSY.System) diff --git a/test/includes.jl b/test/includes.jl index 7f44b5db6..c2c49d6aa 100644 --- a/test/includes.jl +++ b/test/includes.jl @@ -56,7 +56,6 @@ include("test_utils/mbc_system_utils.jl") include("test_utils/mbc_simulation_utils.jl") include("test_utils/events_simulation_utils.jl") include("test_utils/iec_simulation_utils.jl") -include("test_utils/klu_ci_debug.jl") # TEMPORARY CI diagnostic — remove before merge ENV["RUNNING_PSI_TESTS"] = "true" ENV["SIENNA_RANDOM_SEED"] = 1234 # Set a fixed seed for reproducibility in tests diff --git a/test/test_ac_transmission_security_constrained_models.jl b/test/test_ac_transmission_security_constrained_models.jl index 4c2de4350..a0ea8cbc9 100644 --- a/test/test_ac_transmission_security_constrained_models.jl +++ b/test/test_ac_transmission_security_constrained_models.jl @@ -4,10 +4,6 @@ # `comment out unfeasible test` (commit 5fe9232bc) was a real modeling issue, # not KLU instability — re-commented and left as a follow-up. -# TEMPORARY CI diagnostic (PR #1619) — trace every KLU solve in this file so the -# last line before a SIGSEGV pinpoints the crashing call. Remove before merge. -PowerNetworkMatrices.KLUWrapper.PNM_KLU_CI_TRACE[] = true - @testset "Security Constrained branch formulation Network DC-PF with VirtualPTDF + auto-MODF" begin # Guards against regressions on the threaded Woodbury code path: combining # VirtualPTDF with MODF contingency solves has shown KLU-solver instability, @@ -1253,6 +1249,3 @@ end nodal = PSI.get_expression(container, PSI.ActivePowerBalance(), PSY.ACBus) @test size(nodal.data, 1) == length(modf_retained) end - -# TEMPORARY CI diagnostic (PR #1619) — disable KLU solve tracing. Remove before merge. -PowerNetworkMatrices.KLUWrapper.PNM_KLU_CI_TRACE[] = false diff --git a/test/test_utils/klu_ci_debug.jl b/test/test_utils/klu_ci_debug.jl deleted file mode 100644 index 281f1cd4c..000000000 --- a/test/test_utils/klu_ci_debug.jl +++ /dev/null @@ -1,185 +0,0 @@ -# TEMPORARY CI DIAGNOSTIC — remove before merge. -# -# PR #1619 fails on Linux/KLU with `ArgumentError: KLU klu_solve failed: invalid -# argument` (KLU_INVALID) followed by SIGSEGV, in the security-constrained MODF -# Woodbury build. The crash does not reproduce on macOS (heap-layout luck), so we -# capture the deterministic Linux signal directly in CI. -# -# CI builds PSI against the *registered* PowerNetworkMatrices 0.22.0, so we cannot -# edit PNM source for CI. Instead we redefine `KLUWrapper.solve!` at runtime to: -# - print the full klu_solve precondition snapshot when the ccall returns FALSE -# (which {Numeric, Symbolic, B} is NULL? what is ldim / nrhs / status?), then -# throw exactly as the original does, and -# - optionally trace every solve's dimensions (gated by PNM_KLU_CI_TRACE[]) so -# the LAST line before a SIGSEGV pinpoints the crashing call. -# -# All output goes to stderr with an explicit flush so it survives a segfault. - -@eval PowerNetworkMatrices.KLUWrapper begin - # Toggled true only around the crashing testset (see - # test_ac_transmission_security_constrained_models.jl) to bound trace volume. - const PNM_KLU_CI_TRACE = Base.Ref(false) - - function solve!( - cache::KLULinSolveCache{Tv, Ti}, - B::StridedVecOrMat{Tv}, - ) where {Tv, Ti} - is_factored(cache) || error("KLULinSolveCache: not factored yet.") - n = _dim(cache) - size(B, 1) == Int(n) || throw(DimensionMismatch( - "size(B, 1) = $(size(B, 1)), cache n = $(Int(n))", - )) - stride(B, 1) == 1 || throw(ArgumentError( - "B must have unit stride in the first dimension.", - )) - nrhs = size(B, 2) - nrhs == 0 && return B - if PNM_KLU_CI_TRACE[] - # Read the handles' internal dimension directly to detect a corrupted - # or mismatched factorization before the (possibly fatal) ccall. - # `n` lives at byte offset 0 in {klu_numeric, klu_l_numeric} and byte - # offset 40 (4 doubles + 1 ptr) in {klu_symbolic, klu_l_symbolic}; its - # width follows Ti (Int32 for the `klu_*` family, Int64/SuiteSparse_long - # for `klu_l_*`). Read with the matching width so Int32 caches don't - # report `n`+`nblocks` packed as one Int64. - if Ti === Int64 - sym_n = cache.symbolic == C_NULL ? -1 : - Int(Base.unsafe_load(Base.reinterpret(Base.Ptr{Int64}, cache.symbolic), 6)) - num_n = cache.numeric == C_NULL ? -1 : - Int(Base.unsafe_load(Base.reinterpret(Base.Ptr{Int64}, cache.numeric), 1)) - else - # Int32 family: symbolic n at byte 40 → Int32 index 11; numeric n at byte 0 → index 1. - sym_n = cache.symbolic == C_NULL ? -1 : - Int(Base.unsafe_load(Base.reinterpret(Base.Ptr{Int32}, cache.symbolic), 11)) - num_n = cache.numeric == C_NULL ? -1 : - Int(Base.unsafe_load(Base.reinterpret(Base.Ptr{Int32}, cache.numeric), 1)) - end - Base.println( - Base.stderr, - "[KLU-CI-TRACE] pre-solve Ti=$(Ti) n(ldim)=$(Int(n)) nrhs=$(Int(nrhs)) ", - "Symbolic_n=$(sym_n) Numeric_n=$(num_n) ", - "sizeof_common=$(Base.sizeof(eltype(typeof(cache.common)))) ", - "status=$(Int(cache.common[].status)) ", - "sym_null=$(cache.symbolic == C_NULL) ", - "num_null=$(cache.numeric == C_NULL) ", - "cache_id=$(objectid(cache)) ", - "sym=$(UInt(cache.symbolic)) num=$(UInt(cache.numeric)) ", - "B=$(UInt(pointer(B))) size(B)=$(size(B))", - ) - Base.flush(Base.stderr) - end - ok = _solve_call( - Tv, Ti, cache.symbolic, cache.numeric, n, nrhs, pointer(B), cache.common, - ) - if ok == 0 - Base.println( - Base.stderr, - "[KLU-CI-FAIL] klu_solve returned FALSE: ", - "status=$(Int(cache.common[].status)) ", - "n(ldim)=$(Int(n)) nrhs=$(Int(nrhs)) ", - "sym_null=$(cache.symbolic == C_NULL) ", - "num_null=$(cache.numeric == C_NULL) ", - "B_null=$(pointer(B) == C_NULL) ", - "cache_id=$(objectid(cache)) ", - "sym=$(UInt(cache.symbolic)) num=$(UInt(cache.numeric)) ", - "size(B)=$(size(B)) tid=$(Threads.threadid())", - ) - Base.flush(Base.stderr) - klu_throw(cache.common[], "klu_solve") - end - return B - end -end - -# Instrument the handle-free path to identify the `Numeric_n=0` mechanism: -# - If a `[KLU-FREE] cache_id=X` (especially `from_finalizer=true`) appears and a -# `[KLU-CI-TRACE] ... cache_id=X` solve runs on the same X → the numeric was -# FREED out from under an in-use cache (use-after-free; GC/lifetime bug). -# - If `Numeric_n` flips 4→0 for some cache_id with NO `[KLU-FREE]` for it in -# between → the struct was OVERWRITTEN/zeroed without a free (heap overrun). -# The backtrace shows whether the free came from the GC finalizer or an explicit -# call (`symbolic_factor!` / `_recover_factorization!`). Body mirrors the original -# `_free_klu_handles!` exactly; logging is gated by the same trace flag. -@eval PowerNetworkMatrices.KLUWrapper begin - function _free_klu_handles!(cache::KLULinSolveCache{Tv, Ti}) where {Tv, Ti} - if PNM_KLU_CI_TRACE[] && (cache.numeric != C_NULL || cache.symbolic != C_NULL) - bt = Base.stacktrace(Base.backtrace()) - from_finalizer = Base.any( - fr -> Base.occursin( - r"finaliz|run_finalizer|jl_gc|_gc_|gc.c", - Base.string(fr.func) * "|" * Base.string(fr.file), - ), - bt, - ) - Base.println( - Base.stderr, - "[KLU-FREE] cache_id=$(objectid(cache)) ", - "numeric=$(UInt(cache.numeric)) symbolic=$(UInt(cache.symbolic)) ", - "from_finalizer=$(from_finalizer) tid=$(Threads.threadid()) ", - "nframes=$(length(bt))", - ) - for fr in bt[1:min(14, length(bt))] - Base.println(Base.stderr, " @ $(fr.func) $(fr.file):$(fr.line)") - end - Base.flush(Base.stderr) - end - if cache.numeric != C_NULL - num_ref = Base.Ref(cache.numeric) - _free_numeric!(Tv, Ti, num_ref, cache.common) - cache.numeric = num_ref[] - end - if cache.symbolic != C_NULL - sym_ref = Base.Ref(cache.symbolic) - _free_symbolic!(Ti, sym_ref, cache.common) - cache.symbolic = sym_ref[] - end - return nothing - end -end - -@info "[KLU-CI-DEBUG] KLUWrapper.solve! + _free_klu_handles! instrumented for CI diagnostics." - -# --- Parallel AA (Apple Accelerate) solve-path instrumentation --------------- -# The AA backend crashes in the same Woodbury/PTDF solve path on macOS. AA is -# macOS-only, so this block only redefines on Apple. It dumps the opaque factor's -# internal state (status, dimensions, the numeric/symbolic factor pointers, and -# the required-vs-allocated solve-workspace bytes) right before the libSparse -# ccall — so a run under `MallocPreScribble` pinpoints which field is corrupt. -# Shares the same trace toggle as the KLU path. -if Sys.isapple() - @eval PowerNetworkMatrices.AccelerateWrapper begin - function solve!(cache::AAFactorCache, b::StridedVector{Cdouble}) - is_factored(cache) || error("AAFactorCache: not factored yet.") - n = cache.n - length(b) == n || throw(DimensionMismatch( - "length(b) = $(length(b)), cache n = $(n)", - )) - stride(b, 1) == 1 || throw(ArgumentError("b must have unit stride.")) - ws = _ensure_solve_workspace!(cache, 1) - if Main.PowerNetworkMatrices.KLUWrapper.PNM_KLU_CI_TRACE[] - f = cache.numeric - sf = f.symbolicFactorization - req = _solve_workspace_bytes(f, 1) - Base.println( - Base.stderr, - "[AA-CI-TRACE] pre-solve n=$(n) len_b=$(length(b)) ", - "fac_status=$(Int(f.status)) sym_status=$(Int(sf.status)) ", - "sym_rows=$(sf.rowCount) sym_cols=$(sf.columnCount) ", - "numFac=$(UInt(f.numericFactorization)) ", - "symFac=$(UInt(sf.factorization)) ", - "userFactorStorage=$(f.userFactorStorage) ", - "wsStatic=$(f.solveWorkspaceRequiredStatic) ", - "wsPerRHS=$(f.solveWorkspaceRequiredPerRHS) ", - "ws_bytes_req=$(req) ws_have_bytes=$(length(cache.solve_workspace) * 8) ", - "ws_ptr=$(UInt(ws)) nnz=$(cache.nnz) ", - "colStarts_end=$(cache.columnStarts[end]) ", - "n_rowIdx=$(length(cache.rowIndices))", - ) - Base.flush(Base.stderr) - end - GC.@preserve cache _sparse_solve_vector_ws!(cache.numeric, _dense_vector(b), ws) - return b - end - end - @info "[AA-CI-DEBUG] AccelerateWrapper.solve! instrumented for AA diagnostics." -end From f7febf1333071e58391f3c58bdd500097656a422 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Sun, 31 May 2026 16:26:40 +0100 Subject: [PATCH 10/18] monre Symbol Subtitutions --- src/devices_models/devices/common/add_to_expression.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/devices_models/devices/common/add_to_expression.jl b/src/devices_models/devices/common/add_to_expression.jl index a4518efb1..362e89e1a 100644 --- a/src/devices_models/devices/common/add_to_expression.jl +++ b/src/devices_models/devices/common/add_to_expression.jl @@ -2075,7 +2075,7 @@ function _reduced_entry_in_interface( branch_names = PSY.get_name.(reduction_entry) throw( ArgumentError( - "An interface is specified with only part of a double-circuit that has been reduced. + "An interface is specified with only part of a double-circuit that has been reduced. Branches: $(branch_names[in_interface]) are in the interface and branches: $(branch_names[.!in_interface]) are not. Modify the data to include all of or none of the parallel segements.", ), @@ -2197,7 +2197,7 @@ function add_to_expression!( ) where {U <: VariableType, V <: PSY.Reserve, W <: AbstractReservesFormulation} contributing_devices_map = get_contributing_devices_map(model) for (device_type, devices) in contributing_devices_map - device_model = get(devices_template, Symbol(device_type), nothing) + device_model = get(devices_template, nameof(device_type), nothing) isnothing(device_model) && continue expression_type = get_expression_type_for_reserve(U(), device_type, V) add_to_expression!(container, expression_type, U, devices, model) From 4c6014c46baa4ec3dd9ed729771b85fccea65861 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Sun, 31 May 2026 17:00:26 +0100 Subject: [PATCH 11/18] fix issue with the Symbol calls --- src/core/device_model.jl | 2 +- src/operation/problem_template.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core/device_model.jl b/src/core/device_model.jl index 5d9f7c3ed..4b71cc144 100644 --- a/src/core/device_model.jl +++ b/src/core/device_model.jl @@ -149,7 +149,7 @@ function _set_model!( dict::Dict, model::DeviceModel{D, B}, ) where {D <: PSY.Device, B <: AbstractDeviceFormulation} - key = Symbol(D) + key = Symbol(IS.strip_module_name(D)) if haskey(dict, key) @warn "Overwriting $(D) existing model" end diff --git a/src/operation/problem_template.jl b/src/operation/problem_template.jl index 3fbac9596..32c5bde00 100644 --- a/src/operation/problem_template.jl +++ b/src/operation/problem_template.jl @@ -67,9 +67,9 @@ end function get_model(template::ProblemTemplate, ::Type{T}) where {T <: PSY.Device} if T <: PSY.Branch - return get(template.branches, Symbol(T), nothing) + return get(template.branches, Symbol(IS.strip_module_name(T)), nothing) elseif T <: PSY.Device - return get(template.devices, Symbol(T), nothing) + return get(template.devices, Symbol(IS.strip_module_name(T)), nothing) else error("Component $T not present in the template") end From ab6aa09837c7d6bf7ffde7d288d6dd674ff84f54 Mon Sep 17 00:00:00 2001 From: m-bossart Date: Mon, 1 Jun 2026 15:56:42 -0700 Subject: [PATCH 12/18] populate branch maps before deepcopy --- src/core/network_model.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/core/network_model.jl b/src/core/network_model.jl index c0f185259..5bb071194 100644 --- a/src/core/network_model.jl +++ b/src/core/network_model.jl @@ -451,13 +451,8 @@ function instantiate_network_model!( if isempty(model.subnetworks) model.subnetworks = _make_subnetworks_from_subnetwork_axes(ybus) end + PNM.populate_branch_maps_by_type!(PNM.get_network_reduction_data(ybus), _get_filters(branch_models)) model.network_reduction = deepcopy(PNM.get_network_reduction_data(ybus)) - #if !isempty(model.network_reductionget_net_reduction_data) - # TODO: Network reimplement this when it becomes necessary. We don't have any - # reductions that are incompatible right now. - # check_network_reduction_compatibility(T) - #end - PNM.populate_branch_maps_by_type!(model.network_reduction, _get_filters(branch_models)) empty!(model.reduced_branch_tracker) set_number_of_steps!(model.reduced_branch_tracker, number_of_steps) return @@ -663,6 +658,9 @@ function instantiate_network_model!( () -> _build_virtual_ptdf(model, sys, irreducible_buses), ) _validate_ptdf_reduction_flags(model) + PNM.populate_branch_maps_by_type!( + PNM.get_network_reduction_data(model.PTDF_matrix), _get_filters(branch_models), + ) model.network_reduction = deepcopy(PNM.get_network_reduction_data(model.PTDF_matrix)) _set_subnetworks_from_ptdf!(model, sys) if _template_uses_security_constrained_branch(branch_models) @@ -678,7 +676,6 @@ function instantiate_network_model!( _ensure_ptdf_modf_consistent!(model, sys) _consolidate_device_model_outages_with_modf!(branch_models, get_MODF_matrix(model)) end - PNM.populate_branch_maps_by_type!(model.network_reduction, _get_filters(branch_models)) empty!(model.reduced_branch_tracker) set_number_of_steps!(model.reduced_branch_tracker, number_of_steps) return From 994713655b4b4920304b799b4ae5e5726f153e07 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Tue, 2 Jun 2026 00:41:49 +0100 Subject: [PATCH 13/18] fix reg device --- src/PowerSimulations.jl | 3 - src/core/network_model.jl | 5 +- .../regulationdevice_constructor.jl | 142 -------- .../devices/regulation_device.jl | 249 -------------- src/operation/template_validation.jl | 5 + src/services_models/agc.jl | 312 ------------------ src/utils/powersystems_utils.jl | 13 - test/test_problem_template.jl | 128 ++++++- 8 files changed, 136 insertions(+), 721 deletions(-) delete mode 100644 src/devices_models/device_constructors/regulationdevice_constructor.jl delete mode 100644 src/devices_models/devices/regulation_device.jl delete mode 100644 src/services_models/agc.jl diff --git a/src/PowerSimulations.jl b/src/PowerSimulations.jl index 72df0a6e7..f27fe38c1 100644 --- a/src/PowerSimulations.jl +++ b/src/PowerSimulations.jl @@ -687,10 +687,8 @@ include("devices_models/devices/TwoTerminalDC_branches.jl") include("devices_models/devices/HVDCsystems.jl") include("devices_models/devices/source.jl") include("devices_models/devices/reactivepower_device.jl") -#include("devices_models/devices/regulation_device.jl") # Services Models -#include("services_models/agc.jl") include("services_models/reserves.jl") include("services_models/reserve_group.jl") include("services_models/transmission_interface.jl") @@ -716,7 +714,6 @@ include("devices_models/device_constructors/renewablegeneration_constructor.jl") include("devices_models/device_constructors/load_constructor.jl") include("devices_models/device_constructors/source_constructor.jl") include("devices_models/device_constructors/reactivepowerdevice_constructor.jl") -#include("devices_models/device_constructors/regulationdevice_constructor.jl") # Network constructors include("network_models/network_constructor.jl") diff --git a/src/core/network_model.jl b/src/core/network_model.jl index 5bb071194..e5dd4c056 100644 --- a/src/core/network_model.jl +++ b/src/core/network_model.jl @@ -451,7 +451,10 @@ function instantiate_network_model!( if isempty(model.subnetworks) model.subnetworks = _make_subnetworks_from_subnetwork_axes(ybus) end - PNM.populate_branch_maps_by_type!(PNM.get_network_reduction_data(ybus), _get_filters(branch_models)) + PNM.populate_branch_maps_by_type!( + PNM.get_network_reduction_data(ybus), + _get_filters(branch_models), + ) model.network_reduction = deepcopy(PNM.get_network_reduction_data(ybus)) empty!(model.reduced_branch_tracker) set_number_of_steps!(model.reduced_branch_tracker, number_of_steps) diff --git a/src/devices_models/device_constructors/regulationdevice_constructor.jl b/src/devices_models/device_constructors/regulationdevice_constructor.jl deleted file mode 100644 index da3231a2b..000000000 --- a/src/devices_models/device_constructors/regulationdevice_constructor.jl +++ /dev/null @@ -1,142 +0,0 @@ -function construct_device!( - container::OptimizationContainer, - sys::PSY.System, - ::ArgumentConstructStage, - model::DeviceModel{PSY.RegulationDevice{T}, U}, - network_model::NetworkModel{AreaBalancePowerModel}, -) where {T <: PSY.StaticInjection, U <: DeviceLimitedRegulation} - devices = get_available_components(get_component_type(model), sys) - add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model) - - add_to_expression!( - container, - ActivePowerBalance, - ActivePowerTimeSeriesParameter, - devices, - model, - network_model, - ) - - add_variables!(container, DeltaActivePowerUpVariable, devices, U()) - add_variables!(container, DeltaActivePowerDownVariable, devices, U()) - add_variables!(container, AdditionalDeltaActivePowerUpVariable, devices, U()) - add_variables!(container, AdditionalDeltaActivePowerDownVariable, devices, U()) - return -end - -function construct_device!( - container::OptimizationContainer, - sys::PSY.System, - ::ModelConstructStage, - model::DeviceModel{PSY.RegulationDevice{T}, DeviceLimitedRegulation}, - network_model::NetworkModel{AreaBalancePowerModel}, -) where {T <: PSY.StaticInjection} - devices = get_available_components(get_component_type(model), sys) - - add_constraints!( - container, - RegulationLimitsConstraint, - DeltaActivePowerUpVariable, - devices, - model, - network_model, - ) - - add_constraints!(container, RampLimitConstraint, devices, model, network_model) - add_constraints!( - container, - ParticipationAssignmentConstraint, - devices, - model, - network_model, - ) - objective_function!(container, devices, model) - return -end - -function construct_device!( - container::OptimizationContainer, - sys::PSY.System, - ::ArgumentConstructStage, - model::DeviceModel{PSY.RegulationDevice{T}, U}, - network_model::NetworkModel{AreaBalancePowerModel}, -) where {T <: PSY.StaticInjection, U <: ReserveLimitedRegulation} - devices = get_available_components(get_component_type(model), sys) - add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model) - - add_to_expression!( - container, - ActivePowerBalance, - ActivePowerTimeSeriesParameter, - devices, - model, - network_model, - ) - - add_variables!(container, DeltaActivePowerUpVariable, devices, U()) - add_variables!(container, DeltaActivePowerDownVariable, devices, U()) - add_variables!(container, AdditionalDeltaActivePowerUpVariable, devices, U()) - add_variables!(container, AdditionalDeltaActivePowerDownVariable, devices, U()) - return -end - -function construct_device!( - container::OptimizationContainer, - sys::PSY.System, - ::ModelConstructStage, - model::DeviceModel{PSY.RegulationDevice{T}, ReserveLimitedRegulation}, - network_model::NetworkModel{AreaBalancePowerModel}, -) where {T <: PSY.StaticInjection} - devices = get_available_components(get_component_type(model), sys) - - add_constraints!( - container, - RegulationLimitsConstraint, - DeltaActivePowerUpVariable, - devices, - model, - network_model, - ) - - add_constraints!( - container, - ParticipationAssignmentConstraint, - devices, - model, - AreaBalancePowerModel, - ) - objective_function!(container, devices, model) - return -end - -function construct_device!( - container::OptimizationContainer, - sys::PSY.System, - ::ArgumentConstructStage, - model::DeviceModel{PSY.RegulationDevice{<:PSY.StaticInjection}, FixedOutput}, - network_model::NetworkModel{AreaBalancePowerModel}, -) - devices = get_available_components(get_component_type(model), sys) - add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model) - - add_to_expression!( - container, - ActivePowerBalance, - ActivePowerTimeSeriesParameter, - devices, - model, - network_model, - ) - return -end - -function construct_device!( - ::OptimizationContainer, - ::PSY.System, - ::ModelConstructStage, - ::DeviceModel{PSY.RegulationDevice{<:PSY.StaticInjection}, FixedOutput}, - network_model::NetworkModel{AreaBalancePowerModel}, -) - # There is no-op under FixedOutput formulation - return -end diff --git a/src/devices_models/devices/regulation_device.jl b/src/devices_models/devices/regulation_device.jl deleted file mode 100644 index dbc4dc948..000000000 --- a/src/devices_models/devices/regulation_device.jl +++ /dev/null @@ -1,249 +0,0 @@ -#! format: off -get_variable_multiplier(_, ::Type{PSY.RegulationDevice{PSY.ThermalStandard}}, ::DeviceLimitedRegulation) = NaN -############################ DeltaActivePowerUpVariable, RegulationDevice ########################### - -get_variable_binary(::DeltaActivePowerUpVariable, ::Type{<:PSY.RegulationDevice}, ::AbstractRegulationFormulation) = false -get_variable_lower_bound(::DeltaActivePowerUpVariable, ::PSY.RegulationDevice, ::AbstractRegulationFormulation) = 0.0 - -############################ DeltaActivePowerDownVariable, RegulationDevice ########################### - -get_variable_binary(::DeltaActivePowerDownVariable, ::Type{<:PSY.RegulationDevice}, ::AbstractRegulationFormulation) = false -get_variable_lower_bound(::DeltaActivePowerDownVariable, ::PSY.RegulationDevice, ::AbstractRegulationFormulation) = 0.0 - -############################ AdditionalDeltaActivePowerUpVariable, RegulationDevice ########################### - -get_variable_binary(::AdditionalDeltaActivePowerUpVariable, ::Type{<:PSY.RegulationDevice}, ::AbstractRegulationFormulation) = false -get_variable_lower_bound(::AdditionalDeltaActivePowerUpVariable, ::PSY.RegulationDevice, ::AbstractRegulationFormulation) = 0.0 - -############################ AdditionalDeltaActivePowerDownVariable, RegulationDevice ########################### - -get_variable_binary(::AdditionalDeltaActivePowerDownVariable, ::Type{<:PSY.RegulationDevice}, ::AbstractRegulationFormulation) = false -get_variable_lower_bound(::AdditionalDeltaActivePowerDownVariable, ::PSY.RegulationDevice, ::AbstractRegulationFormulation) = 0.0 - -get_multiplier_value(::ActivePowerTimeSeriesParameter, d::PSY.RegulationDevice, _) = PSY.get_max_active_power(d) - -########################Objective Function################################################## -proportional_cost(::PSY.OperationalCost, ::AdditionalDeltaActivePowerUpVariable, d::PSY.RegulationDevice, ::AbstractRegulationFormulation) = isapprox(PSY.get_participation_factor(d).up, 0.0; atol=1e-2) ? SERVICES_SLACK_COST : 1 / PSY.get_participation_factor(d).up -proportional_cost(::PSY.OperationalCost, ::AdditionalDeltaActivePowerDownVariable, d::PSY.RegulationDevice, ::AbstractRegulationFormulation) = isapprox(PSY.get_participation_factor(d).dn, 0.0; atol=1e-2) ? SERVICES_SLACK_COST : 1 / PSY.get_participation_factor(d).dn - -objective_function_multiplier(::VariableType, ::AbstractRegulationFormulation)=OBJECTIVE_FUNCTION_POSITIVE - -#! format: on - -function get_default_time_series_names( - ::Type{<:PSY.RegulationDevice{T}}, - ::Type{<:AbstractRegulationFormulation}, -) where {T <: PSY.StaticInjection} - return Dict{Type{<:TimeSeriesParameter}, String}( - ActivePowerTimeSeriesParameter => "max_active_power", - ) -end - -function get_default_attributes( - ::Type{<:PSY.RegulationDevice{T}}, - ::Type{<:AbstractRegulationFormulation}, -) where {T <: PSY.StaticInjection} - return Dict{String, Any}() -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{S}, - ::Type{DeltaActivePowerUpVariable}, - devices::IS.FlattenIteratorWrapper{T}, - ::DeviceModel{T, DeviceLimitedRegulation}, - ::NetworkModel{AreaBalancePowerModel}, -) where { - S <: RegulationLimitsConstraint, - T <: PSY.RegulationDevice{U}, -} where {U <: PSY.StaticInjection} - var_up = get_variable(container, DeltaActivePowerUpVariable(), T) - var_dn = get_variable(container, DeltaActivePowerDownVariable(), T) - base_points_param = get_parameter(container, ActivePowerTimeSeriesParameter(), T) - multiplier = get_multiplier_array(base_points_param) - - names = [PSY.get_name(g) for g in devices] - time_steps = get_time_steps(container) - - container_up = - add_constraints_container!(container, S(), U, names, time_steps; meta = "up") - container_dn = - add_constraints_container!(container, S(), U, names, time_steps; meta = "dn") - - for d in devices - name = PSY.get_name(d) - limits = PSY.get_active_power_limits(d) - param = get_parameter_column_refs(base_points_param, name) - for t in time_steps - container_up[name, t] = JuMP.@constraint( - container.JuMPmodel, - var_up[name, t] <= limits.max - param[t] * multiplier[name, t] - ) - container_dn[name, t] = JuMP.@constraint( - container.JuMPmodel, - var_dn[name, t] <= param[t] * multiplier[name, t] - limits.min - ) - end - end - return -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{S}, - ::Type{DeltaActivePowerUpVariable}, - devices::IS.FlattenIteratorWrapper{T}, - ::DeviceModel{T, ReserveLimitedRegulation}, - ::NetworkModel{AreaBalancePowerModel}, -) where { - S <: RegulationLimitsConstraint, - T <: PSY.RegulationDevice{U}, -} where {U <: PSY.StaticInjection} - var_up = get_variable(container, DeltaActivePowerUpVariable(), T) - var_dn = get_variable(container, DeltaActivePowerDownVariable(), T) - - names = [PSY.get_name(g) for g in devices] - time_steps = get_time_steps(container) - - container_up = - add_constraints_container!(container, S(), U, names, time_steps; meta = "up") - container_dn = - add_constraints_container!(container, S(), U, names, time_steps; meta = "dn") - - for d in devices - name = PSY.get_name(d) - limit_up = PSY.get_reserve_limit_up(d) - limit_dn = PSY.get_reserve_limit_dn(d) - for t in time_steps - container_up[name, t] = - JuMP.@constraint(container.JuMPmodel, var_up[name, t] <= limit_up) - container_dn[name, t] = - JuMP.@constraint(container.JuMPmodel, var_dn[name, t] <= limit_dn) - end - end - return -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{S}, - devices::IS.FlattenIteratorWrapper{T}, - ::DeviceModel{T, DeviceLimitedRegulation}, - ::NetworkModel{AreaBalancePowerModel}, -) where { - S <: RampLimitConstraint, - T <: PSY.RegulationDevice{U}, -} where {U <: PSY.StaticInjection} - R_up = get_variable(container, DeltaActivePowerUpVariable(), T) - R_dn = get_variable(container, DeltaActivePowerDownVariable(), T) - - resolution = Dates.value(Dates.Minute(get_resolution(container))) - names = [PSY.get_name(g) for g in devices] - time_steps = get_time_steps(container) - - container_up = - add_constraints_container!(container, S(), U, names, time_steps; meta = "up") - container_dn = - add_constraints_container!(container, S(), U, names, time_steps; meta = "dn") - - for d in devices - ramp_limits = PSY.get_ramp_limits(d) - isnothing(ramp_limits) && continue - scaling_factor = resolution * SECONDS_IN_MINUTE - name = PSY.get_name(d) - for t in time_steps - container_up[name, t] = JuMP.@constraint( - container.JuMPmodel, - R_up[name, t] <= ramp_limits.up * scaling_factor - ) - container_dn[name, t] = JuMP.@constraint( - container.JuMPmodel, - R_dn[name, t] <= ramp_limits.down * scaling_factor - ) - end - end - return -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{S}, - devices::IS.FlattenIteratorWrapper{T}, - ::DeviceModel{T, <:AbstractRegulationFormulation}, - ::NetworkModel{AreaBalancePowerModel}, -) where { - S <: ParticipationAssignmentConstraint, - T <: PSY.RegulationDevice{U}, -} where {U <: PSY.StaticInjection} - time_steps = get_time_steps(container) - R_up = get_variable(container, DeltaActivePowerUpVariable(), T) - R_dn = get_variable(container, DeltaActivePowerDownVariable(), T) - R_up_emergency = get_variable(container, AdditionalDeltaActivePowerUpVariable(), T) - R_dn_emergency = get_variable(container, AdditionalDeltaActivePowerUpVariable(), T) - area_reserve_up = get_variable(container, DeltaActivePowerUpVariable(), PSY.AGC) - area_reserve_dn = get_variable(container, DeltaActivePowerDownVariable(), PSY.AGC) - - component_names = PSY.get_name.(devices) - participation_assignment_up = add_constraints_container!( - container, - S(), - T, - component_names, - time_steps; - meta = "up", - ) - participation_assignment_dn = add_constraints_container!( - container, - S(), - T, - component_names, - time_steps; - meta = "dn", - ) - - expr_up = get_expression(container, EmergencyUp(), PSY.Area) - expr_dn = get_expression(container, EmergencyDown(), PSY.Area) - for d in devices - name = PSY.get_name(d) - services = PSY.get_services(d) - if length(services) > 1 - device_agc = (a for a in PSY.get_services(d) if isa(a, PSY.AGC)) - agc_name = PSY.get_name.(device_agc)[1] - area_name = PSY.get_name.(PSY.get_area.(device_agc))[1] - else - device_agc = first(services) - agc_name = PSY.get_name(device_agc) - area_name = PSY.get_name(PSY.get_area(device_agc)) - end - p_factor = PSY.get_participation_factor(d) - for t in time_steps - participation_assignment_up[name, t] = JuMP.@constraint( - container.JuMPmodel, - R_up[name, t] == - (p_factor.up * area_reserve_up[agc_name, t]) + R_up_emergency[name, t] - ) - participation_assignment_dn[name, t] = JuMP.@constraint( - container.JuMPmodel, - R_dn[name, t] == - (p_factor.dn * area_reserve_dn[agc_name, t]) + R_dn_emergency[name, t] - ) - JuMP.add_to_expression!(expr_up[area_name, t], -1.0, R_up_emergency[name, t]) - JuMP.add_to_expression!(expr_dn[area_name, t], -1.0, R_dn_emergency[name, t]) - end - end - - return -end - -function objective_function!( - container::OptimizationContainer, - devs::IS.FlattenIteratorWrapper{T}, - ::DeviceModel{T, U}, -) where { - T <: PSY.RegulationDevice{V}, - U <: AbstractRegulationFormulation, -} where {V <: PSY.StaticInjection} - add_proportional_cost!(container, AdditionalDeltaActivePowerUpVariable(), devs, U()) - add_proportional_cost!(container, AdditionalDeltaActivePowerDownVariable(), devs, U()) - return -end diff --git a/src/operation/template_validation.jl b/src/operation/template_validation.jl index f8250a8fa..ab4f501e3 100644 --- a/src/operation/template_validation.jl +++ b/src/operation/template_validation.jl @@ -207,6 +207,11 @@ function validate_template_impl!(model::OperationModel) model_has_branch_filters = false branch_keys_to_delete = Symbol[] + # Rebuilt fresh on every validation pass so it always matches the branch + # types actually present in the template; otherwise stale entries from an + # earlier build (e.g. a branch type later pruned for an empty device cache) + # leak into the PTDF branch_models lookups and raise a KeyError. + empty!(network_model.modeled_ac_branch_types) for (k, device_model) in model.template.branches make_device_cache!(device_model, system, get_check_components(settings)) if isempty(get_device_cache(device_model)) diff --git a/src/services_models/agc.jl b/src/services_models/agc.jl deleted file mode 100644 index 5437cae12..000000000 --- a/src/services_models/agc.jl +++ /dev/null @@ -1,312 +0,0 @@ -#! format: off -get_variable_multiplier(_, ::Type{<:PSY.AGC}, ::AbstractAGCFormulation) = NaN -########################## ActivePowerVariable, AGC ########################### - -########################## SteadyStateFrequencyDeviation ################################## -get_variable_binary(::SteadyStateFrequencyDeviation, ::Type{<:PSY.AGC}, ::AbstractAGCFormulation) = false - -get_variable_binary(::ActivePowerVariable, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = false -########################## SmoothACE, AggregationTopology ########################### - -get_variable_binary(::SmoothACE, ::Type{<:PSY.AggregationTopology}, ::AbstractAGCFormulation) = false -get_variable_binary(::SmoothACE, ::Type{<:PSY.AGC}, ::AbstractAGCFormulation) = false - -########################## DeltaActivePowerUpVariable, AGC ########################### - -get_variable_binary(::DeltaActivePowerUpVariable, ::Type{<:PSY.AGC}, ::AbstractAGCFormulation) = false -get_variable_lower_bound(::DeltaActivePowerUpVariable, ::PSY.AGC, ::AbstractAGCFormulation) = 0.0 - -########################## DeltaActivePowerDownVariable, AGC ########################### - -get_variable_binary(::DeltaActivePowerDownVariable, ::Type{<:PSY.AGC}, ::AbstractAGCFormulation) = false -get_variable_lower_bound(::DeltaActivePowerDownVariable, ::PSY.AGC, ::AbstractAGCFormulation) = 0.0 - -########################## AdditionalDeltaPowerUpVariable, Area ########################### - -get_variable_binary(::AdditionalDeltaActivePowerUpVariable, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = false -get_variable_lower_bound(::AdditionalDeltaActivePowerUpVariable, ::PSY.Area, ::AbstractAGCFormulation) = 0.0 - -########################## AdditionalDeltaPowerDownVariable, Area ########################### - -get_variable_binary(::AdditionalDeltaActivePowerDownVariable, ::Type{<:PSY.Area}, ::AbstractAGCFormulation) = false -get_variable_lower_bound(::AdditionalDeltaActivePowerDownVariable, ::PSY.Area, ::AbstractAGCFormulation) = 0.0 - -########################## AreaMismatchVariable, AGC ########################### -get_variable_binary(::AreaMismatchVariable, ::Type{<:PSY.AGC}, ::AbstractAGCFormulation) = false - -########################## LiftVariable, Area ########################### -get_variable_binary(::LiftVariable, ::Type{<:PSY.AGC}, ::AbstractAGCFormulation) = false -get_variable_lower_bound(::LiftVariable, ::PSY.AGC, ::AbstractAGCFormulation) = 0.0 - -initial_condition_default(::AreaControlError, d::PSY.AGC, ::AbstractAGCFormulation) = PSY.get_initial_ace(d) -initial_condition_variable(::AreaControlError, d::PSY.AGC, ::AbstractAGCFormulation) = AreaMismatchVariable() - -get_variable_multiplier(::SteadyStateFrequencyDeviation, d::PSY.AGC, ::AbstractAGCFormulation) = -10 * PSY.get_bias(d) - -#! format: on - -function get_default_time_series_names( - ::Type{PSY.AGC}, - ::Type{<:AbstractAGCFormulation}, -) - return Dict{Type{<:TimeSeriesParameter}, String}() -end - -function get_default_attributes( - ::Type{PSY.AGC}, - ::Type{<:AbstractAGCFormulation}, -) - return Dict{String, Any}("aggregated_service_model" => false) -end - -""" -Steady State deviation of the frequency -""" -function add_variables!( - container::OptimizationContainer, - ::Type{T}, -) where {T <: SteadyStateFrequencyDeviation} - time_steps = get_time_steps(container) - variable = add_variable_container!(container, T(), PSY.AGC, time_steps) - for t in time_steps - variable[t] = JuMP.@variable(container.JuMPmodel, base_name = "ΔF_{$(t)}") - end -end - -########################## Initial Condition ########################### - -function _get_variable_initial_value( - d::PSY.Component, - key::InitialConditionKey, - ::AbstractAGCFormulation, - ::Nothing, -) - return _get_ace_error(d, key) -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{T}, - ::Type{LiftVariable}, - agcs::IS.FlattenIteratorWrapper{U}, - ::ServiceModel{PSY.AGC, V}, -) where {T <: AbsoluteValueConstraint, U <: PSY.AGC, V <: PIDSmoothACE} - time_steps = get_time_steps(container) - agc_names = PSY.get_name.(agcs) - container_lb = - add_constraints_container!(container, T(), U, agc_names, time_steps; meta = "lb") - container_ub = - add_constraints_container!(container, T(), U, agc_names, time_steps; meta = "ub") - mismatch = get_variable(container, AreaMismatchVariable(), U) - z = get_variable(container, LiftVariable(), U) - jump_model = get_jump_model(container) - - for t in time_steps, a in agc_names - container_lb[a, t] = JuMP.@constraint(jump_model, mismatch[a, t] <= z[a, t]) - container_ub[a, t] = JuMP.@constraint(jump_model, -1 * mismatch[a, t] <= z[a, t]) - end - return -end - -""" -Expression for the power deviation given deviation in the frequency. This expression allows -updating the response of the frequency depending on commitment decisions -""" -function add_constraints!( - container::OptimizationContainer, - ::Type{T}, - ::Type{SteadyStateFrequencyDeviation}, - agcs::IS.FlattenIteratorWrapper{U}, - ::ServiceModel{PSY.AGC, V}, - sys::PSY.System, -) where {T <: FrequencyResponseConstraint, U <: PSY.AGC, V <: PIDSmoothACE} - time_steps = get_time_steps(container) - agc_names = PSY.get_name.(agcs) - - frequency_response = 0.0 - for agc in agcs - area = PSY.get_area(agc) - frequency_response += PSY.get_load_response(area) - end - - for g in PSY.get_components(PSY.get_available, PSY.RegulationDevice, sys) - d = PSY.get_droop(g) - response = 1 / d - frequency_response += response - end - - IS.@assert_op frequency_response >= 0.0 - - # This value is the one updated later in simulation based on the UC result - inv_frequency_response = 1 / frequency_response - - area_balance = get_variable(container, ActivePowerVariable(), PSY.Area) - frequency = get_variable(container, SteadyStateFrequencyDeviation(), U) - R_up = get_variable(container, DeltaActivePowerUpVariable(), U) - R_dn = get_variable(container, DeltaActivePowerDownVariable(), U) - R_up_emergency = - get_variable(container, AdditionalDeltaActivePowerUpVariable(), PSY.Area) - R_dn_emergency = - get_variable(container, AdditionalDeltaActivePowerUpVariable(), PSY.Area) - - const_container = add_constraints_container!(container, T(), PSY.System, time_steps) - - for t in time_steps - system_balance = sum(area_balance.data[:, t]) - for agc in agcs - a = PSY.get_name(agc) - area_name = PSY.get_name(PSY.get_area(agc)) - JuMP.add_to_expression!(system_balance, R_up[a, t]) - JuMP.add_to_expression!(system_balance, -1.0, R_dn[a, t]) - JuMP.add_to_expression!(system_balance, R_up_emergency[area_name, t]) - JuMP.add_to_expression!(system_balance, -1.0, R_dn_emergency[area_name, t]) - end - const_container[t] = JuMP.@constraint( - container.JuMPmodel, - frequency[t] == -inv_frequency_response * system_balance - ) - end - return -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{T}, - ::Type{SteadyStateFrequencyDeviation}, - agcs::IS.FlattenIteratorWrapper{U}, - model::ServiceModel{PSY.AGC, V}, - sys::PSY.System, -) where {T <: SACEPIDAreaConstraint, U <: PSY.AGC, V <: PIDSmoothACE} - services = get_available_components(model, sys) - time_steps = get_time_steps(container) - agc_names = PSY.get_name.(services) - area_names = [PSY.get_name(PSY.get_area(s)) for s in services] - RAW_ACE = get_expression(container, RawACE(), U) - SACE = get_variable(container, SmoothACE(), U) - SACE_pid = add_constraints_container!( - container, - SACEPIDAreaConstraint(), - U, - agc_names, - time_steps, - ) - - jump_model = get_jump_model(container) - for (ix, service) in enumerate(services) - kp = PSY.get_K_p(service) - ki = PSY.get_K_i(service) - kd = PSY.get_K_d(service) - Δt = convert(Dates.Second, get_resolution(container)).value - a = PSY.get_name(service) - for t in time_steps - if t == 1 - ACE_ini = get_initial_condition(container, AreaControlError(), PSY.AGC)[ix] - ace_exp = get_value(ACE_ini) + kp * ((1 + Δt / (kp / ki)) * (RAW_ACE[a, t])) - SACE_pid[a, t] = JuMP.@constraint(jump_model, SACE[a, t] == ace_exp) - continue - end - SACE_pid[a, t] = JuMP.@constraint( - jump_model, - SACE[a, t] == - SACE[a, t - 1] + - kp * ( - (1 + Δt / (kp / ki) + (kd / kp) / Δt) * (RAW_ACE[a, t]) + - (-1 - 2 * (kd / kp) / Δt) * (RAW_ACE[a, t - 1]) - ) - ) - end - end - return -end - -function add_constraints!( - container::OptimizationContainer, - ::Type{T}, - ::Type{SmoothACE}, - agcs::IS.FlattenIteratorWrapper{U}, - ::ServiceModel{PSY.AGC, V}, - sys::PSY.System, -) where {T <: BalanceAuxConstraint, U <: PSY.AGC, V <: PIDSmoothACE} - time_steps = get_time_steps(container) - agc_names = PSY.get_name.(agcs) - aux_equation = add_constraints_container!( - container, - BalanceAuxConstraint(), - PSY.System, - agc_names, - time_steps, - ) - area_mismatch = get_variable(container, AreaMismatchVariable(), PSY.AGC) - SACE = get_variable(container, SmoothACE(), PSY.AGC) - R_up = get_variable(container, DeltaActivePowerUpVariable(), PSY.AGC) - R_dn = get_variable(container, DeltaActivePowerDownVariable(), PSY.AGC) - R_up_emergency = - get_variable(container, AdditionalDeltaActivePowerUpVariable(), PSY.Area) - R_dn_emergency = - get_variable(container, AdditionalDeltaActivePowerUpVariable(), PSY.Area) - - for t in time_steps - for agc in agcs - a = PSY.get_name(agc) - area_name = PSY.get_name(PSY.get_area(agc)) - aux_equation[a, t] = JuMP.@constraint( - container.JuMPmodel, - -1 * SACE[a, t] == - (R_up[a, t] - R_dn[a, t]) + - (R_up_emergency[area_name, t] - R_dn_emergency[area_name, t]) + - area_mismatch[a, t] - ) - end - end - return -end - -function objective_function!( - container::OptimizationContainer, - agcs::IS.FlattenIteratorWrapper{T}, - ::ServiceModel{<:PSY.AGC, U}, -) where {T <: PSY.AGC, U <: PIDSmoothACE} - add_proportional_cost!(container, LiftVariable(), agcs, U()) - return -end - -# Defined here so we can dispatch on PIDSmoothACE -function add_feedforward_arguments!( - container::OptimizationContainer, - model::ServiceModel{PSY.AGC, PIDSmoothACE}, - areas::IS.FlattenIteratorWrapper{PSY.AGC}, -) - for ff in get_feedforwards(model) - @debug "arguments" ff V _group = LOG_GROUP_FEEDFORWARDS_CONSTRUCTION - add_feedforward_arguments!(container, model, areas, ff) - end - return -end - -function add_feedforward_constraints!( - container::OptimizationContainer, - model::ServiceModel{PSY.AGC, PIDSmoothACE}, - areas::IS.FlattenIteratorWrapper{PSY.AGC}, -) - for ff in get_feedforwards(model) - @debug "arguments" ff V _group = LOG_GROUP_FEEDFORWARDS_CONSTRUCTION - add_feedforward_constraints!(container, model, areas, ff) - end - return -end - -function add_proportional_cost!( - container::OptimizationContainer, - ::U, - agcs::IS.FlattenIteratorWrapper{T}, - ::PIDSmoothACE, -) where {T <: PSY.AGC, U <: LiftVariable} - lift_variable = get_variable(container, U(), T) - for index in Iterators.product(axes(lift_variable)...) - add_to_objective_invariant_expression!( - container, - SERVICES_SLACK_COST * lift_variable[index...], - ) - end - return -end diff --git a/src/utils/powersystems_utils.jl b/src/utils/powersystems_utils.jl index df8562fcc..b9d24985c 100644 --- a/src/utils/powersystems_utils.jl +++ b/src/utils/powersystems_utils.jl @@ -91,19 +91,6 @@ function get_available_components( ) end -#= -function get_available_components( - ::Type{PSY.RegulationDevice{T}}, - sys::PSY.System, -) where {T <: PSY.Component} - return PSY.get_components( - x -> (PSY.get_available(x) && PSY.has_service(x, PSY.AGC)), - PSY.RegulationDevice{T}, - sys, - ) -end -=# - make_system_filename(sys::PSY.System) = make_system_filename(IS.get_uuid(sys)) make_system_filename(sys_uuid::Union{Base.UUID, AbstractString}) = "system-$(sys_uuid).json" diff --git a/test/test_problem_template.jl b/test/test_problem_template.jl index 7c7bd0d57..6dfb2e60e 100644 --- a/test/test_problem_template.jl +++ b/test/test_problem_template.jl @@ -15,7 +15,7 @@ end set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment) @test_logs (:warn, "Overwriting ThermalStandard existing model") set_device_model!( template, - DeviceModel(ThermalStandard, ThermalBasicUnitCommitment), + DeviceModel(ThermalStandard, ThermalBasicUnitCommitment) ) @test PSI.get_formulation(template.devices[:ThermalStandard]) == ThermalBasicUnitCommitment @@ -39,3 +39,129 @@ end @test !isempty(ed_template.branches) @test !isempty(ed_template.services) end + +@testset "Template model keys are namespace-stable (issue #1615)" begin + template = ProblemTemplate() + set_device_model!(template, MonitoredLine, StaticBranchBounds) + set_device_model!(template, Line, StaticBranch) + + # Branch keys must be the bare, module-qualifier-free symbols that the build-time + # network lookup (`nameof(branch_type)` in `instantiate_network_model!`) uses. + # `Symbol(D)` qualifies the module path (e.g. `Symbol("PowerSystems.MonitoredLine")`) + # in some namespace contexts, which produced `KeyError(:MonitoredLine)` during build!. + branch_keys = keys(PSI.get_branch_models(template)) + @test Set(branch_keys) == Set([nameof(MonitoredLine), nameof(Line)]) + @test all(k -> !occursin(".", string(k)), branch_keys) + @test PSI.get_model(template, MonitoredLine) isa PSI.DeviceModel + @test PSI.get_model(template, Line) isa PSI.DeviceModel + + # Core invariant that prevents the KeyError regardless of namespace visibility: + # the key the setter stores equals the one the network model looks up. + for D in (MonitoredLine, Line) + @test Symbol(IS.strip_module_name(D)) == nameof(D) + end +end + +@testset "_copy_template_for_build preserves the caller's template across build!" begin + # `_copy_template_for_build` shares PTDF/MODF matrices by reference (their + # libklu/libSparse factorizations are unsafe to deepcopy) and deepcopies the + # device/branch/service dicts so `validate_template_impl!` and + # `_build_device_model_outages!` mutate the build copy without leaking back. + # The invariants below catch a future regression in either half: a missing + # `finally` restore would null the caller's matrices, and dropping the + # device-dict deepcopy would leak auto-discovered outages onto the caller. + sys = PSB.build_system(PSITestSystems, "c_sys5") + lines = collect(PSY.get_components(PSY.Line, sys)) + @assert length(lines) >= 2 + outage = PSY.GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = [lines[2]] + ) + PSY.add_supplemental_attribute!(sys, lines[1], outage) + + ptdf = PNM.PTDF(sys) + modf = PNM.VirtualMODF(sys) + template = get_thermal_dispatch_template_network( + NetworkModel(PTDFPowerModel; PTDF_matrix = ptdf, MODF_matrix = modf), + ) + set_device_model!(template, DeviceModel(PSY.Line, SecurityConstrainedStaticBranch)) + + caller_network = PSI.get_network_model(template) + caller_line_model = PSI.get_model(template, PSY.Line) + caller_device_keys = Set(keys(PSI.get_device_models(template))) + caller_branch_keys = Set(keys(PSI.get_branch_models(template))) + @assert isempty(PSI.get_outages(caller_line_model)) + + model = DecisionModel(template, sys) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + PSI.ModelBuildStatus.BUILT + + # Caller's matrices survived the try/finally restore (identity, not equality). + @test PSI.get_PTDF_matrix(caller_network) === ptdf + @test PSI.get_MODF_matrix(caller_network) === modf + # Caller's device/branch dicts were not mutated by validate_template_impl!. + @test Set(keys(PSI.get_device_models(template))) == caller_device_keys + @test Set(keys(PSI.get_branch_models(template))) == caller_branch_keys + # Auto-discovered outages stayed on the build copy and did not leak onto the + # caller's DeviceModel. + @test isempty(PSI.get_outages(caller_line_model)) + built_line_model = PSI.get_model(PSI.get_template(model), PSY.Line) + @test !isempty(PSI.get_outages(built_line_model)) + + # Build copy aliases the matrices rather than deepcopying their solver caches. + built_network = PSI.get_network_model(PSI.get_template(model)) + @test PSI.get_PTDF_matrix(built_network) === ptdf + @test PSI.get_MODF_matrix(built_network) === modf + + # Re-building from the same template must succeed — catches a broken + # `finally` that nulls the caller's matrices on the first pass. + model2 = DecisionModel(template, sys) + @test build!(model2; output_dir = mktempdir(; cleanup = true)) == + PSI.ModelBuildStatus.BUILT +end + +@testset "modeled_ac_branch_types is rebuilt per validation pass (issue #1621)" begin + # `validate_template_impl!` populates `network_model.modeled_ac_branch_types` + # from the branch types whose device cache is non-empty. Re-building the same + # model after a branch type is pruned (here HVDC, made unavailable) must NOT + # leave a stale entry in that vector: the PTDF `instantiate_network_model!` + # path indexes `branch_models[nameof(branch_type)]` for every modeled type, + # so a stale `:TwoTerminalGenericHVDCLine` raised `KeyError` on the rebuild. + sys = PSB.build_system(PSITestSystems, "c_sys14_dc"; add_forecasts = true) + ptdf = PNM.PTDF(sys) + + template = ProblemTemplate( + NetworkModel(PTDFPowerModel; use_slacks = false, PTDF_matrix = ptdf), + ) + set_device_model!(template, ThermalStandard, ThermalBasicDispatch) + set_device_model!(template, PowerLoad, StaticPowerLoad) + set_device_model!(template, Line, StaticBranch) + set_device_model!(template, TwoTerminalGenericHVDCLine, HVDCTwoTerminalLossless) + + out = mktempdir(; cleanup = true) + model = DecisionModel( + template, + sys; + name = "issue_1621", + horizon = Hour(1), + initial_time = DateTime(2024, 1, 1, 0) + ) + + # First build models the HVDC line, so its type is recorded. + @test build!(model; output_dir = out) == PSI.ModelBuildStatus.BUILT + built_network = PSI.get_network_model(PSI.get_template(model)) + @test TwoTerminalGenericHVDCLine in built_network.modeled_ac_branch_types + + # Prune the HVDC: its device cache becomes empty on the next validation pass, + # so it is deleted from the template's branches. + for hvdc in PSY.get_components(TwoTerminalGenericHVDCLine, sys) + PSY.set_available!(hvdc, false) + end + + # Rebuilding the same model must succeed (previously KeyError), and the stale + # HVDC entry must be gone from the rebuilt vector. + @test build!(model; output_dir = out) == PSI.ModelBuildStatus.BUILT + rebuilt_network = PSI.get_network_model(PSI.get_template(model)) + @test TwoTerminalGenericHVDCLine ∉ rebuilt_network.modeled_ac_branch_types +end From c5a5df4600d7c39d5a16e71c6039100817e57336 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Tue, 2 Jun 2026 00:59:57 +0100 Subject: [PATCH 14/18] formattet --- test/test_problem_template.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_problem_template.jl b/test/test_problem_template.jl index 6dfb2e60e..a406522d8 100644 --- a/test/test_problem_template.jl +++ b/test/test_problem_template.jl @@ -15,7 +15,7 @@ end set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment) @test_logs (:warn, "Overwriting ThermalStandard existing model") set_device_model!( template, - DeviceModel(ThermalStandard, ThermalBasicUnitCommitment) + DeviceModel(ThermalStandard, ThermalBasicUnitCommitment), ) @test PSI.get_formulation(template.devices[:ThermalStandard]) == ThermalBasicUnitCommitment @@ -76,7 +76,7 @@ end outage = PSY.GeometricDistributionForcedOutage(; mean_time_to_recovery = 10, outage_transition_probability = 0.9999, - monitored_components = [lines[2]] + monitored_components = [lines[2]], ) PSY.add_supplemental_attribute!(sys, lines[1], outage) @@ -145,7 +145,7 @@ end sys; name = "issue_1621", horizon = Hour(1), - initial_time = DateTime(2024, 1, 1, 0) + initial_time = DateTime(2024, 1, 1, 0), ) # First build models the HVDC line, so its type is recorded. From a70fec7ab387a099e92709e5c6a0feb220db8fb4 Mon Sep 17 00:00:00 2001 From: m-bossart Date: Tue, 2 Jun 2026 09:14:44 -0700 Subject: [PATCH 15/18] populate maps and deepcopy after reconciliation --- src/core/network_model.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/core/network_model.jl b/src/core/network_model.jl index e5dd4c056..81577693c 100644 --- a/src/core/network_model.jl +++ b/src/core/network_model.jl @@ -626,8 +626,6 @@ function _ensure_ptdf_modf_consistent!( modf_retained = collect(_retained_buses(PNM.get_network_reduction_data(get_MODF_matrix(model)))) model.PTDF_matrix = _build_virtual_ptdf(model, sys, modf_retained) - model.network_reduction = deepcopy(PNM.get_network_reduction_data(model.PTDF_matrix)) - _set_subnetworks_from_ptdf!(model, sys) _reductions_match(get_PTDF_matrix(model), get_MODF_matrix(model)) && return np = length(_retained_buses(PNM.get_network_reduction_data(get_PTDF_matrix(model)))) nm = length(_retained_buses(PNM.get_network_reduction_data(get_MODF_matrix(model)))) @@ -661,11 +659,6 @@ function instantiate_network_model!( () -> _build_virtual_ptdf(model, sys, irreducible_buses), ) _validate_ptdf_reduction_flags(model) - PNM.populate_branch_maps_by_type!( - PNM.get_network_reduction_data(model.PTDF_matrix), _get_filters(branch_models), - ) - model.network_reduction = deepcopy(PNM.get_network_reduction_data(model.PTDF_matrix)) - _set_subnetworks_from_ptdf!(model, sys) if _template_uses_security_constrained_branch(branch_models) model.MODF_matrix = _resolve_reduced_matrix( get_MODF_matrix(model), @@ -679,6 +672,11 @@ function instantiate_network_model!( _ensure_ptdf_modf_consistent!(model, sys) _consolidate_device_model_outages_with_modf!(branch_models, get_MODF_matrix(model)) end + PNM.populate_branch_maps_by_type!( + PNM.get_network_reduction_data(model.PTDF_matrix), _get_filters(branch_models), + ) + model.network_reduction = deepcopy(PNM.get_network_reduction_data(model.PTDF_matrix)) + _set_subnetworks_from_ptdf!(model, sys) empty!(model.reduced_branch_tracker) set_number_of_steps!(model.reduced_branch_tracker, number_of_steps) return From d531158f27f0a4c51045ddddee047610ef53a5c0 Mon Sep 17 00:00:00 2001 From: m-bossart Date: Tue, 2 Jun 2026 12:38:38 -0700 Subject: [PATCH 16/18] add types --- src/core/network_model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/network_model.jl b/src/core/network_model.jl index 81577693c..4d6478f2d 100644 --- a/src/core/network_model.jl +++ b/src/core/network_model.jl @@ -213,7 +213,7 @@ bus/arc numbering and indexes `modf[arc, outage_spec]`, so the PTDF (which dimensions the nodal balance) and the MODF must agree on the entire map, not just the retained-bus key set. """ -function _reductions_match(a, b) +function _reductions_match(a::PNM.PowerNetworkMatrix, b::PNM.VirtualMODF) return PNM.get_bus_reduction_map(PNM.get_network_reduction_data(a)) == PNM.get_bus_reduction_map(PNM.get_network_reduction_data(b)) end From 24dbed4f41b2d594b29fbb794e50491465707d1a Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Thu, 4 Jun 2026 16:02:22 +0100 Subject: [PATCH 17/18] fix typo --- src/devices_models/devices/common/add_to_expression.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/devices_models/devices/common/add_to_expression.jl b/src/devices_models/devices/common/add_to_expression.jl index 362e89e1a..9b02c91bd 100644 --- a/src/devices_models/devices/common/add_to_expression.jl +++ b/src/devices_models/devices/common/add_to_expression.jl @@ -2077,7 +2077,7 @@ function _reduced_entry_in_interface( ArgumentError( "An interface is specified with only part of a double-circuit that has been reduced. Branches: $(branch_names[in_interface]) are in the interface and branches: $(branch_names[.!in_interface]) are not. - Modify the data to include all of or none of the parallel segements.", + Modify the data to include all of or none of the parallel segments.", ), ) end From 699532631b04cd2a54cc589e230bdea8cb560702 Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Fri, 5 Jun 2026 11:32:20 +0100 Subject: [PATCH 18/18] remove subnetworks kwarg --- src/core/network_model.jl | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/core/network_model.jl b/src/core/network_model.jl index 4d6478f2d..a046d76c1 100644 --- a/src/core/network_model.jl +++ b/src/core/network_model.jl @@ -35,9 +35,6 @@ Establishes the NetworkModel for a given PowerModels formulation type. Enable radial branch reduction when building network matrices. - `reduce_degree_two_branches::Bool` = false Enable degree-two branch reduction when building network matrices. -- `subnetworks::Dict{Int, Set{Int}}` = Dict() - Optional mapping of reference bus → set of mapped buses. If not provided, - subnetworks are inferred from PTDF/VirtualPTDF or discovered from the system. - `duals::Vector{DataType}` = Vector{DataType}() Constraint types for which duals should be recorded. - `power_flow_evaluation::Union{PFS.PowerFlowEvaluationModel, Vector{PFS.PowerFlowEvaluationModel}}` @@ -54,8 +51,6 @@ Establishes the NetworkModel for a given PowerModels formulation type. ptdf = PNM.VirtualPTDF(system) nw = NetworkModel(PTDFPowerModel; PTDF_matrix = ptdf, reduce_radial_branches = true, power_flow_evaluation = PFS.PowerFlowEvaluationModel()) - -nw2 = NetworkModel(CopperPlatePowerModel; subnetworks = Dict(1 => Set([1,2,3]))) """ mutable struct NetworkModel{T <: PM.AbstractPowerModel} use_slacks::Bool @@ -80,7 +75,6 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel} MODF_matrix = nothing, reduce_radial_branches = false, reduce_degree_two_branches = false, - subnetworks = Dict{Int, Set{Int}}(), duals = Vector{DataType}(), power_flow_evaluation::Union{ PFS.PowerFlowEvaluationModel, @@ -93,7 +87,7 @@ mutable struct NetworkModel{T <: PM.AbstractPowerModel} use_slacks, PTDF_matrix, MODF_matrix, - subnetworks, + Dict{Int, Set{Int}}(), Dict{PSY.ACBus, Int}(), duals, PNM.NetworkReductionData(), @@ -597,7 +591,9 @@ function _set_subnetworks_from_ptdf!( model::NetworkModel{<:AbstractPTDFModel}, sys::PSY.System, ) - model.subnetworks = _make_subnetworks_from_subnetwork_axes(model.PTDF_matrix) + if isempty(model.subnetworks) + model.subnetworks = _make_subnetworks_from_subnetwork_axes(model.PTDF_matrix) + end if length(model.subnetworks) > 1 @debug "System Contains Multiple Subnetworks. Assigning buses to subnetworks." _assign_subnetworks_to_buses(model, sys)