diff --git a/Project.toml b/Project.toml index fa73c03..87bb412 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/g1-port", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} PowerNetworkMatrices = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"} [extensions] diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index 41bdde2..4073d75 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -307,6 +307,7 @@ include("network_models/network_constructor.jl") # Services Models include("services_models/service_slacks.jl") include("services_models/reserves.jl") +include("services_models/static_injection_security_constrained_models.jl") include("services_models/reserve_group.jl") # include("services_models/agc.jl") # TODO: needs _get_ace_error include("services_models/transmission_interface.jl") @@ -536,6 +537,8 @@ export FlowActivePowerSlackUpperBound export FlowActivePowerSlackLowerBound export PostContingencyFlowActivePowerSlackUpperBound export PostContingencyFlowActivePowerSlackLowerBound +export PostContingencyActivePowerChangeVariable +export PostContingencyActivePowerReserveDeploymentVariable export FlowActivePowerFromToVariable export FlowActivePowerToFromVariable export FlowReactivePowerFromToVariable @@ -724,6 +727,11 @@ export FlowRateConstraint export FlowRateConstraintFromTo export FlowRateConstraintToFrom export PostContingencyFlowRateConstraint +export PostContingencyGenerationBalanceConstraint +export PostContingencyActivePowerGenerationLimitsConstraint +export PostContingencyCopperPlateBalanceConstraint +export PostContingencyActivePowerVariableLimitsConstraint +export PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint export FlowLimitConstraint export FlowLimitFromToConstraint export FlowLimitToFromConstraint @@ -777,6 +785,11 @@ export FuelConsumptionExpression export ActivePowerRangeExpressionLB export ActivePowerRangeExpressionUB export PostContingencyBranchFlow +export PostContingencyAreaInterchangeFlow +export PostContingencyActivePowerGeneration +export PostContingencyActivePowerBalance +export PostContingencyNodalActivePowerDeployment +export PostContingencyAreaActivePowerDeployment export NetActivePower export DCCurrentBalance export ComponentReserveUpBalanceExpression @@ -874,6 +887,8 @@ export RangeReserve export StepwiseCostReserve export RampReserve export NonSpinningReserve +export SecurityConstrainedContingencyReserve +export SecurityConstrainedRampReserve export ConstantMaxInterfaceFlow export VariableMaxInterfaceFlow diff --git a/src/ac_transmission_models/security_constrained_branch.jl b/src/ac_transmission_models/security_constrained_branch.jl index cc1a0ca..077eabf 100644 --- a/src/ac_transmission_models/security_constrained_branch.jl +++ b/src/ac_transmission_models/security_constrained_branch.jl @@ -205,9 +205,9 @@ function _find_shared_post_contingency_constraint_sources( src_ub = nothing for (key, cc) in get_constraints(container) _is_shared_post_contingency_source(key, cc, target, T, V) || continue - if key.meta == "lb" + if key.meta == POST_CONTINGENCY_LB_META src_lb = cc - elseif key.meta == "ub" + elseif key.meta == POST_CONTINGENCY_UB_META src_ub = cc end !isnothing(src_lb) && !isnothing(src_ub) && break @@ -266,50 +266,6 @@ function _has_other_v_container( return false end -""" -Pre-allocate a `SparseAxisArray` keyed by -`(outage_id::String, monitored_name::String, t::Int)` holding `JuMP.AffExpr` -zeros for every entry produced by `_resolve_monitored_arcs`. The pre-fill is -required so the parallel PTDF expression build below cannot race on Dict resize. -""" -function _add_post_contingency_sparse_expression!( - container::OptimizationContainer, - ::Type{T}, - ::Type{V}, - resolved::Vector{ - Pair{Base.UUID, Vector{Tuple{DataType, String, Tuple{Int, Int}, String}}}, - }, - time_steps::UnitRange{Int}, -) where {T <: PostContingencyExpressions, V <: PSY.ACTransmission} - contents = Dict{Tuple{String, String, Int}, JuMP.AffExpr}() - for (uuid, entries) in resolved - outage_id = string(uuid) - for (_, name, _, _) in entries, t in time_steps - contents[(outage_id, name, t)] = zero(JuMP.AffExpr) - end - end - expr_container = SparseAxisArray(contents) - IOM._assign_container!(container.expressions, ExpressionKey(T, V), expr_container) - return expr_container -end - -""" -Register an empty `SparseAxisArray` keyed by -`(outage_id::String, monitored_name::String, t::Int)` for the given constraint -type and meta tag. -""" -function _add_post_contingency_sparse_constraints!( - container::OptimizationContainer, - ::Type{T}, - ::Type{V}; - meta::String, -) where {T <: ConstraintType, V <: PSY.ACTransmission} - cons_container = - SparseAxisArray(Dict{Tuple{String, String, Int}, JuMP.ConstraintRef}()) - IOM._assign_container!(container.constraints, ConstraintKey(T, V, meta), cons_container) - return cons_container -end - """ For each outage in `device_model.outages`, resolve every monitored component (across every monitored type) to its arc in the active reduction graph — using @@ -412,15 +368,36 @@ function add_constraints!( resolved = _resolve_monitored_arcs(device_model, net_reduction_data) - con_lb = _add_post_contingency_sparse_constraints!(container, T, V; meta = "lb") - con_ub = _add_post_contingency_sparse_constraints!(container, T, V; meta = "ub") + # Sparse `(outage_id, name, t)` constraint containers; the empty axes only fix + # the key tuple type and the build loop fills the ragged entries by assignment. + con_lb = IOM.add_constraints_container!( + container, T, V, String[], String[], time_steps; + sparse = true, meta = POST_CONTINGENCY_LB_META, + ) + con_ub = IOM.add_constraints_container!( + container, T, V, String[], String[], time_steps; + sparse = true, meta = POST_CONTINGENCY_UB_META, + ) use_slacks = get_use_slacks(device_model) - # Local relaxation-slack containers keyed by `(outage_id, name, t)`. Built - # here (not via `add_variables!`) because the post-contingency axes are only - # known after `_resolve_monitored_arcs`; registered after the loop iff used. - slack_ub = SparseAxisArray(Dict{Tuple{String, String, Int}, JuMP.VariableRef}()) - slack_lb = SparseAxisArray(Dict{Tuple{String, String, Int}, JuMP.VariableRef}()) + # Relaxation-slack containers keyed by `(outage_id, name, t)`, registered only + # when slacks are enabled. Built here (not via `add_variables!`) because the + # post-contingency keys are ragged; the axes only fix the key tuple type and + # the build loop assigns the resolved entries. + slack_ub, slack_lb = if use_slacks + ( + IOM.add_variable_container!( + container, PostContingencyFlowActivePowerSlackUpperBound, V, + String[], String[], time_steps; sparse = true, + ), + IOM.add_variable_container!( + container, PostContingencyFlowActivePowerSlackLowerBound, V, + String[], String[], time_steps; sparse = true, + ), + ) + else + (nothing, nothing) + end expressions = get_expression(container, PostContingencyBranchFlow, V) jump_model = get_jump_model(container) @@ -542,20 +519,6 @@ function add_constraints!( end end - if !isempty(slack_ub.data) - IOM._assign_container!( - container.variables, - VariableKey(PostContingencyFlowActivePowerSlackUpperBound, V), - slack_ub, - ) - end - if !isempty(slack_lb.data) - IOM._assign_container!( - container.variables, - VariableKey(PostContingencyFlowActivePowerSlackLowerBound, V), - slack_lb, - ) - end return end @@ -599,8 +562,10 @@ function add_post_contingency_flow_expressions!( net_reduction_data = network_model.network_reduction resolved = _resolve_monitored_arcs(model, net_reduction_data) - expression_container = _add_post_contingency_sparse_expression!( - container, T, V, resolved, time_steps, + # Empty sparse `(outage_id, name, t)` expression container; the axes only fix + # the key tuple type and the build loop fills the ragged entries by assignment. + expression_container = IOM.add_expression_container!( + container, T, V, String[], String[], time_steps; sparse = true, ) fresh_resolved = _copy_existing_post_contingency_expressions!( @@ -713,10 +678,8 @@ function add_post_contingency_flow_expressions!( time_steps = get_time_steps(container) resolved = _resolve_monitored_arcs(model, network_model.network_reduction) - expression_container = - SparseAxisArray(Dict{Tuple{String, String, Int}, JuMP.AffExpr}()) - IOM._assign_container!( - container.expressions, ExpressionKey(T, V), expression_container, + expression_container = IOM.add_expression_container!( + container, T, V, String[], String[], time_steps; sparse = true, ) has_other_v = _has_other_v_container(IOM.get_expressions(container), T, V) diff --git a/src/core/constraints.jl b/src/core/constraints.jl index 7fdab97..472b218 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -244,6 +244,46 @@ is the branch emergency rating (system base / per-unit). """ struct PostContingencyFlowRateConstraint <: PostContingencyConstraintType end +""" +Constraint that closes the per-contingency post-contingency generation balance: +the total active power shed by the outaged generators equals the total reserve +deployment across all contributing generators. + +```math +\\sum_{g \\in \\mathcal{G}_c} p_{g,t} = \\sum_{g \\in \\mathcal{G}} \\Delta rsv_{c,g,t}, +\\quad \\forall c \\in \\mathcal{C},\\ \\forall t +``` +""" +struct PostContingencyGenerationBalanceConstraint <: PostContingencyConstraintType end + +""" +Constraint on the post-contingency active power generation expression +(`PostContingencyActivePowerGeneration`) of each contributing generator +under each outage. Outaged generators are pinned to zero; all other +generators are bounded by their `active_power_limits`. + +```math +P^{\\text{min}}_g \\le p_{g,t} + \\Delta rsv_{c,g,t} \\le P^{\\text{max}}_g, +\\quad \\forall c \\in \\mathcal{C},\\ \\forall g \\notin \\mathcal{G}_c,\\ \\forall t +``` +""" +struct PostContingencyActivePowerGenerationLimitsConstraint <: + PostContingencyConstraintType end + +""" +Constraint that closes the per-area post-contingency power balance for the +`AreaBalancePowerModel` network representation, summing the +`PostContingencyAreaActivePowerDeployment` with the pre-contingency +`ActivePowerBalance` expression for each area. + +```math +\\sum_{g \\in \\mathcal{A}}(\\Delta rsv_{c,g,t} + p_{g,t}\\mathbb{1}_{g \\in \\mathcal{G}_c}) + +\\text{Bal}^{\\text{pre}}_{a,t} = 0,\\quad \\forall a \\in \\mathcal{A},\\ \\forall c,\\ \\forall t +``` +""" +struct PostContingencyCopperPlateBalanceConstraint <: + PostContingencyConstraintType end + """ Struct to create the constraint for branch flow rate limits from the 'from' bus to the 'to' bus. For more information check [Branch Formulations](@ref PowerSystems.Branch-Formulations). @@ -454,6 +494,31 @@ S_k^{\\max}`` for ``k \\in \\{f, t\\}`` via one of two shapes, selected by the struct HVDCVSCApparentPowerLimitConstraint <: ConstraintType end abstract type PowerVariableLimitsConstraint <: ConstraintType end + +abstract type PostContingencyVariableLimitsConstraint <: PowerVariableLimitsConstraint end + +""" +Limits each generator's post-contingency active power to its range. + +```math +P^\\text{min} \\le p_t + \\Delta p_{c, t} \\le P^\\text{max}, +\\quad \\forall c \\in \\mathcal{C},\\ \\forall t +``` +""" +struct PostContingencyActivePowerVariableLimitsConstraint <: + PostContingencyVariableLimitsConstraint end + +""" +Caps reserve deployed under a contingency by the reserve procured pre-contingency. + +```math +\\Delta rsv_{r, c, t} \\le rsv_{r, c, t}, +\\quad \\forall r \\in \\mathcal{R},\\ \\forall c \\in \\mathcal{C},\\ \\forall t +``` +""" +struct PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint <: + PostContingencyVariableLimitsConstraint end + """ Struct to create the constraint to limit active power input expressions. For more information check [Device Formulations](@ref formulation_intro). diff --git a/src/core/definitions.jl b/src/core/definitions.jl index 54d8f1e..df61b22 100644 --- a/src/core/definitions.jl +++ b/src/core/definitions.jl @@ -62,6 +62,11 @@ const UNSET_INI_TIME = Dates.DateTime(0) const ABSOLUTE_TOLERANCE = 1.0e-3 const BALANCE_SLACK_COST = 1e6 const CONSTRAINT_VIOLATION_SLACK_COST = 2e5 +const POST_CONTINGENCY_CONSTRAINT_VIOLATION_SLACK_COST = 1e5 +# Meta tags distinguishing the lower/upper post-contingency flow-rate constraint +# containers; shared so registration and the reuse-lookup never drift. +const POST_CONTINGENCY_LB_META = "lb" +const POST_CONTINGENCY_UB_META = "ub" const SERVICES_SLACK_COST = 1e5 const COST_EPSILON = 1e-3 const PTDF_ZERO_TOL = 1e-9 diff --git a/src/core/expressions.jl b/src/core/expressions.jl index a12d112..3739d9d 100644 --- a/src/core/expressions.jl +++ b/src/core/expressions.jl @@ -6,7 +6,11 @@ struct ComponentReserveUpBalanceExpression <: ExpressionType end struct ComponentReserveDownBalanceExpression <: ExpressionType end struct InterfaceTotalFlow <: ExpressionType end struct PTDFBranchFlow <: ExpressionType end +abstract type PostContingencySystemBalanceExpressions <: SystemBalanceExpressions end +struct PostContingencyActivePowerBalance <: PostContingencySystemBalanceExpressions end +struct PostContingencyAreaInterchangeFlow <: PostContingencyExpressions end struct PostContingencyNodalActivePowerDeployment <: PostContingencyExpressions end +struct PostContingencyAreaActivePowerDeployment <: PostContingencyExpressions end struct RealizedShiftedLoad <: ExpressionType end ################################################################################# @@ -118,6 +122,7 @@ struct StorageReserveBalanceExpression{D, S, Sd} <: # Method extensions for output writing should_write_resulting_value(::Type{InterfaceTotalFlow}) = true should_write_resulting_value(::Type{PTDFBranchFlow}) = true +should_write_resulting_value(::Type{PostContingencyAreaInterchangeFlow}) = true should_write_resulting_value(::Type{RealizedShiftedLoad}) = true should_write_resulting_value(::Type{HydroServedReserveUpExpression}) = true @@ -133,5 +138,7 @@ should_write_resulting_value( # Method extensions for unit conversion convert_output_to_natural_units(::Type{InterfaceTotalFlow}) = true convert_output_to_natural_units(::Type{PostContingencyBranchFlow}) = true +convert_output_to_natural_units(::Type{PostContingencyAreaInterchangeFlow}) = true +convert_output_to_natural_units(::Type{PostContingencyActivePowerGeneration}) = true convert_output_to_natural_units(::Type{PTDFBranchFlow}) = true convert_output_to_natural_units(::Type{RealizedShiftedLoad}) = true diff --git a/src/core/formulations.jl b/src/core/formulations.jl index 76a001f..77a14a6 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -332,6 +332,36 @@ struct RampReserve <: AbstractReservesFormulation end Struct to add non spinning reserve requirements larger than specified requirement """ struct NonSpinningReserve <: AbstractReservesFormulation end + +abstract type AbstractSecurityConstrainedReservesFormulation <: AbstractReservesFormulation end + +""" +Security-constrained contingency reserve formulation: deploys reserves under +each G-1 outage scoped to the reserve `PSY.Service`. The set of contingencies a +service responds to is the `PSY.Outage` supplemental attributes attached to that +service via `add_supplemental_attribute!(sys, service, outage)`; template +validation mirrors those attachments into `service_model.outages`. +A `RequirementTimeSeriesParameter` is optional: when present the requirement / +ramp / participation stack is built; when absent, per-generator +`PostContingencyActivePowerGeneration` limits are applied instead. +Post-contingency branch-flow constraints are added only for the monitored +components listed on each outage's `monitored_components`. + +See also `SecurityConstrainedRampReserve`. +""" +struct SecurityConstrainedContingencyReserve <: + AbstractSecurityConstrainedReservesFormulation end + +""" +Security-constrained ramp reserve formulation: like `RampReserve` for the +pre-contingency requirement/ramp/participation constraints, plus the same +G-1 post-contingency deployment + monitored-branch flow constraints as +`SecurityConstrainedContingencyReserve`. + +See also `SecurityConstrainedContingencyReserve`. +""" +struct SecurityConstrainedRampReserve <: + AbstractSecurityConstrainedReservesFormulation end """ Struct to add a constant maximum transmission flow for specified interface """ diff --git a/src/core/variables.jl b/src/core/variables.jl index e3f8e1d..9aa137c 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -9,6 +9,7 @@ abstract type MultiStartVariable <: VariableType end abstract type AbstractACActivePowerFlow <: VariableType end abstract type AbstractACReactivePowerFlow <: VariableType end +abstract type AbstractContingencyVariableType <: VariableType end # AbstractPiecewiseLinearBlockOffer: moved into IOM """ @@ -205,18 +206,44 @@ Docs abbreviation: ``f^\\text{sl,lo}`` struct FlowActivePowerSlackLowerBound <: AbstractACActivePowerFlow end """ -Struct to dispatch the creation of post-contingency active power flow upper bound slack variables. Relaxes the post-contingency (N-1) emergency-rate upper-bound constraint when `use_slacks = true`. +Post-contingency active power change of a generator under a contingency. + +Docs abbreviation: ``\\Delta p_{g,c}`` +""" +struct PostContingencyActivePowerChangeVariable <: AbstractContingencyVariableType end + +""" +Reserve deployed by a generator under a contingency. + +Docs abbreviation: ``\\Delta rsv_{r,g,c}`` +""" +struct PostContingencyActivePowerReserveDeploymentVariable <: + AbstractContingencyVariableType end + +""" +Abstract supertype for non-negative slack variables that absorb infeasibility +in post-contingency branch-flow inequalities of the security-constrained +reserve formulations. +""" +abstract type AbstractContingencySlackVariableType <: VariableType end + +""" +Non-negative slack relaxing the post-contingency (N-1) emergency-rate upper-bound +constraint when `use_slacks = true`. Docs abbreviation: ``f^\\text{sl,up,N-1}`` """ -struct PostContingencyFlowActivePowerSlackUpperBound <: VariableType end +struct PostContingencyFlowActivePowerSlackUpperBound <: + AbstractContingencySlackVariableType end """ -Struct to dispatch the creation of post-contingency active power flow lower bound slack variables. Relaxes the post-contingency (N-1) emergency-rate lower-bound constraint when `use_slacks = true`. +Non-negative slack relaxing the post-contingency (N-1) emergency-rate lower-bound +constraint when `use_slacks = true`. Docs abbreviation: ``f^\\text{sl,lo,N-1}`` """ -struct PostContingencyFlowActivePowerSlackLowerBound <: VariableType end +struct PostContingencyFlowActivePowerSlackLowerBound <: + AbstractContingencySlackVariableType end """ Struct to dispatch the creation of Phase Shifters Variables @@ -679,6 +706,11 @@ convert_output_to_natural_units(::Type{ActivePowerOutVariable}) = true convert_output_to_natural_units(::Type{EnergyVariable}) = true convert_output_to_natural_units(::Type{ReactivePowerVariable}) = true convert_output_to_natural_units(::Type{ActivePowerReserveVariable}) = true +convert_output_to_natural_units(::Type{PostContingencyActivePowerChangeVariable}) = true +convert_output_to_natural_units( + ::Type{PostContingencyActivePowerReserveDeploymentVariable}, +) = + true convert_output_to_natural_units(::Type{ServiceRequirementVariable}) = true convert_output_to_natural_units(::Type{RateofChangeConstraintSlackUp}) = true convert_output_to_natural_units(::Type{RateofChangeConstraintSlackDown}) = true diff --git a/src/operation/template_validation.jl b/src/operation/template_validation.jl index 5b4ec31..af4cd4a 100644 --- a/src/operation/template_validation.jl +++ b/src/operation/template_validation.jl @@ -88,6 +88,7 @@ function validate_template_impl!(model::IOM.AbstractOptimizationModel) _check_branch_rating_time_series_formulation!(template.branches, system) validate_network_model(network_model, unmodeled_branch_types, model_has_branch_filters) _build_device_model_outages!(template, system) + _build_service_model_outages!(template, system) return end @@ -294,6 +295,112 @@ function _sc_branch_models(template::IOM.AbstractProblemTemplate) ] end +""" +Populate `service_model.outages` for every security-constrained (SC) reserve +`ServiceModel` in the template, in a single pass over each SC service model. + +The user opts a service into responding to a given contingency by attaching +the outage supplemental attribute to the `PSY.Service` instance directly +(`add_supplemental_attribute!(sys, service, outage)`). That attachment is the +sole selection mechanism: a service claims exactly the outages attached to it, +regardless of whether the outaged component is among the service's contributing +devices. + +`PlannedOutage`s attached to the service are still gated by the +`"include_planned_outages"` attribute on the SC `ServiceModel` (default +`false`); other `Outage` subtypes are always claimed. +""" +function _build_service_model_outages!( + template::IOM.AbstractProblemTemplate, + sys::PSY.System, +) + sc_service_models = _sc_reserve_service_models(template) + template_sc_service_keys = Set{Tuple{DataType, String}}( + (get_component_type(m), get_service_name(m)) for m in sc_service_models + ) + # Unconditional (before the early return below) so orphan attachments are + # still reported when no SC service models are registered. + _warn_outages_attached_to_unmodeled_services(sys, template_sc_service_keys) + + isempty(sc_service_models) && return + + modeled_types = Set{DataType}(get_component_types(template)) + uncovered_types = Dict{DataType, Set{Base.UUID}}() + + for m in sc_service_models + empty!(get_outages(m)) + D = get_component_type(m) + service_name = get_service_name(m) + service = PSY.get_component(D, sys, service_name) + if service === nothing + @warn "ServiceModel{$D, $(get_formulation(m))} (service_name=\ + $(service_name)) is in the template but no matching service \ + exists in the system; it will not contribute any \ + post-contingency constraints." _group = + IOM.LOG_GROUP_MODELS_VALIDATION + continue + end + + for outage in PSY.get_supplemental_attributes(PSY.Outage, service) + outage_uuid = IS.get_uuid(outage) + if isempty(PSY.get_monitored_components(outage)) + @warn "Outage $(outage_uuid) ($(typeof(outage))) attached to \ + service $(service_name) has empty monitored_components; \ + no post-contingency variables or constraints will be \ + created for this outage." _group = + IOM.LOG_GROUP_MODELS_VALIDATION + continue + end + _service_skips_outage(outage, m) && continue + + per_type, uncovered = _monitored_components_by_modeled_type( + outage, outage_uuid, sys, modeled_types, + ) + for comp_type in uncovered + push!(get!(uncovered_types, comp_type, Set{Base.UUID}()), outage_uuid) + end + isempty(per_type) && continue + get_outages(m)[outage_uuid] = per_type + end + end + + _warn_uncovered_monitored_types(uncovered_types) + return +end + +# SC reserve service models in the template. +function _sc_reserve_service_models(template::IOM.AbstractProblemTemplate) + return ServiceModel[ + m for m in values(get_service_models(template)) if + get_formulation(m) <: AbstractSecurityConstrainedReservesFormulation + ] +end + +_service_skips_outage(::PSY.Outage, ::ServiceModel) = false +_service_skips_outage(::PSY.PlannedOutage, m::ServiceModel) = + get_attribute(m, "include_planned_outages") !== true + +function _warn_outages_attached_to_unmodeled_services( + sys::PSY.System, + template_sc_service_keys::Set{Tuple{DataType, String}}, +) + for service in PSY.get_components(PSY.Service, sys) + attached_outages = PSY.get_supplemental_attributes(PSY.Outage, service) + isempty(attached_outages) && continue + key = (typeof(service), PSY.get_name(service)) + key in template_sc_service_keys && continue + for outage in attached_outages + @warn "Outage $(IS.get_uuid(outage)) is attached to service \ + $(PSY.get_name(service)) ($(typeof(service))) but the \ + template does not include a security-constrained \ + ServiceModel for it; the outage will not contribute any \ + post-contingency reserve constraints." _group = + IOM.LOG_GROUP_MODELS_VALIDATION + end + end + return +end + # Per SC-model component type, the user's explicit outage-UUID allow-list from # the constructor kwarg: a non-empty set restricts auto-discovery to those # UUIDs; an empty set means auto-discover all. Clears `m.outages` so the main @@ -327,7 +434,11 @@ function _monitored_components_by_modeled_type( continue end comp_type = typeof(component) - if comp_type <: PSY.ACTransmission && comp_type in modeled_types + # `PSY.AreaInterchange` is admitted alongside `PSY.ACTransmission` so the + # AreaBalance service-side builder can pick it up; anything else is routed + # to the skip path (it would `KeyError` in the arc-resolution maps). + admissible = (comp_type <: PSY.ACTransmission || comp_type <: PSY.AreaInterchange) + if admissible && comp_type in modeled_types push!(get!(per_type, comp_type, Set{String}()), PSY.get_name(component)) else push!(uncovered, comp_type) diff --git a/src/services_models/reserves.jl b/src/services_models/reserves.jl index 1340b15..9f47590 100644 --- a/src/services_models/reserves.jl +++ b/src/services_models/reserves.jl @@ -68,6 +68,18 @@ function get_default_time_series_names( ) end +# The returned name (`"requirement"`) is the exact `PSY` time-series name the +# requirement must be stored under on the service; the requirement is optional +# for security-constrained formulations (see `SecurityConstrainedContingencyReserve`). +function get_default_time_series_names( + ::Type{<:PSY.Reserve}, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) + return Dict{Type{<:TimeSeriesParameter}, String}( + RequirementTimeSeriesParameter => "requirement", + ) +end + function get_default_time_series_names( ::Type{<:PSY.ReserveNonSpinning}, ::Type{NonSpinningReserve}, diff --git a/src/services_models/static_injection_security_constrained_models.jl b/src/services_models/static_injection_security_constrained_models.jl new file mode 100644 index 0000000..6d3dc01 --- /dev/null +++ b/src/services_models/static_injection_security_constrained_models.jl @@ -0,0 +1,1642 @@ +# ---------------------------------------------------------------------------- +# Security-constrained reserve service formulations (G-1 with reserve +# deployment + monitored-branch post-contingency flow constraints). +# +# Sparse + monitored counterpart to the legacy dense-per-branch implementation +# in `older_static_injection_security_constrained_models.jl`. Mirrors the +# device-side ac_transmission_security_constrained_models.jl: post-contingency +# flow expressions/constraints (and optional slacks) live in +# `SparseAxisArray`s keyed by `(outage_id::String, monitored_name::String, +# t::Int)`, scoped to the monitored components carried by each outage in +# `service_model.outages[uuid]::Dict{DataType, Set{String}}` (populated by +# `_build_service_model_outages!`). +# +# Per-outage reserve-deployment variables and the per-outage power-balance / +# nodal-deployment / area-deployment expressions remain dense over the +# contributing devices and the modeled bus/area axes — those are independent +# of which branches are monitored. +# ---------------------------------------------------------------------------- + +#! format: off +get_variable_upper_bound( + ::Type{PostContingencyFlowActivePowerSlackUpperBound}, + ::PSY.ACTransmission, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) = nothing +get_variable_lower_bound( + ::Type{PostContingencyFlowActivePowerSlackUpperBound}, + ::PSY.ACTransmission, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) = 0.0 +get_variable_upper_bound( + ::Type{PostContingencyFlowActivePowerSlackLowerBound}, + ::PSY.ACTransmission, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) = nothing +get_variable_lower_bound( + ::Type{PostContingencyFlowActivePowerSlackLowerBound}, + ::PSY.ACTransmission, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) = 0.0 + +# Reserve-deployment variable bounds/binary and signed multipliers for the +# post-contingency expressions. Type-form (POM convention); reads are +# system-base (`PSY.SU`). +get_variable_binary( + ::Type{PostContingencyActivePowerReserveDeploymentVariable}, + ::Type{<:PSY.Reserve}, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) = false +get_variable_upper_bound( + ::Type{PostContingencyActivePowerReserveDeploymentVariable}, + r::PSY.Reserve, + d::PSY.Device, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) = PSY.get_max_active_power(d, PSY.SU) +get_variable_lower_bound( + ::Type{PostContingencyActivePowerReserveDeploymentVariable}, + ::PSY.Reserve, + ::PSY.Device, + ::Type, +) = 0.0 +get_variable_multiplier( + ::Type{<:AbstractContingencyVariableType}, + ::Type{<:PSY.Reserve{PSY.ReserveDown}}, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) = -1.0 +get_variable_multiplier( + ::Type{<:AbstractContingencyVariableType}, + ::Type{<:PSY.Reserve{PSY.ReserveUp}}, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) = 1.0 +get_variable_multiplier( + ::Type{<:VariableType}, + ::Type{<:PSY.Generator}, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) = -1.0 +#! format: on + +# `true` when the formulation always requires a reserve-requirement time series +# (the pre-contingency requirement/ramp/participation stack). Contingency +# reserves make it optional; ramp reserves always require it. +requires_requirement_ts(::Type{SecurityConstrainedContingencyReserve}) = false +requires_requirement_ts(::Type{SecurityConstrainedRampReserve}) = true + +# `true` when `service` carries the specific reserve-requirement time series +# mapped by `RequirementTimeSeriesParameter` — not merely *some* time series. +function _has_requirement_ts( + container::OptimizationContainer, + model::ServiceModel, + service::PSY.AbstractReserve, +) + ts_names = get_time_series_names(model) + haskey(ts_names, RequirementTimeSeriesParameter) || return false + return PSY.has_time_series( + service, + get_default_time_series_type(container), + ts_names[RequirementTimeSeriesParameter], + ) +end + +# ---------------------------------------------------------------------------- +# Helpers: monitored-arc resolution + sparse container scaffolding +# ---------------------------------------------------------------------------- + +""" +Resolve every monitored component in `service_model.outages` to a container +name and arc tuple in the active network reduction. Mirrors the device-side +`_resolve_monitored_arcs` but operates on a `ServiceModel`. Outages whose +monitored component types are not modeled in the network are skipped. + +Returns +`Vector{Pair{UUID, Vector{Tuple{DataType, String, Tuple{Int,Int}, String}}}}` +where each inner tuple is `(monitored_type, container_name, arc, +reduction_kind)`. Outages are sorted by UUID for deterministic axes. +""" +function _resolve_service_monitored_arcs( + service_model::ServiceModel, + net_reduction_data::PNM.NetworkReductionData, +) + name_to_arc_maps = PNM.get_name_to_arc_maps(net_reduction_data) + component_to_reduction_maps = + PNM.get_component_to_reduction_name_map(net_reduction_data) + resolved = + Pair{Base.UUID, Vector{Tuple{DataType, String, Tuple{Int, Int}, String}}}[] + for (uuid, per_type) in get_outages(service_model) + kept = Tuple{DataType, String, Tuple{Int, Int}, String}[] + for (T, names) in per_type + haskey(name_to_arc_maps, T) || continue + name_to_arc = name_to_arc_maps[T] + component_to_reduction = + get(component_to_reduction_maps, T, Dict{String, String}()) + seen = Set{Tuple{Int, Int}}() + for name in sort!(collect(names)) + if haskey(name_to_arc, name) + container_name = name + elseif haskey(component_to_reduction, name) + container_name = component_to_reduction[name] + else + error( + "Monitored component \"$name\" (type $T) for outage $uuid is " * + "absent from both the network-reduction name-to-arc map and " * + "the component-to-reduction map. Verify the component exists " * + "in the system and is modeled with a branch formulation that " * + "produces a PTDFBranchFlow expression.", + ) + end + arc, reduction_kind = name_to_arc[container_name] + arc in seen && continue + push!(seen, arc) + push!(kept, (T, container_name, arc, reduction_kind)) + end + end + isempty(kept) && continue + push!(resolved, uuid => kept) + end + sort!(resolved; by = first) + return resolved +end + +""" +Resolve the outages claimed by `service_model.outages` to the `PSY.Outage` +supplemental attribute objects attached to the reserve service in `sys`. +Returned vector is sorted by UUID for deterministic axes. + +The element type is `PSY.Outage` (rather than `PSY.UnplannedOutage`) so that +the `"include_planned_outages"` opt-in honored by +`_build_service_model_outages!` in `template_validation.jl` can flow through +without a `MethodError` when a `PSY.PlannedOutage` is claimed. + +This is the service-side counterpart to iterating `get_outages(device_model)` +on the AC-branch side: outages are attached to the reserve service (and, +typically, also to the outaged generator), and resolution requires a UUID +lookup against the system. Callers use the resolved objects to query +`PSY.get_associated_components(sys, outage; component_type = PSY.Generator)` +and pin the outaged generator's deployment variable to zero. +""" +function _service_outages(sys::PSY.System, service_model::ServiceModel) + outage_uuids = sort!(collect(keys(get_outages(service_model)))) + return PSY.Outage[ + PSY.get_supplemental_attribute(sys, uuid) for uuid in outage_uuids + ] +end + +""" +Sparse slack variable container keyed by +`(outage_id::String, monitored_name::String, t::Int)`, penalized in the objective +at `POST_CONTINGENCY_CONSTRAINT_VIOLATION_SLACK_COST`. +""" +function add_post_contingency_slack_variables!( + container::OptimizationContainer, + ::Type{T}, + service::R, + service_name::String, + resolved::Vector{ + Pair{Base.UUID, Vector{Tuple{DataType, String, Tuple{Int, Int}, String}}}, + }, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) where {T <: AbstractContingencySlackVariableType, R <: PSY.AbstractReserve} + time_steps = get_time_steps(container) + jump_model = get_jump_model(container) + slack_container = IOM.add_variable_container!( + container, T, R, String[], String[], time_steps; + sparse = true, meta = service_name, + ) + for (uuid, entries) in resolved + outage_id = string(uuid) + for (_, name, _, _) in entries + for t in time_steps + v = JuMP.@variable( + jump_model, + base_name = "$(T)_$(R)_$(service_name)_{$(outage_id), $(name), $(t)}", + lower_bound = 0.0, + start = 0.0, + ) + slack_container[outage_id, name, t] = v + add_to_objective_invariant_expression!( + container, + v * POST_CONTINGENCY_CONSTRAINT_VIOLATION_SLACK_COST, + ) + end + end + end + return slack_container +end + +# Build the upper/lower post-contingency flow slack containers for formulation +# `F`, or `(nothing, nothing)` when slacks are disabled. Shared by the PTDF and +# AreaBalance flow-rate constraints. +function _make_post_contingency_slacks( + container::OptimizationContainer, + service::R, + service_name::String, + resolved::Vector{ + Pair{Base.UUID, Vector{Tuple{DataType, String, Tuple{Int, Int}, String}}}, + }, + ::Type{F}, + use_slacks::Bool, +) where {R <: PSY.AbstractReserve, F <: AbstractSecurityConstrainedReservesFormulation} + use_slacks || return (nothing, nothing) + slack_ub = add_post_contingency_slack_variables!( + container, PostContingencyFlowActivePowerSlackUpperBound, + service, service_name, resolved, F, + ) + slack_lb = add_post_contingency_slack_variables!( + container, PostContingencyFlowActivePowerSlackLowerBound, + service, service_name, resolved, F, + ) + return (slack_ub, slack_lb) +end + +# Emit the `[lb, ub]` post-contingency flow inequalities for one +# `(outage_id, name, t)`. Slacks-disabled and slacks-enabled paths dispatch on the +# slack arguments rather than branching at run time. +function _add_post_contingency_flow_rate_constraint!( + jump_model, + con_ub, + con_lb, + post_cont_flow, + ::Nothing, + ::Nothing, + outage_id::String, + name::String, + t::Int, + lb, + ub, +) + con_ub[outage_id, name, t] = + JuMP.@constraint(jump_model, post_cont_flow[outage_id, name, t] <= ub) + con_lb[outage_id, name, t] = + JuMP.@constraint(jump_model, post_cont_flow[outage_id, name, t] >= lb) + return +end + +function _add_post_contingency_flow_rate_constraint!( + jump_model, + con_ub, + con_lb, + post_cont_flow, + slack_ub, + slack_lb, + outage_id::String, + name::String, + t::Int, + lb, + ub, +) + con_ub[outage_id, name, t] = JuMP.@constraint( + jump_model, + post_cont_flow[outage_id, name, t] - slack_ub[outage_id, name, t] <= ub, + ) + con_lb[outage_id, name, t] = JuMP.@constraint( + jump_model, + post_cont_flow[outage_id, name, t] + slack_lb[outage_id, name, t] >= lb, + ) + return +end + +# ---------------------------------------------------------------------------- +# Reserve deployment variable per (outage, contributing device, t) +# ---------------------------------------------------------------------------- + +""" +Add the per-contingency reserve-deployment scenario variables `Δrsv[c, g, t]` — +how much each contributing generator `g` redeploys under contingency `c` at time +`t`. These are separate from the base dispatch `p[g, t]`: G-1 is corrective in +formulation (each contingency gets its own deployment) but preventive in purpose +(the deployment must be backed by reserve procured pre-contingency). The outaged +generator deploys nothing under its own contingency, so its variable is pinned to +zero; every other variable takes its bounds/warm start from the `get_variable_*` +accessors. +""" +function add_variables!( + container::OptimizationContainer, + sys::PSY.System, + variable_type::Type{T}, + service::R, + service_model::ServiceModel{R, <:AbstractSecurityConstrainedReservesFormulation}, + contributing_devices::Vector{V}, + ::Type{F}, +) where { + T <: AbstractContingencyVariableType, + R <: PSY.AbstractReserve, + V <: PSY.StaticInjection, + F <: AbstractSecurityConstrainedReservesFormulation, +} + @assert !isempty(contributing_devices) + time_steps = get_time_steps(container) + binary = get_variable_binary(variable_type, R, F) + service_name = PSY.get_name(service) + + # Outages claimed by this service live on `service_model.outages` (UUID + # keys). Resolve them to the supplemental attribute objects so we can + # query the associated outaged generators and pin their deployment + # variables to zero under their own contingency. + associated_outages = _service_outages(sys, service_model) + outage_ids = string.(IS.get_uuid.(associated_outages)) + + variable = lazy_container_addition!(container, variable_type, + R, + outage_ids, + [PSY.get_name(d) for d in contributing_devices], + time_steps; + meta = service_name, + ) + + for outage in associated_outages + outage_id = string(IS.get_uuid(outage)) + associated_devices = + PSY.get_associated_components(sys, outage; component_type = PSY.Generator) + for device in contributing_devices + name = PSY.get_name(device) + device_outaged = device in associated_devices + # Bounds and warm start are constant across time steps. + ub = get_variable_upper_bound(variable_type, service, device, F) + lb = get_variable_lower_bound(variable_type, service, device, F) + init = get_variable_warm_start_value(variable_type, device, F) + for t in time_steps + v = JuMP.@variable( + get_jump_model(container), + base_name = "$(T)_$(R)_$(service_name)_{$(outage_id), $(name), $(t)}", + binary = binary, + ) + variable[outage_id, name, t] = v + if device_outaged + # The outaged generator cannot deploy reserves for its own + # contingency; force the variable to zero. + JuMP.set_upper_bound(v, 0.0) + JuMP.set_lower_bound(v, 0.0) + JuMP.set_start_value(v, 0.0) + continue + end + ub === nothing || JuMP.set_upper_bound(v, ub) + (lb === nothing || binary) || JuMP.set_lower_bound(v, lb) + init === nothing || JuMP.set_start_value(v, init) + end + end + end + return +end + +# ---------------------------------------------------------------------------- +# Post-contingency power-balance, nodal-deployment, area-deployment +# expressions. Reserve-deployment contributions and generator-outage +# contributions are added in separate dispatches to avoid `isa` checks. +# ---------------------------------------------------------------------------- + +""" +Add the healthy generators' reserve **deployments** to the per-outage system +power-balance expression: `+Σ_g mult·Δrsv[c, g, t]`. Paired with the +pre-contingency method below (which subtracts the lost unit's `p[g_outaged, t]`), +the `PostContingencyGenerationBalanceConstraint` then forces the deployed reserve +to exactly replace the outaged generation. The outaged generator deploys nothing, +so it is skipped rather than added with a zero coefficient. +""" +function add_to_expression!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + ::Type{U}, + contributing_devices::Union{IS.FlattenIteratorWrapper{V}, Vector{V}}, + service::R, + service_model::ServiceModel{R, F}, + ::NetworkModel{<:PM.AbstractPowerModel}, +) where { + T <: PostContingencyActivePowerBalance, + U <: AbstractContingencyVariableType, + V <: PSY.Generator, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + associated_outages = _service_outages(sys, service_model) + expression = lazy_container_addition!(container, T, + R, + string.(IS.get_uuid.(associated_outages)), + time_steps; + meta = service_name, + ) + reserve_deployment_variable = get_variable(container, U, R, service_name) + mult_default = get_variable_multiplier(U, R, F) + for outage in associated_outages + associated_devices = + PSY.get_associated_components(sys, outage; component_type = PSY.Generator) + outage_id = string(IS.get_uuid(outage)) + for device in contributing_devices + device in associated_devices && continue + name = PSY.get_name(device) + for t in time_steps + add_proportional_to_jump_expression!( + expression[outage_id, t], + reserve_deployment_variable[outage_id, name, t], + mult_default, + ) + end + end + end + return +end + +""" +Add the outaged generators' **pre-contingency** output to the per-outage system +power-balance expression: `−p[g_outaged, t]` (the `Generator` multiplier is +`-1.0`). This is the generation that the healthy units' deployments (above) must +make up. Iterates the system's `(generator, outage)` attribute pairs, keeping only +those whose outage this service responds to. +""" +function add_to_expression!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + ::Type{U}, + attribute_device_map::Vector{ + NamedTuple{(:component, :supplemental_attribute), Tuple{V, PSY.Outage}}, + }, + service::R, + service_model::ServiceModel{R, F}, + ::NetworkModel{<:PM.AbstractPowerModel}, +) where { + T <: PostContingencyActivePowerBalance, + U <: VariableType, + V <: PSY.Generator, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + associated_outages = Set(_service_outages(sys, service_model)) + expression = get_expression(container, T, R, service_name) + for (d, outage) in attribute_device_map + outage in associated_outages || continue + outage_id = string(IS.get_uuid(outage)) + name = PSY.get_name(d) + variable = get_variable(container, U, typeof(d)) + mult = get_variable_multiplier(U, typeof(d), F) + for t in time_steps + add_proportional_to_jump_expression!( + expression[outage_id, t], + variable[name, t], + mult, + ) + end + end + return +end + +""" +PTDF path: place the healthy generators' reserve **deployments** at their buses, +building the per-outage nodal net-injection change `Δinj[c, b, t]`. This nodal +deployment vector is what the monitored-branch flow expressions multiply by the +PTDF rows to check post-contingency line loading. The outaged generator deploys +nothing, so it is skipped. +""" +function add_to_expression!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + ::Type{U}, + contributing_devices::Union{IS.FlattenIteratorWrapper{V}, Vector{V}}, + service::R, + service_model::ServiceModel{R, F}, + network_model::NetworkModel{N}, +) where { + T <: PostContingencyNodalActivePowerDeployment, + U <: AbstractContingencyVariableType, + V <: PSY.Generator, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, + N <: AbstractPTDFModel, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + associated_outages = _service_outages(sys, service_model) + ptdf = get_PTDF_matrix(network_model) + bus_numbers = PNM.get_bus_axis(ptdf) + expression = lazy_container_addition!(container, T, + R, + string.(IS.get_uuid.(associated_outages)), + bus_numbers, + time_steps; + meta = service_name, + ) + reserve_deployment_variable = get_variable(container, U, R, service_name) + mult_default = get_variable_multiplier(U, R, F) + network_reduction = get_network_reduction(network_model) + for outage in associated_outages + associated_devices = + PSY.get_associated_components(sys, outage; component_type = PSY.Generator) + outage_id = string(IS.get_uuid(outage)) + for device in contributing_devices + device in associated_devices && continue + name = PSY.get_name(device) + bus_number = + PNM.get_mapped_bus_number(network_reduction, PSY.get_bus(device)) + for t in time_steps + add_proportional_to_jump_expression!( + expression[outage_id, bus_number, t], + reserve_deployment_variable[outage_id, name, t], + mult_default, + ) + end + end + end + return +end + +""" +PTDF path: place the outaged generators' lost **pre-contingency** output +`−p[g_outaged, t]` at their buses, completing the per-outage nodal net-injection +change. Iterates the system's `(generator, outage)` attribute pairs, keeping only +the outages this service responds to. +""" +function add_to_expression!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + ::Type{U}, + attribute_device_map::Vector{ + NamedTuple{(:component, :supplemental_attribute), Tuple{V, PSY.Outage}}, + }, + service::R, + service_model::ServiceModel{R, F}, + network_model::NetworkModel{N}, +) where { + T <: PostContingencyNodalActivePowerDeployment, + U <: VariableType, + V <: PSY.Generator, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, + N <: AbstractPTDFModel, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + associated_outages = Set(_service_outages(sys, service_model)) + expression = get_expression(container, T, R, service_name) + network_reduction = get_network_reduction(network_model) + for (device, outage) in attribute_device_map + outage in associated_outages || continue + outage_id = string(IS.get_uuid(outage)) + name = PSY.get_name(device) + variable = get_variable(container, U, typeof(device)) + mult = get_variable_multiplier(U, typeof(device), F) + bus_number = PNM.get_mapped_bus_number(network_reduction, PSY.get_bus(device)) + for t in time_steps + add_proportional_to_jump_expression!( + expression[outage_id, bus_number, t], + variable[name, t], + mult, + ) + end + end + return +end + +""" +AreaBalance path: aggregate the healthy generators' reserve **deployments** into +their areas, building the per-outage area net-injection change. Paired with the +pre-contingency method below, this feeds the `PostContingencyCopperPlateBalance` +constraint per area. The outaged generator deploys nothing, so it is skipped. +""" +function add_to_expression!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + ::Type{U}, + contributing_devices::Union{IS.FlattenIteratorWrapper{V}, Vector{V}}, + service::R, + service_model::ServiceModel{R, F}, + ::NetworkModel{<:AreaBalancePowerModel}, +) where { + T <: PostContingencyAreaActivePowerDeployment, + U <: AbstractContingencyVariableType, + V <: PSY.Generator, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + associated_outages = _service_outages(sys, service_model) + area_names = PSY.get_name.(PSY.get_components(PSY.Area, sys)) + expression = lazy_container_addition!(container, T, + R, + string.(IS.get_uuid.(associated_outages)), + area_names, + time_steps; + meta = service_name, + ) + reserve_deployment_variable = get_variable(container, U, R, service_name) + mult_default = get_variable_multiplier(U, R, F) + for outage in associated_outages + associated_devices = + PSY.get_associated_components(sys, outage; component_type = PSY.Generator) + outage_id = string(IS.get_uuid(outage)) + for device in contributing_devices + device in associated_devices && continue + name = PSY.get_name(device) + area_name = PSY.get_name(PSY.get_area(PSY.get_bus(device))) + for t in time_steps + add_proportional_to_jump_expression!( + expression[outage_id, area_name, t], + reserve_deployment_variable[outage_id, name, t], + mult_default, + ) + end + end + end + return +end + +""" +AreaBalance path: aggregate the outaged generators' lost **pre-contingency** +output `−p[g_outaged, t]` into their areas, completing the per-outage area +net-injection change. Iterates the system's `(generator, outage)` attribute pairs, +keeping only the outages this service responds to. +""" +function add_to_expression!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + ::Type{U}, + service::R, + service_model::ServiceModel{R, F}, + ::NetworkModel{<:AreaBalancePowerModel}, +) where { + T <: PostContingencyAreaActivePowerDeployment, + U <: VariableType, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + attribute_device_map = PSY.get_component_supplemental_attribute_pairs( + PSY.Generator, + PSY.Outage, + sys, + ) + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + associated_outages = Set(_service_outages(sys, service_model)) + expression = get_expression(container, T, R, service_name) + for (device, outage) in attribute_device_map + outage in associated_outages || continue + outage_id = string(IS.get_uuid(outage)) + name = PSY.get_name(device) + variable = get_variable(container, U, typeof(device)) + mult = get_variable_multiplier(U, typeof(device), F) + area_name = PSY.get_name(PSY.get_area(PSY.get_bus(device))) + for t in time_steps + add_proportional_to_jump_expression!( + expression[outage_id, area_name, t], + variable[name, t], + mult, + ) + end + end + return +end + +""" +Build the per-(outage, generator, t) post-contingency generation expression +`p[g, t] + Δrsv[c, g, t]` used when the service has no reserve-requirement time +series. With no requirement to bound the deployment, the post-contingency output +itself is bounded directly by the generator's min/max via +`PostContingencyActivePowerGenerationLimitsConstraint` — this is where +`p[g, t] + Δrsv ≤ Pᵐᵃˣ` is enforced. The outaged generator contributes only its +(pinned-to-zero) deployment, so its post-contingency output is forced to zero. +`PostContingencyActivePowerGeneration` is dense over the contributing devices. +""" +function add_to_expression!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + contributing_devices::Union{IS.FlattenIteratorWrapper{V}, Vector{V}}, + service::R, + service_model::ServiceModel{R, F}, + ::NetworkModel{<:PM.AbstractActivePowerModel}, +) where { + T <: PostContingencyActivePowerGeneration, + V <: PSY.Generator, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + associated_outages = _service_outages(sys, service_model) + expression = add_expression_container!(container, T, + R, + string.(IS.get_uuid.(associated_outages)), + PSY.get_name.(contributing_devices), + time_steps; + meta = service_name, + ) + reserve_deployment_variable = + get_variable(container, PostContingencyActivePowerReserveDeploymentVariable, + R, + service_name, + ) + for device in contributing_devices + gen_var = get_variable(container, ActivePowerVariable, typeof(device)) + gen_name = PSY.get_name(device) + for outage in associated_outages + associated_devices = PSY.get_associated_components( + sys, outage; component_type = PSY.Generator, + ) + outage_id = string(IS.get_uuid(outage)) + gen_outaged = device in associated_devices + for t in time_steps + add_proportional_to_jump_expression!( + expression[outage_id, gen_name, t], + reserve_deployment_variable[outage_id, gen_name, t], + 1.0, + ) + gen_outaged && continue + add_proportional_to_jump_expression!( + expression[outage_id, gen_name, t], + gen_var[gen_name, t], + 1.0, + ) + end + end + end + return +end + +# ---------------------------------------------------------------------------- +# Sparse-monitored post-contingency flow expression (PTDF only): +# flow[c, ℓ, t] = pre_flow[ℓ, t] + Σ_b PTDF[ℓ, b] * deployment[c, b, t] +# Only built for monitored components carried by the service-claimed outages +# in `service_model.outages`. The branch type is taken from the monitored +# tuple so the correct `PTDFBranchFlow` container is consulted per component. +# ---------------------------------------------------------------------------- + +""" +Build the post-contingency monitored-branch flow expression +`flow[c, ℓ, t] = pre_flow[ℓ, t] + Σ_b PTDF[ℓ, b]·Δinj[c, b, t]` — the line loading +after the reserve redeployment of contingency `c`. The companion +`PostContingencyFlowRateConstraint` then bounds this by the branch's emergency +rating, i.e. the redeployment must keep monitored branches within their G-1 +ratings. Only the monitored components carried by the service-claimed outages are +instantiated, in a sparse `(outage_id, name, t)` container. +""" +function add_post_contingency_flow_expressions!( + container::OptimizationContainer, + ::Type{T}, + service::R, + service_model::ServiceModel{R, F}, + network_model::NetworkModel{N}, +) where { + T <: PostContingencyBranchFlow, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, + N <: AbstractPTDFModel, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + net_reduction_data = network_model.network_reduction + resolved = _resolve_service_monitored_arcs(service_model, net_reduction_data) + # Empty sparse `(outage_id, name, t)` expression container; the build loop + # fills the ragged entries by assignment. + expression_container = IOM.add_expression_container!( + container, T, R, String[], String[], time_steps; + sparse = true, meta = service_name, + ) + isempty(resolved) && return expression_container + + nodal_deployment = get_expression(container, PostContingencyNodalActivePowerDeployment, + R, + service_name, + ) + ptdf = get_PTDF_matrix(network_model) + + # Cache pre-contingency flow expressions by monitored type and PTDF rows by + # arc. `ptdf` is a `VirtualPTDF`: `ptdf[arc, :]` solves the row via KLU on + # first access (the internal cache is bounded and may evict), so caching the + # row here avoids repeating that solve for an arc monitored across multiple + # outages/time steps. + pre_flow_cache = Dict{DataType, Any}() + ptdf_col_cache = Dict{Any, Any}() + for (uuid, entries) in resolved + outage_id = string(uuid) + # Positional slice over the bus/time axes; matches `ptdf_col`'s + # positional indexing and avoids keyed lookup mismatches when the + # nodal expression's bus axis is a subset of the PTDF column space + # (e.g. AreaPTDFPowerModel). + post_cont_expr = nodal_deployment[outage_id, :, :].data + for (entry_type, name, arc, _) in entries + pre_flow = get!(pre_flow_cache, entry_type) do + get_expression(container, PTDFBranchFlow, entry_type) + end + ptdf_col = get!(() -> ptdf[arc, :], ptdf_col_cache, arc) + for t in time_steps + acc = JuMP.AffExpr(0.0) + JuMP.add_to_expression!(acc, pre_flow[name, t]) + @inbounds for b in eachindex(ptdf_col) + coef = ptdf_col[b] + abs(coef) < PTDF_ZERO_TOL && continue + JuMP.add_to_expression!(acc, coef, post_cont_expr[b, t]) + end + expression_container[outage_id, name, t] = acc + end + end + end + return expression_container +end + +# ---------------------------------------------------------------------------- +# Post-contingency constraints +# ---------------------------------------------------------------------------- + +""" +Per-outage system-wide generation-balance constraint: the +`PostContingencyActivePowerBalance` expression (sum of reserve deployments +minus the outaged generation) must close to zero for every outage and time. +""" +function add_constraints!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + ::Type{U}, + ::Union{IS.FlattenIteratorWrapper{V}, Vector{V}}, + service::R, + service_model::ServiceModel{R, F}, + ::NetworkModel{<:PM.AbstractPowerModel}, +) where { + T <: PostContingencyGenerationBalanceConstraint, + U <: PostContingencyActivePowerBalance, + V <: PSY.Generator, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + associated_outages = _service_outages(sys, service_model) + expressions = get_expression(container, U, R, service_name) + constraint = add_constraints_container!(container, T, + R, + [string(IS.get_uuid(o)) for o in associated_outages], + time_steps; + meta = service_name, + ) + jump_model = get_jump_model(container) + for outage in associated_outages, t in time_steps + outage_id = string(IS.get_uuid(outage)) + constraint[outage_id, t] = + JuMP.@constraint(jump_model, expressions[outage_id, t] == 0) + end + return +end + +""" +Sparse-monitored post-contingency branch flow inequalities. The container is +keyed by `(outage_id::String, monitored_name::String, t::Int)` and only +entries resolved from `service_model.outages` are populated. Limits use the +monitored branch's emergency rating. Optional non-negative slacks relax the +inequalities at `POST_CONTINGENCY_CONSTRAINT_VIOLATION_SLACK_COST`. +""" +function add_constraints!( + container::OptimizationContainer, + ::Type{T}, + ::Type{U}, + service::R, + service_model::ServiceModel{R, F}, + network_model::NetworkModel{<:AbstractPTDFModel}, +) where { + T <: PostContingencyFlowRateConstraint, + U <: PostContingencyBranchFlow, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + net_reduction_data = network_model.network_reduction + all_branch_maps_by_type = PNM.get_all_branch_maps_by_type(net_reduction_data) + resolved = _resolve_service_monitored_arcs(service_model, net_reduction_data) + + con_lb = IOM.add_constraints_container!( + container, T, R, String[], String[], time_steps; + sparse = true, meta = "$(service_name)_$(POST_CONTINGENCY_LB_META)", + ) + con_ub = IOM.add_constraints_container!( + container, T, R, String[], String[], time_steps; + sparse = true, meta = "$(service_name)_$(POST_CONTINGENCY_UB_META)", + ) + isempty(resolved) && return + + post_cont_flow = get_expression(container, U, R, service_name) + jump_model = get_jump_model(container) + + slack_ub, slack_lb = _make_post_contingency_slacks( + container, service, service_name, resolved, F, get_use_slacks(service_model), + ) + + for (uuid, entries) in resolved + outage_id = string(uuid) + for (entry_type, name, arc, reduction_kind) in entries + reduction_entry = + all_branch_maps_by_type[reduction_kind][entry_type][arc] + limits = get_emergency_min_max_limits( + reduction_entry, T, StaticBranch, + ) + for t in time_steps + _add_post_contingency_flow_rate_constraint!( + jump_model, con_ub, con_lb, post_cont_flow, slack_ub, slack_lb, + outage_id, name, t, limits.min, limits.max, + ) + end + end + end + return +end + +# ---------------------------------------------------------------------------- +# AreaBalance network model: post-contingency AreaInterchange flow expression +# and rate-limit constraints. The expression is keyed by +# `(outage_id::String, area_interchange_name::String, t::Int)` and represents +# the from→to flow after the contingency: +# post_flow = pre_flow + Σ_{g ∈ from} deploy_g - Σ_{g ∈ to} deploy_g +# - sign(outaged_side) * P_outaged +# where `sign(outaged_side)` is `+1` if the outaged generator sits in the +# from-area and `-1` if it sits in the to-area. Deployment variables for the +# outaged generator are pinned to zero, so iterating contributing devices is +# safe. Only the AreaInterchanges named in `service_model.outages[uuid]` are +# instantiated. +# ---------------------------------------------------------------------------- + +""" +Resolve every monitored `PSY.AreaInterchange` carried by +`service_model.outages` to its system component. Mirrors +`_resolve_service_monitored_arcs` but for the AreaBalance path where the +monitored object is an AreaInterchange (not a branch arc) and the only +information needed downstream is the component itself. Outages with no +monitored AreaInterchanges are skipped; the returned vector is sorted by +UUID for deterministic axes. +""" +function _resolve_service_monitored_area_interchanges( + sys::PSY.System, + service_model::ServiceModel, +) + resolved = Pair{Base.UUID, Vector{Tuple{String, PSY.AreaInterchange}}}[] + for (uuid, per_type) in get_outages(service_model) + kept = Tuple{String, PSY.AreaInterchange}[] + names = get(per_type, PSY.AreaInterchange, nothing) + if names !== nothing + for name in sort!(collect(names)) + comp = PSY.get_component(PSY.AreaInterchange, sys, name) + comp === nothing && continue + push!(kept, (name, comp)) + end + end + isempty(kept) && continue + push!(resolved, uuid => kept) + end + sort!(resolved; by = first) + return resolved +end + +""" +Build the post-contingency AreaInterchange flow expression for the +AreaBalance network model. See module-level comment above for the formula. +The container is a `SparseAxisArray` keyed by +`(outage_id, area_interchange_name, t)` registered under +`ExpressionKey(PostContingencyAreaInterchangeFlow, R; meta = service_name)`. +""" +function add_post_contingency_flow_expressions!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + service::R, + service_model::ServiceModel{R, F}, + ::NetworkModel{<:AreaBalancePowerModel}, +) where { + T <: PostContingencyAreaInterchangeFlow, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + resolved = _resolve_service_monitored_area_interchanges(sys, service_model) + + expression_container = IOM.add_expression_container!( + container, T, R, String[], String[], time_steps; + sparse = true, meta = service_name, + ) + isempty(resolved) && return expression_container + + # Baseline flow variable for the AreaInterchange (from→to convention). + flow_var = get_variable(container, FlowActivePowerVariable, PSY.AreaInterchange) + + # Reserve-deployment variable keyed by (outage_id, gen_name, t). + reserve_deployment_variable = + get_variable(container, PostContingencyActivePowerReserveDeploymentVariable, + R, + service_name, + ) + contributing_devices = get_contributing_devices(service_model) + + # Pre-compute area assignment for each contributing device so we can + # apply +1 for from-area and -1 for to-area without `isa` checks. + device_areas = Dict{String, String}() + for device in contributing_devices + device_areas[PSY.get_name(device)] = + PSY.get_name(PSY.get_area(PSY.get_bus(device))) + end + + associated_outages = _service_outages(sys, service_model) + outage_by_uuid = Dict(IS.get_uuid(o) => o for o in associated_outages) + for (uuid, entries) in resolved + outage = outage_by_uuid[uuid] + outage_id = string(uuid) + outaged_gens = PSY.get_associated_components( + sys, outage; component_type = PSY.Generator, + ) + for (name, area_interchange) in entries + from_area = PSY.get_name(PSY.get_from_area(area_interchange)) + to_area = PSY.get_name(PSY.get_to_area(area_interchange)) + for t in time_steps + # Initialize the entry, then accumulate; stored by the write-back below. + expr = zero(JuMP.AffExpr) + JuMP.add_to_expression!(expr, flow_var[name, t]) + for device in contributing_devices + gen_name = PSY.get_name(device) + gen_area = device_areas[gen_name] + coef = if gen_area == from_area + 1.0 + elseif gen_area == to_area + -1.0 + else + 0.0 + end + coef == 0.0 && continue + add_proportional_to_jump_expression!( + expr, + reserve_deployment_variable[outage_id, gen_name, t], + coef, + ) + end + # Subtract the outaged generation contribution on the + # outaged side: +pre-contingency power if in from-area, + # -pre-contingency power if in to-area. + for outaged_gen in outaged_gens + outaged_area = + PSY.get_name(PSY.get_area(PSY.get_bus(outaged_gen))) + coef = if outaged_area == from_area + -1.0 + elseif outaged_area == to_area + 1.0 + else + 0.0 + end + coef == 0.0 && continue + gen_var = + get_variable(container, ActivePowerVariable, typeof(outaged_gen) + ) + add_proportional_to_jump_expression!( + expr, gen_var[PSY.get_name(outaged_gen), t], coef, + ) + end + expression_container[outage_id, name, t] = expr + end + end + end + return expression_container +end + +""" +Per-(outage, area_interchange, t) flow-rate inequalities under the +AreaBalance network model. Limits come from +`PSY.get_flow_limits(area_interchange, PSY.SU)` as `[-from_to, +to_from]`. Optional +non-negative slacks relax the inequalities at +`POST_CONTINGENCY_CONSTRAINT_VIOLATION_SLACK_COST`. +""" +function add_constraints!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + ::Type{U}, + service::R, + service_model::ServiceModel{R, F}, + ::NetworkModel{<:AreaBalancePowerModel}, +) where { + T <: PostContingencyFlowRateConstraint, + U <: PostContingencyAreaInterchangeFlow, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + resolved = _resolve_service_monitored_area_interchanges(sys, service_model) + + con_lb = IOM.add_constraints_container!( + container, T, R, String[], String[], time_steps; + sparse = true, meta = "$(service_name)_$(POST_CONTINGENCY_LB_META)", + ) + con_ub = IOM.add_constraints_container!( + container, T, R, String[], String[], time_steps; + sparse = true, meta = "$(service_name)_$(POST_CONTINGENCY_UB_META)", + ) + isempty(resolved) && return + + post_cont_flow = get_expression(container, U, R, service_name) + jump_model = get_jump_model(container) + + # The shared slack helper expects the PTDF-style resolved tuple shape + # `(type, name, arc, reduction_kind)`; build an equivalent shape with + # placeholder arc/reduction values so the slack containers are keyed by the + # same `(outage_id, name, t)`. + slack_resolved = + Pair{Base.UUID, Vector{Tuple{DataType, String, Tuple{Int, Int}, String}}}[ + uuid => Tuple{DataType, String, Tuple{Int, Int}, String}[ + (PSY.AreaInterchange, name, (0, 0), "") for (name, _) in entries + ] for (uuid, entries) in resolved + ] + slack_ub, slack_lb = _make_post_contingency_slacks( + container, service, service_name, slack_resolved, F, + get_use_slacks(service_model), + ) + + for (uuid, entries) in resolved + outage_id = string(uuid) + for (name, area_interchange) in entries + flow_limits = PSY.get_flow_limits(area_interchange, PSY.SU) + ub = flow_limits.to_from + lb = -1.0 * flow_limits.from_to + for t in time_steps + _add_post_contingency_flow_rate_constraint!( + jump_model, con_ub, con_lb, post_cont_flow, slack_ub, slack_lb, + outage_id, name, t, lb, ub, + ) + end + end + end + return +end + +""" +Per-outage upper bound on the reserve-deployment variable by the +pre-contingency reserve variable. Outaged generators are pinned to zero. +""" +function add_constraints!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + ::Type{X}, + ::Type{U}, + contributing_devices::Union{IS.FlattenIteratorWrapper{V}, Vector{V}}, + service::R, + service_model::ServiceModel{R, F}, + ::NetworkModel{<:PM.AbstractPowerModel}, +) where { + T <: PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + X <: VariableType, + U <: AbstractContingencyVariableType, + V <: PSY.Generator, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + associated_outages = _service_outages(sys, service_model) + constraint = add_constraints_container!(container, T, + R, + [string(IS.get_uuid(o)) for o in associated_outages], + [PSY.get_name(d) for d in contributing_devices], + time_steps; + meta = service_name, + ) + variable = get_variable(container, X, R, service_name) + variable_outage = get_variable(container, U, R, service_name) + jump_model = get_jump_model(container) + for outage in associated_outages + associated_devices = PSY.get_associated_components( + sys, outage; component_type = PSY.Generator, + ) + outage_id = string(IS.get_uuid(outage)) + for device in contributing_devices + name = PSY.get_name(device) + gen_outaged = device in associated_devices + for t in time_steps + if gen_outaged + constraint[outage_id, name, t] = JuMP.@constraint( + jump_model, + variable_outage[outage_id, name, t] == 0.0, + ) + continue + end + constraint[outage_id, name, t] = JuMP.@constraint( + jump_model, + variable_outage[outage_id, name, t] <= variable[name, t], + ) + end + end + end + return +end + +""" +Used when the service has no reserve requirement time series. +""" +function add_constraints!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + contributing_devices::Union{IS.FlattenIteratorWrapper{V}, Vector{V}}, + service::R, + service_model::ServiceModel{R, F}, + ::NetworkModel{<:PM.AbstractActivePowerModel}, +) where { + T <: PostContingencyActivePowerGenerationLimitsConstraint, + V <: PSY.Generator, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + time_steps = get_time_steps(container) + service_name = PSY.get_name(service) + associated_outages = _service_outages(sys, service_model) + con_lb = add_constraints_container!(container, T, + R, + string.(IS.get_uuid.(associated_outages)), + PSY.get_name.(contributing_devices), + time_steps; + meta = "$(service_name)_$(POST_CONTINGENCY_LB_META)", + ) + con_ub = add_constraints_container!(container, T, + R, + string.(IS.get_uuid.(associated_outages)), + PSY.get_name.(contributing_devices), + time_steps; + meta = "$(service_name)_$(POST_CONTINGENCY_UB_META)", + ) + expressions = + get_expression(container, PostContingencyActivePowerGeneration, R, service_name) + jump_model = get_jump_model(container) + for device in contributing_devices + name = PSY.get_name(device) + limits = PSY.get_active_power_limits(device, PSY.SU) + for outage in associated_outages + associated_devices = PSY.get_associated_components( + sys, outage; component_type = PSY.Generator, + ) + outage_id = string(IS.get_uuid(outage)) + gen_outaged = device in associated_devices + for t in time_steps + if gen_outaged + con_ub[outage_id, name, t] = JuMP.@constraint( + jump_model, expressions[outage_id, name, t] == 0.0, + ) + con_lb[outage_id, name, t] = JuMP.@constraint( + jump_model, expressions[outage_id, name, t] == 0.0, + ) + continue + end + con_ub[outage_id, name, t] = JuMP.@constraint( + jump_model, expressions[outage_id, name, t] <= limits.max, + ) + con_lb[outage_id, name, t] = JuMP.@constraint( + jump_model, expressions[outage_id, name, t] >= limits.min, + ) + end + end + end + return +end + +""" +AreaBalance path: per-(outage, area, t) power balance after the contingency. The +per-outage area net-injection change (reserve deployments minus the lost unit's +output) plus the base area balance must close to zero — the area-level analogue of +`PostContingencyGenerationBalanceConstraint`, requiring deployed reserve to +replace the outaged generation within each area. +""" +function add_constraints!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + ::Type{U}, + ::Type{Y}, + service::R, + service_model::ServiceModel{R, F}, + ::NetworkModel{<:AreaBalancePowerModel}, +) where { + T <: PostContingencyCopperPlateBalanceConstraint, + U <: PostContingencyAreaActivePowerDeployment, + Y <: ActivePowerBalance, + R <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + time_steps = get_time_steps(container) + devices = PSY.get_components(PSY.Area, sys) + area_names = PSY.get_name.(devices) + service_name = PSY.get_name(service) + associated_outages = _service_outages(sys, service_model) + con = add_constraints_container!(container, T, + R, + string.(IS.get_uuid.(associated_outages)), + area_names, + time_steps; + meta = service_name, + ) + contingency_expression = get_expression(container, U, R, service_name) + area_expression = get_expression(container, Y, PSY.Area) + jump_model = get_jump_model(container) + for outage in associated_outages + outage_id = string(IS.get_uuid(outage)) + for area in devices + area_name = PSY.get_name(area) + for t in time_steps + con[outage_id, area_name, t] = JuMP.@constraint( + jump_model, + contingency_expression[outage_id, area_name, t] + + area_expression[area_name, t] == 0.0, + ) + end + end + end + return +end + +# ---------------------------------------------------------------------------- +# construct_service! dispatches: argument + model construct stages for +# (SecurityConstrainedContingencyReserve, SecurityConstrainedRampReserve) × +# (PTDF, CopperPlate, AreaBalance). +# ---------------------------------------------------------------------------- + +""" +ArgumentConstructStage for both security-constrained reserve formulations. Adds +the optional pre-contingency reserve requirement/variable stack, then the +per-contingency reserve-deployment scenario variables `Δrsv[c, g, t]` — the +variables that let the model react to each G-1 outage separately from the base +dispatch. Returns early (with a warning) when the service has no outages attached. +""" +function construct_service!( + container::OptimizationContainer, + sys::PSY.System, + ::ArgumentConstructStage, + model::ServiceModel{SR, F}, + devices_template::Dict{Symbol, DeviceModel}, + ::Set{<:DataType}, + ::NetworkModel{<:PM.AbstractActivePowerModel}, +) where {SR <: PSY.AbstractReserve, F <: AbstractSecurityConstrainedReservesFormulation} + name = get_service_name(model) + service = PSY.get_component(SR, sys, name) + !PSY.get_available(service) && return + contributing_devices = get_contributing_devices(model) + + if _has_requirement_ts(container, model, service) || requires_requirement_ts(F) + add_parameters!(container, RequirementTimeSeriesParameter, service, model) + add_service_variables!( + container, + ActivePowerReserveVariable, + service, + contributing_devices, + F, + ) + add_to_expression!(container, ActivePowerReserveVariable, model, devices_template) + end + add_feedforward_arguments!(container, model, service) + + associated_outages = _service_outages(sys, model) + if isempty(associated_outages) + @warn "Service $(SR)('$name'): `service_model.outages` is empty; the \ + security-constrained formulation $(F) will not add any \ + post-contingency variables or constraints." + return + end + + add_variables!( + container, + sys, + PostContingencyActivePowerReserveDeploymentVariable, + service, + model, + contributing_devices, + F, + ) + return +end + +# Shared ModelConstructStage helper for the pre-contingency requirement, +# ramp, participation, objective, feedforward and dual hookups. +function _construct_service_pre_contingency!( + container::OptimizationContainer, + sys::PSY.System, + service::PSY.AbstractReserve, + contributing_devices, + model::ServiceModel{SR, F}, + has_requirement_ts::Bool, +) where { + SR <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + if has_requirement_ts + add_constraints!( + container, + RequirementConstraint, + service, + contributing_devices, + model, + ) + add_constraints!( + container, RampConstraint, service, contributing_devices, model, + ) + add_constraints!( + container, + ParticipationFractionConstraint, + service, + contributing_devices, + model, + ) + add_to_objective_function!(container, service, model) + end + add_feedforward_constraints!(container, model, service) + add_constraint_dual!(container, sys, model) + return +end + +# Shared helper for the post-contingency power-balance + generation-balance +# expression/constraint stack that's common to every network model. +function _construct_service_post_contingency_balance!( + container::OptimizationContainer, + sys::PSY.System, + service::PSY.AbstractReserve, + contributing_devices, + model::ServiceModel{SR, F}, + network_model::NetworkModel, +) where { + SR <: PSY.AbstractReserve, + F <: AbstractSecurityConstrainedReservesFormulation, +} + add_to_expression!( + container, sys, PostContingencyActivePowerBalance, + PostContingencyActivePowerReserveDeploymentVariable, + contributing_devices, service, model, network_model, + ) + attribute_device_map = PSY.get_component_supplemental_attribute_pairs( + PSY.Generator, PSY.Outage, sys, + ) + add_to_expression!( + container, sys, PostContingencyActivePowerBalance, ActivePowerVariable, + attribute_device_map, service, model, network_model, + ) + add_constraints!( + container, sys, PostContingencyGenerationBalanceConstraint, + PostContingencyActivePowerBalance, + contributing_devices, service, model, network_model, + ) + return attribute_device_map +end + +""" +ModelConstructStage for both security-constrained reserve formulations. Wires up +the G-1 corrective model: the pre-contingency reserve stack, the per-outage +power-balance/generation-balance constraints (deployed reserve must replace the +lost unit), the network-specific deployment/flow terms (dispatched on the network +model via `_add_post_contingency_network_terms!`), and the link tying each +deployment to procured reserve — either `Δrsv ≤ reserve[g, t]` when a requirement +time series is present, or the direct post-contingency generation limits otherwise. +""" +function construct_service!( + container::OptimizationContainer, + sys::PSY.System, + ::ModelConstructStage, + model::ServiceModel{SR, F}, + ::Dict{Symbol, DeviceModel}, + ::Set{<:DataType}, + network_model::NetworkModel{ + <:Union{AbstractPTDFModel, CopperPlatePowerModel, AreaBalancePowerModel}, + }, +) where {SR <: PSY.AbstractReserve, F <: AbstractSecurityConstrainedReservesFormulation} + name = get_service_name(model) + service = PSY.get_component(SR, sys, name) + !PSY.get_available(service) && return + contributing_devices = get_contributing_devices(model) + + has_requirement_ts = + requires_requirement_ts(F) || _has_requirement_ts(container, model, service) + _construct_service_pre_contingency!( + container, sys, service, contributing_devices, model, has_requirement_ts, + ) + + associated_outages = _service_outages(sys, model) + isempty(associated_outages) && return + + attribute_device_map = _construct_service_post_contingency_balance!( + container, sys, service, contributing_devices, model, network_model, + ) + _add_post_contingency_network_terms!( + container, sys, service, contributing_devices, model, network_model, + attribute_device_map, + ) + + if has_requirement_ts + add_constraints!( + container, sys, + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + ActivePowerReserveVariable, + PostContingencyActivePowerReserveDeploymentVariable, + contributing_devices, service, model, network_model, + ) + else + add_to_expression!( + container, sys, PostContingencyActivePowerGeneration, + contributing_devices, service, model, network_model, + ) + add_constraints!( + container, sys, PostContingencyActivePowerGenerationLimitsConstraint, + contributing_devices, service, model, network_model, + ) + end + return +end + +# ----- Network-specific post-contingency deployment/flow terms ----- + +""" +CopperPlate: no network-specific post-contingency terms. With no transmission +representation there are no monitored branches to keep within emergency ratings, +so the shared system-wide power-balance stack is the whole G-1 model. +""" +function _add_post_contingency_network_terms!( + ::OptimizationContainer, + ::PSY.System, + ::PSY.AbstractReserve, + contributing_devices, + ::ServiceModel{SR, F}, + ::NetworkModel{<:CopperPlatePowerModel}, + attribute_device_map, +) where {SR <: PSY.AbstractReserve, F <: AbstractSecurityConstrainedReservesFormulation} + return +end + +""" +PTDF (DC): build the per-outage nodal net-injection change (reserve deployments +plus the lost unit's `−p[g_outaged]`), then the monitored-branch flow expressions +and their emergency-rating constraints — enforcing that the G-1 redeployment keeps +monitored branches within their post-contingency ratings. +""" +function _add_post_contingency_network_terms!( + container::OptimizationContainer, + sys::PSY.System, + service::PSY.AbstractReserve, + contributing_devices, + model::ServiceModel{SR, F}, + network_model::NetworkModel{<:AbstractPTDFModel}, + attribute_device_map, +) where {SR <: PSY.AbstractReserve, F <: AbstractSecurityConstrainedReservesFormulation} + add_to_expression!( + container, sys, PostContingencyNodalActivePowerDeployment, + PostContingencyActivePowerReserveDeploymentVariable, + contributing_devices, service, model, network_model, + ) + add_to_expression!( + container, sys, PostContingencyNodalActivePowerDeployment, ActivePowerVariable, + attribute_device_map, service, model, network_model, + ) + add_post_contingency_flow_expressions!( + container, PostContingencyBranchFlow, service, model, network_model, + ) + add_constraints!( + container, PostContingencyFlowRateConstraint, PostContingencyBranchFlow, + service, model, network_model, + ) + return +end + +""" +AreaBalance: build the per-outage area net-injection change and per-area +post-contingency balance, then the monitored `AreaInterchange` flow expressions +and their rate limits — enforcing that the G-1 redeployment keeps inter-area flows +within their post-contingency limits. +""" +function _add_post_contingency_network_terms!( + container::OptimizationContainer, + sys::PSY.System, + service::PSY.AbstractReserve, + contributing_devices, + model::ServiceModel{SR, F}, + network_model::NetworkModel{<:AreaBalancePowerModel}, + attribute_device_map, +) where {SR <: PSY.AbstractReserve, F <: AbstractSecurityConstrainedReservesFormulation} + add_to_expression!( + container, sys, PostContingencyAreaActivePowerDeployment, + PostContingencyActivePowerReserveDeploymentVariable, + contributing_devices, service, model, network_model, + ) + add_to_expression!( + container, sys, PostContingencyAreaActivePowerDeployment, ActivePowerVariable, + service, model, network_model, + ) + add_constraints!( + container, sys, PostContingencyCopperPlateBalanceConstraint, + PostContingencyAreaActivePowerDeployment, ActivePowerBalance, + service, model, network_model, + ) + add_post_contingency_flow_expressions!( + container, sys, PostContingencyAreaInterchangeFlow, + service, model, network_model, + ) + add_constraints!( + container, sys, PostContingencyFlowRateConstraint, + PostContingencyAreaInterchangeFlow, + service, model, network_model, + ) + return +end diff --git a/test/Project.toml b/test/Project.toml index 7a61d08..d470cec 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/g1-port", 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_device_hvdc.jl b/test/test_device_hvdc.jl index 7c82870..428c126 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -192,8 +192,6 @@ function _generate_test_vsc_sys(; g = g, dc_current = 0.0, reactive_power_from = 0.0, - dc_voltage_control_from = true, - ac_voltage_control_from = true, dc_setpoint_from = 0.0, ac_setpoint_from = 1.0, converter_loss_from = QuadraticCurve(loss_a, loss_b, loss_c), @@ -203,8 +201,6 @@ function _generate_test_vsc_sys(; power_factor_weighting_fraction_from = 1.0, voltage_limits_from = (min = 0.95, max = 1.05), reactive_power_to = 0.0, - dc_voltage_control_to = true, - ac_voltage_control_to = true, dc_setpoint_to = 0.0, ac_setpoint_to = 1.0, converter_loss_to = QuadraticCurve(loss_a, loss_b, loss_c), diff --git a/test/test_static_injection_security_constrained_models.jl b/test/test_static_injection_security_constrained_models.jl new file mode 100644 index 0000000..cdf1f6a --- /dev/null +++ b/test/test_static_injection_security_constrained_models.jl @@ -0,0 +1,2591 @@ +function get_outage_total_power_by_step_dict( + sys::PSY.System, + variables::Dict{String, DataFrame}, + var_name::String, + associated_outages::Vector{PSY.UnplannedOutage}; + col_name::String = "name", +) + required_variables = variables[var_name] + total_variable_dict = Dict{String, Vector{Float64}}() + for outage in associated_outages + outage_name = string(IS.get_uuid(outage)) + devices = PSY.get_associated_components( + sys, + outage; + component_type = PSY.Generator, + ) + outage_power_v = sum(devices) do device + filter( + x -> x[col_name] == PSY.get_name(device), + required_variables, + )[!, "value"] + end + total_variable_dict[outage_name] = outage_power_v + end + return total_variable_dict +end + +function get_reserve_total_power_by_step_dict( + variables::Dict{String, DataFrame}, + var_name::String, + associated_outages::Vector{PSY.UnplannedOutage}, + contributing_devices::Union{ + IS.FlattenIteratorWrapper{<:PSY.Generator}, + Vector{<:PSY.Generator}, + }; + col_name::String = "name2", + outage_col_name::String = "name", +) + required_variables = variables[var_name] + total_variable_dict = Dict{String, Vector{Float64}}() + for outage in associated_outages + outage_name = string(IS.get_uuid(outage)) + # Scope to this outage: the deployment variable is keyed by + # (outage_id, device_name, t), so filtering on the device alone would + # mix deployments across outages. + outage_power_v = sum(contributing_devices) do device + device_name = PSY.get_name(device) + filter( + x -> x[outage_col_name] == outage_name && x[col_name] == device_name, + required_variables, + )[ + !, + "value", + ] + end + total_variable_dict[outage_name] = outage_power_v + end + return total_variable_dict +end + +function test_reserves_deployment( + power_outage::Float64, + reserve_deployment::Float64; + tol::Float64 = 1e-3, +) + @test isapprox(power_outage, reserve_deployment, atol = tol) +end + +function compare_outage_power_and_deployed_reserves( + sys::PSY.System, + res::OptimizationProblemOutputs, + service::PSY.VariableReserve; + tolerance::Float64 = 1e-3, +) + variablesdict = read_variables(res) + associated_outages = + collect(PSY.get_supplemental_attributes(PSY.UnplannedOutage, service)) + # Fall back: when no outages are attached to the reserve service, resolve + # from the system's outages instead. + if isempty(associated_outages) + all_outages = collect(PSY.get_supplemental_attributes(PSY.UnplannedOutage, sys)) + associated_outages = all_outages + end + outage_dict = get_outage_total_power_by_step_dict( + sys, + variablesdict, + "ActivePowerVariable__ThermalStandard", + associated_outages; + col_name = "name", + ) + contributing_devices = PSY.get_contributing_devices(sys, service) + service_name = PSY.get_name(service) + reserve_dict = get_reserve_total_power_by_step_dict( + variablesdict, + "PostContingencyActivePowerReserveDeploymentVariable__VariableReserve__ReserveUp__" * + service_name, + associated_outages, + contributing_devices; + col_name = "name2", + ) + for outage in associated_outages + outage_name = string(IS.get_uuid(outage)) + for i in 1:length(outage_dict[outage_name]) + test_reserves_deployment( + outage_dict[outage_name][i], + reserve_dict[outage_name][i]; + tol = tolerance, + ) + end + end +end + +@testset "G-n with Ramp reserve deliverability constraints Dispatch with responding reserves only up, including reduction of parallel circuits" begin + for add_parallel_line in [true, false] + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + if add_parallel_line + l4 = get_component(Line, c_sys5, "4") + add_equivalent_ac_transmission_with_parallel_circuits!(c_sys5, l4, PSY.Line) + end + systems = [c_sys5] + objfuncs = [GAEVF, GQEVF, GQEVF] + constraint_keys = [ + IOM.ConstraintKey( + ActivePowerVariableLimitsConstraint, + PSY.ThermalStandard, + "lb", + ), + IOM.ConstraintKey( + ActivePowerVariableLimitsConstraint, + PSY.ThermalStandard, + "ub", + ), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "lb"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "ub"), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_ub", + ), + IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.System), + #IOM.ConstraintKey(NetworkFlowConstraint, PSY.Line), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + ] + PTDF_ref = IdDict{System, PTDF}( + c_sys5 => PTDF(c_sys5), + ) + test_results = IdDict{System, Vector{Int}}( + c_sys5 => [360, 0, 600, 432, 72], + ) + test_obj_values = IdDict{System, Float64}( + c_sys5 => 329000.0, + ) + components_outages_cases = IdDict{System, Vector{String}}( + c_sys5 => ["Alta"], + ) + for (ix, sys) in enumerate(systems) + gen = get_component(ThermalStandard, sys, "Solitude") + set_ramp_limits!(gen, (up = 0.4 * PSY.SU, down = 0.4 * PSY.SU)) #Increase ramp limits to make the problem feasible + components_outages_names = components_outages_cases[sys] + reserve_up = get_component(VariableReserve{ReserveUp}, sys, "Reserve1") + for component_name in components_outages_names + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + template = get_thermal_dispatch_template_network( + NetworkModel(PTDFPowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1", + )) + + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results[sys]..., + false, + ) + psi_checkobjfun_test(ps_model, objfuncs[ix]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_values[sys], + 10000, + ) + res = OptimizationProblemOutputs(ps_model) + compare_outage_power_and_deployed_reserves( + sys, + res, + reserve_up) + end + end +end + +# Exercises the per-service line-scoping path in +# `_monitored_components_by_modeled_type` and the downstream +# `PostContingencyFlowRateConstraint` build for `PTDFPowerModel` when the +# reserve service monitors a strict subset of the system's AC lines instead of +# every line. The constraint key meta does not change (it remains keyed by +# service name) but the per-outage flow constraint container ends up with +# fewer entries, which lowers the MOI counts compared to the all-lines variant. +@testset "G-n with Ramp reserve deliverability constraints PTDFPowerModel with monitored line subset" begin + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + monitored_line_names = ["1", "2"] + systems = [c_sys5] + objfuncs = [GAEVF, GQEVF, GQEVF] + constraint_keys = [ + IOM.ConstraintKey( + ActivePowerVariableLimitsConstraint, + PSY.ThermalStandard, + "lb", + ), + IOM.ConstraintKey( + ActivePowerVariableLimitsConstraint, + PSY.ThermalStandard, + "ub", + ), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "lb"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "ub"), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_ub", + ), + IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.System), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + ] + PTDF_ref = IdDict{System, PTDF}( + c_sys5 => PTDF(c_sys5), + ) + # Counts are smaller than the all-lines baseline `[360, 0, 600, 432, 72]` + # because only the monitored subset contributes + # `PostContingencyFlowRateConstraint` rows per outage step. + test_results = IdDict{System, Vector{Int}}( + c_sys5 => [360, 0, 504, 336, 72], + ) + test_obj_values = IdDict{System, Float64}( + c_sys5 => 329000.0, + ) + components_outages_cases = IdDict{System, Vector{String}}( + c_sys5 => ["Alta"], + ) + for (ix, sys) in enumerate(systems) + gen = get_component(ThermalStandard, sys, "Solitude") + set_ramp_limits!(gen, (up = 0.4 * PSY.SU, down = 0.4 * PSY.SU)) #Increase ramp limits to make the problem feasible + components_outages_names = components_outages_cases[sys] + reserve_up = get_component(VariableReserve{ReserveUp}, sys, "Reserve1") + monitored_subset = + [get_component(Line, sys, n) for n in monitored_line_names] + for component_name in components_outages_names + # --- Create Outage Data with a hand-picked monitored subset --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = monitored_subset, + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + template = get_thermal_dispatch_template_network( + NetworkModel(PTDFPowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1", + )) + + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results[sys]..., + false, + ) + psi_checkobjfun_test(ps_model, objfuncs[ix]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_values[sys], + 10000, + ) + res = OptimizationProblemOutputs(ps_model) + compare_outage_power_and_deployed_reserves( + sys, + res, + reserve_up) + end +end + +@testset "G-n with contingency reserves deliverability constraints including responding reserves only up, reserve requirement, and reduction of parallel circuits" begin + for add_parallel_line in [true, false] + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + + if add_parallel_line + l4 = get_component(Line, c_sys5, "4") + add_equivalent_ac_transmission_with_parallel_circuits!(c_sys5, l4, PSY.Line) + end + systems = [c_sys5] + objfuncs = [GAEVF, GQEVF, GQEVF] + constraint_keys = [ + IOM.ConstraintKey( + ActivePowerVariableLimitsConstraint, + PSY.ThermalStandard, + "lb", + ), + IOM.ConstraintKey( + ActivePowerVariableLimitsConstraint, + PSY.ThermalStandard, + "ub", + ), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "lb"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "ub"), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_ub", + ), + IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.System), + #IOM.ConstraintKey(NetworkFlowConstraint, PSY.Line), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + ] + PTDF_ref = IdDict{System, PTDF}( + c_sys5 => PTDF(c_sys5), + ) + test_results = IdDict{System, Vector{Int}}( + c_sys5 => [360, 0, 600, 432, 72], + ) + test_obj_values = IdDict{System, Float64}( + c_sys5 => 329000.0, + ) + components_outages_cases = IdDict{System, Vector{String}}( + c_sys5 => ["Alta"], + ) + for (ix, sys) in enumerate(systems) + gen = get_component(ThermalStandard, sys, "Solitude") + set_ramp_limits!(gen, (up = 0.4 * PSY.SU, down = 0.4 * PSY.SU)) #Increase ramp limits to make the problem feasible + reserve_up = get_component(VariableReserve{ReserveUp}, sys, "Reserve1") + + components_outages_names = components_outages_cases[sys] + for component_name in components_outages_names + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + template = get_thermal_dispatch_template_network( + NetworkModel(PTDFPowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1", + )) + + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results[sys]..., + false, + ) + psi_checkobjfun_test(ps_model, objfuncs[ix]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_values[sys], + 10000, + ) + res = OptimizationProblemOutputs(ps_model) + compare_outage_power_and_deployed_reserves( + sys, + res, + reserve_up) + end + end +end + +@testset "G-n with contingency reserves deliverability constraints including responding reserves only up, NO reserve requirement, and reduction of parallel circuits" begin + for add_parallel_line in [true, false] + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + + if add_parallel_line + l4 = get_component(Line, c_sys5, "4") + add_equivalent_ac_transmission_with_parallel_circuits!(c_sys5, l4, PSY.Line) + end + systems = [c_sys5] + objfuncs = [GAEVF, GQEVF, GQEVF] + constraint_keys = [ + IOM.ConstraintKey( + ActivePowerVariableLimitsConstraint, + PSY.ThermalStandard, + "lb", + ), + IOM.ConstraintKey( + ActivePowerVariableLimitsConstraint, + PSY.ThermalStandard, + "ub", + ), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "lb"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "ub"), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_ub", + ), + IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.System), + #IOM.ConstraintKey(NetworkFlowConstraint, PSY.Line), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_lb", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_ub", + ), + ] + PTDF_ref = IdDict{System, PTDF}( + c_sys5 => PTDF(c_sys5), + ) + test_results = IdDict{System, Vector{Int}}( + c_sys5 => [240, 0, 504, 504, 96], + ) + test_obj_values = IdDict{System, Float64}( + c_sys5 => 329000.0, + ) + components_outages_cases = IdDict{System, Vector{String}}( + c_sys5 => ["Alta"], + ) + for (ix, sys) in enumerate(systems) + gen = get_component(ThermalStandard, sys, "Solitude") + set_ramp_limits!(gen, (up = 0.4 * PSY.SU, down = 0.4 * PSY.SU)) #Increase ramp limits to make the problem feasible + reserve_up = get_component(VariableReserve{ReserveUp}, sys, "Reserve1") + remove_time_series!( + sys, + Deterministic, + reserve_up, + "requirement", + ) + components_outages_names = components_outages_cases[sys] + for component_name in components_outages_names + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + template = get_thermal_dispatch_template_network( + NetworkModel(PTDFPowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1", + )) + + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results[sys]..., + false, + ) + psi_checkobjfun_test(ps_model, objfuncs[ix]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_values[sys], + 10000, + ) + res = OptimizationProblemOutputs(ps_model) + compare_outage_power_and_deployed_reserves( + sys, + res, + reserve_up) + end + end +end + +#This test ensures that the security constrained models build even when there are devices without set_device_model!() +@testset "Test if G-n with Ramp reserve deliverability constraints builds when there is a device without set_device_model!()" begin + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + + l4 = get_component(Line, c_sys5, "4") + add_equivalent_ac_transmission_with_parallel_circuits!( + c_sys5, + l4, + PSY.Line, + PSY.MonitoredLine, + ) + remove_component!(c_sys5, l4) + + systems = [c_sys5] + + PTDF_ref = IdDict{System, PTDF}( + c_sys5 => PTDF(c_sys5), + ) + + components_outages_cases = IdDict{System, Vector{String}}( + c_sys5 => ["Alta"], + ) + for (ix, sys) in enumerate(systems) + components_outages_names = components_outages_cases[sys] + for component_name in components_outages_names + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, sys, "Reserve1") + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + + template = + PowerOperationsProblemTemplate( + NetworkModel(PTDFPowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + set_device_model!(template, ThermalStandard, ThermalBasicDispatch) + set_device_model!(template, PowerLoad, StaticPowerLoad) + #set_device_model!(template, MonitoredLine, StaticBranchBounds) + set_device_model!(template, Line, StaticBranch) + set_device_model!(template, Transformer2W, StaticBranch) + set_device_model!(template, TapTransformer, StaticBranch) + set_device_model!(template, TwoTerminalGenericHVDCLine, HVDCTwoTerminalLossless) + + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1", + )) + + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + end +end +@testset "Test SecurityConstrainedContingencyReserve with different BranchFormulations" begin + for line_formulation in [StaticBranch, StaticBranchUnbounded, StaticBranchBounds] + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + l4 = get_component(Line, c_sys5, "4") + add_equivalent_ac_transmission_with_parallel_circuits!( + c_sys5, + l4, + PSY.Line, + PSY.MonitoredLine, + ) + remove_component!(c_sys5, l4) + + systems = [c_sys5] + + PTDF_ref = IdDict{System, PTDF}( + c_sys5 => PTDF(c_sys5), + ) + + components_outages_cases = IdDict{System, Vector{String}}( + c_sys5 => ["Alta"], + ) + for (ix, sys) in enumerate(systems) + components_outages_names = components_outages_cases[sys] + for component_name in components_outages_names + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, sys, "Reserve1") + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + + template = + PowerOperationsProblemTemplate( + NetworkModel(PTDFPowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + set_device_model!(template, ThermalStandard, ThermalBasicDispatch) + set_device_model!(template, PowerLoad, StaticPowerLoad) + #set_device_model!(template, MonitoredLine, StaticBranchBounds) + set_device_model!(template, Line, line_formulation) + set_device_model!(template, Transformer2W, StaticBranch) + set_device_model!(template, TapTransformer, StaticBranch) + set_device_model!(template, TwoTerminalGenericHVDCLine, HVDCTwoTerminalLossless) + + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1", + )) + + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + constraints = ps_model.internal.container.constraints + flow_rate_cons = constraints[IOM.ConstraintKey{ + PostContingencyFlowRateConstraint, + VariableReserve{ReserveUp}, + }( + "Reserve1_lb", + )] + @test length(flow_rate_cons) == 1 * 5 * 24 + end + end +end + +@testset "Test if G-n with Contingency reserve deliverability constraints builds when there is a device without set_device_model!()" begin + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + + l4 = get_component(Line, c_sys5, "4") + add_equivalent_ac_transmission_with_parallel_circuits!( + c_sys5, + l4, + PSY.Line, + PSY.MonitoredLine, + ) + remove_component!(c_sys5, l4) + + systems = [c_sys5] + + PTDF_ref = IdDict{System, PTDF}( + c_sys5 => PTDF(c_sys5), + ) + + components_outages_cases = IdDict{System, Vector{String}}( + c_sys5 => ["Alta"], + ) + for (ix, sys) in enumerate(systems) + components_outages_names = components_outages_cases[sys] + for component_name in components_outages_names + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, sys, "Reserve1") + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + + template = + PowerOperationsProblemTemplate( + NetworkModel(PTDFPowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + set_device_model!(template, ThermalStandard, ThermalBasicDispatch) + set_device_model!(template, PowerLoad, StaticPowerLoad) + #set_device_model!(template, MonitoredLine, StaticBranchBounds) + set_device_model!(template, Line, StaticBranch) + set_device_model!(template, Transformer2W, StaticBranch) + set_device_model!(template, TapTransformer, StaticBranch) + set_device_model!(template, TwoTerminalGenericHVDCLine, HVDCTwoTerminalLossless) + + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1", + )) + + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + end +end + +@testset "G-n with Ramp reserve deliverability constraints UC allowing 2 reserve products to respond" begin + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + systems = [c_sys5] + objfuncs = [GAEVF, GQEVF, GQEVF] + constraint_keys = [ + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "lb"), + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "ub"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "lb"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "ub"), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_ub", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve11_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve11_ub", + ), + IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.System), + #IOM.ConstraintKey(NetworkFlowConstraint, PSY.Line), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve11", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve11", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve11", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve11", + ), + ] + PTDF_ref = IdDict{System, PTDF}( + c_sys5 => PTDF(c_sys5), + ) + test_results = IdDict{System, Vector{Int}}( + c_sys5 => [960, 0, 1296, 600, 240], + ) + test_obj_values = IdDict{System, Float64}( + c_sys5 => 254242.0, + ) + components_outages_cases = IdDict{System, Vector{String}}( + c_sys5 => ["Alta"], + ) + for (ix, sys) in enumerate(systems) + components_outages_names = components_outages_cases[sys] + for component_name in components_outages_names + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, sys, "Reserve1") + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + reserve_up2 = get_component(VariableReserve{ReserveUp}, sys, "Reserve11") + add_supplemental_attribute!(sys, reserve_up2, transition_data) + end + + template = get_thermal_dispatch_template_network( + NetworkModel(PTDFPowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + + set_device_model!( + template, + ThermalStandard, + ThermalStandardUnitCommitment, + ) + + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1", + )) + + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve11", + )) + + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results[sys]..., + true, + ) + psi_checkobjfun_test(ps_model, objfuncs[ix]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_values[sys], + 10000, + ) + end +end + +@testset "G-n with Ramp reserve deliverability constraints with AreaPTDFPowerModel w/wo Reserve Slacks" begin + reserve_slacks = [false, true] + objfuncs = [GAEVF] + constraint_keys = [ + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "lb"), + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "ub"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "lb"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "ub"), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_ub", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_ub", + ), + IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.Area), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + ] + test_results = IdDict{Bool, Vector{Int}}( + reserve_slacks[1] => [744, 0, 1536, 1200, 168], + reserve_slacks[2] => [2040, 0, 1536, 1200, 168], + ) + test_obj_values = IdDict{Bool, Float64}( + reserve_slacks[1] => 497000.0, + reserve_slacks[2] => 497000.0, + ) + components_outages_cases = (["Alta_1", "Alta_2"], ["Reserve1_1", "Reserve1_2"]) + + for reserve_slack in reserve_slacks + sys = PSB.build_system(PSISystems, "two_area_pjm_DA"; add_reserves = true) + transform_single_time_series!(sys, Hour(24), Hour(1)) + + components_outages_names, reserve_names = components_outages_cases + for (component_name, reserve_name) in + zip(components_outages_names, reserve_names) + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + + template = get_thermal_dispatch_template_network( + NetworkModel(AreaPTDFPowerModel; PTDF_matrix = PTDF(sys)), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1_1"; + use_slacks = reserve_slack, + )) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1_2"; + use_slacks = reserve_slack, + )) + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results[reserve_slack]..., + false, + ) + psi_checkobjfun_test(ps_model, objfuncs[1]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_values[reserve_slack], + 10000, + ) + res = OptimizationProblemOutputs(ps_model) + for reserve_name in reserve_names + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + compare_outage_power_and_deployed_reserves( + sys, + res, + reserve_up) + end + end +end + +# Exercises the per-service line-scoping path in +# `_monitored_components_by_modeled_type` for the `AreaPTDFPowerModel` +# network model. Each reserve service monitors a different hand-picked +# subset of AC lines, which keeps the constraint key meta keyed by service +# name but reduces the number of `PostContingencyFlowRateConstraint` rows +# compared to the all-lines baseline. +@testset "G-n with Ramp reserve deliverability constraints with AreaPTDFPowerModel and monitored line subset" begin + objfuncs = [GAEVF] + monitored_line_names_per_service = (["1_1", "2_1"], ["1_2", "2_2"]) + constraint_keys = [ + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "lb"), + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "ub"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "lb"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "ub"), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_ub", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_ub", + ), + IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.Area), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + ] + # Counts are smaller than the all-lines baseline `[744, 0, 1536, 1200, 168]` + # because each service monitors only two AC lines instead of all 13. + test_results = [744, 0, 1008, 672, 168] + test_obj_value = 497000.0 + components_outages_cases = (["Alta_1", "Alta_2"], ["Reserve1_1", "Reserve1_2"]) + + sys = PSB.build_system(PSISystems, "two_area_pjm_DA"; add_reserves = true) + transform_single_time_series!(sys, Hour(24), Hour(1)) + + components_outages_names, reserve_names = components_outages_cases + for (component_name, reserve_name, monitored_names) in zip( + components_outages_names, + reserve_names, + monitored_line_names_per_service, + ) + monitored_subset = [get_component(Line, sys, n) for n in monitored_names] + # --- Create Outage Data with a hand-picked monitored subset --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = monitored_subset, + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + + template = get_thermal_dispatch_template_network( + NetworkModel(AreaPTDFPowerModel; PTDF_matrix = PTDF(sys)), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1_1", + )) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1_2", + )) + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results..., + false, + ) + psi_checkobjfun_test(ps_model, objfuncs[1]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_value, + 10000, + ) + res = OptimizationProblemOutputs(ps_model) + for reserve_name in reserve_names + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + compare_outage_power_and_deployed_reserves( + sys, + res, + reserve_up) + end +end + +@testset "G-n with Contingency reserve deliverability constraints with AreaPTDFPowerModel, reserves only up, reserve requirement" begin + reserve_slacks = [false, true] + objfuncs = [GAEVF] + constraint_keys = [ + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "lb"), + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "ub"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "lb"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "ub"), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_ub", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_ub", + ), + IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.Area), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + ] + test_results = IdDict{Bool, Vector{Int}}( + reserve_slacks[1] => [744, 0, 1536, 1200, 168], + reserve_slacks[2] => [2040, 0, 1536, 1200, 168], + ) + test_obj_values = IdDict{Bool, Float64}( + reserve_slacks[1] => 497000.0, + reserve_slacks[2] => 497000.0, + ) + components_outages_cases = (["Alta_1", "Alta_2"], ["Reserve1_1", "Reserve1_2"]) + + for reserve_slack in reserve_slacks + sys = PSB.build_system(PSISystems, "two_area_pjm_DA"; add_reserves = true) + transform_single_time_series!(sys, Hour(24), Hour(1)) + components_outages_names, reserve_names = components_outages_cases + for (component_name, reserve_name) in zip(components_outages_names, reserve_names) + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + + template = get_thermal_dispatch_template_network( + NetworkModel(AreaPTDFPowerModel; PTDF_matrix = PTDF(sys)), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1_1"; + use_slacks = reserve_slack, + )) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1_2"; + use_slacks = reserve_slack, + )) + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results[reserve_slack]..., + false, + ) + psi_checkobjfun_test(ps_model, objfuncs[1]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_values[reserve_slack], + 10000, + ) + res = OptimizationProblemOutputs(ps_model) + for reserve_name in reserve_names + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + compare_outage_power_and_deployed_reserves( + sys, + res, + reserve_up) + end + end +end + +@testset "G-n with Contingency reserve deliverability constraints with AreaPTDFPowerModel, reserves only up, NO reserve requirement" begin + c_sys5_2area = PSB.build_system(PSISystems, "two_area_pjm_DA") + transform_single_time_series!(c_sys5_2area, Hour(24), Hour(1)) + systems = [c_sys5_2area] + objfuncs = [GAEVF] + constraint_keys = [ + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "lb"), + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "ub"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "lb"), + IOM.ConstraintKey(FlowRateConstraint, PSY.Line, "ub"), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_ub", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_ub", + ), + IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.Area), + #IOM.ConstraintKey(NetworkFlowConstraint, PSY.Line), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_lb", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_ub", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_lb", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_ub", + ), + ] + PTDF_ref = IdDict{System, PTDF}( + c_sys5_2area => PTDF(c_sys5_2area), + ) + test_results = IdDict{System, Vector{Int}}( + c_sys5_2area => [504, 0, 1344, 1344, 216], + ) + test_obj_values = IdDict{System, Float64}( + c_sys5_2area => 497000.0, + ) + components_outages_cases = IdDict{System, Tuple{Vector{String}, Vector{String}}}( + c_sys5_2area => (["Alta_1", "Alta_2"], ["Reserve1_1", "Reserve1_2"]), + ) + for (ix, sys) in enumerate(systems) + components_outages_names, reserve_names = components_outages_cases[sys] + contributing_devices = get_components( + g -> get_name(get_area(get_bus(g))) == "Area1", + ThermalStandard, + sys, + ) + add_reserve_product_without_requirement_time_series!( + sys, + "Reserve1_1", + "Up", + contributing_devices, + ) + contributing_devices = get_components( + g -> get_name(get_area(get_bus(g))) == "Area2", + ThermalStandard, + sys, + ) + add_reserve_product_without_requirement_time_series!( + sys, + "Reserve1_2", + "Up", + contributing_devices, + ) + + for (component_name, reserve_name) in zip(components_outages_names, reserve_names) + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + + template = get_thermal_dispatch_template_network( + NetworkModel(AreaPTDFPowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1_1", + )) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1_2", + )) + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results[sys]..., + false, + ) + psi_checkobjfun_test(ps_model, objfuncs[ix]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_values[sys], + 10000, + ) + res = OptimizationProblemOutputs(ps_model) + for reserve_name in reserve_names + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + compare_outage_power_and_deployed_reserves( + sys, + res, + reserve_up) + end + end +end + +@testset "G-n with Ramp reserve deliverability constraints with CopperPlatePowerModel" begin + c_sys5_2area = PSB.build_system(PSISystems, "two_area_pjm_DA"; add_reserves = true) + transform_single_time_series!(c_sys5_2area, Hour(24), Hour(1)) + systems = [c_sys5_2area] + objfuncs = [GAEVF] + constraint_keys = [ + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "lb"), + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "ub"), IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.System), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + ] + PTDF_ref = IdDict{System, PTDF}( + c_sys5_2area => PTDF(c_sys5_2area), + ) + test_results = IdDict{System, Vector{Int}}( + c_sys5_2area => [720, 0, 624, 288, 120], + ) + test_obj_values = IdDict{System, Float64}( + c_sys5_2area => 497494.48, + ) + components_outages_cases = IdDict{System, Tuple{Vector{String}, Vector{String}}}( + c_sys5_2area => (["Alta_1", "Alta_2"], ["Reserve1_1", "Reserve1_2"]), + ) + for (ix, sys) in enumerate(systems) + components_outages_names, reserve_names = components_outages_cases[sys] + for (component_name, reserve_name) in zip(components_outages_names, reserve_names) + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + + template = get_thermal_dispatch_template_network( + NetworkModel(CopperPlatePowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1_1", + )) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1_2", + )) + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results[sys]..., + false, + ) + psi_checkobjfun_test(ps_model, objfuncs[ix]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_values[sys], + 10000, + ) + res = OptimizationProblemOutputs(ps_model) + for reserve_name in reserve_names + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + compare_outage_power_and_deployed_reserves( + sys, + res, + reserve_up) + end + end +end + +@testset "G-n with Contingency reserve deliverability constraints with CopperPlatePowerModel with Reserve Requirement" begin + c_sys5_2area = PSB.build_system(PSISystems, "two_area_pjm_DA"; add_reserves = true) + transform_single_time_series!(c_sys5_2area, Hour(24), Hour(1)) + systems = [c_sys5_2area] + objfuncs = [GAEVF] + constraint_keys = [ + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "lb"), + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "ub"), IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.System), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + ] + PTDF_ref = IdDict{System, PTDF}( + c_sys5_2area => PTDF(c_sys5_2area), + ) + test_results = IdDict{System, Vector{Int}}( + c_sys5_2area => [720, 0, 624, 288, 120], + ) + test_obj_values = IdDict{System, Float64}( + c_sys5_2area => 497494.48, + ) + components_outages_cases = IdDict{System, Tuple{Vector{String}, Vector{String}}}( + c_sys5_2area => (["Alta_1", "Alta_2"], ["Reserve1_1", "Reserve1_2"]), + ) + for (ix, sys) in enumerate(systems) + components_outages_names, reserve_names = components_outages_cases[sys] + for (component_name, reserve_name) in zip(components_outages_names, reserve_names) + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + + template = get_thermal_dispatch_template_network( + NetworkModel(CopperPlatePowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1_1", + )) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1_2", + )) + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results[sys]..., + false, + ) + psi_checkobjfun_test(ps_model, objfuncs[ix]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_values[sys], + 10000, + ) + res = OptimizationProblemOutputs(ps_model) + for reserve_name in reserve_names + reserve_up = get_component(VariableReserve{ReserveUp}, sys, reserve_name) + compare_outage_power_and_deployed_reserves( + sys, + res, + reserve_up) + end + end +end + +@testset "G-n with Contingency reserve deliverability constraints with CopperPlatePowerModel with NO Reserve Requirement" begin + c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + + systems = [c_sys5] + objfuncs = [GAEVF] + constraint_keys = [ + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "lb"), + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "ub"), IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.System), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_lb", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_ub", + )] + PTDF_ref = IdDict{System, PTDF}( + c_sys5 => PTDF(c_sys5), + ) + test_results = IdDict{System, Vector{Int}}( + c_sys5 => [240, 0, 216, 216, 96], + ) + test_obj_values = IdDict{System, Float64}( + c_sys5 => 329000.0, + ) + components_outages_cases = IdDict{System, Vector{String}}( + c_sys5 => ["Alta"], + ) + for (ix, sys) in enumerate(systems) + reserve_up = get_component(VariableReserve{ReserveUp}, sys, "Reserve1") + remove_time_series!( + sys, + Deterministic, + reserve_up, + "requirement", + ) + + components_outages_names = components_outages_cases[sys] + for component_name in components_outages_names + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, sys, component_name) + add_supplemental_attribute!(sys, component, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + end + + template = get_thermal_dispatch_template_network( + NetworkModel(CopperPlatePowerModel; PTDF_matrix = PTDF_ref[sys]), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1", + )) + + ps_model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + psi_constraint_test(ps_model, constraint_keys) + moi_tests( + ps_model, + test_results[sys]..., + false, + ) + psi_checkobjfun_test(ps_model, objfuncs[ix]) + psi_checksolve_test( + ps_model, + [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL], + test_obj_values[sys], + 10000, + ) + res = OptimizationProblemOutputs(ps_model) + compare_outage_power_and_deployed_reserves( + sys, + res, + reserve_up) + end +end + +@testset "G-n with Ramp reserve deliverability constraints with AreaBalance PowerModel" begin + constraint_keys = [ + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "lb"), + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "ub"), IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.Area), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyCopperPlateBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyCopperPlateBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_ub", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_ub", + ), + ] + + c_sys = PSB.build_system(PSISystems, "two_area_pjm_DA"; add_reserves = true) + transform_single_time_series!(c_sys, Hour(24), Hour(1)) + components_outages_names, reserve_names = + (["Alta_1", "Alta_2"], ["Reserve1_1", "Reserve1_2"]) + + for (component_name, reserve_name) in zip(components_outages_names, reserve_names) + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = vcat( + collect(get_components(ACTransmission, c_sys)), + collect(get_components(AreaInterchange, c_sys)), + ), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, c_sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, c_sys, reserve_name) + add_supplemental_attribute!(c_sys, component, transition_data) + add_supplemental_attribute!(c_sys, reserve_up, transition_data) + end + + template = get_thermal_dispatch_template_network(NetworkModel(AreaBalancePowerModel)) + set_device_model!(template, AreaInterchange, StaticBranch) + + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1_1", + )) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1_2", + )) + + ps_model = + DecisionModel(template, c_sys; resolution = Hour(1), optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + + psi_constraint_test(ps_model, constraint_keys) + + @test solve!(ps_model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + moi_tests(ps_model, 744, 0, 696, 360, 240, false) + + opt_container = IOM.get_optimization_container(ps_model) + copper_plate_constraints = + IOM.get_constraint(opt_container, CopperPlateBalanceConstraint(), PSY.Area) + @test size(copper_plate_constraints) == (2, 24) + + # Re-recorded under psy6: structure matches PSI exactly (moi_tests + constraint + # keys above pass), only the `two_area_pjm_DA` system data drifted across the + # PowerSystems major version, shifting the optimal cost. + psi_checksolve_test(ps_model, [MOI.OPTIMAL], 493287.9128996057, 1) + + results = OptimizationProblemOutputs(ps_model) + interarea_flow = read_variable( + results, + "FlowActivePowerVariable__AreaInterchange"; + table_format = TableFormat.WIDE, + ) + # The values for these tests come from the data + @test all(interarea_flow[!, "1_2"] .<= 150 + 1e-6) + @test all(interarea_flow[!, "1_2"] .>= -150 - 1e-6) + + load = read_parameter( + results, + "ActivePowerTimeSeriesParameter__PowerLoad"; + table_format = TableFormat.WIDE, + ) + thermal_gen = read_variable( + results, + "ActivePowerVariable__ThermalStandard"; + table_format = TableFormat.WIDE, + ) + + zone_1_load = sum(eachcol(load[!, ["Bus4_1", "Bus3_1", "Bus2_1"]])) + zone_1_gen = sum( + eachcol( + thermal_gen[ + !, + ["Solitude_1", "Park City_1", "Sundance_1", "Brighton_1", "Alta_1"], + ], + ), + ) + @test all( + isapprox.( + sum(zone_1_gen .+ zone_1_load .- interarea_flow[!, "1_2"]; dims = 2), + 0.0; + atol = 1e-3, + ), + ) + + zone_2_load = sum(eachcol(load[!, ["Bus4_2", "Bus3_2", "Bus2_2"]])) + zone_2_gen = sum( + eachcol( + thermal_gen[ + !, + ["Solitude_2", "Park City_2", "Sundance_2", "Brighton_2", "Alta_2"], + ], + ), + ) + @test all( + isapprox.( + sum(zone_2_gen .+ zone_2_load .+ interarea_flow[!, "1_2"]; dims = 2), + 0.0; + atol = 1e-3, + ), + ) + + res = OptimizationProblemOutputs(ps_model) + for reserve_name in reserve_names + reserve_up = get_component(VariableReserve{ReserveUp}, c_sys, reserve_name) + compare_outage_power_and_deployed_reserves( + c_sys, + res, + reserve_up) + end +end + +@testset "G-n with Contingency reserve deliverability constraints with AreaBalancePowerModel with Reserve Requirement" begin + constraint_keys = [ + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "lb"), + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "ub"), IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.Area), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RequirementConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + RampConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyActivePowerReserveDeploymentVariableLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyCopperPlateBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyCopperPlateBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_ub", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_ub", + ), + ] + + c_sys = PSB.build_system(PSISystems, "two_area_pjm_DA"; add_reserves = true) + transform_single_time_series!(c_sys, Hour(24), Hour(1)) + components_outages_names, reserve_names = + (["Alta_1", "Alta_2"], ["Reserve1_1", "Reserve1_2"]) + + for (component_name, reserve_name) in zip(components_outages_names, reserve_names) + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = vcat( + collect(get_components(ACTransmission, c_sys)), + collect(get_components(AreaInterchange, c_sys)), + ), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, c_sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, c_sys, reserve_name) + add_supplemental_attribute!(c_sys, component, transition_data) + add_supplemental_attribute!(c_sys, reserve_up, transition_data) + end + + template = get_thermal_dispatch_template_network(NetworkModel(AreaBalancePowerModel)) + set_device_model!(template, AreaInterchange, StaticBranch) + + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1_1", + )) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1_2", + )) + + ps_model = + DecisionModel(template, c_sys; resolution = Hour(1), optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + + psi_constraint_test(ps_model, constraint_keys) + + @test solve!(ps_model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + moi_tests(ps_model, 744, 0, 696, 360, 240, false) + + opt_container = IOM.get_optimization_container(ps_model) + copper_plate_constraints = + IOM.get_constraint(opt_container, CopperPlateBalanceConstraint(), PSY.Area) + @test size(copper_plate_constraints) == (2, 24) + + # Re-recorded under psy6: structure matches PSI exactly (moi_tests + constraint + # keys above pass), only the `two_area_pjm_DA` system data drifted across the + # PowerSystems major version, shifting the optimal cost. + psi_checksolve_test(ps_model, [MOI.OPTIMAL], 493287.9128996057, 1) + + results = OptimizationProblemOutputs(ps_model) + interarea_flow = read_variable( + results, + "FlowActivePowerVariable__AreaInterchange"; + table_format = TableFormat.WIDE, + ) + # The values for these tests come from the data + @test all(interarea_flow[!, "1_2"] .<= 150 + 1e-6) + @test all(interarea_flow[!, "1_2"] .>= -150 - 1e-6) + + load = read_parameter( + results, + "ActivePowerTimeSeriesParameter__PowerLoad"; + table_format = TableFormat.WIDE, + ) + thermal_gen = read_variable( + results, + "ActivePowerVariable__ThermalStandard"; + table_format = TableFormat.WIDE, + ) + + zone_1_load = sum(eachcol(load[!, ["Bus4_1", "Bus3_1", "Bus2_1"]])) + zone_1_gen = sum( + eachcol( + thermal_gen[ + !, + ["Solitude_1", "Park City_1", "Sundance_1", "Brighton_1", "Alta_1"], + ], + ), + ) + @test all( + isapprox.( + sum(zone_1_gen .+ zone_1_load .- interarea_flow[!, "1_2"]; dims = 2), + 0.0; + atol = 1e-3, + ), + ) + + zone_2_load = sum(eachcol(load[!, ["Bus4_2", "Bus3_2", "Bus2_2"]])) + zone_2_gen = sum( + eachcol( + thermal_gen[ + !, + ["Solitude_2", "Park City_2", "Sundance_2", "Brighton_2", "Alta_2"], + ], + ), + ) + @test all( + isapprox.( + sum(zone_2_gen .+ zone_2_load .+ interarea_flow[!, "1_2"]; dims = 2), + 0.0; + atol = 1e-3, + ), + ) + + res = OptimizationProblemOutputs(ps_model) + for reserve_name in reserve_names + reserve_up = get_component(VariableReserve{ReserveUp}, c_sys, reserve_name) + compare_outage_power_and_deployed_reserves( + c_sys, + res, + reserve_up) + end +end + +@testset "G-n with Contingency reserve deliverability constraints with AreaBalancePowerModel with NO Reserve Requirement" begin + constraint_keys = [ + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "lb"), + IOM.ConstraintKey(ActivePowerVariableLimitsConstraint, PSY.ThermalStandard, "ub"), IOM.ConstraintKey(CopperPlateBalanceConstraint, PSY.Area), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_lb", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_ub", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_lb", + ), + IOM.ConstraintKey( + PostContingencyActivePowerGenerationLimitsConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_ub", + ), + IOM.ConstraintKey( + PostContingencyCopperPlateBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + IOM.ConstraintKey( + PostContingencyCopperPlateBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1_ub", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_lb", + ), + IOM.ConstraintKey( + PostContingencyFlowRateConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2_ub", + ), + ] + + c_sys = PSB.build_system(PSISystems, "two_area_pjm_DA"; add_reserves = true) + + reserve_up = get_component(VariableReserve{ReserveUp}, c_sys, "Reserve1_1") + remove_time_series!( + c_sys, + SingleTimeSeries, + reserve_up, + "requirement", + ) + reserve_up = get_component(VariableReserve{ReserveUp}, c_sys, "Reserve1_2") + remove_time_series!( + c_sys, + SingleTimeSeries, + reserve_up, + "requirement", + ) + + transform_single_time_series!(c_sys, Hour(24), Hour(1)) + components_outages_names, reserve_names = + (["Alta_1", "Alta_2"], ["Reserve1_1", "Reserve1_2"]) + + for (component_name, reserve_name) in zip(components_outages_names, reserve_names) + # --- Create Outage Data --- + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = vcat( + collect(get_components(ACTransmission, c_sys)), + collect(get_components(AreaInterchange, c_sys)), + ), + ) + # --- Add Outage Supplemental attribute to device and services that should respond --- + component = get_component(ThermalStandard, c_sys, component_name) + reserve_up = get_component(VariableReserve{ReserveUp}, c_sys, reserve_name) + add_supplemental_attribute!(c_sys, component, transition_data) + add_supplemental_attribute!(c_sys, reserve_up, transition_data) + end + + template = get_thermal_dispatch_template_network(NetworkModel(AreaBalancePowerModel)) + set_device_model!(template, AreaInterchange, StaticBranch) + + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1_1", + )) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedContingencyReserve, + "Reserve1_2", + )) + + ps_model = + DecisionModel(template, c_sys; resolution = Hour(1), optimizer = HiGHS_optimizer) + + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + + psi_constraint_test(ps_model, constraint_keys) + + @test solve!(ps_model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + + moi_tests(ps_model, 504, 0, 504, 504, 288, false) + + opt_container = IOM.get_optimization_container(ps_model) + copper_plate_constraints = + IOM.get_constraint(opt_container, CopperPlateBalanceConstraint(), PSY.Area) + @test size(copper_plate_constraints) == (2, 24) + + psi_checksolve_test(ps_model, [MOI.OPTIMAL], 482055.7647083302, 1) + + results = OptimizationProblemOutputs(ps_model) + interarea_flow = read_variable( + results, + "FlowActivePowerVariable__AreaInterchange"; + table_format = TableFormat.WIDE, + ) + # The values for these tests come from the data + @test all(interarea_flow[!, "1_2"] .<= 150) + @test all(interarea_flow[!, "1_2"] .>= -150) + + load = read_parameter( + results, + "ActivePowerTimeSeriesParameter__PowerLoad"; + table_format = TableFormat.WIDE, + ) + thermal_gen = read_variable( + results, + "ActivePowerVariable__ThermalStandard"; + table_format = TableFormat.WIDE, + ) + + zone_1_load = sum(eachcol(load[!, ["Bus4_1", "Bus3_1", "Bus2_1"]])) + zone_1_gen = sum( + eachcol( + thermal_gen[ + !, + ["Solitude_1", "Park City_1", "Sundance_1", "Brighton_1", "Alta_1"], + ], + ), + ) + @test all( + isapprox.( + sum(zone_1_gen .+ zone_1_load .- interarea_flow[!, "1_2"]; dims = 2), + 0.0; + atol = 1e-3, + ), + ) + + zone_2_load = sum(eachcol(load[!, ["Bus4_2", "Bus3_2", "Bus2_2"]])) + zone_2_gen = sum( + eachcol( + thermal_gen[ + !, + ["Solitude_2", "Park City_2", "Sundance_2", "Brighton_2", "Alta_2"], + ], + ), + ) + @test all( + isapprox.( + sum(zone_2_gen .+ zone_2_load .+ interarea_flow[!, "1_2"]; dims = 2), + 0.0; + atol = 1e-3, + ), + ) + + res = OptimizationProblemOutputs(ps_model) + for reserve_name in reserve_names + reserve_up = get_component(VariableReserve{ReserveUp}, c_sys, reserve_name) + compare_outage_power_and_deployed_reserves( + c_sys, + res, + reserve_up) + end +end + +# Regression test for per-service outage scoping under the +# attachment-as-the-rule contract: a security-constrained reserve service +# responds to exactly the outages attached to it via +# `add_supplemental_attribute!(sys, service, outage)`. Generator attachment +# is required for the post-contingency build (so the outaged generator can +# be pinned to zero deployment), but it is the *service* attachment that +# selects which `ServiceModel` claims the outage. Membership in the +# service's contributing-devices set is irrelevant to the selection. +@testset "SC reserve outage attachment scopes responding services" begin + sys = PSB.build_system(PSISystems, "two_area_pjm_DA"; add_reserves = true) + transform_single_time_series!(sys, Hour(24), Hour(1)) + + reserve1 = get_component(VariableReserve{ReserveUp}, sys, "Reserve1_1") + + # Attach an UnplannedOutage to a single Area1 generator and to the + # reserve that should respond. Reserve1_2 is intentionally NOT attached. + alta1 = get_component(ThermalStandard, sys, "Alta_1") + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + add_supplemental_attribute!(sys, alta1, transition_data) + add_supplemental_attribute!(sys, reserve1, transition_data) + outage_uuid = IS.get_uuid(transition_data) + + template = get_thermal_dispatch_template_network( + NetworkModel(AreaPTDFPowerModel; PTDF_matrix = PTDF(sys)), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1_1", + )) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1_2", + )) + + # --- Unit-level: attachment scoping populates only the responding ServiceModel --- + POM._build_service_model_outages!(template, sys) + + services = IOM.get_service_models(template) + sm1 = services[("Reserve1_1", Symbol(VariableReserve{ReserveUp}))] + sm2 = services[("Reserve1_2", Symbol(VariableReserve{ReserveUp}))] + @test haskey(sm1.outages, outage_uuid) + @test !haskey(sm2.outages, outage_uuid) + + # --- Build-level: post-contingency constraints fire only on Reserve1_1 --- + ps_model = + DecisionModel(template, sys; resolution = Hour(1), optimizer = HiGHS_optimizer) + @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + + container = IOM.get_optimization_container(ps_model) + @test IOM.has_container_key( + container, + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ) + @test !IOM.has_container_key( + container, + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_2", + ) + cons_resp = IOM.get_constraint( + container, + IOM.ConstraintKey( + PostContingencyGenerationBalanceConstraint, + PSY.VariableReserve{ReserveUp}, + "Reserve1_1", + ), + ) + @test size(cons_resp) == (1, 24) +end + +# Regression test for the single-reserve case: when only one SC reserve +# service is in the template, an outage attached to both the outaged +# generator and the reserve must end up in that ServiceModel.outages dict. +@testset "SC reserve outage attachment covers single-reserve case" begin + sys = PSB.build_system(PSITestSystems, "c_sys5_uc"; add_reserves = true) + + transition_data = GeometricDistributionForcedOutage(; + mean_time_to_recovery = 10, + outage_transition_probability = 0.9999, + monitored_components = collect(get_components(ACTransmission, sys)), + ) + alta = get_component(ThermalStandard, sys, "Alta") + reserve_up = get_component(VariableReserve{ReserveUp}, sys, "Reserve1") + add_supplemental_attribute!(sys, alta, transition_data) + add_supplemental_attribute!(sys, reserve_up, transition_data) + outage_uuid = IS.get_uuid(transition_data) + + template = get_thermal_dispatch_template_network( + NetworkModel(PTDFPowerModel; PTDF_matrix = PTDF(sys)), + ) + set_service_model!(template, + ServiceModel( + VariableReserve{ReserveUp}, + SecurityConstrainedRampReserve, + "Reserve1", + )) + + POM._build_service_model_outages!(template, sys) + + services = IOM.get_service_models(template) + sm = services[("Reserve1", Symbol(VariableReserve{ReserveUp}))] + @test haskey(sm.outages, outage_uuid) +end