From 7871cc286ba2ad46ef2851fb162f91a6381ad27c Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Mon, 15 Jun 2026 13:28:36 -0400 Subject: [PATCH 01/11] Add time-varying ORDC (ReserveDemandTimeSeriesCurve) support Port of PowerSimulations.jl#1629 (time-varying ORDC), adapted to the psy6 data model where a time-varying ORDC is a `ReserveDemandTimeSeriesCurve` whose `variable` is a `CostCurve{TimeSeriesPiecewiseIncrementalCurve}` (rather than a `ReserveDemandCurve` carrying a bare `TimeSeriesKey`). - add_parameters.jl: service-level piecewise slope/breakpoint cost-parameter population for `ReserveDemandTimeSeriesCurve` (calc_additional_axes, get_max_tranches, _get_time_series_name, _add_parameters!, _param_to_vars). - reserves.jl / market_bid_overrides.jl: drop the "not supported" error and route ORDC cost through the decremental block-offer path, reading per-timestep data from the parameter arrays (time-varying) or the curve (static); time-varying cost goes to the variant objective expression. - services_constructor.jl: route time-series-backed ORDC to parameter setup. - Broaden shared reserve / production-cost-expression methods to Union{ReserveDemandCurve, ReserveDemandTimeSeriesCurve}; update hybrid/storage disambiguators to keep dispatch unambiguous. - test/test_services_constructor.jl: build, and build & solve, time-varying ORDC. - Fix a stray `end` in hybrid_systems.jl that blocked precompilation. - Point the IOM source at the ac/ordc branch (companion changes live there). Co-Authored-By: Claude Opus 4.8 (1M context) --- Project.toml | 2 +- src/common_models/add_parameters.jl | 164 ++++++++++++++++++-- src/common_models/add_to_expression.jl | 10 +- src/common_models/market_bid_overrides.jl | 90 ++++++++--- src/energy_storage_models/storage_models.jl | 2 +- src/hybrid_system_models/hybrid_systems.jl | 5 +- src/services_models/reserves.jl | 62 ++++++-- src/services_models/services_constructor.jl | 10 ++ test/Project.toml | 2 +- test/test_services_constructor.jl | 144 +++++++++++++++++ 10 files changed, 437 insertions(+), 54 deletions(-) create mode 100644 test/test_services_constructor.jl diff --git a/Project.toml b/Project.toml index ed92579..bf40441 100644 --- a/Project.toml +++ b/Project.toml @@ -27,7 +27,7 @@ PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6" [sources] InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl"} PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} -InfrastructureOptimizationModels = {rev = "main", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} +InfrastructureOptimizationModels = {rev = "ac/ordc", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} PowerNetworkMatrices = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"} [extensions] diff --git a/src/common_models/add_parameters.jl b/src/common_models/add_parameters.jl index e962e41..15fbc1c 100644 --- a/src/common_models/add_parameters.jl +++ b/src/common_models/add_parameters.jl @@ -45,6 +45,29 @@ function add_parameters!( return end +# Time-varying ORDC: piecewise slope/breakpoint cost parameters for a +# `ReserveDemandTimeSeriesCurve` service. +function add_parameters!( + container::OptimizationContainer, + ::Type{T}, + service::U, + model::ServiceModel{U, V}, +) where { + T <: Union{ + AbstractPiecewiseLinearSlopeParameter, + AbstractPiecewiseLinearBreakpointParameter, + }, + U <: PSY.ReserveDemandTimeSeriesCurve, + V <: AbstractServiceFormulation, +} + if get_rebuild_model(get_settings(container)) && + has_container_key(container, T, U, PSY.get_name(service)) + return + end + _add_parameters!(container, T, service, model) + return +end + function add_branch_parameters!( container::OptimizationContainer, ::Type{T}, @@ -210,7 +233,7 @@ function _add_time_series_parameters!( # parameter axis, so iterating `initial_values` in order matches the container's # first axis row-by-row. for (i, (ts_uuid, raw_ts_vals)) in enumerate(initial_values) - ts_vals = _unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) + ts_vals = IOM._unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) @assert all(_size_wrapper.(ts_vals) .== (length.(additional_axes),)) for step in time_steps @@ -330,7 +353,7 @@ function _add_time_series_parameters!( interval = ts_interval, ) ts_vals = - _unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) + IOM._unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) @assert all(_size_wrapper.(ts_vals) .== (length.(additional_axes),)) end # `_resolve_branch_multiplier` is `DeviceModel`-aware: it applies the @@ -379,6 +402,19 @@ _get_time_series_name( ) where {T <: ParameterType} = get_time_series_names(model)[T] +# Time-varying ORDC: the slope/breakpoint parameters read from the curve time +# series referenced by the service's `variable` field. +_get_time_series_name( + ::Type{ + <:Union{ + DecrementalPiecewiseLinearSlopeParameter, + DecrementalPiecewiseLinearBreakpointParameter, + }, + }, + service::PSY.ReserveDemandTimeSeriesCurve, + ::ServiceModel, +) = IS.get_name(IS.get_time_series_key(PSY.get_variable(service))) + # The fact that we're seeing these parameters means that we should # have a time-varying MBC/IEC, so the `get_time_series_key` call should be valid. @@ -467,6 +503,15 @@ _param_to_vars( ::Type{<:AbstractDeviceFormulation}, ) = (PiecewiseLinearBlockDecrementalOffer,) +# Time-varying ORDC: the reserve demand curve is offered as a decremental block. +_param_to_vars( + ::Union{ + Type{DecrementalPiecewiseLinearSlopeParameter}, + Type{DecrementalPiecewiseLinearBreakpointParameter}, + }, + ::Type{<:AbstractServiceFormulation}, +) = + (PiecewiseLinearBlockDecrementalOffer,) ################################################################################# # calc_additional_axes — default implementations (no additional axes) @@ -495,16 +540,43 @@ calc_additional_axes( } where {D <: PSY.Service} = () ################################################################################# -# _unwrap_for_param — default implementation (identity) +# Piecewise linear parameter helpers (time-varying ORDC) ################################################################################# -_unwrap_for_param(::ParameterType, ts_elem, expected_axs) = ts_elem +function get_max_tranches(component::PSY.Component, piecewise_ts::IS.TimeSeriesKey) + data = PSY.get_data(PSY.get_time_series(component, piecewise_ts)) + return IOM._get_max_tranches(data) +end -################################################################################# -# Piecewise linear parameter helpers -# NOTE: _unwrap_for_param overloads, get_max_tranches, make_tranche_axis, and -# lookup_additional_axes belong in PSI (multi-timestep update path), not POM. -################################################################################# +function calc_additional_axes( + ::OptimizationContainer, + ::Type{P}, + service::U, + ::ServiceModel{U, W}, +) where { + P <: AbstractPiecewiseLinearSlopeParameter, + U <: PSY.ReserveDemandTimeSeriesCurve, + W <: AbstractServiceFormulation, +} + ts_key = IS.get_time_series_key(PSY.get_variable(service)) + max_tranches = get_max_tranches(service, ts_key) + return (IOM.make_tranche_axis(max_tranches),) +end + +function calc_additional_axes( + ::OptimizationContainer, + ::Type{P}, + service::U, + ::ServiceModel{U, W}, +) where { + P <: AbstractPiecewiseLinearBreakpointParameter, + U <: PSY.ReserveDemandTimeSeriesCurve, + W <: AbstractServiceFormulation, +} + ts_key = IS.get_time_series_key(PSY.get_variable(service)) + max_tranches = get_max_tranches(service, ts_key) + return (IOM.make_tranche_axis(max_tranches + 1),) # one more breakpoint than tranches +end ################################################################################# # _add_parameters! for ObjectiveFunctionParameter @@ -575,7 +647,7 @@ function _add_parameters!( ts_name; interval = ts_interval, ) - ts_vals = _unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) + ts_vals = IOM._unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) @assert all(_size_wrapper.(ts_vals) .== (length.(additional_axes),)) # PWL/cost-function path: parameter values from `_unwrap_for_param` are # tuples-of-floats, while the multiplier itself is a scalar Float64. @@ -591,6 +663,78 @@ function _add_parameters!( return end +################################################################################# +# _add_parameters! for time-varying ORDC slope/breakpoint cost parameters +# +# Mirrors the device `ObjectiveFunctionParameter` path above but for a single +# `ReserveDemandTimeSeriesCurve` service whose cost curve is a time series. Allocates the +# 3-D `(name, tranche, time)` container and fills it from the curve time series. +################################################################################# + +function _add_parameters!( + container::OptimizationContainer, + ::Type{T}, + service::U, + model::ServiceModel{U, V}, +) where { + T <: Union{ + AbstractPiecewiseLinearSlopeParameter, + AbstractPiecewiseLinearBreakpointParameter, + }, + U <: PSY.ReserveDemandTimeSeriesCurve, + V <: AbstractServiceFormulation, +} + ts_type = get_default_time_series_type(container) + if !(ts_type <: Union{PSY.AbstractDeterministic, PSY.StaticTimeSeries}) + error( + "add_parameters! for $T is not compatible with $ts_type", + ) + end + ts_name = _get_time_series_name(T, service, model) + time_steps = get_time_steps(container) + name = PSY.get_name(service) + model_interval = get_interval(get_settings(container)) + ts_interval = model_interval + + additional_axes = calc_additional_axes(container, T, service, model) + param_container = add_param_container!( + container, + T, + U, + _param_to_vars(T, V), + SOSStatusVariable.NO_VARIABLE, + false, + _get_expected_time_series_eltype(T), + [name], + additional_axes..., + time_steps; + meta = name, + ) + # NOTE: ObjectiveFunctionParameter containers carry `CostFunctionAttributes`, which + # (unlike `TimeSeriesAttributes`) have no `set_subsystem!` method — mirror the device + # ObjectiveFunctionParameter path and skip it. + + jump_model = get_jump_model(container) + param_instance = T() + multiplier = get_multiplier_value(T, service, V) + raw_ts_vals = IOM.get_time_series_initial_values!( + container, + ts_type, + service, + ts_name; + interval = ts_interval, + ) + ts_vals = IOM._unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) + @assert all(_size_wrapper.(ts_vals) .== (length.(additional_axes),)) + parent_mult = IOM.get_multiplier_array_data(param_container) + IOM._set_multiplier_at!(parent_mult, Float64(multiplier), 1) + parent_param = IOM.get_parameter_array_data(param_container) + for t in time_steps + IOM._set_parameter_at!(parent_param, jump_model, ts_vals[t], 1, t) + end + return +end + ################################################################################# # _add_parameters! for ServiceModel TimeSeriesParameter ################################################################################# diff --git a/src/common_models/add_to_expression.jl b/src/common_models/add_to_expression.jl index eb203a3..adbef21 100644 --- a/src/common_models/add_to_expression.jl +++ b/src/common_models/add_to_expression.jl @@ -2548,7 +2548,10 @@ function add_cost_to_expression!( cost_expression::JuMP.AbstractJuMPScalar, component::T, time_period::Int, -) where {S <: CostExpressions, T <: PSY.ReserveDemandCurve} +) where { + S <: CostExpressions, + T <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, +} if has_container_key(container, S, T, PSY.get_name(component)) device_cost_expression = get_expression(container, S, T, PSY.get_name(component)) component_name = PSY.get_name(component) @@ -2871,7 +2874,10 @@ function add_to_expression!( cost_expression::JuMP.AbstractJuMPScalar, component::T, time_period::Int, -) where {S <: CostExpressions, T <: PSY.ReserveDemandCurve} +) where { + S <: CostExpressions, + T <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, +} if has_container_key(container, S, T, PSY.get_name(component)) device_cost_expression = get_expression(container, S, T, PSY.get_name(component)) diff --git a/src/common_models/market_bid_overrides.jl b/src/common_models/market_bid_overrides.jl index 4c5a69a..0558946 100644 --- a/src/common_models/market_bid_overrides.jl +++ b/src/common_models/market_bid_overrides.jl @@ -293,7 +293,9 @@ _vom_offer_direction(::Type{<:AbstractControllablePowerLoadFormulation}) = ################################################################################# """ -PWL block offer constraints for ORDC (ReserveDemandCurve). +PWL block offer constraints for ORDC (ReserveDemandCurve). ORDC is offered as a +decremental block (the demand curve is a willingness-to-pay), so it uses the +decremental block-offer variable/constraint family. """ function add_pwl_constraint_delta!( container::OptimizationContainer, @@ -302,12 +304,15 @@ function add_pwl_constraint_delta!( break_points::Vector{Float64}, pwl_vars::Vector{JuMP.VariableRef}, period::Int, -) where {T <: PSY.ReserveDemandCurve, U <: ServiceRequirementVariable} +) where { + T <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, + U <: ServiceRequirementVariable, +} name = PSY.get_name(component) variables = get_variable(container, U, T, name) const_container = lazy_container_addition!( container, - PiecewiseLinearBlockIncrementalOfferConstraint, + PiecewiseLinearBlockDecrementalOfferConstraint, T, axes(variables)...; meta = name, @@ -325,38 +330,83 @@ function add_pwl_constraint_delta!( end """ -PWL cost terms for StepwiseCostReserve (AbstractServiceFormulation). +Return `(breakpoints, slopes)` for an ORDC service at time step `t`. For a static +`ReserveDemandCurve` the data is read directly from the curve; for a time-varying +`ReserveDemandTimeSeriesCurve` it is read from the pre-allocated decremental +slope/breakpoint parameter arrays. Mirrors IOM's device-side `_get_raw_pwl_data`. +""" +function _get_reserve_pwl_data( + container::OptimizationContainer, + component::T, + variable_cost::PSY.CostCurve, + t::Int, +) where {T <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}} + base_power = get_model_base_power(container) + device_base_power = PSY.get_base_power(component, PSY.NU) + # `get_power_units` returns the curve's `IS.AbstractUnitSystem` instance for both the + # static and time-series curves (the parameter arrays store data in those same units). + unit_system = PSY.get_power_units(variable_cost) + + if !IS.is_time_series_backed(variable_cost) + cost_component = PSY.get_function_data(PSY.get_value_curve(variable_cost)) + breakpoint_cost_component = IS.get_x_coords(cost_component) + slope_cost_component = IS.get_y_coords(cost_component) + else + name = PSY.get_name(component) + SlopeParam = DecrementalPiecewiseLinearSlopeParameter + slope_arr = get_parameter_array(container, SlopeParam, T, name) + slope_mult = get_parameter_multiplier_array(container, SlopeParam, T, name) + @assert size(slope_arr) == size(slope_mult) + slope_cost_component = + [slope_arr[name, s, t] * slope_mult[name, s, t] for s in axes(slope_arr)[2]] + + BreakpointParam = DecrementalPiecewiseLinearBreakpointParameter + bp_arr = get_parameter_array(container, BreakpointParam, T, name) + bp_mult = get_parameter_multiplier_array(container, BreakpointParam, T, name) + @assert size(bp_arr) == size(bp_mult) + breakpoint_cost_component = + [bp_arr[name, p, t] * bp_mult[name, p, t] for p in axes(bp_arr)[2]] + + @assert length(slope_cost_component) == length(breakpoint_cost_component) - 1 + end + + breakpoints, slopes = get_piecewise_curve_per_system_unit( + breakpoint_cost_component, + slope_cost_component, + unit_system, + base_power, + device_base_power, + ) + return breakpoints, slopes +end + +""" +PWL cost terms for StepwiseCostReserve (AbstractServiceFormulation). Handles both +static (`ReserveDemandCurve`) and time-varying (`ReserveDemandTimeSeriesCurve`) +reserve demand curves. """ function add_pwl_term_delta!( container::OptimizationContainer, component::T, - cost_data::PSY.CostCurve{PSY.PiecewiseIncrementalCurve}, + cost_data::PSY.CostCurve, ::Type{U}, ::Type{V}, -) where {T <: PSY.Component, U <: VariableType, V <: AbstractServiceFormulation} +) where { + T <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, + U <: VariableType, + V <: AbstractServiceFormulation, +} multiplier = objective_function_multiplier(U, V) resolution = get_resolution(container) dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR - base_power = get_model_base_power(container) - value_curve = PSY.get_value_curve(cost_data) - power_units = PSY.get_power_units(cost_data) - cost_component = PSY.get_function_data(value_curve) - device_base_power = PSY.get_base_power(component, PSY.NU) - data = get_piecewise_curve_per_system_unit( - cost_component, - power_units, - base_power, - device_base_power, - ) name = PSY.get_name(component) time_steps = get_time_steps(container) pwl_cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) - slopes = IS.get_y_coords(data) - break_points = PSY.get_x_coords(data) for t in time_steps + break_points, slopes = _get_reserve_pwl_data(container, component, cost_data, t) pwl_vars = add_pwl_variables_delta!( container, - PiecewiseLinearBlockIncrementalOffer, + PiecewiseLinearBlockDecrementalOffer, T, name, t, diff --git a/src/energy_storage_models/storage_models.jl b/src/energy_storage_models/storage_models.jl index beb6ff2..cf0e7d1 100644 --- a/src/energy_storage_models/storage_models.jl +++ b/src/energy_storage_models/storage_models.jl @@ -52,7 +52,7 @@ end function get_variable_upper_bound(::Type{ActivePowerReserveVariable}, r::PSY.Reserve, d::PSY.Storage, ::Type{<:AbstractReservesFormulation}) return PSY.get_max_output_fraction(r) * (PSY.get_output_active_power_limits(d, PSY.SU).max + PSY.get_input_active_power_limits(d, PSY.SU).max) end -function get_variable_upper_bound(::Type{ActivePowerReserveVariable}, r::PSY.ReserveDemandCurve, d::PSY.Storage, ::Type{<:AbstractReservesFormulation}) +function get_variable_upper_bound(::Type{ActivePowerReserveVariable}, r::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, d::PSY.Storage, ::Type{<:AbstractReservesFormulation}) return PSY.get_max_output_fraction(r) * (PSY.get_output_active_power_limits(d, PSY.SU).max + PSY.get_input_active_power_limits(d, PSY.SU).max) end diff --git a/src/hybrid_system_models/hybrid_systems.jl b/src/hybrid_system_models/hybrid_systems.jl index c61fa56..1878adf 100644 --- a/src/hybrid_system_models/hybrid_systems.jl +++ b/src/hybrid_system_models/hybrid_systems.jl @@ -450,10 +450,10 @@ function get_variable_upper_bound( ) end -# Disambiguate against the generic ReserveDemandCurve method in services_models/reserves.jl. +# Disambiguate against the generic ORDC method in services_models/reserves.jl. function get_variable_upper_bound( ::Type{ActivePowerReserveVariable}, - r::PSY.ReserveDemandCurve, + r::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, d::PSY.HybridSystem, ::Type{<:AbstractReservesFormulation}, ) @@ -1616,7 +1616,6 @@ function _emit_coverage_constraint!( ) end end - end return end diff --git a/src/services_models/reserves.jl b/src/services_models/reserves.jl index 3705c83..27246c6 100644 --- a/src/services_models/reserves.jl +++ b/src/services_models/reserves.jl @@ -7,7 +7,7 @@ get_variable_binary(::Type{ActivePowerReserveVariable}, ::Type{<:PSY.Reserve}, : function get_variable_upper_bound(::Type{ActivePowerReserveVariable}, r::PSY.Reserve, d::PSY.Device, ::Type{<:AbstractReservesFormulation}) return PSY.get_max_output_fraction(r) * PSY.get_max_active_power(d, PSY.SU) end -get_variable_upper_bound(::Type{ActivePowerReserveVariable}, r::PSY.ReserveDemandCurve, d::PSY.Device, ::Type{<:AbstractReservesFormulation}) = PSY.get_max_active_power(d, PSY.SU) +get_variable_upper_bound(::Type{ActivePowerReserveVariable}, r::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, d::PSY.Device, ::Type{<:AbstractReservesFormulation}) = PSY.get_max_active_power(d, PSY.SU) get_variable_lower_bound(::Type{ActivePowerReserveVariable}, ::PSY.Reserve, ::PSY.Device, ::Type) = 0.0 ############################### ActivePowerReserveVariable, ReserveNonSpinning ######################################### @@ -19,9 +19,9 @@ get_variable_lower_bound(::Type{ActivePowerReserveVariable}, ::PSY.ReserveNonSpi ############################### ServiceRequirementVariable, ReserveDemandCurve ################################ -get_variable_binary(::Type{ServiceRequirementVariable}, ::Type{<:PSY.ReserveDemandCurve}, ::Type{<:AbstractReservesFormulation}) = false -get_variable_upper_bound(::Type{ServiceRequirementVariable}, ::PSY.ReserveDemandCurve, d::PSY.Component, ::Type{<:AbstractReservesFormulation}) = PSY.get_max_active_power(d, PSY.SU) -get_variable_lower_bound(::Type{ServiceRequirementVariable}, ::PSY.ReserveDemandCurve, ::PSY.Component, ::Type{<:AbstractReservesFormulation}) = 0.0 +get_variable_binary(::Type{ServiceRequirementVariable}, ::Type{<:Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}}, ::Type{<:AbstractReservesFormulation}) = false +get_variable_upper_bound(::Type{ServiceRequirementVariable}, ::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, d::PSY.Component, ::Type{<:AbstractReservesFormulation}) = PSY.get_max_active_power(d, PSY.SU) +get_variable_lower_bound(::Type{ServiceRequirementVariable}, ::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, ::PSY.Component, ::Type{<:AbstractReservesFormulation}) = 0.0 # `VariableReserve` stores `requirement` as a dimensionless factor that scales its # requirement and needs an explicit unit system (PS6 made the reserve requirement @@ -35,7 +35,9 @@ get_parameter_multiplier(::Type{<:VariableValueParameter}, d::Type{<:PSY.Abstrac get_initial_parameter_value(::Type{<:VariableValueParameter}, d::Type{<:PSY.AbstractReserve}, ::Type{<:AbstractReservesFormulation}) = 0.0 objective_function_multiplier(::Type{ServiceRequirementVariable}, ::Type{StepwiseCostReserve}) = -1.0 -uses_compact_power(::PSY.ReserveDemandCurve, ::StepwiseCostReserve)=false +uses_compact_power(::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, ::StepwiseCostReserve)=false +get_multiplier_value(::Type{<:AbstractPiecewiseLinearBreakpointParameter}, ::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, ::Type{<:AbstractReservesFormulation}) = 1.0 +get_multiplier_value(::Type{<:AbstractPiecewiseLinearSlopeParameter}, ::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, ::Type{<:AbstractReservesFormulation}) = 1.0 #! format: on function get_initial_conditions_service_model( @@ -101,7 +103,7 @@ function add_reserve_variables!( formulation, ) where { T <: ServiceRequirementVariable, - D <: PSY.ReserveDemandCurve, + D <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, } time_steps = get_time_steps(container) service_name = PSY.get_name(service) @@ -312,7 +314,7 @@ function add_constraints!( ::U, ::ServiceModel{SR, StepwiseCostReserve}, ) where { - SR <: PSY.ReserveDemandCurve, + SR <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, } where {D <: PSY.Component} time_steps = get_time_steps(container) @@ -496,9 +498,12 @@ end function add_to_objective_function!( container::OptimizationContainer, - service::PSY.ReserveDemandCurve{T}, - ::ServiceModel{PSY.ReserveDemandCurve{T}, SR}, -) where {T <: PSY.ReserveDirection, SR <: StepwiseCostReserve} + service::S, + ::ServiceModel{S, SR}, +) where { + S <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, + SR <: StepwiseCostReserve, +} add_reserves_variable_cost!(container, ServiceRequirementVariable, service, SR) return end @@ -509,7 +514,11 @@ function add_reserves_variable_cost!( ::Type{U}, service::T, ::Type{V}, -) where {T <: PSY.ReserveDemandCurve, U <: VariableType, V <: StepwiseCostReserve} +) where { + T <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, + U <: VariableType, + V <: StepwiseCostReserve, +} _add_reserves_variable_cost_to_objective!(container, U, service, V) return end @@ -528,14 +537,13 @@ function _add_reserves_variable_cost_to_objective!( variable_cost = PSY.get_variable(component) if variable_cost isa Nothing error("ReserveDemandCurve $(component.name) does not have cost data.") - elseif typeof(variable_cost) <: PSY.TimeSeriesKey - error( - "Timeseries curve for ReserveDemandCurve $(component.name) is not supported yet.", - ) end pwl_cost_expressions = add_pwl_term_delta!(container, component, variable_cost, T, U) + # A time-series-backed curve changes across simulation steps, so its cost goes + # into the variant objective expression; a static curve into the invariant one. + is_t_variant = IS.is_time_series_backed(variable_cost) for t in time_steps add_to_expression!( container, @@ -544,7 +552,29 @@ function _add_reserves_variable_cost_to_objective!( component, t, ) - add_to_objective_invariant_expression!(container, pwl_cost_expressions[t]) + if is_t_variant + IOM.add_to_objective_variant_expression!(container, pwl_cost_expressions[t]) + else + add_to_objective_invariant_expression!(container, pwl_cost_expressions[t]) + end + end + return +end + +""" +Add the decremental piecewise slope/breakpoint cost parameters for a time-varying +ORDC (`ReserveDemandTimeSeriesCurve`) service. +""" +function process_stepwise_cost_reserve_parameters!( + container::OptimizationContainer, + model::ServiceModel, + service::D, +) where {D <: PSY.ReserveDemandTimeSeriesCurve} + for param in ( + DecrementalPiecewiseLinearBreakpointParameter, + DecrementalPiecewiseLinearSlopeParameter, + ) + add_parameters!(container, param, service, model) end return end diff --git a/src/services_models/services_constructor.jl b/src/services_models/services_constructor.jl index 08507ea..aa20a4e 100644 --- a/src/services_models/services_constructor.jl +++ b/src/services_models/services_constructor.jl @@ -193,6 +193,15 @@ function construct_service!( return end +# Add the time-varying ORDC cost parameters only when the demand curve is backed by +# a time series (a `ReserveDemandTimeSeriesCurve`, whose `variable` is a +# time-series-backed `CostCurve`). Static `ReserveDemandCurve` services are a no-op. +function _maybe_process_stepwise(container, model, service) + IS.is_time_series_backed(PSY.get_variable(service)) && + process_stepwise_cost_reserve_parameters!(container, model, service) + return +end + function construct_service!( container::OptimizationContainer, sys::PSY.System, @@ -206,6 +215,7 @@ function construct_service!( service = PSY.get_component(SR, sys, name) !PSY.get_available(service) && return contributing_devices = get_contributing_devices(model) + _maybe_process_stepwise(container, model, service) add_reserve_variables!( container, ServiceRequirementVariable, diff --git a/test/Project.toml b/test/Project.toml index 7a61d08..85136d9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -34,7 +34,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [sources] -InfrastructureOptimizationModels = {rev = "main", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} +InfrastructureOptimizationModels = {rev = "ac/ordc", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl"} PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} PowerFlows = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerFlows.jl"} diff --git a/test/test_services_constructor.jl b/test/test_services_constructor.jl new file mode 100644 index 0000000..4a50e90 --- /dev/null +++ b/test/test_services_constructor.jl @@ -0,0 +1,144 @@ +# Time-varying ORDC. Adapted from PowerSimulations.jl PR #1629 to the psy6 data model, +# where a time-varying ORDC is a `ReserveDemandTimeSeriesCurve` whose `variable` is a +# `CostCurve{TimeSeriesPiecewiseIncrementalCurve}` (rather than a `ReserveDemandCurve` +# carrying a bare `TimeSeriesKey`). The multi-step Simulation scenario from that PR is +# omitted: the Simulation framework is not yet available in the IOM/POM split. + +# Build a time-varying ORDC (`ReserveDemandTimeSeriesCurve`) from the system's existing +# static ORDC baseline curve, backing it with a deterministic cost-curve forecast. +function _add_ts_ordc!( + sys, + name::String, + static_ordc; + incrs_x = (0.0, 0.0, 0.0), + incrs_y = (0.0, 0.0, 0.0), + create_extra_tranches = false, +) + baseline_curve = PSY.get_variable(static_ordc) + power_units = PSY.get_power_units(baseline_curve) + fd = PSY.get_function_data(PSY.get_value_curve(baseline_curve)) + + # Construct with a stub TS curve so the component can be added; the real forecast is + # attached and set below. + stub = stub_ts_offer_curve(; power_units = power_units) + ordc_ts = ReserveDemandTimeSeriesCurve{ReserveUp}( + stub, + name, + true, + PSY.get_time_frame(static_ordc), + ) + add_service!(sys, ordc_ts, get_components(ThermalStandard, sys)) + + pwl_ts = make_deterministic_ts( + sys, + "variable_cost", + fd, + incrs_x, + incrs_y; + override_min_x = 0.0, + override_max_x = last(get_x_coords(fd)), + create_extra_tranches = create_extra_tranches, + ) + pwl_key = add_time_series!(sys, ordc_ts, pwl_ts) + PSY.set_variable!(ordc_ts, PSY.make_market_bid_ts_curve(pwl_key, nothing, power_units)) + return ordc_ts +end + +@testset "Test ORDC time series (build)" begin + c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + static_ordc = first(get_components(PSY.ReserveDemandCurve, c_sys5_uc)) + _add_ts_ordc!(c_sys5_uc, "ORDC_TS", static_ordc) + + template = get_thermal_standard_uc_template() + set_service_model!( + template, + ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve1"), + ) + set_service_model!( + template, + ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve11"), + ) + set_service_model!( + template, + ServiceModel(VariableReserve{ReserveDown}, RangeReserve, "Reserve2"), + ) + set_service_model!( + template, + ServiceModel( + ReserveDemandTimeSeriesCurve{ReserveUp}, + StepwiseCostReserve, + "ORDC_TS", + ), + ) + model = DecisionModel( + template, + c_sys5_uc; + store_variable_names = true, + optimizer = HiGHS_optimizer, + ) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT +end + +@testset "Test ORDC time series (build & solve)" begin + c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + static_ordc = first(get_components(PSY.ReserveDemandCurve, c_sys5_uc)) + # Two time-varying ORDCs with different per-timestep tranche counts to exercise the + # tranche-axis padding and per-service (meta-keyed) parameter containers. + _add_ts_ordc!( + c_sys5_uc, + "ORDC_TS1", + static_ordc; + incrs_x = (0.03, 0.13, 0.07), + incrs_y = (0.03, 0.13, 0.07), + create_extra_tranches = true, + ) + _add_ts_ordc!( + c_sys5_uc, + "ORDC_TS2", + static_ordc; + incrs_x = (0.03, 0.13, 0.07), + incrs_y = (0.02, 0.14, 0.08), + create_extra_tranches = true, + ) + + template = PowerOperationsProblemTemplate( + NetworkModel(CopperPlatePowerModel; use_slacks = true), + ) + set_device_model!(template, ThermalStandard, ThermalDispatchNoMin) + set_device_model!(template, PowerLoad, StaticPowerLoad) + set_service_model!( + template, + ServiceModel(VariableReserve{ReserveUp}, RangeReserve, "Reserve1"), + ) + set_service_model!( + template, + ServiceModel(VariableReserve{ReserveDown}, RangeReserve, "Reserve2"), + ) + set_service_model!( + template, + ServiceModel( + ReserveDemandTimeSeriesCurve{ReserveUp}, + StepwiseCostReserve, + "ORDC_TS1", + ), + ) + set_service_model!( + template, + ServiceModel( + ReserveDemandTimeSeriesCurve{ReserveUp}, + StepwiseCostReserve, + "ORDC_TS2", + ), + ) + model = DecisionModel( + template, + c_sys5_uc; + name = "UC", + store_variable_names = true, + optimizer = HiGHS_optimizer, + ) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED +end From e4f86ac8faf1a3846384a330733d5dfbe4f82a2e Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Mon, 15 Jun 2026 15:21:10 -0400 Subject: [PATCH 02/11] Clean up time-varying ORDC reserve cost path - Track the IOM helper renames: IOM._get_max_tranches/_unwrap_for_param -> IOM.get_max_tranches/unwrap_for_param at the call sites. - Use the semantic IOM.is_time_variant alias instead of raw IS.is_time_series_backed when splitting the reserve cost into the variant vs invariant objective expression. - Make _maybe_process_stepwise dispatch on PSY.ReserveDemandTimeSeriesCurve (with a no-op fallback) instead of unconditionally inspecting get_variable(service), so a misconfigured service surfaces in the later, more specific cost-data validation rather than throwing during construction. - Tidy ORDC comments/error message. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/common_models/add_parameters.jl | 27 +++++---------------- src/common_models/market_bid_overrides.jl | 9 +------ src/services_models/reserves.jl | 4 +-- src/services_models/services_constructor.jl | 16 ++++++------ 4 files changed, 17 insertions(+), 39 deletions(-) diff --git a/src/common_models/add_parameters.jl b/src/common_models/add_parameters.jl index 15fbc1c..d1086e3 100644 --- a/src/common_models/add_parameters.jl +++ b/src/common_models/add_parameters.jl @@ -45,8 +45,6 @@ function add_parameters!( return end -# Time-varying ORDC: piecewise slope/breakpoint cost parameters for a -# `ReserveDemandTimeSeriesCurve` service. function add_parameters!( container::OptimizationContainer, ::Type{T}, @@ -229,11 +227,8 @@ function _add_time_series_parameters!( jump_model = get_jump_model(container) param_instance = T() parent_param = IOM.get_parameter_array_data(param_container) - # `collect(keys(initial_values))` was passed to `add_param_container!` above as the - # parameter axis, so iterating `initial_values` in order matches the container's - # first axis row-by-row. for (i, (ts_uuid, raw_ts_vals)) in enumerate(initial_values) - ts_vals = IOM._unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) + ts_vals = IOM.unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) @assert all(_size_wrapper.(ts_vals) .== (length.(additional_axes),)) for step in time_steps @@ -242,8 +237,6 @@ function _add_time_series_parameters!( end parent_mult = IOM.get_multiplier_array_data(param_container) - # `devices_with_time_series` was used to build `device_names`, which is the - # multiplier array's first axis, so enumeration index `i` lines up. for (i, device) in enumerate(devices_with_time_series) multiplier = get_multiplier_value(T, device, W) device_name = PSY.get_name(device) @@ -353,7 +346,7 @@ function _add_time_series_parameters!( interval = ts_interval, ) ts_vals = - IOM._unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) + IOM.unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) @assert all(_size_wrapper.(ts_vals) .== (length.(additional_axes),)) end # `_resolve_branch_multiplier` is `DeviceModel`-aware: it applies the @@ -503,7 +496,6 @@ _param_to_vars( ::Type{<:AbstractDeviceFormulation}, ) = (PiecewiseLinearBlockDecrementalOffer,) -# Time-varying ORDC: the reserve demand curve is offered as a decremental block. _param_to_vars( ::Union{ Type{DecrementalPiecewiseLinearSlopeParameter}, @@ -545,7 +537,7 @@ calc_additional_axes( function get_max_tranches(component::PSY.Component, piecewise_ts::IS.TimeSeriesKey) data = PSY.get_data(PSY.get_time_series(component, piecewise_ts)) - return IOM._get_max_tranches(data) + return IOM.get_max_tranches(data) end function calc_additional_axes( @@ -647,9 +639,9 @@ function _add_parameters!( ts_name; interval = ts_interval, ) - ts_vals = IOM._unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) + ts_vals = IOM.unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) @assert all(_size_wrapper.(ts_vals) .== (length.(additional_axes),)) - # PWL/cost-function path: parameter values from `_unwrap_for_param` are + # PWL/cost-function path: parameter values from `unwrap_for_param` are # tuples-of-floats, while the multiplier itself is a scalar Float64. IOM._set_multiplier_at!( parent_mult, @@ -665,10 +657,6 @@ end ################################################################################# # _add_parameters! for time-varying ORDC slope/breakpoint cost parameters -# -# Mirrors the device `ObjectiveFunctionParameter` path above but for a single -# `ReserveDemandTimeSeriesCurve` service whose cost curve is a time series. Allocates the -# 3-D `(name, tranche, time)` container and fills it from the curve time series. ################################################################################# function _add_parameters!( @@ -710,9 +698,6 @@ function _add_parameters!( time_steps; meta = name, ) - # NOTE: ObjectiveFunctionParameter containers carry `CostFunctionAttributes`, which - # (unlike `TimeSeriesAttributes`) have no `set_subsystem!` method — mirror the device - # ObjectiveFunctionParameter path and skip it. jump_model = get_jump_model(container) param_instance = T() @@ -724,7 +709,7 @@ function _add_parameters!( ts_name; interval = ts_interval, ) - ts_vals = IOM._unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) + ts_vals = IOM.unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) @assert all(_size_wrapper.(ts_vals) .== (length.(additional_axes),)) parent_mult = IOM.get_multiplier_array_data(param_container) IOM._set_multiplier_at!(parent_mult, Float64(multiplier), 1) diff --git a/src/common_models/market_bid_overrides.jl b/src/common_models/market_bid_overrides.jl index 0558946..28ee73c 100644 --- a/src/common_models/market_bid_overrides.jl +++ b/src/common_models/market_bid_overrides.jl @@ -293,7 +293,7 @@ _vom_offer_direction(::Type{<:AbstractControllablePowerLoadFormulation}) = ################################################################################# """ -PWL block offer constraints for ORDC (ReserveDemandCurve). ORDC is offered as a +PWL block offer constraints for ORDC. ORDC is offered as a decremental block (the demand curve is a willingness-to-pay), so it uses the decremental block-offer variable/constraint family. """ @@ -343,8 +343,6 @@ function _get_reserve_pwl_data( ) where {T <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}} base_power = get_model_base_power(container) device_base_power = PSY.get_base_power(component, PSY.NU) - # `get_power_units` returns the curve's `IS.AbstractUnitSystem` instance for both the - # static and time-series curves (the parameter arrays store data in those same units). unit_system = PSY.get_power_units(variable_cost) if !IS.is_time_series_backed(variable_cost) @@ -380,11 +378,6 @@ function _get_reserve_pwl_data( return breakpoints, slopes end -""" -PWL cost terms for StepwiseCostReserve (AbstractServiceFormulation). Handles both -static (`ReserveDemandCurve`) and time-varying (`ReserveDemandTimeSeriesCurve`) -reserve demand curves. -""" function add_pwl_term_delta!( container::OptimizationContainer, component::T, diff --git a/src/services_models/reserves.jl b/src/services_models/reserves.jl index 27246c6..e8d2e84 100644 --- a/src/services_models/reserves.jl +++ b/src/services_models/reserves.jl @@ -536,14 +536,14 @@ function _add_reserves_variable_cost_to_objective!( # FIXME clashes with name of a function...ick. variable_cost = PSY.get_variable(component) if variable_cost isa Nothing - error("ReserveDemandCurve $(component.name) does not have cost data.") + error("ORDC curve $(PSY.get_name(component)) does not have cost data.") end pwl_cost_expressions = add_pwl_term_delta!(container, component, variable_cost, T, U) # A time-series-backed curve changes across simulation steps, so its cost goes # into the variant objective expression; a static curve into the invariant one. - is_t_variant = IS.is_time_series_backed(variable_cost) + is_t_variant = is_time_variant(variable_cost) for t in time_steps add_to_expression!( container, diff --git a/src/services_models/services_constructor.jl b/src/services_models/services_constructor.jl index aa20a4e..6e8623d 100644 --- a/src/services_models/services_constructor.jl +++ b/src/services_models/services_constructor.jl @@ -193,14 +193,14 @@ function construct_service!( return end -# Add the time-varying ORDC cost parameters only when the demand curve is backed by -# a time series (a `ReserveDemandTimeSeriesCurve`, whose `variable` is a -# time-series-backed `CostCurve`). Static `ReserveDemandCurve` services are a no-op. -function _maybe_process_stepwise(container, model, service) - IS.is_time_series_backed(PSY.get_variable(service)) && - process_stepwise_cost_reserve_parameters!(container, model, service) - return -end +# Time-varying ORDC parameters are added only for ReserveDemandTimeSeriesCurve; +# static ReserveDemandCurve is a no-op. Dispatch on the service type rather than +# inspecting `get_variable(service)` so a misconfigured service surfaces in the +# later, more specific validation in `_add_reserves_variable_cost_to_objective!` +# instead of throwing here. +_maybe_process_stepwise(container, model, service::PSY.ReserveDemandTimeSeriesCurve) = + process_stepwise_cost_reserve_parameters!(container, model, service) +_maybe_process_stepwise(container, model, service) = nothing function construct_service!( container::OptimizationContainer, From 47b65c3a4c8845445c7975bfa8c4b04ddd933e30 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Mon, 15 Jun 2026 17:08:52 -0400 Subject: [PATCH 03/11] Adapt native branch admittance/3W code to PNM psy6 API MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The native DCP/ACP branch builders and the PowerModels translator called PNM functions that do not exist in the pinned PowerNetworkMatrices psy6 rev: branch_admittance and winding_admittance / three_winding_arcs (the latter two never landed in PNM), and the single-arg get_equivalent_physical_branch_parameters (psy6 added a required NetworkReductionData argument). - AC_branches.jl: add POM-local π-admittance adapters (_branch_admittance, _segment_admittance, _winding_admittance, _three_winding_arcs) that rebuild the (g, b, g_fr, b_fr, g_to, b_to, tap, shift) NamedTuples the Ohm's-law and 3W builders consume, sourced from PSY accessors and psy6's reduction-aware EquivalentBranch / ThreeWindingTransformerWinding API. - pm_translator.jl: thread net_reduction_data into the four series/parallel get_branch_to_pm methods for the new two-arg signature; route the call site by segment type; fix a stale double_circuit reference in a BranchesSeries method. - test_native_dcp_acp_models.jl: migrate the admittance primitive tests to the new POM helpers; rewrite the winding unit test against a real ThreeWindingTransformerWinding. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/ac_transmission_models/AC_branches.jl | 131 +++++++++++++++++++--- src/network_models/pm_translator.jl | 56 ++++++--- test/test_native_dcp_acp_models.jl | 124 +++++++++++++++----- 3 files changed, 253 insertions(+), 58 deletions(-) diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index 7381f0d..3eede3e 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -1411,6 +1411,57 @@ function _reverse_admittance(adm) ) end +# ── Native π-admittance helpers ────────────────────────────────────────────── +# The native DCP/ACP branch builders need each branch's π-model admittance as a +# `(g, b, g_fr, b_fr, g_to, b_to, tap, shift)` NamedTuple (PowerModels convention: +# `g + im*b == 1/(r + im*x)` is the series admittance, `*_fr`/`*_to` are the from/to +# shunts, `tap`/`shift` the transformer ratio/phase). These reproduce, from PSY data and +# PNM's reduction-aware `EquivalentBranch`, exactly what `pm_translator.jl`'s +# `get_branch_to_pm` extracts per branch type — the single source of truth for this +# conversion. (They replace the former `PNM.branch_admittance`, which never landed in PNM.) +_branch_tap(::PSY.ACTransmission) = 1.0 +_branch_tap(b::PSY.TapTransformer) = PSY.get_tap(b) +_branch_tap(b::PSY.PhaseShiftingTransformer) = PSY.get_tap(b) +_branch_shift(::PSY.ACTransmission) = 0.0 +_branch_shift(b::PSY.PhaseShiftingTransformer) = PSY.get_α(b) + +# Lines (and other shunt-on-both-ends branches): π-shunt split from/to. +function _branch_admittance(b::PSY.ACTransmission) + ys = 1.0 / (PSY.get_r(b, PSY.SU) + PSY.get_x(b, PSY.SU) * im) + b_sh = PSY.get_b(b, PSY.SU) + return ( + g = real(ys), b = imag(ys), + g_fr = 0.0, b_fr = b_sh.from, g_to = 0.0, b_to = b_sh.to, + tap = 1.0, shift = 0.0, + ) +end + +# Two-winding transformers: magnetizing shunt allocated to the primary (from) side only. +function _branch_admittance( + b::Union{PSY.Transformer2W, PSY.TapTransformer, PSY.PhaseShiftingTransformer}, +) + ys = 1.0 / (PSY.get_r(b, PSY.SU) + PSY.get_x(b, PSY.SU) * im) + yt = PSY.get_primary_shunt(b, PSY.SU) + return ( + g = real(ys), b = imag(ys), + g_fr = real(yt), b_fr = imag(yt), g_to = 0.0, b_to = 0.0, + tap = _branch_tap(b), shift = _branch_shift(b), + ) +end + +# Reduction-aggregated arc (series chain / parallel group): use PNM's reduction-aware +# equivalent π-parameters. Series/parallel equivalents of lines carry `tap == 1`. +function _segment_admittance(segment, nr::PNM.NetworkReductionData) + eb = PNM.get_equivalent_physical_branch_parameters(segment, nr) + ys = 1.0 / (PNM.get_equivalent_r(eb) + PNM.get_equivalent_x(eb) * im) + return ( + g = real(ys), b = imag(ys), + g_fr = PNM.get_equivalent_g_from(eb), b_fr = PNM.get_equivalent_b_from(eb), + g_to = PNM.get_equivalent_g_to(eb), b_to = PNM.get_equivalent_b_to(eb), + tap = PNM.get_equivalent_tap(eb), shift = PNM.get_equivalent_shift(eb), + ) +end + # Reduction-aware admittance for the retained arc (from_no -> to_no). Returns the PNM # series/parallel equivalent π-tuple oriented from->to when the arc was aggregated by # reduction, or `nothing` when the arc is direct (caller falls back to the branch's own). @@ -1420,13 +1471,13 @@ function _reduced_arc_admittance(nr::PNM.NetworkReductionData, from_no::Int, to_ arc = (from_no, to_no) rev = (to_no, from_no) if haskey(series_map, arc) - return PNM.branch_admittance(series_map[arc], nr) + return _segment_admittance(series_map[arc], nr) elseif haskey(series_map, rev) - return _reverse_admittance(PNM.branch_admittance(series_map[rev], nr)) + return _reverse_admittance(_segment_admittance(series_map[rev], nr)) elseif haskey(parallel_map, arc) - return PNM.branch_admittance(parallel_map[arc], nr) + return _segment_admittance(parallel_map[arc], nr) elseif haskey(parallel_map, rev) - return _reverse_admittance(PNM.branch_admittance(parallel_map[rev], nr)) + return _reverse_admittance(_segment_admittance(parallel_map[rev], nr)) end return nothing end @@ -1436,9 +1487,9 @@ end # series/parallel-aggregated arc. function _resolve_branch_admittance(network_model, branch, from_no::Int, to_no::Int) nr = get_network_reduction(network_model) - isempty(nr) && return PNM.branch_admittance(branch) + isempty(nr) && return _branch_admittance(branch) eq = _reduced_arc_admittance(nr, from_no, to_no) - return eq === nothing ? PNM.branch_admittance(branch) : eq + return eq === nothing ? _branch_admittance(branch) : eq end # One-pass per-branch network geometry for the native DCP/ACP builders. Computes each @@ -1455,7 +1506,7 @@ function _branch_geometries(sys::PSY.System, network_model, devices) collapsed = fr.number == to.number adm = if collapsed - PNM.branch_admittance(d) + _branch_admittance(d) else _resolve_branch_admittance(network_model, d, fr.number, to.number) end @@ -1972,12 +2023,58 @@ end # storage 2D (name × time) without inventing a new container shape. ################################################################################ +# ── Native three-winding helpers ───────────────────────────────────────────── +# Decompose a Transformer3W into its three wye-model windings via PNM's +# `ThreeWindingTransformerWinding`, exposing the per-winding data the native 3W builders +# need: the star-point arc (for reduction-aware bus mapping), the winding rating, and the +# winding object (for admittance). `suffix` reproduces PNM's per-winding name component +# (`"_winding_"`) so variable/constraint names stay stable. +# (Replaces the former `PNM.three_winding_arcs` / `PNM.winding_admittance`, which never +# landed in PNM; psy6 exposes the winding decomposition through this type instead.) +function _three_winding_arcs(d::PSY.Transformer3W) + star_arcs = ( + PSY.get_primary_star_arc(d), + PSY.get_secondary_star_arc(d), + PSY.get_tertiary_star_arc(d), + ) + out = NamedTuple[] + for i in 1:3 + w = PNM.ThreeWindingTransformerWinding(d, i) + push!( + out, + ( + suffix = "winding_$i", + arc = star_arcs[i], + rating = PNM.get_equivalent_rating(w), + winding = w, + ), + ) + end + return out +end + +# No phase shift for windings; only the (optional) winding tap matters for the ACP π-model. +_winding_tap(::PNM.ThreeWindingTransformerWinding) = 1.0 +_winding_tap(w::PNM.ThreeWindingTransformerWinding{PSY.PhaseShiftingTransformer3W}) = + PNM.get_equivalent_tap(w) + +# Per-winding π-admittance NamedTuple (shunt on the primary side only, no phase shift). +function _winding_admittance(w::PNM.ThreeWindingTransformerWinding) + ys = 1.0 / (PNM.get_equivalent_r(w) + PNM.get_equivalent_x(w) * im) + b_sh = PNM.get_equivalent_b(w) + return ( + g = real(ys), b = imag(ys), + g_fr = 0.0, b_fr = b_sh.from, g_to = 0.0, b_to = b_sh.to, + tap = _winding_tap(w), + ) +end + "Build the list of per-winding variable names for a set of Transformer3W devices." function _three_winding_var_names(devices) names = String[] for d in devices dname = PSY.get_name(d) - for w in PNM.three_winding_arcs(d) + for w in _three_winding_arcs(d) push!(names, dname * "_" * w.suffix) end end @@ -2046,7 +2143,7 @@ function add_to_expression!( time_steps = get_time_steps(container) for d in devices dname = PSY.get_name(d) - for w in PNM.three_winding_arcs(d) + for w in _three_winding_arcs(d) wname = dname * "_" * w.suffix from_no = PNM.get_mapped_bus_number(network_reduction, PSY.get_from(w.arc)) to_no = PNM.get_mapped_bus_number(network_reduction, PSY.get_to(w.arc)) @@ -2082,7 +2179,7 @@ for (E, V, isfrom) in ( time_steps = get_time_steps(container) for d in devices dname = PSY.get_name(d) - for w in PNM.three_winding_arcs(d) + for w in _three_winding_arcs(d) wname = dname * "_" * w.suffix terminal_bus_obj = $isfrom ? PSY.get_from(w.arc) : PSY.get_to(w.arc) bus_no = PNM.get_mapped_bus_number(network_reduction, terminal_bus_obj) @@ -2116,9 +2213,9 @@ function add_constraints!( for d in devices dname = PSY.get_name(d) - for w in PNM.three_winding_arcs(d) + for w in _three_winding_arcs(d) wname = dname * "_" * w.suffix - adm = PNM.winding_admittance(w) + adm = _winding_admittance(w.winding) fr = _retained_bus(number_to_name, network_model, PSY.get_from(w.arc)) to = _retained_bus(number_to_name, network_model, PSY.get_to(w.arc)) fr.number == to.number && continue @@ -2173,9 +2270,9 @@ function add_constraints!( for d in devices dname = PSY.get_name(d) - for w in PNM.three_winding_arcs(d) + for w in _three_winding_arcs(d) wname = dname * "_" * w.suffix - adm = PNM.winding_admittance(w) + adm = _winding_admittance(w.winding) g, b, g_fr, b_fr, g_to, b_to, tm = adm.g, adm.b, adm.g_fr, adm.b_fr, adm.g_to, adm.b_to, adm.tap fr = _retained_bus(number_to_name, network_model, PSY.get_from(w.arc)) @@ -2246,7 +2343,7 @@ function add_constraints!( ) for d in devices dname = PSY.get_name(d) - for w in PNM.three_winding_arcs(d) + for w in _three_winding_arcs(d) wname = dname * "_" * w.suffix for t in time_steps cons_lb[wname, t] = JuMP.@constraint( @@ -2278,7 +2375,7 @@ function add_constraints!( ) for d in devices dname = PSY.get_name(d) - for w in PNM.three_winding_arcs(d) + for w in _three_winding_arcs(d) wname = dname * "_" * w.suffix r2 = w.rating^2 for t in time_steps @@ -2308,7 +2405,7 @@ function add_constraints!( ) for d in devices dname = PSY.get_name(d) - for w in PNM.three_winding_arcs(d) + for w in _three_winding_arcs(d) wname = dname * "_" * w.suffix r2 = w.rating^2 for t in time_steps diff --git a/src/network_models/pm_translator.jl b/src/network_models/pm_translator.jl index 6e25218..69a1aa6 100644 --- a/src/network_models/pm_translator.jl +++ b/src/network_models/pm_translator.jl @@ -439,8 +439,13 @@ function get_branch_to_pm( T::Type{<:AbstractBranchFormulation}, U::Type{<:AbstractPowerModel}, device_model::DeviceModel, + network_reduction_data::PNM.NetworkReductionData, ) - equivalent_branch = PNM.get_equivalent_physical_branch_parameters(double_circuit) + equivalent_branch = + PNM.get_equivalent_physical_branch_parameters( + double_circuit, + network_reduction_data, + ) rating = branch_rating(double_circuit, device_model) PM_branch = Dict{String, Any}( "br_r" => PNM.get_equivalent_r(equivalent_branch), @@ -473,8 +478,13 @@ function get_branch_to_pm( T::Type{StaticBranchUnbounded}, U::Type{<:AbstractPowerModel}, device_model::DeviceModel, + network_reduction_data::PNM.NetworkReductionData, ) - equivalent_branch = PNM.get_equivalent_physical_branch_parameters(double_circuit) + equivalent_branch = + PNM.get_equivalent_physical_branch_parameters( + double_circuit, + network_reduction_data, + ) rating = branch_rating(double_circuit, device_model) PM_branch = Dict{String, Any}( "br_r" => PNM.get_equivalent_r(equivalent_branch), @@ -507,8 +517,10 @@ function get_branch_to_pm( T::Type{<:AbstractBranchFormulation}, U::Type{<:AbstractPowerModel}, ::DeviceModel, + network_reduction_data::PNM.NetworkReductionData, ) - equivalent_branch = PNM.get_equivalent_physical_branch_parameters(series_chain) + equivalent_branch = + PNM.get_equivalent_physical_branch_parameters(series_chain, network_reduction_data) PM_branch = Dict{String, Any}( "br_r" => PNM.get_equivalent_r(equivalent_branch), "shift" => 0.0, @@ -536,8 +548,10 @@ function get_branch_to_pm( T::Type{StaticBranchUnbounded}, U::Type{<:AbstractPowerModel}, ::DeviceModel, + network_reduction_data::PNM.NetworkReductionData, ) - equivalent_branch = PNM.get_equivalent_physical_branch_parameters(series_chain) + equivalent_branch = + PNM.get_equivalent_physical_branch_parameters(series_chain, network_reduction_data) PM_branch = Dict{String, Any}( "br_r" => PNM.get_equivalent_r(equivalent_branch), "shift" => 0.0, @@ -546,7 +560,7 @@ function get_branch_to_pm( "g_fr" => 0.0, "b_fr" => PNM.get_equivalent_b_from(equivalent_branch), "f_bus" => arc_tuple[1], - "br_status" => Float64(PNM.get_equivalent_available(double_circuit)), + "br_status" => Float64(PNM.get_equivalent_available(series_chain)), "t_bus" => arc_tuple[2], "b_to" => PNM.get_equivalent_b_to(equivalent_branch), "index" => ix, @@ -733,14 +747,30 @@ function get_branches_to_pm( for (_, (arc_tuple, reduction)) in name_to_arc_map arc_tuple ∈ modeled_arc_tuples && continue # This is the PowerModels equivalent of the branch and constraint tracker. reduction_entry = all_branch_maps_by_type[reduction][comp_type][arc_tuple] - PM_branches["$(ix)"] = get_branch_to_pm( - ix, - arc_tuple, - reduction_entry, - get_formulation(device_model), - S, - device_model, - ) + # Reduction-aggregated arcs (series/parallel segments) need the reduction data + # to resolve their equivalent π-parameters; direct branches do not. + PM_branches["$(ix)"] = + if reduction_entry isa + Union{PNM.AbstractBranchesParallel, PNM.BranchesSeries} + get_branch_to_pm( + ix, + arc_tuple, + reduction_entry, + get_formulation(device_model), + S, + device_model, + net_reduction_data, + ) + else + get_branch_to_pm( + ix, + arc_tuple, + reduction_entry, + get_formulation(device_model), + S, + device_model, + ) + end if PM_branches["$(ix)"]["br_status"] == true f = PM_branches["$(ix)"]["f_bus"] t = PM_branches["$(ix)"]["t_bus"] diff --git a/test/test_native_dcp_acp_models.jl b/test/test_native_dcp_acp_models.jl index c3e471c..3025dab 100644 --- a/test/test_native_dcp_acp_models.jl +++ b/test/test_native_dcp_acp_models.jl @@ -152,7 +152,7 @@ end @testset "branch_admittance primitives" begin sys = PSB.build_system(PSITestSystems, "c_sys5") line = first(PSY.get_components(PSY.Line, sys)) - a = PNM.branch_admittance(line) + a = PowerOperationsModels._branch_admittance(line) r, x = PSY.get_r(line, PSY.SU), PSY.get_x(line, PSY.SU) y = inv(complex(r, x)) @test a.g ≈ real(y) @@ -196,7 +196,7 @@ import PowerNetworkMatrices as PNM (from_no, to_no), chain = first(series_map) resolved = PowerOperationsModels._reduced_arc_admittance(nr, from_no, to_no) @test resolved !== nothing - expected = PNM.branch_admittance(chain, nr) + expected = PowerOperationsModels._segment_admittance(chain, nr) @test isapprox(resolved.b, expected.b; atol = 1e-9) # Non-triviality: the series equivalent is the MERGED admittance of the chain, so @@ -205,7 +205,7 @@ import PowerNetworkMatrices as PNM # for the reduced arc, which is wrong. members = collect(chain) @test length(members) >= 2 - member_b = PNM.branch_admittance(members[1]).b + member_b = PowerOperationsModels._branch_admittance(members[1]).b @test !isapprox(resolved.b, member_b; rtol = 1e-3) # Reversed-orientation arc exercises the `_reverse_admittance` path: series b is @@ -224,32 +224,100 @@ import PowerNetworkMatrices as PNM end @testset "Transformer3W _winding_admittance star-arc decomposition" begin - # Unit test the per-winding admittance helper: for a winding with series - # impedance R + jX the helper must return the admittance 1/(R + jX), i.e. - # g = R / (R^2 + X^2), b = -X / (R^2 + X^2) - # with zero shunts, zero phase shift, and the winding's tap passed through. - R = 0.01 - X = 0.1 - w = ( - suffix = "winding_1", - arc = nothing, - r = R, - x = X, - rating = 1.0, - tap = 1.0, + # Unit test the per-winding admittance helper against a real PNM + # `ThreeWindingTransformerWinding`: for a winding with series impedance R + jX the + # helper must return the series admittance 1/(R + jX), the winding's PNM shunt on the + # from/to sides, no phase shift, and (here) a unit tap. R/X are read back through PNM + # so the assertion is robust to per-unit base conversions. + sys = PSB.build_system(PSITestSystems, "c_sys5_ml") + busD = PSY.get_component(PSY.ACBus, sys, "nodeD") + star_bus = PSY.ACBus(; + number = 103, + name = "Star_Bus_T3W", + available = true, + bustype = PSY.ACBusTypes.PQ, + angle = 0.0, + magnitude = 1.0, + voltage_limits = (min = 0.95, max = 1.05), + base_voltage = 230.0, + area = PSY.get_area(busD), + load_zone = PSY.get_load_zone(busD), ) - adm = PNM.winding_admittance(w) - denom = R^2 + X^2 - @test isapprox(adm.g, R / denom; atol = 1e-12) - @test isapprox(adm.b, -X / denom; atol = 1e-12) - @test adm.g_fr == 0.0 - @test adm.b_fr == 0.0 - @test adm.g_to == 0.0 - @test adm.b_to == 0.0 - @test adm.tap == w.tap - @test adm.shift == 0.0 - # Cross-check against direct complex inversion. - y = inv(complex(R, X)) + PSY.add_component!(sys, star_bus) + sec_bus = PSY.ACBus(; + number = 101, + name = "Bus3WT_1", + available = true, + bustype = PSY.ACBusTypes.PQ, + angle = 0.0, + magnitude = 1.0, + voltage_limits = (min = 0.95, max = 1.05), + base_voltage = 230.0, + area = PSY.get_area(busD), + load_zone = PSY.get_load_zone(busD), + ) + PSY.add_component!(sys, sec_bus) + ter_bus = PSY.ACBus(; + number = 102, + name = "Bus3WT_2", + available = true, + bustype = PSY.ACBusTypes.PQ, + angle = 0.0, + magnitude = 1.0, + voltage_limits = (min = 0.95, max = 1.05), + base_voltage = 230.0, + area = PSY.get_area(busD), + load_zone = PSY.get_load_zone(busD), + ) + PSY.add_component!(sys, ter_bus) + transformer3w = PSY.Transformer3W(; + name = "Transformer3W_busD", + available = true, + primary_star_arc = PSY.Arc(; from = busD, to = star_bus), + secondary_star_arc = PSY.Arc(; from = sec_bus, to = star_bus), + tertiary_star_arc = PSY.Arc(; from = ter_bus, to = star_bus), + star_bus = star_bus, + active_power_flow_primary = 0.0, + reactive_power_flow_primary = 0.0, + active_power_flow_secondary = 0.0, + reactive_power_flow_secondary = 0.0, + active_power_flow_tertiary = 0.0, + reactive_power_flow_tertiary = 0.0, + r_primary = 0.01, + x_primary = 0.1, + r_secondary = 0.01, + x_secondary = 0.1, + r_tertiary = 0.01, + x_tertiary = 0.1, + r_12 = 0.01, + x_12 = 0.1, + r_23 = 0.01, + x_23 = 0.1, + r_13 = 0.01, + x_13 = 0.1, + base_power_12 = 100.0, + base_power_23 = 100.0, + base_power_13 = 100.0, + rating = nothing, + rating_primary = 1.0, + rating_secondary = 1.0, + rating_tertiary = 0.5, + ) + PSY.add_component!(sys, transformer3w) + + w = PNM.ThreeWindingTransformerWinding(transformer3w, 1) + adm = PowerOperationsModels._winding_admittance(w) + + r = PNM.get_equivalent_r(w) + x = PNM.get_equivalent_x(w) + y = inv(complex(r, x)) @test isapprox(adm.g, real(y); atol = 1e-12) @test isapprox(adm.b, imag(y); atol = 1e-12) + + b_sh = PNM.get_equivalent_b(w) + @test adm.g_fr == 0.0 + @test adm.b_fr == b_sh.from + @test adm.g_to == 0.0 + @test adm.b_to == b_sh.to + @test adm.tap == 1.0 end From f672cd43f7c0a083e5ca482eb0dfff668cedaf59 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Mon, 15 Jun 2026 19:48:44 -0400 Subject: [PATCH 04/11] Replace isa dispatch in pm_translator; de-flake SC MODF/PTDF tests - pm_translator: collapse the runtime `isa` branch that selected the reduction-aware `get_branch_to_pm` overload into multiple dispatch, via a 7-arg forwarding method that drops the network-reduction data for direct branches (the parallel/series methods remain strictly more specific). - test_ac_transmission_security_constrained_models: compare post-contingency MODF/PTDF expressions up to factorization round-off (`_affexpr_approx_equal`) instead of exact `isequal_canonical`. Apple Accelerate's sparse factorization (PNM's macOS default) is not bit-reproducible across builds, so two independently built matrices differ by ~1e-15; the test's intentionally-fresh ground-truth matrix surfaced this as intermittent failures. - Minor: comment cleanups (AC_branches, reserves, services_constructor) and "repin once IOM merged" TODOs in Project.toml/test/Project.toml. Co-Authored-By: Claude Opus 4.8 (1M context) --- Project.toml | 1 + src/ac_transmission_models/AC_branches.jl | 10 +--- src/network_models/pm_translator.jl | 48 +++++++++---------- src/services_models/reserves.jl | 2 - src/services_models/services_constructor.jl | 5 -- test/Project.toml | 1 + ...ransmission_security_constrained_models.jl | 21 ++++++-- 7 files changed, 46 insertions(+), 42 deletions(-) diff --git a/Project.toml b/Project.toml index bf40441..04b526f 100644 --- a/Project.toml +++ b/Project.toml @@ -27,6 +27,7 @@ PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6" [sources] InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl"} PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} +# todo: repin once iom merged InfrastructureOptimizationModels = {rev = "ac/ordc", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} PowerNetworkMatrices = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"} diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index 3eede3e..d53a24a 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -1415,10 +1415,7 @@ end # The native DCP/ACP branch builders need each branch's π-model admittance as a # `(g, b, g_fr, b_fr, g_to, b_to, tap, shift)` NamedTuple (PowerModels convention: # `g + im*b == 1/(r + im*x)` is the series admittance, `*_fr`/`*_to` are the from/to -# shunts, `tap`/`shift` the transformer ratio/phase). These reproduce, from PSY data and -# PNM's reduction-aware `EquivalentBranch`, exactly what `pm_translator.jl`'s -# `get_branch_to_pm` extracts per branch type — the single source of truth for this -# conversion. (They replace the former `PNM.branch_admittance`, which never landed in PNM.) +# shunts, `tap`/`shift` the transformer ratio/phase). _branch_tap(::PSY.ACTransmission) = 1.0 _branch_tap(b::PSY.TapTransformer) = PSY.get_tap(b) _branch_tap(b::PSY.PhaseShiftingTransformer) = PSY.get_tap(b) @@ -2027,10 +2024,7 @@ end # Decompose a Transformer3W into its three wye-model windings via PNM's # `ThreeWindingTransformerWinding`, exposing the per-winding data the native 3W builders # need: the star-point arc (for reduction-aware bus mapping), the winding rating, and the -# winding object (for admittance). `suffix` reproduces PNM's per-winding name component -# (`"_winding_"`) so variable/constraint names stay stable. -# (Replaces the former `PNM.three_winding_arcs` / `PNM.winding_admittance`, which never -# landed in PNM; psy6 exposes the winding decomposition through this type instead.) +# winding object (for admittance). function _three_winding_arcs(d::PSY.Transformer3W) star_arcs = ( PSY.get_primary_star_arc(d), diff --git a/src/network_models/pm_translator.jl b/src/network_models/pm_translator.jl index 69a1aa6..4cdc7dc 100644 --- a/src/network_models/pm_translator.jl +++ b/src/network_models/pm_translator.jl @@ -572,6 +572,18 @@ function get_branch_to_pm( return PM_branch end +function get_branch_to_pm( + ix::Int, + arc_tuple::Tuple{Int, Int}, + branch, + formulation::Type{<:AbstractBranchFormulation}, + network::Type{<:AbstractPowerModel}, + device_model::DeviceModel, + ::PNM.NetworkReductionData, +) + return get_branch_to_pm(ix, arc_tuple, branch, formulation, network, device_model) +end + function get_branch_to_pm( ix::Int, branch::PSY.TwoTerminalGenericHVDCLine, @@ -747,30 +759,18 @@ function get_branches_to_pm( for (_, (arc_tuple, reduction)) in name_to_arc_map arc_tuple ∈ modeled_arc_tuples && continue # This is the PowerModels equivalent of the branch and constraint tracker. reduction_entry = all_branch_maps_by_type[reduction][comp_type][arc_tuple] - # Reduction-aggregated arcs (series/parallel segments) need the reduction data - # to resolve their equivalent π-parameters; direct branches do not. - PM_branches["$(ix)"] = - if reduction_entry isa - Union{PNM.AbstractBranchesParallel, PNM.BranchesSeries} - get_branch_to_pm( - ix, - arc_tuple, - reduction_entry, - get_formulation(device_model), - S, - device_model, - net_reduction_data, - ) - else - get_branch_to_pm( - ix, - arc_tuple, - reduction_entry, - get_formulation(device_model), - S, - device_model, - ) - end + # Reduction-aggregated arcs (series/parallel segments) need the reduction + # data to resolve their equivalent π-parameters; direct branches dispatch + # to a method that ignores it. + PM_branches["$(ix)"] = get_branch_to_pm( + ix, + arc_tuple, + reduction_entry, + get_formulation(device_model), + S, + device_model, + net_reduction_data, + ) if PM_branches["$(ix)"]["br_status"] == true f = PM_branches["$(ix)"]["f_bus"] t = PM_branches["$(ix)"]["t_bus"] diff --git a/src/services_models/reserves.jl b/src/services_models/reserves.jl index e8d2e84..31becd7 100644 --- a/src/services_models/reserves.jl +++ b/src/services_models/reserves.jl @@ -541,8 +541,6 @@ function _add_reserves_variable_cost_to_objective!( pwl_cost_expressions = add_pwl_term_delta!(container, component, variable_cost, T, U) - # A time-series-backed curve changes across simulation steps, so its cost goes - # into the variant objective expression; a static curve into the invariant one. is_t_variant = is_time_variant(variable_cost) for t in time_steps add_to_expression!( diff --git a/src/services_models/services_constructor.jl b/src/services_models/services_constructor.jl index 6e8623d..b6a81bf 100644 --- a/src/services_models/services_constructor.jl +++ b/src/services_models/services_constructor.jl @@ -193,11 +193,6 @@ function construct_service!( return end -# Time-varying ORDC parameters are added only for ReserveDemandTimeSeriesCurve; -# static ReserveDemandCurve is a no-op. Dispatch on the service type rather than -# inspecting `get_variable(service)` so a misconfigured service surfaces in the -# later, more specific validation in `_add_reserves_variable_cost_to_objective!` -# instead of throwing here. _maybe_process_stepwise(container, model, service::PSY.ReserveDemandTimeSeriesCurve) = process_stepwise_cost_reserve_parameters!(container, model, service) _maybe_process_stepwise(container, model, service) = nothing diff --git a/test/Project.toml b/test/Project.toml index 85136d9..4c4c216 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -34,6 +34,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [sources] +# TODO repin once IOM merged InfrastructureOptimizationModels = {rev = "ac/ordc", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl"} PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} diff --git a/test/test_ac_transmission_security_constrained_models.jl b/test/test_ac_transmission_security_constrained_models.jl index 8395e41..cc2876b 100644 --- a/test/test_ac_transmission_security_constrained_models.jl +++ b/test/test_ac_transmission_security_constrained_models.jl @@ -12,6 +12,21 @@ # auto-populated by `_build_device_model_outages!` during template validation — # no manual stand-ins. +# The ground-truth columns below are read from a freshly-built VirtualMODF/ +# VirtualPTDF (intentionally independent of the production matrix). On macOS the +# default sparse backend is Apple Accelerate, whose factorization is internally +# multithreaded and NOT bit-reproducible across builds, so two independent +# factorizations of the identical ABA matrix differ by ~1e-15. That noise is +# benign (a single production matrix is built once and reused, so it is +# self-consistent), but it breaks exact `isequal_canonical` on the re-derived +# coefficients. Compare the affine expressions up to that round-off instead: the +# tolerance is ~1e7× the observed noise yet far below any real coefficient bug. +function _affexpr_approx_equal(actual, expected; atol = 1e-8) + d = actual - expected # matching terms cancel to ~1e-15 + coeff_resid = maximum((abs(c) for (c, _) in JuMP.linear_terms(d)); init = 0.0) + return coeff_resid ≤ atol && abs(JuMP.constant(d)) ≤ atol +end + @testset "Post-contingency expressions and constraints match MODF ground truth" begin c_sys5 = PSB.build_system(PSITestSystems, "c_sys5") all_branches = collect(get_components(PSY.ACTransmission, c_sys5)) @@ -103,7 +118,7 @@ JuMP.add_to_expression!(expected, modf_col[i], nodal_balance[i, t]) end actual = pcbf[outage_id_str, name, t] - @test JuMP.isequal_canonical(actual, expected) + @test _affexpr_approx_equal(actual, expected) # --- Constraint RHS equality: emergency rating in PER-UNIT. --- # The post-contingency flow expression carries an affine constant @@ -221,7 +236,7 @@ end JuMP.add_to_expression!(expected, ptdf_col[i], nodal_balance[i, t]) end actual = pbf[name, t] - @test JuMP.isequal_canonical(actual, expected) + @test _affexpr_approx_equal(actual, expected) end end end @@ -320,7 +335,7 @@ end JuMP.add_to_expression!(expected, modf_col[i], nodal_balance[i, t]) end actual = pcbf[outage_id_str, name, t] - @test JuMP.isequal_canonical(actual, expected) + @test _affexpr_approx_equal(actual, expected) end end @test n_checked >= 1 From 732ba8fa5dd90e654bdcd2d76f89bc865dc0b7cc Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Mon, 15 Jun 2026 20:02:24 -0400 Subject: [PATCH 05/11] better error catching --- src/services_models/reserves.jl | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/src/services_models/reserves.jl b/src/services_models/reserves.jl index 31becd7..9993508 100644 --- a/src/services_models/reserves.jl +++ b/src/services_models/reserves.jl @@ -536,12 +536,30 @@ function _add_reserves_variable_cost_to_objective!( # FIXME clashes with name of a function...ick. variable_cost = PSY.get_variable(component) if variable_cost isa Nothing - error("ORDC curve $(PSY.get_name(component)) does not have cost data.") + error("ORDC curve $(component_name) does not have cost data.") + elseif !(variable_cost isa PSY.CostCurve) + error( + "ORDC curve $(component_name) has cost data of type $(typeof(variable_cost)), \ + but a `PSY.CostCurve` is required for the StepwiseCostReserve formulation.", + ) + end + + # A time-series-backed cost varies across simulation steps and is read from + # per-timestep parameter arrays, which are only populated for + # `ReserveDemandTimeSeriesCurve` (see `process_stepwise_cost_reserve_parameters!`). + # Reject a time-series-backed cost on any other reserve type up front, rather + # than failing later with a missing-parameter error. + is_t_variant = is_time_variant(variable_cost) + if is_t_variant && !(component isa PSY.ReserveDemandTimeSeriesCurve) + error( + "ORDC curve $(component_name) of type $(typeof(component)) has a \ + time-series-backed cost; a `PSY.ReserveDemandTimeSeriesCurve` is required \ + for time-varying ORDC cost.", + ) end pwl_cost_expressions = add_pwl_term_delta!(container, component, variable_cost, T, U) - is_t_variant = is_time_variant(variable_cost) for t in time_steps add_to_expression!( container, From d82738a5679fc1f3109825d5074666cbb347cdd0 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 16 Jun 2026 12:46:37 -0400 Subject: [PATCH 06/11] Unify ORDC reserve PWL cost path with IOM's OfferDirection abstraction The StepwiseCostReserve cost path (over ReserveDemandCurve / ReserveDemandTimeSeriesCurve) duplicated logic IOM already provides generically. Route it through the OfferDirection machinery instead of hardcoding Decremental*: - reserves.jl: add the `_reserve_offer_direction` trait (decremental); derive the stepwise param loop from it. - market_bid_plumbing.jl: add `get_offer_curves` for reserve curves; thread `meta` through `_get_pwl_data` and the static `_get_raw_pwl_data`. - market_bid_overrides.jl: delete `_get_reserve_pwl_data` and the reserve `add_pwl_constraint_delta!`; slim `add_pwl_term_delta!` to call the IOM helpers with `dir` + `meta = name` (both static and time-series curves now flow through `IOM._get_pwl_data`). - add_parameters.jl: share the device ObjectiveFunctionParameter body via `_add_objective_function_parameters!(...; meta)`; device + ORDC `_add_parameters!` become thin wrappers; `_calc_additional_axes` dispatches the collection-vs-single `calc_additional_axes` call by model type. Depends on the meta-aware IOM delta-PWL helpers (IOM ac/ordc d43dd12). Co-Authored-By: Claude Opus 4.8 (1M context) --- src/common_models/add_parameters.jl | 112 ++++++++++---------- src/common_models/market_bid_overrides.jl | 122 ++++++---------------- src/common_models/market_bid_plumbing.jl | 20 +++- src/services_models/reserves.jl | 11 +- 4 files changed, 114 insertions(+), 151 deletions(-) diff --git a/src/common_models/add_parameters.jl b/src/common_models/add_parameters.jl index d1086e3..67131e2 100644 --- a/src/common_models/add_parameters.jl +++ b/src/common_models/add_parameters.jl @@ -570,20 +570,42 @@ function calc_additional_axes( return (IOM.make_tranche_axis(max_tranches + 1),) # one more breakpoint than tranches end +# The shared `_add_objective_function_parameters!` helper holds the active components +# as a vector, but `calc_additional_axes` is called differently per model: devices have +# collection forms, while a service is passed singly (its concrete two-parameter type +# can't unify with the ServiceModel's partially-applied component type through an +# invariant `Vector`, and its forms take one service). Dispatch on the model type. +_calc_additional_axes( + container::OptimizationContainer, + ::Type{P}, + active, + model::DeviceModel, +) where {P <: ParameterType} = calc_additional_axes(container, P, active, model) +_calc_additional_axes( + container::OptimizationContainer, + ::Type{P}, + active, + model::ServiceModel, +) where {P <: ParameterType} = calc_additional_axes(container, P, only(active), model) + ################################################################################# # _add_parameters! for ObjectiveFunctionParameter ################################################################################# -function _add_parameters!( +# Shared body for ObjectiveFunctionParameter cost params (slope/breakpoint, +# cost-at-min, startup/shutdown). Devices are batched into one names-axis container +# (`meta` empty); ORDC services pass a 1-element collection with `meta = name` +# because services are constructed one `ServiceModel` at a time (one container per +# service, not a shared names-axis batch). `W` is the device/service formulation. +function _add_objective_function_parameters!( container::OptimizationContainer, - param::Type{T}, - devices::U, - model::DeviceModel{D, W}, -) where { - T <: ObjectiveFunctionParameter, - U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, - W <: AbstractDeviceFormulation, -} where {D <: PSY.Component} + ::Type{T}, + devices::V, + model, + ::Type{W}; + meta = IOM.CONTAINER_KEY_EMPTY_META, +) where {T <: ObjectiveFunctionParameter, V, W} + D = eltype(devices) ts_type = get_default_time_series_type(container) if !(ts_type <: Union{PSY.AbstractDeterministic, PSY.StaticTimeSeries}) error( @@ -610,7 +632,7 @@ function _add_parameters!( end jump_model = get_jump_model(container) - additional_axes = calc_additional_axes(container, param, active_devices, model) + additional_axes = _calc_additional_axes(container, T, active_devices, model) param_container = add_param_container!( container, T, @@ -621,7 +643,8 @@ function _add_parameters!( _get_expected_time_series_eltype(T), device_names, additional_axes..., - time_steps, + time_steps; + meta = meta, ) param_instance = T() @@ -655,8 +678,29 @@ function _add_parameters!( return end +# Device entry: batch all components of a type into one names-axis container +# (`meta` defaults to empty, reproducing the prior key exactly). +function _add_parameters!( + container::OptimizationContainer, + ::Type{T}, + devices::U, + model::DeviceModel{D, W}, +) where { + T <: ObjectiveFunctionParameter, + U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, + W <: AbstractDeviceFormulation, +} where {D <: PSY.Component} + _add_objective_function_parameters!(container, T, devices, model, W) + return +end + ################################################################################# # _add_parameters! for time-varying ORDC slope/breakpoint cost parameters +# +# Same cost-parameter machinery as the device path, but per-service: pass the +# single service as a 1-element collection and key the container by `meta = name` +# (services are constructed one ServiceModel at a time, so they can't share a +# names-axis batch the way devices do). ################################################################################# function _add_parameters!( @@ -672,51 +716,9 @@ function _add_parameters!( U <: PSY.ReserveDemandTimeSeriesCurve, V <: AbstractServiceFormulation, } - ts_type = get_default_time_series_type(container) - if !(ts_type <: Union{PSY.AbstractDeterministic, PSY.StaticTimeSeries}) - error( - "add_parameters! for $T is not compatible with $ts_type", - ) - end - ts_name = _get_time_series_name(T, service, model) - time_steps = get_time_steps(container) - name = PSY.get_name(service) - model_interval = get_interval(get_settings(container)) - ts_interval = model_interval - - additional_axes = calc_additional_axes(container, T, service, model) - param_container = add_param_container!( - container, - T, - U, - _param_to_vars(T, V), - SOSStatusVariable.NO_VARIABLE, - false, - _get_expected_time_series_eltype(T), - [name], - additional_axes..., - time_steps; - meta = name, - ) - - jump_model = get_jump_model(container) - param_instance = T() - multiplier = get_multiplier_value(T, service, V) - raw_ts_vals = IOM.get_time_series_initial_values!( - container, - ts_type, - service, - ts_name; - interval = ts_interval, + _add_objective_function_parameters!( + container, T, [service], model, V; meta = PSY.get_name(service), ) - ts_vals = IOM.unwrap_for_param.((param_instance,), raw_ts_vals, (additional_axes,)) - @assert all(_size_wrapper.(ts_vals) .== (length.(additional_axes),)) - parent_mult = IOM.get_multiplier_array_data(param_container) - IOM._set_multiplier_at!(parent_mult, Float64(multiplier), 1) - parent_param = IOM.get_parameter_array_data(param_container) - for t in time_steps - IOM._set_parameter_at!(parent_param, jump_model, ts_vals[t], 1, t) - end return end diff --git a/src/common_models/market_bid_overrides.jl b/src/common_models/market_bid_overrides.jl index 28ee73c..acafe7c 100644 --- a/src/common_models/market_bid_overrides.jl +++ b/src/common_models/market_bid_overrides.jl @@ -293,95 +293,27 @@ _vom_offer_direction(::Type{<:AbstractControllablePowerLoadFormulation}) = ################################################################################# """ -PWL block offer constraints for ORDC. ORDC is offered as a -decremental block (the demand curve is a willingness-to-pay), so it uses the -decremental block-offer variable/constraint family. +Add the delta PWL objective terms for an ORDC service (`StepwiseCostReserve` over +`ReserveDemandCurve` / `ReserveDemandTimeSeriesCurve`). The demand curve is a +willingness-to-pay (concave) decremental offer, so the whole path routes through +IOM's generic `OfferDirection` machinery via `_reserve_offer_direction`: + +* `IOM._get_pwl_data` fetches the curve (`get_offer_curves` → `PSY.get_variable`) + and dispatches the static vs. time-series read internally. +* `IOM._block_offer_var` / `IOM._block_offer_constraint` pick the block-offer + variable/constraint family. +* `IOM.add_pwl_constraint_delta!` builds the block-sum constraint. + +All containers are keyed by `meta = name` because services are built one at a +time (one container per service, not a shared names-axis batch). + +`cost_data` is accepted to match the caller's signature but is re-resolved inside +`IOM._get_pwl_data`; it is the same `PSY.get_variable(component)` curve. """ -function add_pwl_constraint_delta!( - container::OptimizationContainer, - component::T, - ::Type{U}, - break_points::Vector{Float64}, - pwl_vars::Vector{JuMP.VariableRef}, - period::Int, -) where { - T <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, - U <: ServiceRequirementVariable, -} - name = PSY.get_name(component) - variables = get_variable(container, U, T, name) - const_container = lazy_container_addition!( - container, - PiecewiseLinearBlockDecrementalOfferConstraint, - T, - axes(variables)...; - meta = name, - ) - add_pwl_block_offer_constraints!( - get_jump_model(container), - const_container, - name, - period, - variables[name, period], - pwl_vars, - break_points, - ) - return -end - -""" -Return `(breakpoints, slopes)` for an ORDC service at time step `t`. For a static -`ReserveDemandCurve` the data is read directly from the curve; for a time-varying -`ReserveDemandTimeSeriesCurve` it is read from the pre-allocated decremental -slope/breakpoint parameter arrays. Mirrors IOM's device-side `_get_raw_pwl_data`. -""" -function _get_reserve_pwl_data( - container::OptimizationContainer, - component::T, - variable_cost::PSY.CostCurve, - t::Int, -) where {T <: Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}} - base_power = get_model_base_power(container) - device_base_power = PSY.get_base_power(component, PSY.NU) - unit_system = PSY.get_power_units(variable_cost) - - if !IS.is_time_series_backed(variable_cost) - cost_component = PSY.get_function_data(PSY.get_value_curve(variable_cost)) - breakpoint_cost_component = IS.get_x_coords(cost_component) - slope_cost_component = IS.get_y_coords(cost_component) - else - name = PSY.get_name(component) - SlopeParam = DecrementalPiecewiseLinearSlopeParameter - slope_arr = get_parameter_array(container, SlopeParam, T, name) - slope_mult = get_parameter_multiplier_array(container, SlopeParam, T, name) - @assert size(slope_arr) == size(slope_mult) - slope_cost_component = - [slope_arr[name, s, t] * slope_mult[name, s, t] for s in axes(slope_arr)[2]] - - BreakpointParam = DecrementalPiecewiseLinearBreakpointParameter - bp_arr = get_parameter_array(container, BreakpointParam, T, name) - bp_mult = get_parameter_multiplier_array(container, BreakpointParam, T, name) - @assert size(bp_arr) == size(bp_mult) - breakpoint_cost_component = - [bp_arr[name, p, t] * bp_mult[name, p, t] for p in axes(bp_arr)[2]] - - @assert length(slope_cost_component) == length(breakpoint_cost_component) - 1 - end - - breakpoints, slopes = get_piecewise_curve_per_system_unit( - breakpoint_cost_component, - slope_cost_component, - unit_system, - base_power, - device_base_power, - ) - return breakpoints, slopes -end - function add_pwl_term_delta!( container::OptimizationContainer, component::T, - cost_data::PSY.CostCurve, + ::PSY.CostCurve, ::Type{U}, ::Type{V}, ) where { @@ -389,6 +321,10 @@ function add_pwl_term_delta!( U <: VariableType, V <: AbstractServiceFormulation, } + dir = _reserve_offer_direction(component) + # objective_function_multiplier(U, V) == IOM._objective_sign(dir) for the + # decremental ORDC offer (OBJECTIVE_FUNCTION_NEGATIVE); kept as the service-side + # multiplier API. multiplier = objective_function_multiplier(U, V) resolution = get_resolution(container) dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR @@ -396,17 +332,27 @@ function add_pwl_term_delta!( time_steps = get_time_steps(container) pwl_cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) for t in time_steps - break_points, slopes = _get_reserve_pwl_data(container, component, cost_data, t) + break_points, slopes = IOM._get_pwl_data(dir, container, component, t; meta = name) pwl_vars = add_pwl_variables_delta!( container, - PiecewiseLinearBlockDecrementalOffer, + IOM._block_offer_var(dir), T, name, t, length(slopes); upper_bound = Inf, ) - add_pwl_constraint_delta!(container, component, U, break_points, pwl_vars, t) + add_pwl_constraint_delta!( + container, + component, + U, + V, + break_points, + pwl_vars, + t, + IOM._block_offer_constraint(dir); + meta = name, + ) pwl_cost_expressions[t] = get_pwl_cost_expression_delta(pwl_vars, slopes, multiplier * dt) end diff --git a/src/common_models/market_bid_plumbing.jl b/src/common_models/market_bid_plumbing.jl index a09e23e..a50b882 100644 --- a/src/common_models/market_bid_plumbing.jl +++ b/src/common_models/market_bid_plumbing.jl @@ -104,6 +104,15 @@ get_offer_curves(::IOM.DecrementalOffer, op_cost::PSY.OfferCurveCost) = get_offer_curves(::IOM.IncrementalOffer, op_cost::PSY.OfferCurveCost) = get_output_offer_curves(op_cost) +# direction and ORDC reserve service: the demand curve is carried on the service's +# `variable` field (a CostCurve, static or time-series-backed) rather than split +# across incremental/decremental sides. Direction is irrelevant to the lookup; the +# service-side direction trait (`_reserve_offer_direction`) is decremental. +get_offer_curves( + ::IOM.OfferDirection, + service::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, +) = PSY.get_variable(service) + ################################################################################# # Section 3: _get_parameter_field Dispatch Table # Maps parameter types to PSY getter functions. @@ -380,12 +389,13 @@ function IOM._get_pwl_data( dir::IOM.OfferDirection, container::OptimizationContainer, component::T, - time::Int, + time::Int; + meta = IOM.CONTAINER_KEY_EMPTY_META, ) where {T <: IS.InfrastructureSystemsComponent} name = IS.get_name(component) cost_data = get_offer_curves(dir, component) breakpoint_cost_component, slope_cost_component, unit_system = - IOM._get_raw_pwl_data(dir, container, T, name, cost_data, time) + IOM._get_raw_pwl_data(dir, container, T, name, cost_data, time; meta = meta) breakpoints, slopes = IOM.get_piecewise_curve_per_system_unit( breakpoint_cost_component, @@ -397,14 +407,16 @@ function IOM._get_pwl_data( return breakpoints, slopes end -# static curve: read directly from the cost curve +# static curve: read directly from the cost curve. `meta` is accepted (and +# ignored) so the kwarg threaded by `_get_pwl_data` resolves for either branch. function IOM._get_raw_pwl_data( ::IOM.OfferDirection, ::OptimizationContainer, ::Type{<:IS.InfrastructureSystemsComponent}, ::String, cost_data::IS.CostCurve{IS.PiecewiseIncrementalCurve}, - ::Int, + ::Int; + meta = IOM.CONTAINER_KEY_EMPTY_META, ) cost_component = IS.get_function_data(IS.get_value_curve(cost_data)) return IS.get_x_coords(cost_component), diff --git a/src/services_models/reserves.jl b/src/services_models/reserves.jl index 9993508..1340b15 100644 --- a/src/services_models/reserves.jl +++ b/src/services_models/reserves.jl @@ -38,6 +38,11 @@ objective_function_multiplier(::Type{ServiceRequirementVariable}, ::Type{Stepwis uses_compact_power(::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, ::StepwiseCostReserve)=false get_multiplier_value(::Type{<:AbstractPiecewiseLinearBreakpointParameter}, ::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, ::Type{<:AbstractReservesFormulation}) = 1.0 get_multiplier_value(::Type{<:AbstractPiecewiseLinearSlopeParameter}, ::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}, ::Type{<:AbstractReservesFormulation}) = 1.0 +# ORDC demand curves are willingness-to-pay (concave), i.e. a decremental offer. +# Routes the reserve PWL cost path through IOM's OfferDirection dispatch; making +# this incremental is a one-line change here. Mirrors `_onvar_offer_direction` / +# `_vom_offer_direction` in market_bid_overrides.jl. +_reserve_offer_direction(::Union{PSY.ReserveDemandCurve, PSY.ReserveDemandTimeSeriesCurve}) = IOM.DecrementalOffer() #! format: on function get_initial_conditions_service_model( @@ -586,10 +591,8 @@ function process_stepwise_cost_reserve_parameters!( model::ServiceModel, service::D, ) where {D <: PSY.ReserveDemandTimeSeriesCurve} - for param in ( - DecrementalPiecewiseLinearBreakpointParameter, - DecrementalPiecewiseLinearSlopeParameter, - ) + dir = _reserve_offer_direction(service) + for param in (IOM._breakpoint_param(dir), IOM._slope_param(dir)) add_parameters!(container, param, service, model) end return From 43992cf7d7ff55b4db46a23dcc3ad5f20d173221 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 16 Jun 2026 15:23:26 -0400 Subject: [PATCH 07/11] Port PSI#1633: model_all_branches + reduction pruning, SC slack aliasing Ports the still-missing pieces of PowerSimulations.jl PR #1633 into POM, adapted to the POM/IOM split (the core per-contingency SC slack feature was already present on this branch). - MonitoredLine `model_all_branches` attribute (AC_branches.jl): pins both endpoint buses of every covered MonitoredLine irreducible so zero-impedance lines survive the network reduction. Forwarded `irreducible_buses` to the no-reduction VirtualPTDF path so the pin actually takes effect there. - Network-reduction pruning/warnings (instantiate_network_model.jl): drop and warn about branch types fully merged away by the reduction, and warn about individual MonitoredLines partially merged away while their type survives. - Multi-type outage post-contingency slack aliasing (security_constrained_branch.jl): alias shared slack refs into the reusing branch type's container so has_container_key/get_variable stay consistent regardless of build order. - Tests for the SC slack reach, multi-type slack aliasing, and the MonitoredLine reduction retain/prune/warn behavior. Companion IOM change (compute_conflict! hardening + covariant check_conflict_status) lands on InfrastructureOptimizationModels.jl ac/ordc. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/ac_transmission_models/AC_branches.jl | 31 +++ .../security_constrained_branch.jl | 58 ++++- .../instantiate_network_model.jl | 132 ++++++++++- ...ransmission_security_constrained_models.jl | 222 ++++++++++++++++++ test/test_device_branch_constructors.jl | 113 +++++++++ 5 files changed, 552 insertions(+), 4 deletions(-) diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index d53a24a..22011a0 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -103,6 +103,37 @@ function get_default_attributes( ) end +""" +`MonitoredLine` DeviceModel attribute. When `true`, both endpoint buses of every +monitored line are pinned irreducible so zero-impedance lines survive the network +reduction. Defaults to `false` (such lines are reduced away and not modeled). For +the "base case flowgate" use case. +""" +const MODEL_ALL_BRANCHES_KEY = "model_all_branches" + +# Specialize the generic `ACTransmission` defaults for `MonitoredLine` to add +# `MODEL_ALL_BRANCHES_KEY` (default `false`) alongside the inherited keys. +function get_default_attributes( + ::Type{PSY.MonitoredLine}, + ::Type{V}, +) where {V <: AbstractBranchFormulation} + return Dict{String, Any}( + PARALLEL_BRANCH_MAX_RATING_KEY => "single_element_contingency", + MODEL_ALL_BRANCHES_KEY => false, + ) +end + +function get_default_attributes( + ::Type{PSY.MonitoredLine}, + ::Type{V}, +) where {V <: AbstractSecurityConstrainedStaticBranch} + return Dict{String, Any}( + PARALLEL_BRANCH_MAX_RATING_KEY => "single_element_contingency", + "include_planned_outages" => false, + MODEL_ALL_BRANCHES_KEY => false, + ) +end + # Resolve the per-DeviceModel attribute to one of the explicit PNM rating functions. # `MixedBranchesParallel` ignores the attribute and always uses the plain sum, since # the constituent branches may carry different DeviceModel preferences and there is diff --git a/src/ac_transmission_models/security_constrained_branch.jl b/src/ac_transmission_models/security_constrained_branch.jl index 1e50b92..8d96268 100644 --- a/src/ac_transmission_models/security_constrained_branch.jl +++ b/src/ac_transmission_models/security_constrained_branch.jl @@ -215,6 +215,40 @@ function _find_shared_post_contingency_constraint_sources( return src_lb, src_ub end +""" +Locate the post-contingency slack variable containers that a shared constraint for +`(outage_id, name, t)` already references, registered by the source DeviceModel under +a component type other than `V`. Returns `(ub_source, lb_source)`; either slot is +`nothing` when the shared constraint was built without that slack (the source model +had `use_slacks = false`). Lets a reusing model alias those refs into its own slack +container so `get_variable`/`has_container_key` stay consistent with the constraints +it reuses, regardless of branch-model build order. +""" +function _find_shared_post_contingency_slack_sources( + container::OptimizationContainer, + ::Type{V}, + outage_id::String, + name::String, + t::Int, +) where {V <: PSY.ACTransmission} + target = (outage_id, name, t) + src_ub = nothing + src_lb = nothing + for (key, vc) in IOM.get_variables(container) + get_component_type(key) === V && continue + entry = get_entry_type(key) + if entry === PostContingencyFlowActivePowerSlackUpperBound && + haskey(vc.data, target) + src_ub = vc + elseif entry === PostContingencyFlowActivePowerSlackLowerBound && + haskey(vc.data, target) + src_lb = vc + end + !isnothing(src_ub) && !isnothing(src_lb) && break + end + return src_ub, src_lb +end + """ Fast-path precheck: returns `true` iff any container of entry type `T` exists under a component type other than `V`. @@ -404,14 +438,32 @@ function add_constraints!( container, T, V, outage_id, name, first(time_steps), ) if !isnothing(src_lb) && !isnothing(src_ub) - # Reuse the first claimer's constraint refs verbatim; its - # slacks were already created and penalized, so do NOT add - # new ones here. + # Reuse the shared constraint refs verbatim. The shared + # constraint already references the source model's slack (if + # any); alias those slack refs into this model's container so + # `get_variable`/`has_container_key` for `V` stay consistent + # with the constraints it reuses. + src_slack_ub = nothing + src_slack_lb = nothing + if use_slacks + src_slack_ub, src_slack_lb = + _find_shared_post_contingency_slack_sources( + container, V, outage_id, name, first(time_steps), + ) + end for t in time_steps con_ub[outage_id, name, t] = src_ub.data[(outage_id, name, t)] con_lb[outage_id, name, t] = src_lb.data[(outage_id, name, t)] + isnothing(src_slack_ub) || ( + slack_ub[outage_id, name, t] = + src_slack_ub.data[(outage_id, name, t)] + ) + isnothing(src_slack_lb) || ( + slack_lb[outage_id, name, t] = + src_slack_lb.data[(outage_id, name, t)] + ) end continue end diff --git a/src/network_models/instantiate_network_model.jl b/src/network_models/instantiate_network_model.jl index 68b485d..e046830 100644 --- a/src/network_models/instantiate_network_model.jl +++ b/src/network_models/instantiate_network_model.jl @@ -173,9 +173,121 @@ function _get_irreducible_buses_due_to_monitored_components( end end _add_outage_monitored_irreducible_buses!(irreducible_buses, sys, branch_models) + # `model_all_branches` MonitoredLine models pin their lines so zero-impedance + # ones survive the reduction instead of being merged away. + _add_model_all_branches_irreducible_buses!(irreducible_buses, branch_models) return collect(irreducible_buses) end +# Pin both endpoint buses of every branch a `model_all_branches` MonitoredLine model +# covers. Dispatch on the model type so it is a no-op for other branch types. +function _add_model_all_branches_irreducible_buses!( + irreducible_buses::Set{Int64}, + branch_models::BranchModelContainer, +) + for m in values(branch_models) + _pin_model_all_branches!(irreducible_buses, m) + end + return +end + +_pin_model_all_branches!(::Set{Int64}, ::DeviceModel) = nothing + +function _pin_model_all_branches!( + irreducible_buses::Set{Int64}, + m::DeviceModel{PSY.MonitoredLine}, +) + get_attribute(m, MODEL_ALL_BRANCHES_KEY) === true || return + # The device cache is the modeled set (available + filter_function). + for branch in get_device_cache(m) + _push_component_buses!(irreducible_buses, branch) + end + return +end + +# Drop (and warn about) any branch type whose components were all merged away by the +# reduction — e.g. a lone zero-impedance monitored line. Such a type has no surviving +# arc in `name_to_arc_maps`, so building its flow vars/constraints would fail. Absence +# from the map alone is not enough: types that never use it (e.g. HVDC) are also +# absent, so we prune only when an endpoint bus was actually removed by the reduction. +# Uncommon; `model_all_branches` keeps such lines instead. +function _prune_fully_reduced_branch_models!( + network_model::NetworkModel, + branch_models::BranchModelContainer, +) + merged_buses = Set{Int64}() + for removed in values(PNM.get_bus_reduction_map(network_model.network_reduction)) + union!(merged_buses, removed) + end + isempty(merged_buses) && return + name_to_arc_maps = PNM.get_name_to_arc_maps(network_model.network_reduction) + pruned = DataType[] + for branch_type in network_model.modeled_branch_types + branch_type <: PSY.ACTransmission || continue + haskey(branch_models, nameof(branch_type)) || continue + survived = get(name_to_arc_maps, branch_type, nothing) + isnothing(survived) || isempty(survived) || continue + buses = Set{Int64}() + for component in get_device_cache(branch_models[nameof(branch_type)]) + _push_component_buses!(buses, component) + end + isdisjoint(buses, merged_buses) && continue + push!(pruned, branch_type) + end + for branch_type in pruned + @warn "All components of branch type $(branch_type) were merged away by the " * + "network reduction (e.g. a zero-impedance branch merge). The " * + "$(branch_type) DeviceModel is dropped from the template and will not " * + "be modeled. Use the `model_all_branches` attribute on a MonitoredLine " * + "model to retain such branches through the reduction." + delete!(branch_models, nameof(branch_type)) + filter!(!=(branch_type), network_model.modeled_branch_types) + end + return +end + +# Warn about individual monitored lines the reduction merged away while their type +# still has surviving members. The whole-type prune above misses this partial case, +# so without a message the dropped line is silently unmodeled. Suggest +# `model_all_branches` to retain it. +function _warn_partially_reduced_monitored_lines!( + network_model::NetworkModel, + branch_models::BranchModelContainer, +) + removed_arcs = PNM.get_removed_arcs(network_model.network_reduction) + isempty(removed_arcs) && return + for m in values(branch_models) + _warn_reduced_monitored_lines!(removed_arcs, m) + end + return +end + +_warn_reduced_monitored_lines!(removed_arcs::Set{Tuple{Int, Int}}, ::DeviceModel) = nothing + +function _warn_reduced_monitored_lines!( + removed_arcs::Set{Tuple{Int, Int}}, + m::DeviceModel{PSY.MonitoredLine}, +) + dropped = [ + PSY.get_name(ml) for ml in get_device_cache(m) if + _branch_arc_removed(ml, removed_arcs) + ] + isempty(dropped) && return + @warn "MonitoredLine(s) $(dropped) were merged away by the network reduction " * + "(near-zero impedance) and will not be modeled or monitored, though other " * + "MonitoredLines remain. Set the `model_all_branches` attribute on the " * + "MonitoredLine DeviceModel to force all monitored lines to be modeled " * + "through the reduction." + return +end + +function _branch_arc_removed(branch::PSY.Branch, removed_arcs) + arc = PSY.get_arc(branch) + from = PSY.get_number(PSY.get_from(arc)) + to = PSY.get_number(PSY.get_to(arc)) + return (from, to) in removed_arcs || (to, from) in removed_arcs +end + function _get_unmodeled_branch_types( branch_models::BranchModelContainer, sys::PSY.System, @@ -261,6 +373,11 @@ function IOM.instantiate_network_model!( model.network_reduction, IOM._get_filters(branch_models), ) + # After the reduction is known and the branch maps populated, before the + # device constructors run: drop branch types fully merged away (else their + # flow vars/constraints would fail to build) and warn about partial drops. + _prune_fully_reduced_branch_models!(model, branch_models) + _warn_partially_reduced_monitored_lines!(model, branch_models) _reset_reduced_branch_tracker!(model, number_of_steps) return end @@ -367,7 +484,15 @@ function IOM.instantiate_network_model!( irreducible_buses = irreducible_buses, ) else - ptdf = PNM.VirtualPTDF(sys; tol = PTDF_ZERO_TOL) + # No radial/degree-two reduction requested, but irreducible buses may + # still be pinned (e.g. `model_all_branches` MonitoredLines, outages, + # DLRs). Forward them so the base zero-impedance branch coalescing + # cannot merge a pinned bus away. + ptdf = PNM.VirtualPTDF( + sys; + tol = PTDF_ZERO_TOL, + irreducible_buses = irreducible_buses, + ) end model.PTDF_matrix = ptdf model.network_reduction = deepcopy(ptdf.network_reduction_data) @@ -418,6 +543,11 @@ function IOM.instantiate_network_model!( model.network_reduction, IOM._get_filters(branch_models), ) + # After the reduction is known and the branch maps populated, before the + # device constructors run: drop branch types fully merged away (else their + # flow vars/constraints would fail to build) and warn about partial drops. + _prune_fully_reduced_branch_models!(model, branch_models) + _warn_partially_reduced_monitored_lines!(model, branch_models) _reset_reduced_branch_tracker!(model, number_of_steps) return end diff --git a/test/test_ac_transmission_security_constrained_models.jl b/test/test_ac_transmission_security_constrained_models.jl index cc2876b..0922bd2 100644 --- a/test/test_ac_transmission_security_constrained_models.jl +++ b/test/test_ac_transmission_security_constrained_models.jl @@ -1260,3 +1260,225 @@ end @test lim.max ≈ rb @test lim.min ≈ -rb end + +# Attach a single outage (monitoring every branch) to every branch of the system, +# so every branch is both an outaged and a monitored component. +function _attach_all_branch_outages!(sys) + all_branches = collect(PSY.get_components(PSY.ACTransmission, sys)) + for branch in all_branches + outage = PSY.GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = all_branches, + ) + PSY.add_supplemental_attribute!(sys, branch, outage) + end + return sys +end + +# Every monitored post-contingency rate constraint references the slack +# variable for its `(outage_id, name, t)`. A constraint without a slack term +# means `use_slacks` was honored only for the pre-contingency flow limits. +function _post_contingency_constraints_have_slack(container, slack_refs) + has_slack(ref) = + any(haskey(JuMP.constraint_object(ref).func.terms, v) for v in slack_refs) + all_have = true + n = 0 + for meta in ("lb", "ub") + cons = IOM.get_constraints(container)[IOM.ConstraintKey( + POM.PostContingencyFlowRateConstraint, PSY.Line, meta, + )] + for ref in values(cons.data) + n += 1 + all_have &= has_slack(ref) + end + end + return n, all_have +end + +@testset "Security Constrained branch slacks reach post-contingency flow constraints" begin + c_sys5 = _attach_all_branch_outages!(PSB.build_system(PSITestSystems, "c_sys5")) + + function _build_sc(use_slacks) + template = get_thermal_dispatch_template_network( + NetworkModel(PTDFPowerModel; PTDF_matrix = PNM.PTDF(c_sys5)), + ) + set_device_model!( + template, + DeviceModel( + PSY.Line, + POM.SecurityConstrainedStaticBranch; + use_slacks = use_slacks, + ), + ) + set_device_model!( + template, + DeviceModel( + PSY.Transformer2W, + POM.SecurityConstrainedStaticBranch; + use_slacks = use_slacks, + ), + ) + set_device_model!( + template, + DeviceModel( + PSY.TapTransformer, + POM.SecurityConstrainedStaticBranch; + use_slacks = use_slacks, + ), + ) + model = DecisionModel(template, c_sys5; optimizer = HiGHS_optimizer) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + return model + end + + # use_slacks = true: per-contingency slack variables exist and are wired + # into every post-contingency rate inequality. + model = _build_sc(true) + container = IOM.get_optimization_container(model) + @test IOM.has_container_key( + container, + POM.PostContingencyFlowActivePowerSlackUpperBound, + PSY.Line, + ) + @test IOM.has_container_key( + container, + POM.PostContingencyFlowActivePowerSlackLowerBound, + PSY.Line, + ) + slack_ub = + IOM.get_variable( + container, + POM.PostContingencyFlowActivePowerSlackUpperBound, + PSY.Line, + ) + slack_lb = + IOM.get_variable( + container, + POM.PostContingencyFlowActivePowerSlackLowerBound, + PSY.Line, + ) + @test !isempty(slack_ub.data) + @test !isempty(slack_lb.data) + slack_refs = Set{JuMP.VariableRef}() + union!(slack_refs, values(slack_ub.data)) + union!(slack_refs, values(slack_lb.data)) + n, all_have = _post_contingency_constraints_have_slack(container, slack_refs) + @test n > 0 + @test all_have + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + # use_slacks = false: no post-contingency slack variables, hard limits. + model_hard = _build_sc(false) + container_hard = IOM.get_optimization_container(model_hard) + @test !IOM.has_container_key( + container_hard, + POM.PostContingencyFlowActivePowerSlackUpperBound, + PSY.Line, + ) + @test !IOM.has_container_key( + container_hard, + POM.PostContingencyFlowActivePowerSlackLowerBound, + PSY.Line, + ) +end + +# An outage attached to components of more than one branch type makes multiple SC +# DeviceModels claim the same monitored `(outage, name, t)` constraints. The second +# model to build reuses the first's constraint refs; the slack-container aliasing +# registers those slacks under the reusing type too, so `has_container_key` / +# `get_variable` stay consistent regardless of the Dict-ordered build order. +@testset "Multi-type outage post-contingency slacks reach every reusing branch type" begin + c_sys14 = PSB.build_system(PSITestSystems, "c_sys14") + all_branches = collect(PSY.get_components(PSY.ACTransmission, c_sys14)) + line = PSY.get_component(PSY.Line, c_sys14, "Line1") + transformer = first(PSY.get_components(PSY.Transformer2W, c_sys14)) + # One outage attached to both a Line and a Transformer2W: a multi-type outage. + outage = PSY.GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = all_branches, + ) + PSY.add_supplemental_attribute!(c_sys14, line, outage) + PSY.add_supplemental_attribute!(c_sys14, transformer, outage) + + template = get_thermal_dispatch_template_network( + NetworkModel( + PTDFPowerModel; + PTDF_matrix = PNM.VirtualPTDF(c_sys14), + MODF_matrix = PNM.VirtualMODF(c_sys14), + ), + ) + set_device_model!( + template, + DeviceModel(PSY.Line, POM.SecurityConstrainedStaticBranch; use_slacks = true), + ) + set_device_model!( + template, + DeviceModel( + PSY.Transformer2W, + POM.SecurityConstrainedStaticBranch; + use_slacks = true, + ), + ) + set_device_model!( + template, + DeviceModel( + PSY.TapTransformer, + POM.SecurityConstrainedStaticBranch; + use_slacks = true, + ), + ) + model = DecisionModel(template, c_sys14; optimizer = HiGHS_optimizer) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + container = IOM.get_optimization_container(model) + + # Both reusing branch types own non-empty slack containers: the second to build + # aliases the first's refs instead of registering an empty container. + slack_refs = Set{JuMP.VariableRef}() + for V in (PSY.Line, PSY.Transformer2W), + S in ( + POM.PostContingencyFlowActivePowerSlackUpperBound, + POM.PostContingencyFlowActivePowerSlackLowerBound, + ) + + @test IOM.has_container_key(container, S, V) + slacks = IOM.get_variable(container, S, V) + @test !isempty(slacks.data) + union!(slack_refs, values(slacks.data)) + end + + # The aliased slacks are shared between the two component-type containers. + line_ub = IOM.get_variable( + container, + POM.PostContingencyFlowActivePowerSlackUpperBound, + PSY.Line, + ) + transformer_ub = IOM.get_variable( + container, + POM.PostContingencyFlowActivePowerSlackUpperBound, + PSY.Transformer2W, + ) + @test !isempty(intersect(Set(values(line_ub.data)), Set(values(transformer_ub.data)))) + + # Every post-contingency rate constraint, including the reused ones stored under + # the second component type, references a registered slack. + has_slack(ref) = + any(haskey(JuMP.constraint_object(ref).func.terms, v) for v in slack_refs) + n = 0 + all_have = true + for V in (PSY.Line, PSY.Transformer2W), meta in ("lb", "ub") + cons = IOM.get_constraints(container)[IOM.ConstraintKey( + POM.PostContingencyFlowRateConstraint, V, meta, + )] + for ref in values(cons.data) + n += 1 + all_have &= has_slack(ref) + end + end + @test n > 0 + @test all_have + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED +end diff --git a/test/test_device_branch_constructors.jl b/test/test_device_branch_constructors.jl index bce422b..74c96b5 100644 --- a/test/test_device_branch_constructors.jl +++ b/test/test_device_branch_constructors.jl @@ -897,3 +897,116 @@ end IOM.ModelBuildStatus.BUILT @test solve!(model_ac) == IOM.RunStatus.SUCCESSFULLY_FINALIZED end + +# A bus is "merged away" by the reduction when it appears in a value set of the bus +# reduction map (coalesced into a representative key). The retained representatives +# are the keys; `get_removed_buses` is unrelated to zero-impedance coalescing here. +_bus_merged_away(nrd, b) = any(b in s for s in values(PNM.get_bus_reduction_map(nrd))) + +# A zero-impedance `MonitoredLine` is merged away by the reduction. With +# `model_all_branches = true` its buses are pinned so it survives and is modeled; +# with the default `false` it is reduced away and its (sole-of-type) DeviceModel is +# pruned. Both build; only the flow-rate constraint differs. +@testset "MonitoredLine model_all_branches retains zero-impedance branch" begin + function _build_zib_monitored_line(model_all_branches) + sys = PSB.build_system(PSITestSystems, "c_sys5_ml") + # Force MonitoredLine "1" to be zero-impedance: r == 0 and a tiny reactance + # push it above the zero-impedance threshold, so the reduction merges its + # endpoints unless they are pinned irreducible. + ml = PSY.get_component(MonitoredLine, sys, "1") + PSY.set_r!(ml, 0.0 * PSY.SU) + PSY.set_x!(ml, 1e-5 * PSY.SU) + # No `PTDF_matrix` provided, so a VirtualPTDF is built and the reduction runs. + template = get_thermal_dispatch_template_network(NetworkModel(PTDFPowerModel)) + set_device_model!( + template, + DeviceModel( + MonitoredLine, + StaticBranch; + attributes = Dict{String, Any}( + "model_all_branches" => model_all_branches, + ), + ), + ) + model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + status = build!(model; output_dir = mktempdir(; cleanup = true)) + return model, ml, status + end + + # The attribute defaults to false. + default_model = DeviceModel(MonitoredLine, StaticBranch) + @test POM.get_attribute(default_model, "model_all_branches") == false + + # true: line retained, build succeeds, buses not merged, line modeled. + model, ml, status = _build_zib_monitored_line(true) + @test status == IOM.ModelBuildStatus.BUILT + arc = PSY.get_arc(ml) + from_bus = PSY.get_number(PSY.get_from(arc)) + to_bus = PSY.get_number(PSY.get_to(arc)) + nm = IOM.get_network_model(IOM.get_template(model)) + nrd = IOM.get_PTDF_matrix(nm).network_reduction_data + @test !_bus_merged_away(nrd, from_bus) + @test !_bus_merged_away(nrd, to_bus) + @test haskey(IOM.get_branch_models(IOM.get_template(model)), :MonitoredLine) + container = IOM.get_optimization_container(model) + @test IOM.has_container_key(container, FlowRateConstraint, MonitoredLine, "ub") + + # false (default): line reduced away, its sole-of-type DeviceModel pruned, + # build succeeds with no MonitoredLine flow-rate constraint. + model_default, ml_default, status_default = _build_zib_monitored_line(false) + @test status_default == IOM.ModelBuildStatus.BUILT + nm_d = IOM.get_network_model(IOM.get_template(model_default)) + nrd_d = IOM.get_PTDF_matrix(nm_d).network_reduction_data + @test _bus_merged_away(nrd_d, PSY.get_number(PSY.get_to(PSY.get_arc(ml_default)))) + @test !haskey(IOM.get_branch_models(IOM.get_template(model_default)), :MonitoredLine) + container_default = IOM.get_optimization_container(model_default) + @test !IOM.has_container_key( + container_default, + FlowRateConstraint, + MonitoredLine, + "ub", + ) +end + +# Partial reduction: with multiple monitored lines and `model_all_branches = false`, +# a single near-zero-impedance monitored line is merged away while the type survives. +# Whereas a fully-reduced type is pruned, here the reduced line is silently unmodeled, +# so the build must still succeed and emit an actionable warning that names the line +# and points the user at `model_all_branches`. +@testset "MonitoredLine partial reduction warns and drops only the reduced line" begin + sys = PSB.build_system(PSITestSystems, "c_sys5_ml") + # MonitoredLine "1" forced near-zero impedance so the reduction merges it away. + ml = PSY.get_component(MonitoredLine, sys, "1") + PSY.set_r!(ml, 0.0 * PSY.SU) + PSY.set_x!(ml, 1e-5 * PSY.SU) + # A second MonitoredLine (converted from a healthy Line) keeps the type non-empty. + line = first(PSY.get_components(Line, sys)) + survivor = PSY.get_name(line) + PSY.convert_component!( + sys, + line, + MonitoredLine; + flow_limits = (from_to = 1.0, to_from = 1.0), + ) + + template = get_thermal_dispatch_template_network(NetworkModel(PTDFPowerModel)) + set_device_model!(template, DeviceModel(MonitoredLine, StaticBranch)) + model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + output_dir = mktempdir(; cleanup = true) + @test build!(model; output_dir = output_dir) == IOM.ModelBuildStatus.BUILT + + # The surviving monitored line is modeled; the merged-away one is dropped. + container = IOM.get_optimization_container(model) + constraint_names = axes( + IOM.get_constraints(container)[IOM.ConstraintKey( + FlowRateConstraint, MonitoredLine, "ub", + )], + )[1] + @test survivor in constraint_names + @test !("1" in constraint_names) + + # The drop is reported with an actionable warning naming the line. + log_contents = read(joinpath(output_dir, "operation_problem.log"), String) + @test occursin("MonitoredLine(s) [\"1\"]", log_contents) + @test occursin("model_all_branches", log_contents) +end From bb1a4ab46db8ac267223799db759dd98b3415d1c Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 16 Jun 2026 18:48:39 -0400 Subject: [PATCH 08/11] moved stuff to PNM; copilot review --- Project.toml | 2 +- src/ac_transmission_models/AC_branches.jl | 210 +++--------------- .../security_constrained_branch.jl | 30 ++- .../instantiate_network_model.jl | 6 +- test/Project.toml | 2 +- test/test_native_dcp_acp_models.jl | 20 +- 6 files changed, 60 insertions(+), 210 deletions(-) diff --git a/Project.toml b/Project.toml index 04b526f..63a2fa6 100644 --- a/Project.toml +++ b/Project.toml @@ -29,7 +29,7 @@ InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/ PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} # todo: repin once iom merged InfrastructureOptimizationModels = {rev = "ac/ordc", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} -PowerNetworkMatrices = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"} +PowerNetworkMatrices = {rev = "ac/port-branch-admittance", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"} [extensions] PowerFlowsExt = "PowerFlows" diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index 22011a0..11f6d0d 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -1424,107 +1424,22 @@ function add_to_objective_function!( return end -# Flip a π-admittance tuple to the opposite orientation (from<->to). Reduced equivalents -# may be keyed by an arc whose orientation is reversed vs the surviving branch's retained -# from->to; reorient so coefficients match (from_bus, to_bus). g/b are symmetric; from/to -# shunts swap; phase shift negates. Reduced line equivalents have tap == 1. -function _reverse_admittance(adm) - @assert adm.tap == 1.0 "Cannot reorient a reduced arc with a non-unit tap ($(adm.tap))." - return ( - g = adm.g, - b = adm.b, - g_fr = adm.g_to, - b_fr = adm.b_to, - g_to = adm.g_fr, - b_to = adm.b_fr, - tap = adm.tap, - shift = -adm.shift, - ) -end - -# ── Native π-admittance helpers ────────────────────────────────────────────── -# The native DCP/ACP branch builders need each branch's π-model admittance as a -# `(g, b, g_fr, b_fr, g_to, b_to, tap, shift)` NamedTuple (PowerModels convention: -# `g + im*b == 1/(r + im*x)` is the series admittance, `*_fr`/`*_to` are the from/to -# shunts, `tap`/`shift` the transformer ratio/phase). -_branch_tap(::PSY.ACTransmission) = 1.0 -_branch_tap(b::PSY.TapTransformer) = PSY.get_tap(b) -_branch_tap(b::PSY.PhaseShiftingTransformer) = PSY.get_tap(b) -_branch_shift(::PSY.ACTransmission) = 0.0 -_branch_shift(b::PSY.PhaseShiftingTransformer) = PSY.get_α(b) - -# Lines (and other shunt-on-both-ends branches): π-shunt split from/to. -function _branch_admittance(b::PSY.ACTransmission) - ys = 1.0 / (PSY.get_r(b, PSY.SU) + PSY.get_x(b, PSY.SU) * im) - b_sh = PSY.get_b(b, PSY.SU) - return ( - g = real(ys), b = imag(ys), - g_fr = 0.0, b_fr = b_sh.from, g_to = 0.0, b_to = b_sh.to, - tap = 1.0, shift = 0.0, - ) -end - -# Two-winding transformers: magnetizing shunt allocated to the primary (from) side only. -function _branch_admittance( - b::Union{PSY.Transformer2W, PSY.TapTransformer, PSY.PhaseShiftingTransformer}, -) - ys = 1.0 / (PSY.get_r(b, PSY.SU) + PSY.get_x(b, PSY.SU) * im) - yt = PSY.get_primary_shunt(b, PSY.SU) - return ( - g = real(ys), b = imag(ys), - g_fr = real(yt), b_fr = imag(yt), g_to = 0.0, b_to = 0.0, - tap = _branch_tap(b), shift = _branch_shift(b), - ) -end - -# Reduction-aggregated arc (series chain / parallel group): use PNM's reduction-aware -# equivalent π-parameters. Series/parallel equivalents of lines carry `tap == 1`. -function _segment_admittance(segment, nr::PNM.NetworkReductionData) - eb = PNM.get_equivalent_physical_branch_parameters(segment, nr) - ys = 1.0 / (PNM.get_equivalent_r(eb) + PNM.get_equivalent_x(eb) * im) - return ( - g = real(ys), b = imag(ys), - g_fr = PNM.get_equivalent_g_from(eb), b_fr = PNM.get_equivalent_b_from(eb), - g_to = PNM.get_equivalent_g_to(eb), b_to = PNM.get_equivalent_b_to(eb), - tap = PNM.get_equivalent_tap(eb), shift = PNM.get_equivalent_shift(eb), - ) -end - -# Reduction-aware admittance for the retained arc (from_no -> to_no). Returns the PNM -# series/parallel equivalent π-tuple oriented from->to when the arc was aggregated by -# reduction, or `nothing` when the arc is direct (caller falls back to the branch's own). -function _reduced_arc_admittance(nr::PNM.NetworkReductionData, from_no::Int, to_no::Int) - series_map = PNM.get_series_branch_map(nr) - parallel_map = PNM.get_parallel_branch_map(nr) - arc = (from_no, to_no) - rev = (to_no, from_no) - if haskey(series_map, arc) - return _segment_admittance(series_map[arc], nr) - elseif haskey(series_map, rev) - return _reverse_admittance(_segment_admittance(series_map[rev], nr)) - elseif haskey(parallel_map, arc) - return _segment_admittance(parallel_map[arc], nr) - elseif haskey(parallel_map, rev) - return _reverse_admittance(_segment_admittance(parallel_map[rev], nr)) - end - return nothing -end - # Admittance for `branch`'s ohm's law given its retained endpoints: the branch's own -# π-parameters for a direct/un-reduced arc, or PNM's reduction-aware equivalent for a -# series/parallel-aggregated arc. +# π-parameters (`PNM.branch_admittance`) for a direct/un-reduced arc, or PNM's +# reduction-aware equivalent (`PNM.reduced_arc_admittance`) for a series/parallel-aggregated +# arc. The π-admittance primitives themselves now live in PNM (`BranchAdmittance.jl`). function _resolve_branch_admittance(network_model, branch, from_no::Int, to_no::Int) nr = get_network_reduction(network_model) - isempty(nr) && return _branch_admittance(branch) - eq = _reduced_arc_admittance(nr, from_no, to_no) - return eq === nothing ? _branch_admittance(branch) : eq + isempty(nr) && return PNM.branch_admittance(branch) + eq = PNM.reduced_arc_admittance(nr, from_no, to_no) + return eq === nothing ? PNM.branch_admittance(branch) : eq end # One-pass per-branch network geometry for the native DCP/ACP builders. Computes each # branch's retained endpoints and reduction-aware admittance ONCE so the ohm's-law and # angle-limit builders don't each rebuild `number_to_name` and re-map endpoints. -# Each entry is a NamedTuple value-bag (same idiom as `_retained_bus` / `branch_admittance` -# in this file); `collapsed` marks branches whose endpoints fold into one retained bus. +# Each entry is a NamedTuple value-bag (same idiom as `_retained_bus`); `collapsed` marks +# branches whose endpoints fold into one retained bus. function _branch_geometries(sys::PSY.System, network_model, devices) number_to_name = _retained_number_to_name(sys, network_model) geoms = NamedTuple[] @@ -1534,7 +1449,7 @@ function _branch_geometries(sys::PSY.System, network_model, devices) collapsed = fr.number == to.number adm = if collapsed - _branch_admittance(d) + PNM.branch_admittance(d) else _resolve_branch_admittance(network_model, d, fr.number, to.number) end @@ -1554,45 +1469,9 @@ function _branch_geometries(sys::PSY.System, network_model, devices) return geoms end -""" - branch_flow_limits(branch) -> NamedTuple - -Returns directional flow limits in per-unit MVA: `(from_to::Float64, to_from::Float64)`. -For symmetric branches both fields equal `PSY.get_rating(branch)`. -""" -function branch_flow_limits end - -function branch_flow_limits(b::PSY.MonitoredLine) - fl = PSY.get_flow_limits(b, PSY.SU) - return (from_to = fl.from_to, to_from = fl.to_from) -end - -function branch_flow_limits( - b::Union{ - PSY.Line, - PSY.Transformer2W, - PSY.TapTransformer, - PSY.PhaseShiftingTransformer, - }, -) - r = PSY.get_rating(b, PSY.SU) - return (from_to = r, to_from = r) -end - -function branch_flow_limits(b::PNM.BranchesParallel) - r = PNM.get_equivalent_rating(b) - return (from_to = r, to_from = r) -end - -function branch_flow_limits(b::PNM.BranchesSeries) - r = PNM.get_equivalent_rating(b) - return (from_to = r, to_from = r) -end - -function branch_flow_limits(w::PNM.ThreeWindingTransformerWinding) - r = PNM.get_equivalent_rating(w) - return (from_to = r, to_from = r) -end +# `branch_flow_limits(branch)` (directional per-unit MVA flow limits, including the +# `BranchesParallel`/`BranchesSeries`/`ThreeWindingTransformerWinding` reduction wrappers) +# now lives in PNM (`BranchAdmittance.jl`); call `PNM.branch_flow_limits`. ################################## Native ACP apparent-power rate constraints ############### @@ -2052,54 +1931,17 @@ end ################################################################################ # ── Native three-winding helpers ───────────────────────────────────────────── -# Decompose a Transformer3W into its three wye-model windings via PNM's -# `ThreeWindingTransformerWinding`, exposing the per-winding data the native 3W builders -# need: the star-point arc (for reduction-aware bus mapping), the winding rating, and the -# winding object (for admittance). -function _three_winding_arcs(d::PSY.Transformer3W) - star_arcs = ( - PSY.get_primary_star_arc(d), - PSY.get_secondary_star_arc(d), - PSY.get_tertiary_star_arc(d), - ) - out = NamedTuple[] - for i in 1:3 - w = PNM.ThreeWindingTransformerWinding(d, i) - push!( - out, - ( - suffix = "winding_$i", - arc = star_arcs[i], - rating = PNM.get_equivalent_rating(w), - winding = w, - ), - ) - end - return out -end - -# No phase shift for windings; only the (optional) winding tap matters for the ACP π-model. -_winding_tap(::PNM.ThreeWindingTransformerWinding) = 1.0 -_winding_tap(w::PNM.ThreeWindingTransformerWinding{PSY.PhaseShiftingTransformer3W}) = - PNM.get_equivalent_tap(w) - -# Per-winding π-admittance NamedTuple (shunt on the primary side only, no phase shift). -function _winding_admittance(w::PNM.ThreeWindingTransformerWinding) - ys = 1.0 / (PNM.get_equivalent_r(w) + PNM.get_equivalent_x(w) * im) - b_sh = PNM.get_equivalent_b(w) - return ( - g = real(ys), b = imag(ys), - g_fr = 0.0, b_fr = b_sh.from, g_to = 0.0, b_to = b_sh.to, - tap = _winding_tap(w), - ) -end +# The per-winding decomposition (`PNM.three_winding_arcs`, returning a `suffix`/`arc`/ +# `rating`/`winding` NamedTuple per wye-model winding) and the per-winding π-admittance +# (`PNM.winding_admittance`) now live in PNM (`BranchAdmittance.jl`); the native 3W builders +# below call them directly. "Build the list of per-winding variable names for a set of Transformer3W devices." function _three_winding_var_names(devices) names = String[] for d in devices dname = PSY.get_name(d) - for w in _three_winding_arcs(d) + for w in PNM.three_winding_arcs(d) push!(names, dname * "_" * w.suffix) end end @@ -2168,7 +2010,7 @@ function add_to_expression!( time_steps = get_time_steps(container) for d in devices dname = PSY.get_name(d) - for w in _three_winding_arcs(d) + for w in PNM.three_winding_arcs(d) wname = dname * "_" * w.suffix from_no = PNM.get_mapped_bus_number(network_reduction, PSY.get_from(w.arc)) to_no = PNM.get_mapped_bus_number(network_reduction, PSY.get_to(w.arc)) @@ -2204,7 +2046,7 @@ for (E, V, isfrom) in ( time_steps = get_time_steps(container) for d in devices dname = PSY.get_name(d) - for w in _three_winding_arcs(d) + for w in PNM.three_winding_arcs(d) wname = dname * "_" * w.suffix terminal_bus_obj = $isfrom ? PSY.get_from(w.arc) : PSY.get_to(w.arc) bus_no = PNM.get_mapped_bus_number(network_reduction, terminal_bus_obj) @@ -2238,9 +2080,9 @@ function add_constraints!( for d in devices dname = PSY.get_name(d) - for w in _three_winding_arcs(d) + for w in PNM.three_winding_arcs(d) wname = dname * "_" * w.suffix - adm = _winding_admittance(w.winding) + adm = PNM.winding_admittance(w.winding) fr = _retained_bus(number_to_name, network_model, PSY.get_from(w.arc)) to = _retained_bus(number_to_name, network_model, PSY.get_to(w.arc)) fr.number == to.number && continue @@ -2295,9 +2137,9 @@ function add_constraints!( for d in devices dname = PSY.get_name(d) - for w in _three_winding_arcs(d) + for w in PNM.three_winding_arcs(d) wname = dname * "_" * w.suffix - adm = _winding_admittance(w.winding) + adm = PNM.winding_admittance(w.winding) g, b, g_fr, b_fr, g_to, b_to, tm = adm.g, adm.b, adm.g_fr, adm.b_fr, adm.g_to, adm.b_to, adm.tap fr = _retained_bus(number_to_name, network_model, PSY.get_from(w.arc)) @@ -2368,7 +2210,7 @@ function add_constraints!( ) for d in devices dname = PSY.get_name(d) - for w in _three_winding_arcs(d) + for w in PNM.three_winding_arcs(d) wname = dname * "_" * w.suffix for t in time_steps cons_lb[wname, t] = JuMP.@constraint( @@ -2400,7 +2242,7 @@ function add_constraints!( ) for d in devices dname = PSY.get_name(d) - for w in _three_winding_arcs(d) + for w in PNM.three_winding_arcs(d) wname = dname * "_" * w.suffix r2 = w.rating^2 for t in time_steps @@ -2430,7 +2272,7 @@ function add_constraints!( ) for d in devices dname = PSY.get_name(d) - for w in _three_winding_arcs(d) + for w in PNM.three_winding_arcs(d) wname = dname * "_" * w.suffix r2 = w.rating^2 for t in time_steps diff --git a/src/ac_transmission_models/security_constrained_branch.jl b/src/ac_transmission_models/security_constrained_branch.jl index 8d96268..cc1a0ca 100644 --- a/src/ac_transmission_models/security_constrained_branch.jl +++ b/src/ac_transmission_models/security_constrained_branch.jl @@ -438,18 +438,24 @@ function add_constraints!( container, T, V, outage_id, name, first(time_steps), ) if !isnothing(src_lb) && !isnothing(src_ub) - # Reuse the shared constraint refs verbatim. The shared - # constraint already references the source model's slack (if - # any); alias those slack refs into this model's container so - # `get_variable`/`has_container_key` for `V` stay consistent - # with the constraints it reuses. - src_slack_ub = nothing - src_slack_lb = nothing - if use_slacks - src_slack_ub, src_slack_lb = - _find_shared_post_contingency_slack_sources( - container, V, outage_id, name, first(time_steps), - ) + src_slack_ub, src_slack_lb = + _find_shared_post_contingency_slack_sources( + container, V, outage_id, name, first(time_steps), + ) + source_has_slacks = + !isnothing(src_slack_ub) || !isnothing(src_slack_lb) + if source_has_slacks != use_slacks + error( + "Inconsistent `use_slacks` among security-constrained " * + "branch models that share post-contingency constraints " * + "for outage $outage_id, monitored component \"$name\": " * + "the shared constraints were built " * + (source_has_slacks ? "with" : "without") * + " slacks, but the `$V` branch model sets " * + "`use_slacks = $use_slacks`. Configure `use_slacks` " * + "identically across all security-constrained branch " * + "`DeviceModel`s whose outage/monitored sets overlap.", + ) end for t in time_steps con_ub[outage_id, name, t] = diff --git a/src/network_models/instantiate_network_model.jl b/src/network_models/instantiate_network_model.jl index e046830..f5fad03 100644 --- a/src/network_models/instantiate_network_model.jl +++ b/src/network_models/instantiate_network_model.jl @@ -235,11 +235,13 @@ function _prune_fully_reduced_branch_models!( push!(pruned, branch_type) end for branch_type in pruned + hint = branch_type === PSY.MonitoredLine ? + " Use the `model_all_branches` attribute on the MonitoredLine DeviceModel to retain such lines through the reduction." : + " Consider adjusting the network-reduction settings/tolerance to avoid merging all branches of this type." @warn "All components of branch type $(branch_type) were merged away by the " * "network reduction (e.g. a zero-impedance branch merge). The " * "$(branch_type) DeviceModel is dropped from the template and will not " * - "be modeled. Use the `model_all_branches` attribute on a MonitoredLine " * - "model to retain such branches through the reduction." + "be modeled.$hint" delete!(branch_models, nameof(branch_type)) filter!(!=(branch_type), network_model.modeled_branch_types) end diff --git a/test/Project.toml b/test/Project.toml index 4c4c216..4efdb03 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -39,7 +39,7 @@ InfrastructureOptimizationModels = {rev = "ac/ordc", url = "https://github.com/S InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl"} PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} PowerFlows = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerFlows.jl"} -PowerNetworkMatrices = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"} +PowerNetworkMatrices = {rev = "ac/port-branch-admittance", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"} PowerSystemCaseBuilder = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystemCaseBuilder.jl"} [compat] diff --git a/test/test_native_dcp_acp_models.jl b/test/test_native_dcp_acp_models.jl index 3025dab..5e10f81 100644 --- a/test/test_native_dcp_acp_models.jl +++ b/test/test_native_dcp_acp_models.jl @@ -1,3 +1,5 @@ +import PowerNetworkMatrices as PNM + # build! routes log records emitted inside build_model! to the model's internal # file logger (operation_problem.log), so they are NOT visible to @test_logs at # the call site. To assert (or rule out) the REF-less pin warning we read it. @@ -152,7 +154,7 @@ end @testset "branch_admittance primitives" begin sys = PSB.build_system(PSITestSystems, "c_sys5") line = first(PSY.get_components(PSY.Line, sys)) - a = PowerOperationsModels._branch_admittance(line) + a = PNM.branch_admittance(line) r, x = PSY.get_r(line, PSY.SU), PSY.get_x(line, PSY.SU) y = inv(complex(r, x)) @test a.g ≈ real(y) @@ -164,14 +166,12 @@ end @testset "branch_flow_limits MonitoredLine" begin sys = PSB.build_system(PSITestSystems, "c_sys5_ml") ml = first(PSY.get_components(PSY.MonitoredLine, sys)) - fl = PowerOperationsModels.branch_flow_limits(ml) + fl = PNM.branch_flow_limits(ml) psy_fl = PSY.get_flow_limits(ml, PSY.SU) @test fl.from_to == psy_fl.from_to @test fl.to_from == psy_fl.to_from end -import PowerNetworkMatrices as PNM - @testset "reduced arc admittance uses PNM series equivalent, not original branch" begin # `case11_network_reductions` is purpose-built to produce series arcs under the # radial + degree-two reduction (c_sys5/c_sys14 yield no reducible chains). @@ -194,9 +194,9 @@ import PowerNetworkMatrices as PNM @test !isempty(series_map) # degree-2 reduction produces series arcs (from_no, to_no), chain = first(series_map) - resolved = PowerOperationsModels._reduced_arc_admittance(nr, from_no, to_no) + resolved = PNM.reduced_arc_admittance(nr, from_no, to_no) @test resolved !== nothing - expected = PowerOperationsModels._segment_admittance(chain, nr) + expected = PNM.branch_admittance(chain, nr) @test isapprox(resolved.b, expected.b; atol = 1e-9) # Non-triviality: the series equivalent is the MERGED admittance of the chain, so @@ -205,13 +205,13 @@ import PowerNetworkMatrices as PNM # for the reduced arc, which is wrong. members = collect(chain) @test length(members) >= 2 - member_b = PowerOperationsModels._branch_admittance(members[1]).b + member_b = PNM.branch_admittance(members[1]).b @test !isapprox(resolved.b, member_b; rtol = 1e-3) # Reversed-orientation arc exercises the `_reverse_admittance` path: series b is # symmetric, from/to shunts swap, and any phase shift negates. if !haskey(series_map, (to_no, from_no)) - reversed = PowerOperationsModels._reduced_arc_admittance(nr, to_no, from_no) + reversed = PNM.reduced_arc_admittance(nr, to_no, from_no) @test reversed !== nothing @test isapprox(reversed.b, resolved.b; atol = 1e-9) @test isapprox(reversed.b_fr, resolved.b_to; atol = 1e-9) @@ -220,7 +220,7 @@ import PowerNetworkMatrices as PNM # A direct (un-reduced) arc resolves to `nothing` — the caller falls back to the # branch's own admittance. - @test PowerOperationsModels._reduced_arc_admittance(nr, -1, -2) === nothing + @test PNM.reduced_arc_admittance(nr, -1, -2) === nothing end @testset "Transformer3W _winding_admittance star-arc decomposition" begin @@ -306,7 +306,7 @@ end PSY.add_component!(sys, transformer3w) w = PNM.ThreeWindingTransformerWinding(transformer3w, 1) - adm = PowerOperationsModels._winding_admittance(w) + adm = PNM.winding_admittance(w) r = PNM.get_equivalent_r(w) x = PNM.get_equivalent_x(w) From 437f0c7f298c500ee1f58346f3527873565bb61f Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Tue, 16 Jun 2026 18:51:56 -0400 Subject: [PATCH 09/11] formatting --- src/network_models/instantiate_network_model.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/network_models/instantiate_network_model.jl b/src/network_models/instantiate_network_model.jl index f5fad03..c6e2bda 100644 --- a/src/network_models/instantiate_network_model.jl +++ b/src/network_models/instantiate_network_model.jl @@ -235,9 +235,11 @@ function _prune_fully_reduced_branch_models!( push!(pruned, branch_type) end for branch_type in pruned - hint = branch_type === PSY.MonitoredLine ? - " Use the `model_all_branches` attribute on the MonitoredLine DeviceModel to retain such lines through the reduction." : - " Consider adjusting the network-reduction settings/tolerance to avoid merging all branches of this type." + hint = if branch_type === PSY.MonitoredLine + " Use the `model_all_branches` attribute on the MonitoredLine DeviceModel to retain such lines through the reduction." + else + " Consider adjusting the network-reduction settings/tolerance to avoid merging all branches of this type." + end @warn "All components of branch type $(branch_type) were merged away by the " * "network reduction (e.g. a zero-impedance branch merge). The " * "$(branch_type) DeviceModel is dropped from the template and will not " * From 5319afce2949af6705a4fce7ce6702bc845a4db8 Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 17 Jun 2026 10:07:18 -0400 Subject: [PATCH 10/11] remove some comments --- src/ac_transmission_models/AC_branches.jl | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/src/ac_transmission_models/AC_branches.jl b/src/ac_transmission_models/AC_branches.jl index 11f6d0d..2434329 100644 --- a/src/ac_transmission_models/AC_branches.jl +++ b/src/ac_transmission_models/AC_branches.jl @@ -1427,7 +1427,7 @@ end # Admittance for `branch`'s ohm's law given its retained endpoints: the branch's own # π-parameters (`PNM.branch_admittance`) for a direct/un-reduced arc, or PNM's # reduction-aware equivalent (`PNM.reduced_arc_admittance`) for a series/parallel-aggregated -# arc. The π-admittance primitives themselves now live in PNM (`BranchAdmittance.jl`). +# arc. function _resolve_branch_admittance(network_model, branch, from_no::Int, to_no::Int) nr = get_network_reduction(network_model) isempty(nr) && return PNM.branch_admittance(branch) @@ -1469,10 +1469,6 @@ function _branch_geometries(sys::PSY.System, network_model, devices) return geoms end -# `branch_flow_limits(branch)` (directional per-unit MVA flow limits, including the -# `BranchesParallel`/`BranchesSeries`/`ThreeWindingTransformerWinding` reduction wrappers) -# now lives in PNM (`BranchAdmittance.jl`); call `PNM.branch_flow_limits`. - ################################## Native ACP apparent-power rate constraints ############### """ @@ -1930,12 +1926,6 @@ end # storage 2D (name × time) without inventing a new container shape. ################################################################################ -# ── Native three-winding helpers ───────────────────────────────────────────── -# The per-winding decomposition (`PNM.three_winding_arcs`, returning a `suffix`/`arc`/ -# `rating`/`winding` NamedTuple per wye-model winding) and the per-winding π-admittance -# (`PNM.winding_admittance`) now live in PNM (`BranchAdmittance.jl`); the native 3W builders -# below call them directly. - "Build the list of per-winding variable names for a set of Transformer3W devices." function _three_winding_var_names(devices) names = String[] From c94b1708cac0e2c6818c36dd748c7732a05497dc Mon Sep 17 00:00:00 2001 From: Anthony Costarelli Date: Wed, 17 Jun 2026 10:35:09 -0400 Subject: [PATCH 11/11] more comment deletion --- src/common_models/market_bid_overrides.jl | 16 +--------------- 1 file changed, 1 insertion(+), 15 deletions(-) diff --git a/src/common_models/market_bid_overrides.jl b/src/common_models/market_bid_overrides.jl index acafe7c..2d713cf 100644 --- a/src/common_models/market_bid_overrides.jl +++ b/src/common_models/market_bid_overrides.jl @@ -294,21 +294,7 @@ _vom_offer_direction(::Type{<:AbstractControllablePowerLoadFormulation}) = """ Add the delta PWL objective terms for an ORDC service (`StepwiseCostReserve` over -`ReserveDemandCurve` / `ReserveDemandTimeSeriesCurve`). The demand curve is a -willingness-to-pay (concave) decremental offer, so the whole path routes through -IOM's generic `OfferDirection` machinery via `_reserve_offer_direction`: - -* `IOM._get_pwl_data` fetches the curve (`get_offer_curves` → `PSY.get_variable`) - and dispatches the static vs. time-series read internally. -* `IOM._block_offer_var` / `IOM._block_offer_constraint` pick the block-offer - variable/constraint family. -* `IOM.add_pwl_constraint_delta!` builds the block-sum constraint. - -All containers are keyed by `meta = name` because services are built one at a -time (one container per service, not a shared names-axis batch). - -`cost_data` is accepted to match the caller's signature but is re-resolved inside -`IOM._get_pwl_data`; it is the same `PSY.get_variable(component)` curve. +`ReserveDemandCurve` / `ReserveDemandTimeSeriesCurve`). """ function add_pwl_term_delta!( container::OptimizationContainer,