diff --git a/.claude/pom_port_plan.md b/.claude/pom_port_plan.md deleted file mode 100644 index 7fccf0ed..00000000 --- a/.claude/pom_port_plan.md +++ /dev/null @@ -1,183 +0,0 @@ -# POM port plan — PSI → PowerOperationsModels.jl - -Scope: bring POM (`PowerOperationsModels.jl`, the **formulation** library) up to date with PSI -(`PowerSimulations.jl`), covering (1) missing formulation **tests**, (2) post-fork PSI **code** -changes that belong in POM, and (3) the **in-progress G-1** generator-contingency feature. - -Companion file: `iom_port_plan.md` (generic optimization-core changes → InfrastructureOptimizationModels). - -## Topology / ground rules -- Sienna porting is two-headed: **formulation specifics → POM**, **generic optimization core - (optimization_container, dual_processing, settings, serialization, parameters, lifecycle, - objective-function machinery, model-export) → IOM**. Check the right repo before porting. -- POM is a **restructured** port (e.g. PF logic lives in `ext/PowerFlowsExt/`; type-based dispatch - `::Type{F}` vs PSI's instance-based). Verify by **symbol/behavior**, not file path; adapt, don't - copy verbatim. -- Fork baseline ≈ PSI #1503 (POM/IOM forked 2026-01-03). Latest PSI PR swept = **#1640**. -- "Already ported" / "code present" claims below were symbol-verified against the local POM clone. - ---- - -## Workstream G1 — in-progress generator contingency (PSI **open PR #1617**) — TIME-SENSITIVE -PR #1617 (`sm/g-1_monitored_c`, OPEN, not draft, `mergeable_state=blocked` on review; author -SebastianManriqueM) adds the **reserve/service-side** security-constrained contingency layer. -POM already carries the **branch-side** foundation (`src/ac_transmission_models/security_constrained_branch.jl`, -outage population in `src/operation/template_validation.jl`, `SecurityConstrainedStaticBranch`, -`PostContingencyBranchFlow/FlowRateConstraint`, slack vars). So the POM port ≈ the #1617 diff -re-applied onto POM's layout. The three older `g-1` branches are a stale lineage already merged to -PSI `main` (and into POM) — **do not** port from them. - -Port from PSI branch `origin/sm/g-1_monitored_c`: -- **Formulations** (`core/formulations.jl`): `AbstractSecurityConstrainedReservesFormulation`, - `SecurityConstrainedContingencyReserve`, `SecurityConstrainedRampReserve`. -- **Variables** (`core/variables.jl`): `AbstractContingencySlackVariableType`, - `PostContingencyActivePowerReserveDeploymentVariable`; re-parent - `PostContingencyFlowActivePowerSlackUpper/LowerBound` under the slack supertype. -- **Expressions** (`core/expressions.jl`): `PostContingencyActivePowerGeneration`, - `PostContingencyAreaInterchangeFlow`, `PostContingencyAreaActivePowerDeployment` - (+ `should_write_resulting_value`/`convert_result_to_natural_units`). -- **Constraints** (`core/constraints.jl`): `PostContingencyActivePowerGenerationLimitsConstraint`, - `PostContingencyCopperPlateBalanceConstraint`, `PostContingencyGenerationBalanceConstraint`, - `PostContingencyRampConstraint`. -- **Definitions**: `POST_CONTINGENCY_CONSTRAINT_VIOLATION_SLACK_COST = 1e5`. -- **ServiceModel** (`core/service_model.jl`): `outages::Dict{UUID,Dict{DataType,Set{String}}}` field - + `get_outages` (outages attached as `PSY.Outage` supplemental attributes on the `PSY.Service`). -- **Template validation** (`operation/template_validation.jl`): `_build_service_model_outages!`, - `_sc_reserve_service_models`, `_service_skips_outage` (dispatch on `PSY.PlannedOutage`/`PSY.Outage`), - `_warn_outages_attached_to_unmodeled_services`; **widen `_monitored_components_by_modeled_type` - to admit `PSY.AreaInterchange`** (POM's copy currently admits only `PSY.ACTransmission`). -- **Reserves** (`services_models/reserves.jl`): SC `get_default_time_series_names(::Reserve, ::AbstractSecurityConstrainedReservesFormulation)`. -- **New file**: `services_models/static_injection_security_constrained_models.jl` (sparse - post-contingency containers keyed by `(outage_id, monitored_name, t)`; `construct_service!` - paths for CopperPlate / AreaBalance / PTDF / AreaPTDF). Register in the POM module + exports. -- **New test**: `test/test_static_injection_security_constrained_models.jl` (~2,580 lines; both - formulations × CopperPlate/AreaBalance/AreaPTDF/PTDF; slack on/off; parallel-circuit reduction; - mixed branch formulations). - -> Track #1617 to merge; if it changes in review, re-diff before porting. - ---- - -## Workstream M — merged-PR code gaps in POM (symbol-verified ABSENT) -Bugfixes and small features that landed in PSI after the fork and are not yet in POM. Each line: PR — what — where in POM. - -**Bugfixes (do first; small, low-risk):** -- **#1519** — TwoTerminalHVDC `AreaBalancePowerModel` `add_to_expression!` bug: uses - `get_expression(container, T, PSY.ACBus)` (should be `PSY.Area`), undefined `network_reduction`/`W`. - → `src/common_models/add_to_expression.jl:954`. -- **#1527** — `round(Float64, time_limits.up*steps_per_hour, RoundUp)` should drop the `Float64` - positional arg (align with `RoundingMode`). → `src/static_injector_models/thermal_generation.jl:1364`. -- **#1535** — VOM cost missing resolution scaling (issue #1531): add - `dt = Dates.value(resolution)/MILLISECONDS_IN_HOUR` to `_add_vom_cost_to_objective_helper!`. - → `src/common_models/market_bid_plumbing.jl:555`. (Unblocks the VOM-normalization MBC tests.) -- **#1587** — improve parallel-branch interface error message (list offending branch names). - → `src/common_models/add_to_expression.jl:2256` (minor). -- **#1508** — network-reduction/AreaInterchange fix: `modeled_branch_types`→`modeled_ac_branch_types` - and the `device_names_with_branches` undefined-constraint bug (`device_names = PSY.get_name.(devices)`). - → `src/area_interchange.jl` + reduction path. Adapt to POM's `instantiate_network_model.jl`. - -**Small features:** -- **#1549** — Source default TS names: `get_default_time_series_names(::Type{<:PSY.Source}, …)` - returns `Dict()` in POM; should map `ActivePowerOut/InTimeSeriesParameter => "max_active_power_out/in"`. - → `src/static_injector_models/source.jl`. -- **#1573** — Source `FixedOutput` constructor (all POM Source constructors are `ImportExportSourceModel`). - → `src/static_injector_models/source_constructor.jl`. (Unblocks "FixedOutput Source w/ PTDF+TS" test.) -- **#1538** — feedforward arguments for renewables: add `add_feedforward_arguments!` and - `get_parameter_multiplier(::VariableValueParameter, ::PSY.RenewableGen, ::AbstractRenewableFormulation)`. - → `src/static_injector_models/renewablegeneration_constructor.jl`, `renewable_generation.jl`. -- **#1614** — renewable curtailment incentive: flip `objective_function_multiplier(::AbstractRenewableDispatchFormulation)` - from `OBJECTIVE_FUNCTION_NEGATIVE` to `POSITIVE` and feed `curtailment_cost` into the dispatch - objective (not just the reporting `CurtailmentCostExpression`). → `src/static_injector_models/renewable_generation.jl:24`. - (Unblocks "curtailment_cost incentive affects dispatch" test.) -- **#1605** — AreaPTDF + `InterconnectingConverter`: replace - `error("AreaPTDFPowerModel doesn't support InterconnectingConverter")` with the area/DC-bus - injection `add_to_expression!`. → `src/mt_hvdc_models/HVDCsystems.jl:284`. -- **#1566** — headroom-proportional slack **recompute from results**: `_update_headroom_participation_factors!` - / `_accumulate_headroom!` / `get_active_power_limits_for_power_flow` absent. → `ext/PowerFlowsExt/` + network_models. -- **#1612 (POM part)** — split In/Out headroom-proportional slack: data-map is ported, but - `_accumulate_in_out_headroom!` / `_find_paired_out` / `_pf_in_out_discharge_max` absent (old - `# Skip storage devices for now` TODO still present). → `ext/PowerFlowsExt/pf_headroom.jl`. -- **#1622** — route reactive-power TS to AC PF on active-power-only network models: add - `_add_pf_only_time_series_parameters!` + `PF_ONLY_TS_PARAMS_BY_CATEGORY`, thread `template` into - `add_power_flow_data!` (currently `add_power_flow_data!(container, transmission_model, sys)` at - `src/operation/build_problem.jl:157`). → `ext/PowerFlowsExt/pf_input_mapping.jl`. (Unblocks "reactive power on PTDFPowerModel" test.) - -**Large feature — Dynamic Line Ratings (DLR), SPLIT POM+IOM:** -- **#1559 + #1561** — DLR entirely absent in both clones (distinct from the static - branch-rating-time-series feature POM already has). POM part: `DynamicBranchRatingTimeSeriesParameter`, - default `"dynamic_line_ratings"`, DLR network-reduction handling, `get_dynamic_branch_rating_min_max_limits`, - `get_equivalent_dynamic_branch_rating`. IOM part: generic param-type registration (see iom plan). - Substantial; scope as its own effort. - -**Verify (not symbol-pinned — diff-review before deciding):** -- **#1509** (remove old N-1/G-1 code — arch diverged), **#1613** (network bugfixes), **#1619** (pnm/pf - logic; mostly version bumps = PSI-only). Confirm whether anything material is missing. - ---- - -## Workstream T — formulation TEST ports (code present in POM) -These are test-only ports against existing POM code (biggest coverage win, low risk). Land one PR -per file. Annotated where a test is **code-blocked** by a Workstream M/C item. - -- **`test_services_constructor.jl`** (~23 testsets): RangeReserve (thermal dispatch/UC, renewable), - RampReserve, StepwiseCostReserve, reserves+slacks, participation-factor limits, service-level - feedforwards, GroupReserve + errors, ConstantReserve, AGC; Transmission Interface (Constant/Variable - MaxFlow, TS limits, feedforwards, validation-under-reductions, double-circuit, on AreaInterchange, - with AreaBalance, with AreaPTDF). Code all present (incl. AGC, more complete than PSI). -- **`test_network_constructors.jl`** (~23 testsets): PowerModels iteration (All-PowerModels, - CopperPlate, PTDF, NFA, ACR/ACT, DCPLL, Unsupported guard); Area networks (AreaBalance ±slacks - ±TimeSeries, AreaPTDF ±DoubleCircuit ±TimeSeries); HVDC subnetworks; reductions (Ward, radial, - PTDF StaticBranch/StaticBranchBounds, subnetwork, branch-filter edge cases, parallel/series bounds, - full PowerModels×reduction matrix ±slacks). Code present. -- **`test_power_flow_in_the_loop.jl`** Source paths (~6): in/out variables, Source FixedOutput - (⚠ needs #1573), reactive-on-PTDF (⚠ needs #1622), in/out headroom (⚠ needs #1612), LCC HVDC. -- **`test_market_bid_cost.jl`** equivalence subset: no-TS vs constant-TS (Thermal/Renewable/PowerLoad), - RenewableDispatch MBC, Renewable-vs-Thermal compare, time-varying slopes/breakpoints/min-gen/everything - (⚠ tranche-count + concavity need Workstream C), VOM normalization (⚠ needs #1535), 3d results, - heterogeneous TS names, single TS. -- **Small device drops**: MonitoredLine asymmetric flow limits (#1604 test path exists — confirm/port), - renewable curtailment incentive (⚠ needs #1614), FixedOutput Source+PTDF+TS (⚠ needs #1573). - ---- - -## Workstream C — Tier-0 code blockers in POM (port code, then test) -- **Event framework — ported, including the runtime interface.** `src/event_models/` - (`event_model.jl`, `event_traits.jl`, `event_arguments.jl`, `event_constraints.jl`, - `event_runtime.jl`) provides the `EventModel`/`AbstractEventCondition` family, template-level - `set_event_model!` attachment, build-time discovery and time-series validation, event - parameters/constraints for thermal, renewable, load, hydro, pump-turbine, and storage devices - (hydro and storage go beyond the original PSI-parity scope), and the pure runtime functions a - simulation calls between solves (`outage_occurred`, `time_to_recover`, `event_step_values`, - `countdown_trajectory`, `event_parameter_keys`, `required_inputs`/`is_triggered`). PSI #1664 - consumes that interface; PSI keeps the state arrays, clock, and RNG. `test/test_events.jl` and - `test/test_event_runtime.jl` cover construction, traits, attachment, discovery/validation - errors, exclusion from the initialization problem, per-device constraint coefficients, the - runtime arithmetic, and E2E build/solve including forced-zero output under a - `FixedForcedOutage` event. Two IOM-side changes landed for it: the `SupplementalAttribute` - type bound (IOM #162) and the `share_template_references!` template-copy hook (IOM #164). -- **SSS `EnergyTargetFeedforward` — ported** as the storage twin of - `ReservoirTargetFeedforward`; the two storage feedforward testsets are live. -- **MBC variable-tranche-count** and **MBC concavity/convexity validation** — absent; small code adds - that unblock the remaining time-varying-tranche and validation MBC tests. - ---- - -## Out of scope for POM (PSI-only — do NOT port) -Simulation orchestration (`test_simulation_*`, simulation_state/partitions/store/results-IO), -`test_model_emulation`, `test_recorder_events`, `test_print`, `test_jump_utils`, docs/tutorials, -CI, version bumps, org renames. PSI template-library helpers (`_copy_template_for_build`, -"provided templates", namespace-stable keys) are architecture-specific to PSI; `modeled_ac_branch_types` -is tracked in POM as the `branches_modeled` trait (already present). - ---- - -## Suggested execution order -1. **Workstream T (1a services, 1b network)** — pure test ports, existing code, biggest coverage win. -2. **Workstream M bugfixes** (#1519, #1527, #1535, #1587, #1508) — small, unblock MBC/network tests. -3. **Workstream M small features** (#1549, #1573, #1538, #1614, #1605, #1622, #1566, #1612) — then - port the PF-Source / MBC / curtailment tests they unblock. -4. **Workstream G1** (#1617) — track upstream merge; port reserve/service SC layer + its test file. -5. **Workstream C** — event framework, runtime interface, and their tests: **done**; IOM #162 and - #164 are merged and every IOM pin is back on `main`. MBC tranche/concavity → remaining MBC - tests still open. -6. **DLR (#1559/#1561)** and the **verify** items — scope separately. diff --git a/Project.toml b/Project.toml index c6989cea..d9527342 100644 --- a/Project.toml +++ b/Project.toml @@ -35,7 +35,8 @@ PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6" [sources] InfrastructureSystems = {url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl.git", rev = "IS4"} -InfrastructureOptimizationModels = {rev = "main", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} +# TEMPORARY: repin to `rev = "main"` once IOM #171 (ac/service-containers) merges. +InfrastructureOptimizationModels = {rev = "ac/service-containers", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} PowerSystems = {url = "https://github.com/Sienna-Platform/PowerSystems.jl.git", rev = "psy6"} PowerNetworkMatrices = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"} # PSY psy6 depends on these unregistered packages; [sources] of non-root projects are @@ -61,7 +62,7 @@ DocStringExtensions = "~0.8, ~0.9" InfrastructureSystems = "3" InteractiveUtils = "1.11.0" JuMP = "^1.28" -PowerFlows = "^0.24" +PowerFlows = "0.25" PowerNetworkMatrices = "^0.24" PowerSystems = "^5.10" PrettyTables = "3" diff --git a/ext/PowerFlowsExt/pf_headroom.jl b/ext/PowerFlowsExt/pf_headroom.jl index 41855c12..636f60dd 100644 --- a/ext/PowerFlowsExt/pf_headroom.jl +++ b/ext/PowerFlowsExt/pf_headroom.jl @@ -13,7 +13,7 @@ _accumulate_headroom!( ::OptimizationContainerKey{<:ISOPT.ParameterType, <:PSY.Component}, ::Dict{String, Int}, ::Int, - ::Matrix{PSY.ACBusTypes}, + ::Matrix{PSY.ACBusTypes.Value}, ::Vector{Dict{Tuple{DataType, String}, Float64}}, ) = nothing @@ -27,7 +27,7 @@ _accumulate_headroom!( ::OptimizationContainerKey{<:ISOPT.OptimizationKeyType, <:PSY.Storage}, ::Dict{String, Int}, ::Int, - ::Matrix{PSY.ACBusTypes}, + ::Matrix{PSY.ACBusTypes.Value}, ::Vector{Dict{Tuple{DataType, String}, Float64}}, ) = nothing @@ -40,7 +40,7 @@ _accumulate_headroom!( ::OptimizationContainerKey{<:ISOPT.ParameterType, <:PSY.Storage}, ::Dict{String, Int64}, ::Int, - ::Matrix{PSY.ACBusTypes}, + ::Matrix{PSY.ACBusTypes.Value}, ::Vector{Dict{Tuple{DataType, String}, Float64}}, ) = nothing @@ -60,7 +60,7 @@ function _accumulate_headroom!( key::OptimizationContainerKey{<:ISOPT.OptimizationKeyType, U}, component_map::Dict{String, Int}, n_time_steps::Int, - bus_types::Matrix{PSY.ACBusTypes}, + bus_types::Matrix{PSY.ACBusTypes.Value}, computed_gspf::Vector{Dict{Tuple{DataType, String}, Float64}}, ) where {U <: PSY.Component} result = lookup_value(container, key) @@ -129,7 +129,7 @@ function _update_headroom_participation_factors!( PFS.get_computed_gspf(pf_data)::Vector{Dict{Tuple{DataType, String}, Float64}} n_time_steps = length(get_time_steps(container)) - bus_types = PFS.get_bus_type(pf_data)::Matrix{PSY.ACBusTypes} + bus_types = PFS.get_bus_type(pf_data)::Matrix{PSY.ACBusTypes.Value} bus_slack_pf = PFS.get_bus_slack_participation_factors( pf_data, diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index 9089fbab..f6f5b73b 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -346,6 +346,7 @@ include("services_models/reserve_group.jl") # include("services_models/agc.jl") # TODO: needs _get_ace_error include("services_models/transmission_interface.jl") include("services_models/services_constructor.jl") +include("services_models/security_constrained_injectors.jl") # Hybrid System Models (after services_models since they share reserve infrastructure) include("hybrid_system_models/hybrid_systems.jl") @@ -677,6 +678,12 @@ export HVDCReactivePowerToVariable export ShiftUpActivePowerVariable export ShiftDownActivePowerVariable +# G-1 Variables +export PostContingencyDeploymentVariable +export PostContingencyDeviationVariable +export PostGeneratorContingencyFlowSlackUpperBound +export PostGeneratorContingencyFlowSlackLowerBound + ######## Hydro Formulations ######## export HydroDispatchRunOfRiver export HydroCommitmentRunOfRiver @@ -868,6 +875,9 @@ export ShiftDownActivePowerVariableLimitsConstraint export RealizedShiftedLoadMinimumBoundConstraint export NonAnticipativityConstraint export HVDCDCControlConstraint +export PostContingencyBalanceConstraint +export PostContingencyGenerationConstraint +export PostContingencyDeploymentConstraint ################################################################################# # Exports - Expression Types (defined in core/expressions.jl) @@ -900,6 +910,10 @@ export ComponentReserveDownBalanceExpression export InterfaceTotalFlow export PTDFBranchFlow export BThetaBranchFlow +export PostContingencyNodalDeployment +export PostContingencyAreaDeployment +export PostContingencyInterchangeFlow +export PostContingencyTotalDeployment ################################################################################# # Exports - Parameter Types (defined in core/parameters.jl) @@ -1007,6 +1021,10 @@ export VariableMaxInterfaceFlow export ReserveLimitedRegulation export DeviceLimitedRegulation +# G-1 Formulations +export SecurityConstrainedContingencyReserve +export SecurityConstrainedRampReserve + ################################################################################# # Exports - Network Formulation Types (defined in core/network_formulations.jl) ################################################################################# diff --git a/src/common_models/add_to_expression.jl b/src/common_models/add_to_expression.jl index 59acd065..f1e8c0ea 100644 --- a/src/common_models/add_to_expression.jl +++ b/src/common_models/add_to_expression.jl @@ -1971,7 +1971,7 @@ function add_to_expression!( W <: AbstractReservesFormulation, } service_name = PSY.get_name(service) - variable = get_variable(container, U, X) + variable = get_variable(container, U, IOM.ComponentPairKey{V, X}) if !has_container_key(container, T, V) add_expressions!(container, T, devices, model) end @@ -2410,7 +2410,7 @@ function add_to_expression!( W <: AbstractReservesFormulation, } service_name = PSY.get_name(service) - variable = get_variable(container, U, X) + variable = get_variable(container, U, IOM.ComponentPairKey{V, X}) if !has_container_key(container, T, V) add_expressions!(container, T, devices, model) end @@ -2444,7 +2444,7 @@ function add_to_expression!( W <: AbstractReservesFormulation, } service_name = PSY.get_name(service) - variable = get_variable(container, U, X) + variable = get_variable(container, U, IOM.ComponentPairKey{V, X}) if !has_container_key(container, T, V) add_expressions!(container, T, devices, model) end @@ -2478,7 +2478,7 @@ function add_to_expression!( W <: AbstractReservesFormulation, } service_name = PSY.get_name(service) - variable = get_variable(container, U, X) + variable = get_variable(container, U, IOM.ComponentPairKey{V, X}) if !has_container_key(container, T, V) add_expressions!(container, T, devices, model) end diff --git a/src/common_models/converter_control.jl b/src/common_models/converter_control.jl index 3ac476a4..480374bc 100644 --- a/src/common_models/converter_control.jl +++ b/src/common_models/converter_control.jl @@ -18,7 +18,7 @@ end # AC control on one terminal/converter: AC_VOLTAGE pins the regulated bus # VoltageMagnitude; AC_REACTIVE_POWER pins the reactive injection to its setpoint. function _fix_converter_ac_control!( - mode::PSY.VSCACControlModes, + mode::PSY.VSCACControlModes.Value, setpoint::Float64, vm, bus_name::String, @@ -43,7 +43,7 @@ end # magnitude deviation phi = |V| - 1, so AC_VOLTAGE pins phi to setpoint - 1. # AC_REACTIVE_POWER pins the reactive injection to its setpoint, unshifted. function _fix_converter_ac_control_lpacc!( - mode::PSY.VSCACControlModes, + mode::PSY.VSCACControlModes.Value, setpoint::Float64, phi, bus_name::String, @@ -70,7 +70,7 @@ end function _fill_converter_dc_control!( jump_model, con::AbstractArray, - mode::PSY.VSCDCControlModes, + mode::PSY.VSCDCControlModes.Value, setpoint::Float64, droop_gain::Float64, vdc_var, @@ -103,7 +103,7 @@ end # through the RegulatedVoltageMagnitude aux variable by the caller, not here): # AC_REACTIVE_POWER pins the reactive injection; AC_VOLTAGE is a no-op here. function _fix_converter_ac_reactive!( - mode::PSY.VSCACControlModes, + mode::PSY.VSCACControlModes.Value, setpoint::Float64, q_var, name::String, diff --git a/src/core/constraints.jl b/src/core/constraints.jl index 49350d10..3643a8e4 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -1274,3 +1274,33 @@ Single award variable per (device, service): the device's merged offer curve pri provision states (documented approximation). """ struct OfflineReserveBandConstraint <: ConstraintType end + +################################################################################# +# G-1 constraints +################################################################################# + +""" +Post-contingency power balance under outage ``o``. Copper plate and PTDF networks balance the +system: ``\\sum_d \\Delta_{d,o,t} = \\sum_{g \\in G_o} p_{g,t}``. Area balance networks +balance each area ``a`` over the interchanges ``i`` into and out of it: + +``\\Delta P_{a,o,t} + \\sum_{i \\to a} \\Delta f_{i,o,t} - \\sum_{i \\leftarrow a} \\Delta f_{i,o,t} = 0``. + +See [`SecurityConstrainedContingencyReserve`](@ref). +""" +struct PostContingencyBalanceConstraint <: ConstraintType end + +""" +A contributing device's post-contingency output stays within its maximum: +``p_{d,t} + \\Delta_{d,o,t} \\le P^{max}_d``. + +See [`SecurityConstrainedContingencyReserve`](@ref). +""" +struct PostContingencyGenerationConstraint <: ConstraintType end + +""" +A procured reserve deploys at most its award: ``\\delta_{s,d,o,t} \\le r_{s,d,t}``. + +See [`SecurityConstrainedContingencyReserve`](@ref). +""" +struct PostContingencyDeploymentConstraint <: ConstraintType end diff --git a/src/core/expressions.jl b/src/core/expressions.jl index bf4c9b37..45aafb08 100644 --- a/src/core/expressions.jl +++ b/src/core/expressions.jl @@ -15,7 +15,15 @@ variable); `FlowRateConstraint` rows are written directly on it. Reportable as a mirroring `PTDFBranchFlow`. """ struct BThetaBranchFlow <: ExpressionType end -struct PostContingencyNodalActivePowerDeployment <: PostContingencyExpressions end + +""" +Post-contingency change in injection at bus ``n`` under outage ``o``: +``\\Delta P_{n,o,t} = \\sum_{d \\in n} \\Delta_{d,o,t} - \\sum_{g \\in G_o \\cap n} p_{g,t}``. + +Sparse, indexed `(bus number, outage, time)`. See +[`SecurityConstrainedContingencyReserve`](@ref). +""" +struct PostContingencyNodalDeployment <: PostContingencyExpressions end struct RealizedShiftedLoad <: ExpressionType end ################################################################################# @@ -124,6 +132,37 @@ Aggregation of reserve variables allocated to the storage subcomponent of a hybr struct StorageReserveBalanceExpression{D, S, Sd} <: ReserveAggregationExpression{D, S, Sd} end +################################################################################# +# G-1 expressions +################################################################################# + +""" +Post-contingency change in injection in area ``a`` under outage ``o``: +``\\Delta P_{a,o,t} = \\sum_{d \\in a} \\Delta_{d,o,t} - \\sum_{g \\in G_o \\cap a} p_{g,t}``. + +Sparse, indexed `(area, outage, time)`. See [`SecurityConstrainedContingencyReserve`](@ref). +""" +struct PostContingencyAreaDeployment <: ExpressionType end + +""" +Post-contingency flow over monitored area interchange ``i`` under outage ``o``: +``f^o_{i,t} = f_{i,t} + \\Delta f_{i,o,t}``. + +Sparse, indexed `(interchange, outage, time)`, meta `"G1"`. See +[`SecurityConstrainedContingencyReserve`](@ref). +""" +struct PostContingencyInterchangeFlow <: ExpressionType end + +""" +Post-contingency deployment of contributing device ``d`` under outage ``o``, summed across the +security-constrained reserves responding to it: +``\\Delta_{d,o,t} = \\sum_s \\delta_{s,d,o,t}``. + +Sparse, indexed `(device, outage, time)`. Stored because the balance, generation-limit, and +locational-deployment terms all read it. See [`SecurityConstrainedContingencyReserve`](@ref). +""" +struct PostContingencyTotalDeployment <: ExpressionType end + # Method extensions for output writing should_write_resulting_value(::Type{InterfaceTotalFlow}) = true should_write_resulting_value(::Type{PTDFBranchFlow}) = true diff --git a/src/core/formulations.jl b/src/core/formulations.jl index b32b677f..cd4c4bca 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -821,3 +821,72 @@ unit and a standalone copy produce identical objective coefficients). When regularization slacks. """ struct HybridDispatchWithReserves <: AbstractHybridFormulationWithReserves end + +############################### G-1 Formulations ##################################### + +abstract type AbstractSecurityConstrainedReservesFormulation <: AbstractReservesFormulation end + +""" +Security-constrained (G-1) contingency reserve for `PSY.OnlineReserve{PSY.ReserveUp}` and +`PSY.OfflineReserve`: reserves must be deliverable after the loss of any generator outage +they respond to. + +**Sets.** Each reserve ``s`` responds to the `PSY.Outage`s ``O_s`` attached to it with +`add_supplemental_attribute!(sys, service, outage)`. Outage ``o`` takes offline the +generators ``G_o`` associated with it and monitors the components listed in its +`monitored_components`. ``D_s`` are the contributing devices of ``s``. Outaged generators +that are not modeled are skipped with a warning. + +**Procurement.** A reserve with a requirement time series is procured as under +[`RampReserve`](@ref) without the ramp limit: awards ``r_{s,d,t}``, the requirement, +participation-fraction limits, and the reserve cost. A reserve without one is not procured +and only deploys post-contingency. + +**Deployment.** Under each outage ``o \\in O_s``, every contributor ``d \\in D_s \\setminus G_o`` +deploys ``\\delta_{s,d,o,t} \\ge 0`` ([`PostContingencyDeploymentVariable`](@ref)), capped by +its award when ``s`` is procured ([`PostContingencyDeploymentConstraint`](@ref)). A device's +deployment across reserves, ``\\Delta_{d,o,t}`` +([`PostContingencyTotalDeployment`](@ref)), keeps it within its maximum +([`PostContingencyGenerationConstraint`](@ref)). + +**Balance.** Deployment replaces the outaged generation +([`PostContingencyBalanceConstraint`](@ref)): + + - `CopperPlateNetworkModel` and PTDF networks balance the system: + ``\\sum_d \\Delta_{d,o,t} = \\sum_{g \\in G_o} p_{g,t}``. + - `AreaBalanceNetworkModel` balances each area. Every modeled area interchange carries a + deviation ``\\Delta f_{i,o,t}`` ([`PostContingencyDeviationVariable`](@ref)) that moves + deployment between areas. Without an `AreaInterchange` device model each area covers its + own outages. + +**Flow limits.** Monitored components get post-contingency flow limits at their emergency +rating (`PostContingencyFlowRateConstraint`, metas `"G1_lb"`/`"G1_ub"`): + + - PTDF networks: monitored branches carry + ``f^o_{\\ell,t} = f_{\\ell,t} + \\sum_n PTDF_{\\ell,n} \\Delta P_{n,o,t}`` + (`PostContingencyBranchFlow`, meta `"G1"`), with ``\\Delta P_{n,o,t}`` the nodal change + in injection ([`PostContingencyNodalDeployment`](@ref)). Parallel circuits and reduced + branches are limited once, on their reduced entry. + - `AreaBalanceNetworkModel`: monitored interchanges carry + ``f^o_{i,t} = f_{i,t} + \\Delta f_{i,o,t}`` + ([`PostContingencyInterchangeFlow`](@ref)). Deviations exist on every modeled + interchange and enter every area balance; only monitored interchanges are limited. + +Monitored components must be modeled (including by the branch model's `filter_function`), or +template validation fails. With `use_slacks = true` the flow limits are relaxed by +[`PostGeneratorContingencyFlowSlackUpperBound`](@ref) and +[`PostGeneratorContingencyFlowSlackLowerBound`](@ref). Flow limits are shared by every +reserve responding to an outage, so their service models must agree on `use_slacks`. + +See also [`SecurityConstrainedRampReserve`](@ref). +""" +struct SecurityConstrainedContingencyReserve <: + AbstractSecurityConstrainedReservesFormulation end + +""" +Same as [`SecurityConstrainedContingencyReserve`](@ref), except every reserve is procured: +a requirement time series is required, and spinning reserves' awards are also limited by +the contributing devices' ramp rates over the reserve time frame (`RampConstraint`), as in +[`RampReserve`](@ref). +""" +struct SecurityConstrainedRampReserve <: AbstractSecurityConstrainedReservesFormulation end diff --git a/src/core/variables.jl b/src/core/variables.jl index cc6104fd..00defcea 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -841,6 +841,48 @@ Reserve allocated to one side of a hybrid system's storage subcomponent. Paramet struct HybridStorageSubcomponentReserveVariable{Sd <: ReserveSide} <: AbstractHybridReserveVariableType end +################################################################################# +# G-1 Variables +################################################################################# + +""" +Reserve deployed by contributing device ``d`` on security-constrained reserve ``s`` under +outage ``o``: ``\\delta_{s,d,o,t} \\ge 0``. Devices taken offline by ``o`` get none. + +Keyed on `ComponentPairKey{D, S}`, indexed `(service, device, outage, time)`. See +[`SecurityConstrainedContingencyReserve`](@ref). +""" +struct PostContingencyDeploymentVariable <: VariableType end + +""" +Post-contingency change in flow ``\\Delta f_{i,o,t}`` over modeled area interchange ``i`` +under outage ``o``. Free. + +Indexed `(interchange, outage, time)` over every modeled interchange. See +[`SecurityConstrainedContingencyReserve`](@ref). +""" +struct PostContingencyDeviationVariable <: VariableType end + +""" +Relaxes the upper post-contingency flow limit of security-constrained reserves when +`use_slacks = true`: ``f^o_{\\ell,t} - \\sigma^+_{\\ell,o,t} \\le R^{max}_\\ell``, +``\\sigma^+ \\ge 0``. + +A separate type from [`PostContingencyFlowActivePowerSlackUpperBound`](@ref) because variable +containers take no `meta`, and branch security-constrained models key that one by the same +branch type. +""" +struct PostGeneratorContingencyFlowSlackUpperBound <: VariableType end + +""" +Relaxes the lower post-contingency flow limit of security-constrained reserves when +`use_slacks = true`: ``f^o_{\\ell,t} + \\sigma^-_{\\ell,o,t} \\ge R^{min}_\\ell``, +``\\sigma^- \\ge 0``. + +See [`PostGeneratorContingencyFlowSlackUpperBound`](@ref). +""" +struct PostGeneratorContingencyFlowSlackLowerBound <: VariableType end + const MULTI_START_VARIABLES = (HotStartVariable, WarmStartVariable, ColdStartVariable) should_write_resulting_value(::Type{PiecewiseLinearCostVariable}) = false @@ -881,3 +923,5 @@ convert_output_to_natural_units(::Type{HVDCLosses}) = true convert_output_to_natural_units(::Type{InterfaceFlowSlackUp}) = true convert_output_to_natural_units(::Type{InterfaceFlowSlackDown}) = true convert_output_to_natural_units(::Type{ActivePowerPumpVariable}) = true +convert_output_to_natural_units(::Type{PostContingencyDeploymentVariable}) = true +convert_output_to_natural_units(::Type{PostContingencyDeviationVariable}) = true diff --git a/src/energy_storage_models/storage_models.jl b/src/energy_storage_models/storage_models.jl index dc1b7def..56dce775 100644 --- a/src/energy_storage_models/storage_models.jl +++ b/src/energy_storage_models/storage_models.jl @@ -724,7 +724,7 @@ function add_to_expression!( W <: AbstractReservesFormulation, } s_name = PSY.get_name(service) - variable = get_variable(container, U, V) + variable = get_variable(container, U, IOM.ComponentPairKey{UV, V}) for d in devices name = PSY.get_name(d) expression = get_expression(container, T, UV, _service_container_meta(service)) diff --git a/src/hybrid_system_models/hybrid_systems.jl b/src/hybrid_system_models/hybrid_systems.jl index a8329e28..ef06556a 100644 --- a/src/hybrid_system_models/hybrid_systems.jl +++ b/src/hybrid_system_models/hybrid_systems.jl @@ -744,7 +744,7 @@ function add_to_expression!( W <: AbstractReservesFormulation, } s_name = PSY.get_name(service) - variable = get_variable(container, U, V) + variable = get_variable(container, U, IOM.ComponentPairKey{UV, V}) for d in devices name = PSY.get_name(d) expression = get_expression(container, T, UV, _service_container_meta(service)) @@ -2065,7 +2065,7 @@ function add_constraints!( names, time_steps; meta = "$(s_type)_$s_name") # System-level reserve variable for this service, keyed `(service, device, time)`. - sys_reserve = get_variable(container, ActivePowerReserveVariable, s_type) + sys_reserve = _reserve_variable(container, V, s_type) # Per-hybrid reserve variables for this service r_out = get_variable( container, diff --git a/src/operation/decision_model.jl b/src/operation/decision_model.jl index 3a7899fa..19f72b56 100644 --- a/src/operation/decision_model.jl +++ b/src/operation/decision_model.jl @@ -231,8 +231,9 @@ function solve!( # `power_units = :component_base` is what PSY stores internally, so the # write needs no unit conversion and no round-trip ledger — a model's # system carries one only when it was built from a document. - !ispath(sys_dir) && - PSY.to_file(sys, sys_dir; power_units = :component_base) + # TODO: re-enable before merging, once PSY.to_file is fixed. + #!ispath(sys_dir) && + # PSY.to_file(sys, sys_dir; power_units = :component_base) end end @info "\n$(RUN_OPERATION_MODEL_TIMER)\n" diff --git a/src/operation/template_validation.jl b/src/operation/template_validation.jl index c131dc4b..4744de21 100644 --- a/src/operation/template_validation.jl +++ b/src/operation/template_validation.jl @@ -109,6 +109,7 @@ function validate_template_impl!(model::IOM.AbstractOptimizationModel) delete!(template.branches, k) end _check_interface_branches(template, system, network_model) + _check_security_constrained_reserve_monitors(template, system, network_model) _check_security_constrained_three_winding_transformer(template.branches) _check_security_constrained_network(template.branches, network_model) _check_security_constrained_phase_control(template.branches, network_model) diff --git a/src/services_models/reserve_group.jl b/src/services_models/reserve_group.jl index c495c1b5..991f68d2 100644 --- a/src/services_models/reserve_group.jl +++ b/src/services_models/reserve_group.jl @@ -112,10 +112,13 @@ function check_activeservice_variables( ) where {T <: PSY.Service} for service in contributing_services service_name = PSY.get_name(service) - variable = get_variable(container, ActivePowerReserveVariable, typeof(service)) - # The container is keyed `(service_name, device_name, time)` and shared by the whole - # service type, so check for this service's own entries, not just that it exists. - any(k -> k[1] == service_name, keys(variable.data)) || error( + # Containers are keyed `(service_name, device_name, time)` and shared by the whole + # service type, so check for this service's own entries, not just that one exists. + has_entries = any( + variable -> any(k -> k[1] == service_name, keys(variable.data)), + _reserve_variables(container, typeof(service)), + ) + has_entries || error( "The contributing service $service_name has no ActivePowerReserveVariable \ entries; it must be modeled before the group reserve that references it.", ) @@ -196,8 +199,8 @@ end # Collect the group's contributing reserve variables into one bucket per time step, so the # constraint loop above indexes straight in rather than re-scanning per `(group, t)`. Services -# of the same type share one `(service_name, device_name, time)` container, so each container -# is scanned once. +# of the same type share their `(service_name, device_name, time)` containers, so each service +# type is scanned once. function _group_member_variables( container::OptimizationContainer, contributing_services::Vector{<:PSY.Service}, @@ -210,11 +213,27 @@ function _group_member_variables( rtype = typeof(r) rtype in scanned && continue push!(scanned, rtype) - reserve_variable = get_variable(container, ActivePowerReserveVariable, rtype) - for (key, var) in reserve_variable.data - key[1] in member_names || continue - push!(index[key[3]], var) + for reserve_variable in _reserve_variables(container, rtype) + for (key, var) in reserve_variable.data + key[1] in member_names || continue + push!(index[key[3]], var) + end end end return index end + +# Award containers of every contributing device type for service type `SR`. Groups only know +# their member services, not the members' contributing device types. Keys store the service +# type with its unit parameter stripped, so match on that form. +function _reserve_variables( + container::OptimizationContainer, + ::Type{SR}, +) where {SR <: PSY.Service} + S = IOM.canonical_component_type(SR) + return [ + variable for (key, variable) in IOM.get_variables(container) if + IOM.get_entry_type(key) === ActivePowerReserveVariable && + get_component_type(key) <: IOM.ComponentPairKey{<:PSY.Component, S} + ] +end diff --git a/src/services_models/reserve_offers.jl b/src/services_models/reserve_offers.jl index e60dad28..f1a2bfd0 100644 --- a/src/services_models/reserve_offers.jl +++ b/src/services_models/reserve_offers.jl @@ -22,7 +22,7 @@ _cost_offers_reserve(cost::Union{PSY.MarketBidCost, PSY.MarketBidTimeSeriesCost} _cost_offers_reserve(::PSY.OperationalCost, service) = false # Price every contributing device that offers into `service` by its offer curve; returns the set of -# device names so priced (the flat-cost pass skips them). +# `(device type, device name)` so priced (the flat-cost pass skips them). # A group has no contributing devices, so it can carry no per-device offers; with # `GroupReserve <: AbstractReserve` the generic method below would otherwise accept it. # Offers live on the group's contributing services and are priced by their own models. @@ -45,19 +45,19 @@ function add_reserve_offer_costs!( model::ServiceModel{SR, T}, ) where {SR <: PSY.AbstractReserve, T <: AbstractReservesFormulation} service_name = PSY.get_name(service) - award = get_variable(container, ActivePowerReserveVariable, SR) time_steps = get_time_steps(container) resolution = get_resolution(container) dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR base_p = get_model_base_power(container) initial_time = IOM.get_initial_time(container) jump_model = get_jump_model(container) - offered = Set{String}() + offered = Set{Tuple{DataType, String}}() for (device_type, devices) in get_contributing_devices_map(model, service_name) offering = [d for d in devices if _has_reserve_offer(d, service)] isempty(offering) && continue names = [PSY.get_name(d) for d in offering] + award = _reserve_variable(container, device_type, SR) # Block var keyed `(service, device, segment, time)` via # `IOM.sparse_variable_key_type(PiecewiseLinearBlockReserveOffer)`. Segments vary per # device and time, so they are filled sparsely below. @@ -99,7 +99,7 @@ function add_reserve_offer_costs!( add_to_objective_invariant_expression!( container, get_pwl_cost_expression_delta(pwl_vars, slopes, dt)) end - push!(offered, dev_name) + push!(offered, (device_type, dev_name)) end end return offered diff --git a/src/services_models/reserves.jl b/src/services_models/reserves.jl index 4b74a58f..9508c22a 100644 --- a/src/services_models/reserves.jl +++ b/src/services_models/reserves.jl @@ -148,6 +148,15 @@ function get_default_time_series_names( return Dict{Type{<:TimeSeriesParameter}, String}() end +function get_default_time_series_names( + ::Type{<:PSY.AbstractReserve}, + ::Type{<:AbstractSecurityConstrainedReservesFormulation}, +) + return Dict{Type{<:TimeSeriesParameter}, String}( + RequirementTimeSeriesParameter => "requirement", + ) +end + function get_default_attributes( ::Type{<:PSY.AbstractReserve}, ::Type{<:AbstractReservesFormulation}, @@ -188,46 +197,77 @@ function add_reserve_variables!( return end -# Sum the reserve provision of one service across its contributing devices at time `t`, -# reading the service type's sparse container keyed `(service_name, device_name, time)`. +_reserve_variable( + container::OptimizationContainer, + ::Type{D}, + ::Type{SR}, +) where {D <: PSY.Component, SR <: PSY.Service} = + get_variable(container, ActivePowerReserveVariable, IOM.ComponentPairKey{D, SR}) + +# Per time step, the sum of one service's awards across all its contributing device types. function _sum_service_reserves( - reserve_variable::SparseAxisArray, + container::OptimizationContainer, + ::Type{SR}, service_name::String, - contributing_devices::U, - t::Int, + contributing_devices::AbstractDict, extra::Int, -) where { - U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, -} where {D <: PSY.Component} - acc = IOM.get_hinted_aff_expr(length(contributing_devices) + extra) - for d in contributing_devices - JuMP.add_to_expression!(acc, reserve_variable[(service_name, PSY.get_name(d), t)]) +) where {SR <: PSY.Service} + n_terms = sum(length, values(contributing_devices); init = 0) + extra + acc = [IOM.get_hinted_aff_expr(n_terms) for _ in get_time_steps(container)] + for (device_type, devices) in contributing_devices + _sum_service_reserves!( + acc, + _reserve_variable(container, device_type, SR), + service_name, + devices, + ) end return acc end +function _sum_service_reserves!( + acc::Vector{JuMP.AffExpr}, + reserve_variable::SparseAxisArray, + service_name::String, + devices::Vector{D}, +) where {D <: PSY.Component} + for d in devices, t in eachindex(acc) + JuMP.add_to_expression!( + acc[t], + reserve_variable[(service_name, PSY.get_name(d), t)], + ) + end + return +end + ################################## Reserve Requirement Constraint ########################## function add_constraints!( container::OptimizationContainer, T::Type{RequirementConstraint}, service::SR, - contributing_devices::U, + contributing_devices::AbstractDict, model::ServiceModel{SR, V}, -) where { - SR <: PSY.AbstractReserve, - V <: AbstractReservesFormulation, - U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, -} where {D <: PSY.Component} +) where {SR <: PSY.AbstractReserve, V <: AbstractReservesFormulation} time_steps = get_time_steps(container) service_name = PSY.get_name(service) # Dense container keyed `[service_name, time]`, built per type; fill this service's row. constraint = get_constraint(container, T, SR) - reserve_variable = get_variable(container, ActivePowerReserveVariable, SR) use_slacks = get_use_slacks(model) - use_slacks && (slack_vars = get_variable(container, ReserveRequirementSlack, SR)) requirement = _get_requirement(service) jump_model = get_jump_model(container) - extra = use_slacks ? 1 : 0 + resource_expression = _sum_service_reserves( + container, + SR, + service_name, + contributing_devices, + use_slacks ? 1 : 0, + ) + if use_slacks + slack_vars = get_variable(container, ReserveRequirementSlack, SR) + for t in time_steps + JuMP.add_to_expression!(resource_expression[t], slack_vars[service_name, t]) + end + end # A static reserve gets a scalar requirement RHS; a time-varying one scales it by an # attached requirement series (resolved by the model-configured name). @@ -237,20 +277,10 @@ function add_constraints!( get_parameter(container, RequirementTimeSeriesParameter, SR) param = get_parameter_column_refs(param_container, service_name) for t in time_steps - resource_expression = - _sum_service_reserves(reserve_variable, service_name, - contributing_devices, - t, extra) - use_slacks && - JuMP.add_to_expression!( - resource_expression, - slack_vars[service_name, t], - ) - constraint[service_name, t] = - JuMP.@constraint( - jump_model, - resource_expression >= param[t] * requirement - ) + constraint[service_name, t] = JuMP.@constraint( + jump_model, + resource_expression[t] >= param[t] * requirement + ) end else ts_vector = IOM.get_time_series( @@ -260,31 +290,16 @@ function add_constraints!( interval = get_interval(get_settings(container)), ) for t in time_steps - resource_expression = - _sum_service_reserves(reserve_variable, service_name, - contributing_devices, - t, extra) - use_slacks && - JuMP.add_to_expression!( - resource_expression, - slack_vars[service_name, t], - ) constraint[service_name, t] = JuMP.@constraint( jump_model, - resource_expression >= ts_vector[t] * requirement + resource_expression[t] >= ts_vector[t] * requirement ) end end else for t in time_steps - resource_expression = - _sum_service_reserves(reserve_variable, service_name, contributing_devices, - t, - extra) - use_slacks && - JuMP.add_to_expression!(resource_expression, slack_vars[service_name, t]) constraint[service_name, t] = - JuMP.@constraint(jump_model, resource_expression >= requirement) + JuMP.@constraint(jump_model, resource_expression[t] >= requirement) end end return @@ -294,13 +309,9 @@ function add_constraints!( container::OptimizationContainer, T::Type{ParticipationFractionConstraint}, service::SR, - contributing_devices::U, + contributing_devices::Vector{D}, model::ServiceModel{SR, V}, -) where { - SR <: PSY.AbstractReserve, - V <: AbstractReservesFormulation, - U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, -} where {D <: PSY.Device} +) where {SR <: PSY.AbstractReserve, V <: AbstractReservesFormulation, D <: PSY.Component} max_participation_factor = PSY.get_max_participation_factor(service) if max_participation_factor >= 1.0 @@ -310,13 +321,16 @@ function add_constraints!( time_steps = get_time_steps(container) service_name = PSY.get_name(service) # Sparse constraint container keyed `(service_name, device_name, time)`. - cons = lazy_container_addition!(container, T, SR, - [service_name], - [PSY.get_name(d) for d in contributing_devices], - time_steps; + cons = lazy_container_addition!( + container, + T, + IOM.ComponentPairKey{D, SR}, + String[], + String[], + Int[]; sparse = true, ) - var_r = get_variable(container, ActivePowerReserveVariable, SR) + var_r = _reserve_variable(container, D, SR) jump_model = get_jump_model(container) requirement = _get_requirement(service) cap = requirement * max_participation_factor @@ -370,11 +384,11 @@ function add_to_objective_function!( # Devices that submitted a reserve OFFER are priced by their offer curve; the rest keep the # flat DEFAULT_RESERVE_COST. offered = add_reserve_offer_costs!(container, service, model) - contributing_names = - [PSY.get_name(d) for d in get_contributing_devices(model, PSY.get_name(service))] - add_reserves_proportional_cost!( - container, ActivePowerReserveVariable, service, T, contributing_names; - skip_devices = offered) + for devices in values(get_contributing_devices_map(model, PSY.get_name(service))) + add_reserves_proportional_cost!( + container, ActivePowerReserveVariable, service, T, devices; + skip_devices = offered) + end return end @@ -382,32 +396,22 @@ function add_constraints!( container::OptimizationContainer, T::Type{RequirementConstraint}, service::SR, - contributing_devices::U, + contributing_devices::AbstractDict, ::ServiceModel{SR, StepwiseCostReserve}, -) where { - SR <: PSY.AbstractReserve, - U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, -} where {D <: PSY.Component} +) where {SR <: PSY.AbstractReserve} time_steps = get_time_steps(container) service_name = PSY.get_name(service) # Dense container keyed `[service_name, time]`, built per type; fill this service's row. constraint = get_constraint(container, T, SR) - reserve_variable = get_variable(container, ActivePowerReserveVariable, SR) requirement_variable = get_variable(container, ServiceRequirementVariable, SR) jump_model = get_jump_model(container) + resource_expression = + _sum_service_reserves(container, SR, service_name, contributing_devices, 0) for t in time_steps - resource_expression = - _sum_service_reserves( - reserve_variable, - service_name, - contributing_devices, - t, - 0, - ) constraint[service_name, t] = JuMP.@constraint( jump_model, - resource_expression >= requirement_variable[service_name, t] + resource_expression[t] >= requirement_variable[service_name, t] ) end @@ -440,82 +444,48 @@ function _get_ramp_constraint_contributing_devices( return filtered_device end -function add_constraints!( - container::OptimizationContainer, - T::Type{RampConstraint}, - service::SR, - contributing_devices::Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, - ::ServiceModel{SR, V}, -) where { - SR <: PSY.Reserve{PSY.ReserveUp}, - V <: AbstractReservesFormulation, - D <: PSY.Component, -} - ramp_devices = _get_ramp_constraint_contributing_devices(service, contributing_devices) - service_name = PSY.get_name(service) - if !isempty(ramp_devices) - jump_model = get_jump_model(container) - time_steps = get_time_steps(container) - time_frame = PSY.get_time_frame(service) - variable = get_variable(container, ActivePowerReserveVariable, SR) - device_name_set = [PSY.get_name(d) for d in ramp_devices] - con_up = lazy_container_addition!(container, T, - SR, - [service_name], - device_name_set, - time_steps; - sparse = true, - ) - for d in ramp_devices, t in time_steps - name = PSY.get_name(d) - ramp_limits = PSY.get_ramp_limits(d, PSY.SU / u"minute") - con_up[(service_name, name, t)] = JuMP.@constraint( - jump_model, - variable[(service_name, name, t)] <= ramp_limits.up * time_frame - ) - end - else - @warn "Data doesn't contain contributing devices with ramp limits for service $service_name, consider adjusting your formulation" - end - return -end +_directional_ramp_limit(ramp_limits, ::Type{<:PSY.Reserve{PSY.ReserveUp}}) = + ramp_limits.up +_directional_ramp_limit(ramp_limits, ::Type{<:PSY.Reserve{PSY.ReserveDown}}) = + ramp_limits.down function add_constraints!( container::OptimizationContainer, T::Type{RampConstraint}, service::SR, - contributing_devices::Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, + contributing_devices::Vector{D}, ::ServiceModel{SR, V}, ) where { - SR <: PSY.Reserve{PSY.ReserveDown}, + SR <: Union{PSY.Reserve{PSY.ReserveUp}, PSY.Reserve{PSY.ReserveDown}}, V <: AbstractReservesFormulation, D <: PSY.Component, } ramp_devices = _get_ramp_constraint_contributing_devices(service, contributing_devices) service_name = PSY.get_name(service) - if !isempty(ramp_devices) - jump_model = get_jump_model(container) - time_steps = get_time_steps(container) - time_frame = PSY.get_time_frame(service) - variable = get_variable(container, ActivePowerReserveVariable, SR) - device_name_set = [PSY.get_name(d) for d in ramp_devices] - con_down = lazy_container_addition!(container, T, - SR, - [service_name], - device_name_set, - time_steps; - sparse = true, + if isempty(ramp_devices) + @warn "Contributing $(D) devices on service $service_name have no binding ramp limits; no ramp constraints are added for them." + return + end + jump_model = get_jump_model(container) + time_steps = get_time_steps(container) + time_frame = PSY.get_time_frame(service) + variable = _reserve_variable(container, D, SR) + cons = lazy_container_addition!( + container, + T, + IOM.ComponentPairKey{D, SR}, + String[], + String[], + Int[]; + sparse = true, + ) + for d in ramp_devices, t in time_steps + name = PSY.get_name(d) + limit = _directional_ramp_limit(PSY.get_ramp_limits(d, PSY.SU / u"minute"), SR) + cons[(service_name, name, t)] = JuMP.@constraint( + jump_model, + variable[(service_name, name, t)] <= limit * time_frame ) - for d in ramp_devices, t in time_steps - name = PSY.get_name(d) - ramp_limits = PSY.get_ramp_limits(d, PSY.SU / u"minute") - con_down[(service_name, name, t)] = JuMP.@constraint( - jump_model, - variable[(service_name, name, t)] <= ramp_limits.down * time_frame - ) - end - else - @warn "Data doesn't contain contributing devices with ramp limits for service $service_name, consider adjusting your formulation" end return end @@ -524,13 +494,9 @@ function add_constraints!( container::OptimizationContainer, T::Type{ReservePowerConstraint}, service::SR, - contributing_devices::U, + contributing_devices::Vector{D}, ::ServiceModel{SR, V}, -) where { - SR <: PSY.OfflineReserve, - V <: AbstractReservesFormulation, - U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, -} where {D <: PSY.Component} +) where {SR <: PSY.OfflineReserve, V <: AbstractReservesFormulation, D <: PSY.Component} time_steps = get_time_steps(container) resolution = get_resolution(container) if resolution > Dates.Minute(1) @@ -540,62 +506,37 @@ function add_constraints!( minutes_per_period = Dates.value(Dates.Second(resolution)) / 60 end service_name = PSY.get_name(service) - cons = lazy_container_addition!(container, T, - SR, - [service_name], - [PSY.get_name(d) for d in contributing_devices], - time_steps; + cons = lazy_container_addition!( + container, + T, + IOM.ComponentPairKey{D, SR}, + String[], + String[], + Int[]; sparse = true, ) - var_r = get_variable(container, ActivePowerReserveVariable, SR) + var_r = _reserve_variable(container, D, SR) + varstatus = get_variable(container, OnVariable, D) reserve_response_time = PSY.get_time_frame(service) jump_model = get_jump_model(container) for d in contributing_devices - # Function barrier: `contributing_devices` may have an abstract element type, so the - # callee specializes on the concrete types and dispatches once per device rather than - # once per timestep. - varstatus = get_variable(container, OnVariable, typeof(d)) - _add_reserve_power_constraint_device!( - cons, - var_r, - varstatus, - d, - service_name, - reserve_response_time, - minutes_per_period, - jump_model, - time_steps, - ) - end - return -end - -function _add_reserve_power_constraint_device!( - cons, - var_r, - varstatus, - d::D, - service_name::String, - reserve_response_time, - minutes_per_period, - jump_model, - time_steps, -) where {D <: PSY.Component} - name = PSY.get_name(d) - startup_time = PSY.get_time_limits(d).up - ramp_limits = _get_ramp_limits(d) - if reserve_response_time > startup_time - reserve_limit = - PSY.get_active_power_limits(d, PSY.SU).min + - (reserve_response_time - startup_time) * minutes_per_period * ramp_limits.up - else - reserve_limit = 0.0 - end - for t in time_steps - cons[(service_name, name, t)] = JuMP.@constraint( - jump_model, - var_r[(service_name, name, t)] <= (1 - varstatus[name, t]) * reserve_limit - ) + name = PSY.get_name(d) + startup_time = PSY.get_time_limits(d).up + ramp_limits = _get_ramp_limits(d) + if reserve_response_time > startup_time + reserve_limit = + PSY.get_active_power_limits(d, PSY.SU).min + + (reserve_response_time - startup_time) * minutes_per_period * + ramp_limits.up + else + reserve_limit = 0.0 + end + for t in time_steps + cons[(service_name, name, t)] = JuMP.@constraint( + jump_model, + var_r[(service_name, name, t)] <= (1 - varstatus[name, t]) * reserve_limit + ) + end end return end @@ -721,27 +662,26 @@ function process_stepwise_cost_reserve_parameters!( return end +# `skip_devices` are priced by their offer curve in `add_reserve_offer_costs!` instead. function add_reserves_proportional_cost!( container::OptimizationContainer, ::Type{U}, service::T, ::Type{V}, - contributing_names::Vector{String}; - skip_devices = Set{String}(), + contributing_devices::Vector{D}; + skip_devices = Set{Tuple{DataType, String}}(), ) where { T <: PSY.AbstractReserve, U <: ActivePowerReserveVariable, V <: AbstractReservesFormulation, + D <: PSY.Component, } - base_p = get_model_base_power(container) service_name = PSY.get_name(service) - reserve_variable = get_variable(container, U, T) - # Index this service's slice of the `(service, device, time)` container by its contributing - # device names, so each provision is priced once without scanning the whole container. - # `skip_devices` are priced by their offer curve in `add_reserve_offer_costs!` instead. - cost = DEFAULT_RESERVE_COST / base_p - for name in contributing_names - name in skip_devices && continue + reserve_variable = get_variable(container, U, IOM.ComponentPairKey{D, T}) + cost = DEFAULT_RESERVE_COST / get_model_base_power(container) + for d in contributing_devices + name = PSY.get_name(d) + (D, name) in skip_devices && continue for t in get_time_steps(container) add_to_objective_invariant_expression!( container, diff --git a/src/services_models/security_constrained_injectors.jl b/src/services_models/security_constrained_injectors.jl new file mode 100644 index 00000000..3bab9e67 --- /dev/null +++ b/src/services_models/security_constrained_injectors.jl @@ -0,0 +1,880 @@ +const _PER_TYPE = Dict{DataType, Set{String}} +const _OUTAGE_MAP = Dict{Int, _PER_TYPE} + +const _G1_META = "G1" + +_validate_reserve_formulation(::ServiceModel) = false +_validate_reserve_formulation( + ::ServiceModel{<:_RESERVE_UP, <:AbstractSecurityConstrainedReservesFormulation}, +) = true +_validate_reserve_formulation( + ::ServiceModel{<:PSY.AbstractReserve, <:AbstractSecurityConstrainedReservesFormulation}, +) = throw( + IS.ConflictingInputsError( + "Security-constrained formulations currently only support reserve-up services.", + ), +) + +_valid_component_type(::PSY.ACTransmission, ::NetworkModel{<:AbstractPTDFNetworkModel}) = + true +_valid_component_type(::PSY.AreaInterchange, ::NetworkModel{AreaBalanceNetworkModel}) = true +_valid_component_type(::PSY.Component, ::NetworkModel) = false + +""" +Outages attached to security-constrained reserves, by UUID, and whether each one's +post-contingency flow limits are relaxed. + +Flow constraints are shared by every reserve responding to an outage, so the service models +carrying that outage must agree on `use_slacks`. +""" +function _security_constrained_outages( + sys::PSY.System, + services_template::ServicesModelContainer, +) + outages = Dict{Int, PSY.Outage}() + use_slacks = Dict{Int, Bool}() + for model in values(services_template) + _validate_reserve_formulation(model) || continue + model_slacks = get_use_slacks(model) + for service in _services_with_contributors(model, sys) + for outage in PSY.get_supplemental_attributes(PSY.Outage, service) + uuid = IS.get_id(outage) + if get!(use_slacks, uuid, model_slacks) != model_slacks + throw( + IS.ConflictingInputsError( + "Outage $uuid is attached to security-constrained reserves \ + with different `use_slacks` settings; set `use_slacks` \ + consistently across their service models.", + ), + ) + end + outages[uuid] = outage + end + end + end + return outages, use_slacks +end + +""" +Reject monitored components of security-constrained reserve outages that the template does +not model, including those excluded by a branch model's `filter_function`. Post-contingency +flows are built only on modeled components, so monitored components must be a subset of the +modeled ones. +""" +function _check_security_constrained_reserve_monitors( + template::PowerOperationsProblemTemplate, + sys::PSY.System, + network_model::NetworkModel, +) + problems = String[] + checked = Set{Int}() + for model in values(get_service_models(template)) + _validate_reserve_formulation(model) || continue + for service in get_available_components(model, sys), + outage in PSY.get_supplemental_attributes(PSY.Outage, service) + + uuid = IS.get_id(outage) + uuid in checked && continue + push!(checked, uuid) + for component_uuid in PSY.get_monitored_components(outage) + component = IS.get_component(sys, component_uuid) + _valid_component_type(component, network_model) || continue + PSY.get_available(component) || continue + branch_model = get_model(template, typeof(component)) + name = PSY.get_name(component) + if isnothing(branch_model) || + !any(c -> PSY.get_name(c) == name, get_device_cache(branch_model)) + push!( + problems, + "Outage $uuid monitors $(typeof(component)) $name, which the \ + template does not model (no branch model, or excluded by its \ + filter_function).", + ) + end + end + end + end + isempty(problems) || throw(IS.ConflictingInputsError(join(problems, "\n"))) + return +end + +# Outaged generators whose power the model can remove; the rest are skipped with a warning. +function _outaged_generators( + container::OptimizationContainer, + sys::PSY.System, + outage::PSY.Outage, +) + outaged = _PER_TYPE() + for generator in + PSY.get_associated_components(sys, outage; component_type = PSY.Generator) + T = typeof(generator) + name = PSY.get_name(generator) + if has_container_key(container, ActivePowerVariable, T) && + name in axes(get_variable(container, ActivePowerVariable, T), 1) + push!(get!(Set{String}, outaged, T), name) + else + @warn "Generator $name ($T) outaged by outage $(IS.get_id(outage)) is not \ + modeled; it is left out of the post-contingency balance." _group = + LOG_GROUP_SERVICE_CONSTUCTORS + end + end + return outaged +end + +function _monitored_components( + sys::PSY.System, + outages::Dict{Int, PSY.Outage}, + network_model::NetworkModel, +) + monitored_components = _OUTAGE_MAP() + for (uuid, outage) in outages + monitored = _PER_TYPE() + for component_uuid in PSY.get_monitored_components(outage) + component = IS.get_component(sys, component_uuid) + _valid_component_type(component, network_model) || continue + PSY.get_available(component) || continue + typeof(component) in network_model.modeled_branch_types || continue + push!(get!(Set{String}, monitored, typeof(component)), PSY.get_name(component)) + end + monitored_components[uuid] = monitored + end + return monitored_components +end + +# Parallel circuits share one reduced entry, so each entry is constrained once. +function _flow_entries( + network_model::NetworkModel{<:AbstractPTDFNetworkModel}, + ::Type{T}, + names::Set{String}, +) where {T <: PSY.ACTransmission} + reduction_name_map = + PNM.get_component_to_reduction_name_map(get_branch_catalog(network_model), T) + return Set{String}(reduction_name_map[name] for name in names) +end + +_flow_entries( + ::NetworkModel{AreaBalanceNetworkModel}, + ::Type{PSY.AreaInterchange}, + names::Set{String}, +) = names + +function _post_contingency_flow_limits( + ::PSY.System, + network_model::NetworkModel{<:AbstractPTDFNetworkModel}, + ::Type{T}, + entry_name::String, +) where {T <: PSY.ACTransmission} + catalog = get_branch_catalog(network_model) + arc = PNM.get_name_to_arc_map(catalog, T)[entry_name] + return _emergency_flow_limits(PNM.get_reduction_entry(catalog, arc)) +end + +function _post_contingency_flow_limits( + sys::PSY.System, + ::NetworkModel{AreaBalanceNetworkModel}, + ::Type{PSY.AreaInterchange}, + name::String, +) + limits = PSY.get_flow_limits(PSY.get_component(PSY.AreaInterchange, sys, name), PSY.SU) + return (min = -limits.to_from, max = limits.from_to) +end + +################################## Construction ########################################### + +function _construct_post_contingency!( + container::OptimizationContainer, + sys::PSY.System, + ::ArgumentConstructStage, + services_template::ServicesModelContainer, + network_model::NetworkModel, +) + outages, use_slacks = _security_constrained_outages(sys, services_template) + isempty(outages) && return + _add_post_contingency_flow_slacks!( + container, + _monitored_components(sys, outages, network_model), + use_slacks, + network_model, + ) + return +end + +# Every outage gets a balance row per network region. Modeled interchanges carry a deviation +# variable per outage, balanced across areas; monitored interchanges (a subset of the modeled +# ones) and monitored branches additionally get post-contingency flow limits. +function _construct_post_contingency!( + container::OptimizationContainer, + sys::PSY.System, + ::ModelConstructStage, + services_template::ServicesModelContainer, + network_model::NetworkModel, +) + outages, use_slacks = _security_constrained_outages(sys, services_template) + isempty(outages) && return + outaged_generators = _OUTAGE_MAP( + uuid => _outaged_generators(container, sys, outage) for (uuid, outage) in outages + ) + monitored_components = _monitored_components(sys, outages, network_model) + uuids = sort!(collect(keys(outages))) + # AreaInterchange flow variables are created by the branch constructors, which run after the + # services argument stage. + _add_post_contingency_deviation_variables!(container, uuids, network_model) + _add_post_contingency_locational_deployment!( + container, + sys, + outaged_generators, + network_model, + ) + _add_post_contingency_flow!(container, monitored_components, network_model) + _add_post_contingency_balance_constraints!( + container, + sys, + outaged_generators, + uuids, + network_model, + ) + _add_post_contingency_generation_constraints!(container, sys) + _add_post_contingency_flow_constraints!( + container, + sys, + monitored_components, + use_slacks, + network_model, + ) + _add_post_contingency_slack_costs!(container, monitored_components) + return +end + +################################## Deployment ############################################# + +# Each service deploys its own contributors under every outage it responds to, except the +# contributors that outage takes offline; the per-device sum across services is what the +# outage-level constraints see. +function _add_post_contingency_deployment!( + container::OptimizationContainer, + sys::PSY.System, + service::R, + devices::Vector{D}, +) where {R <: PSY.Service, D <: PSY.Component} + jump_model = get_jump_model(container) + service_name = PSY.get_name(service) + # Lazy: every service of type `R` that `D` contributes to shares this container. + deployment = lazy_container_addition!( + container, + PostContingencyDeploymentVariable, + IOM.ComponentPairKey{D, R}, + String[], + String[], + Int[], + Int[]; + sparse = true, + ) + # Lazy: shared by every service `D` contributes to. Keyed `(device name, outage, time)`. + total = lazy_container_addition!( + container, + PostContingencyTotalDeployment, + D, + String[], + Int[], + Int[]; + sparse = true, + ) + for outage in PSY.get_supplemental_attributes(PSY.Outage, service) + uuid = IS.get_id(outage) + outaged = Set{String}( + PSY.get_name(c) for + c in PSY.get_associated_components(sys, outage; component_type = D) + ) + for d in devices + name = PSY.get_name(d) + name in outaged && continue + for t in get_time_steps(container) + var = + deployment[service_name, name, uuid, t] = JuMP.@variable( + jump_model, + base_name = "PostContingencyDeploymentVariable_$(D)_$(R)_{$(service_name), $(name), $(uuid), $(t)}", + lower_bound = 0.0, + ) + JuMP.add_to_expression!( + get!(JuMP.AffExpr, total.data, (name, uuid, t)), + var, + ) + end + end + end + return +end + +# Deployment can only draw on procured reserve. Unprocured reserves are limited by generation +# headroom alone. +function add_constraints!( + container::OptimizationContainer, + ::Type{PostContingencyDeploymentConstraint}, + service::R, + devices::Vector{D}, + ::ServiceModel{R, <:AbstractSecurityConstrainedReservesFormulation}, +) where {R <: _RESERVE_UP, D <: PSY.Component} + key = IOM.ComponentPairKey{D, R} + jump_model = get_jump_model(container) + service_name = PSY.get_name(service) + deployment = get_variable(container, PostContingencyDeploymentVariable, key) + award = _reserve_variable(container, D, R) + cons = lazy_container_addition!( + container, + PostContingencyDeploymentConstraint, + key, + String[], + String[], + Int[], + Int[]; + sparse = true, + ) + uuids = [IS.get_id(o) for o in PSY.get_supplemental_attributes(PSY.Outage, service)] + for d in devices, uuid in uuids, t in get_time_steps(container) + name = PSY.get_name(d) + # Devices the outage takes offline have no deployment. + r = get(deployment.data, (service_name, name, uuid, t), nothing) + isnothing(r) && continue + cons[service_name, name, uuid, t] = + JuMP.@constraint(jump_model, r <= award[service_name, name, t]) + end + return +end + +# Device types are only known from the container keys: services register a total-deployment +# container per contributing device type. +_total_deployments(container::OptimizationContainer) = [ + (get_component_type(key), expr) for + (key, expr) in IOM.get_expressions(container) if + IOM.get_entry_type(key) === PostContingencyTotalDeployment +] + +################################## Variables ############################################## + +function _add_post_contingency_deviation_variables!( + container::OptimizationContainer, + uuids::Vector{Int}, + ::NetworkModel{AreaBalanceNetworkModel}, +) + if !has_container_key(container, FlowActivePowerVariable, PSY.AreaInterchange) + @warn "An AreaBalanceNetworkModel with security-constrained reserves needs PSY.AreaInterchange(s) and DeviceModel{PSY.AreaInterchange} for reserve deployment to cross area boundaries. Otherwise, each area must cover its own outages." _group = + LOG_GROUP_SERVICE_CONSTUCTORS + return + end + flow = get_variable(container, FlowActivePowerVariable, PSY.AreaInterchange) + jump_model = get_jump_model(container) + time_steps = get_time_steps(container) + names = axes(flow, 1) + var = add_variable_container!( + container, + PostContingencyDeviationVariable, + PSY.AreaInterchange, + names, + uuids, + time_steps, + ) + for name in names, uuid in uuids, t in time_steps + var[name, uuid, t] = JuMP.@variable( + jump_model, + base_name = "PostContingencyDeviationVariable_AreaInterchange_{$(name), $(uuid), $(t)}", + ) + end + return +end + +_add_post_contingency_deviation_variables!( + ::OptimizationContainer, + ::Vector{Int}, + ::NetworkModel, +) = nothing + +function _add_post_contingency_flow_slacks!( + container::OptimizationContainer, + monitored_components::_OUTAGE_MAP, + use_slacks::Dict{Int, Bool}, + network_model::NetworkModel{<:Union{AbstractPTDFNetworkModel, AreaBalanceNetworkModel}}, +) + jump_model = get_jump_model(container) + time_steps = get_time_steps(container) + for (uuid, per_type) in monitored_components + use_slacks[uuid] || continue + for (component_type, names) in per_type + # Lazy: slack containers are per component type and shared across outages. + # Keyed `(flow entry, outage, time)`, sparse since outages monitor different entries. + slack_ub = lazy_container_addition!( + container, + PostGeneratorContingencyFlowSlackUpperBound, + component_type, + String[], + Int[], + Int[]; + sparse = true, + ) + slack_lb = lazy_container_addition!( + container, + PostGeneratorContingencyFlowSlackLowerBound, + component_type, + String[], + Int[], + Int[]; + sparse = true, + ) + for entry_name in _flow_entries(network_model, component_type, names), + t in time_steps + + slack_ub[entry_name, uuid, t] = JuMP.@variable( + jump_model, + base_name = "PostGeneratorContingencyFlowSlackUpperBound_$(component_type)_{$(entry_name), $(uuid), $(t)}", + lower_bound = 0.0, + ) + slack_lb[entry_name, uuid, t] = JuMP.@variable( + jump_model, + base_name = "PostGeneratorContingencyFlowSlackLowerBound_$(component_type)_{$(entry_name), $(uuid), $(t)}", + lower_bound = 0.0, + ) + end + end + end + return +end + +_add_post_contingency_flow_slacks!( + ::OptimizationContainer, + ::_OUTAGE_MAP, + ::Dict{Int, Bool}, + ::NetworkModel, +) = nothing + +################################## Expressions ############################################ + +# Keyed `(bus number or area name, outage, time)`, sparse since an outage only touches the +# locations of its deployments and outaged generators. +_add_locational_deployment_container!( + container::OptimizationContainer, + ::NetworkModel{<:AbstractPTDFNetworkModel}, +) = add_expression_container!( + container, + PostContingencyNodalDeployment, + PSY.ACBus, + String[], + Int[], + Int[]; + sparse = true, +) +_add_locational_deployment_container!( + container::OptimizationContainer, + ::NetworkModel{AreaBalanceNetworkModel}, +) = add_expression_container!( + container, + PostContingencyAreaDeployment, + PSY.Area, + String[], + Int[], + Int[]; + sparse = true, +) + +_location_key(component, network_model::NetworkModel{<:AbstractPTDFNetworkModel}) = string( + PNM.get_mapped_bus_number(get_network_reduction(network_model), PSY.get_bus(component)), +) +_location_key(component, ::NetworkModel{AreaBalanceNetworkModel}) = + PSY.get_name(PSY.get_area(PSY.get_bus(component))) + +# [contributing device deployment] minus [outaged generator power] per bus or area +function _add_post_contingency_locational_deployment!( + container::OptimizationContainer, + sys::PSY.System, + outaged_generators::_OUTAGE_MAP, + network_model::NetworkModel{<:Union{AbstractPTDFNetworkModel, AreaBalanceNetworkModel}}, +) + expr = _add_locational_deployment_container!(container, network_model) + + for (device_type, total) in _total_deployments(container) + locations = Dict{String, String}() + for ((name, uuid, t), deployed) in total.data + key = get!(locations, name) do + _location_key(PSY.get_component(device_type, sys, name), network_model) + end + JuMP.add_to_expression!(get!(JuMP.AffExpr, expr.data, (key, uuid, t)), deployed) + end + end + + for (uuid, per_type) in outaged_generators + for (generator_type, names) in per_type + power = get_variable(container, ActivePowerVariable, generator_type) + for name in names + key = _location_key( + PSY.get_component(generator_type, sys, name), + network_model, + ) + for t in get_time_steps(container) + JuMP.add_to_expression!( + get!(JuMP.AffExpr, expr.data, (key, uuid, t)), + -1.0, + power[name, t], + ) + end + end + end + end + return +end + +_add_post_contingency_locational_deployment!( + ::OptimizationContainer, + ::PSY.System, + ::_OUTAGE_MAP, + ::NetworkModel, +) = nothing + +# Monitored names resolve to their reduction entry, the row the pre-contingency +# `PTDFBranchFlow` and the post-contingency expression are keyed by. +function _add_post_contingency_flow!( + container::OptimizationContainer, + monitored_components::_OUTAGE_MAP, + network_model::NetworkModel{<:AbstractPTDFNetworkModel}, +) + time_steps = get_time_steps(container) + catalog = get_branch_catalog(network_model) + ptdf = get_network_matrix(network_model) + bus_axis = PNM.get_bus_axis(ptdf) + + nodal = get_expression(container, PostContingencyNodalDeployment, PSY.ACBus) + buses = Dict{Int, Set{String}}() + for (bus, uuid, t) in keys(nodal.data) + push!(get!(Set{String}, buses, uuid), bus) + end + + # Per reduced entry, its PTDF row keyed by the nodal bus key, dropping entries below + # PTDF_ZERO_TOL. + nonzero_factors = Dict{Tuple{DataType, String}, Dict{String, Float64}}() + for (uuid, per_type) in monitored_components + outage_buses = get(buses, uuid, Set{String}()) + for (line_type, names) in per_type + # Keyed `(flow entry, outage, time)`, sparse since outages monitor different entries. + expr = lazy_container_addition!( + container, + PostContingencyBranchFlow, + line_type, + String[], + Int[], + Int[]; + sparse = true, + meta = _G1_META, + ) + pre_flow = get_expression(container, PTDFBranchFlow, line_type) + arc_map = PNM.get_name_to_arc_map(catalog, line_type) + for entry_name in _flow_entries(network_model, line_type, names) + factors = get!(nonzero_factors, (line_type, entry_name)) do + ptdf_col = ptdf[arc_map[entry_name], :] + Dict{String, Float64}( + string(bus_axis[i]) => ptdf_col[i] for + i in eachindex(ptdf_col) if abs(ptdf_col[i]) > PTDF_ZERO_TOL + ) + end + for t in time_steps + ex = + expr[entry_name, uuid, t] = IOM.get_hinted_aff_expr( + length(JuMP.linear_terms(pre_flow[entry_name, t])) + + length(outage_buses), + ) + JuMP.add_to_expression!(ex, pre_flow[entry_name, t]) + for bus in outage_buses + haskey(factors, bus) || continue + JuMP.add_to_expression!(ex, factors[bus], nodal[bus, uuid, t]) + end + end + end + end + end + return +end + +function _add_post_contingency_flow!( + container::OptimizationContainer, + monitored_components::_OUTAGE_MAP, + ::NetworkModel{AreaBalanceNetworkModel}, +) + has_container_key( + container, + PostContingencyDeviationVariable, + PSY.AreaInterchange, + ) || return + time_steps = get_time_steps(container) + # Keyed `(interchange, outage, time)`, sparse since outages monitor different interchanges. + expr = add_expression_container!( + container, + PostContingencyInterchangeFlow, + PSY.AreaInterchange, + String[], + Int[], + Int[]; + sparse = true, + meta = _G1_META, + ) + flow = get_variable(container, FlowActivePowerVariable, PSY.AreaInterchange) + deviation = get_variable( + container, + PostContingencyDeviationVariable, + PSY.AreaInterchange, + ) + for (uuid, per_type) in monitored_components + for name in get(per_type, PSY.AreaInterchange, Set{String}()), t in time_steps + ex = expr[name, uuid, t] = JuMP.AffExpr(0.0) + JuMP.add_to_expression!(ex, flow[name, t]) + JuMP.add_to_expression!(ex, deviation[name, uuid, t]) + end + end + return +end + +_add_post_contingency_flow!(::OptimizationContainer, ::_OUTAGE_MAP, ::NetworkModel) = + nothing + +################################## Constraints ############################################ + +function _add_post_contingency_balance_constraints!( + container::OptimizationContainer, + ::PSY.System, + outaged_generators::_OUTAGE_MAP, + uuids::Vector{Int}, + ::NetworkModel, +) + time_steps = get_time_steps(container) + jump_model = get_jump_model(container) + cons = add_constraints_container!( + container, + PostContingencyBalanceConstraint, + PSY.System, + uuids, + time_steps, + ) + + balance = Dict((uuid, t) => JuMP.AffExpr(0.0) for uuid in uuids, t in time_steps) + for (uuid, per_type) in outaged_generators + for (generator_type, names) in per_type + power = get_variable(container, ActivePowerVariable, generator_type) + for name in names, t in time_steps + JuMP.add_to_expression!(balance[uuid, t], -1.0, power[name, t]) + end + end + end + for (_, total) in _total_deployments(container) + for ((_, uuid, t), deployed) in total.data + JuMP.add_to_expression!(balance[uuid, t], deployed) + end + end + for ((uuid, t), ex) in balance + cons[uuid, t] = JuMP.@constraint(jump_model, ex == 0.0) + end + return +end + +function _add_post_contingency_balance_constraints!( + container::OptimizationContainer, + sys::PSY.System, + ::_OUTAGE_MAP, + uuids::Vector{Int}, + ::NetworkModel{AreaBalanceNetworkModel}, +) + time_steps = get_time_steps(container) + jump_model = get_jump_model(container) + + # Area name => (sign, interchange name). Without an AreaInterchange model there are no + # deviations, so each area covers its own outages. + interchanges = Dict{String, Vector{Tuple{Float64, String}}}() + if has_container_key( + container, + PostContingencyDeviationVariable, + PSY.AreaInterchange, + ) + deviation = get_variable( + container, + PostContingencyDeviationVariable, + PSY.AreaInterchange, + ) + for name in axes(deviation, 1) + interchange = PSY.get_component(PSY.AreaInterchange, sys, name) + from_area = PSY.get_name(PSY.get_from_area(interchange)) + to_area = PSY.get_name(PSY.get_to_area(interchange)) + push!( + get!(Vector{Tuple{Float64, String}}, interchanges, from_area), + (-1.0, name), + ) + push!(get!(Vector{Tuple{Float64, String}}, interchanges, to_area), (1.0, name)) + end + end + deployment = get_expression(container, PostContingencyAreaDeployment, PSY.Area) + + area_names = PSY.get_name.(PSY.get_components(PSY.Area, sys)) + cons = add_constraints_container!( + container, + PostContingencyBalanceConstraint, + PSY.Area, + area_names, + uuids, + time_steps, + ) + # Every area needs a row, or deviations into an area without deployment are unconstrained. + for area_name in area_names, uuid in uuids, t in time_steps + balance = JuMP.AffExpr(0.0) + JuMP.add_to_expression!( + balance, + get(deployment.data, (area_name, uuid, t), zero(JuMP.AffExpr)), + ) + for (sign, interchange_name) in get(interchanges, area_name, ()) + JuMP.add_to_expression!(balance, sign, deviation[interchange_name, uuid, t]) + end + cons[area_name, uuid, t] = JuMP.@constraint(jump_model, balance == 0.0) + end + return +end + +# Redundant for devices bounded by `PostContingencyDeploymentConstraint`, since their awards +# already sit in the device headroom; kept unconditional so no device can deploy past pmax. +function _add_post_contingency_generation_constraints!( + container::OptimizationContainer, + sys::PSY.System, +) + jump_model = get_jump_model(container) + for (device_type, total) in _total_deployments(container) + cons = add_constraints_container!( + container, + PostContingencyGenerationConstraint, + device_type, + String[], + Int[], + Int[]; + sparse = true, + ) + power = get_variable(container, ActivePowerVariable, device_type) + limits = Dict{String, Float64}() + for ((name, uuid, t), deployed) in total.data + limit = get!(limits, name) do + PSY.get_max_active_power(PSY.get_component(device_type, sys, name), PSY.SU) + end + cons[name, uuid, t] = + JuMP.@constraint(jump_model, power[name, t] + deployed <= limit) + end + end + return +end + +_post_contingency_flow_expression(::NetworkModel{<:AbstractPTDFNetworkModel}) = + PostContingencyBranchFlow +_post_contingency_flow_expression(::NetworkModel{AreaBalanceNetworkModel}) = + PostContingencyInterchangeFlow + +function _add_post_contingency_flow_constraints!( + container::OptimizationContainer, + sys::PSY.System, + monitored_components::_OUTAGE_MAP, + use_slacks::Dict{Int, Bool}, + network_model::NetworkModel{<:Union{AbstractPTDFNetworkModel, AreaBalanceNetworkModel}}, +) + jump_model = get_jump_model(container) + time_steps = get_time_steps(container) + for (uuid, per_type) in monitored_components + slacked = use_slacks[uuid] + for (component_type, names) in per_type + flow = get_expression( + container, + _post_contingency_flow_expression(network_model), + component_type, + _G1_META, + ) + cons_lb = lazy_container_addition!( + container, + PostContingencyFlowRateConstraint, + component_type, + String[], + Int[], + Int[]; + sparse = true, + meta = "$(_G1_META)_lb", + ) + cons_ub = lazy_container_addition!( + container, + PostContingencyFlowRateConstraint, + component_type, + String[], + Int[], + Int[]; + sparse = true, + meta = "$(_G1_META)_ub", + ) + if slacked + slack_ub = get_variable( + container, + PostGeneratorContingencyFlowSlackUpperBound, + component_type, + ) + slack_lb = get_variable( + container, + PostGeneratorContingencyFlowSlackLowerBound, + component_type, + ) + end + for entry_name in _flow_entries(network_model, component_type, names) + lims = _post_contingency_flow_limits( + sys, + network_model, + component_type, + entry_name, + ) + for t in time_steps + f = flow[entry_name, uuid, t] + if slacked + cons_ub[entry_name, uuid, t] = JuMP.@constraint( + jump_model, + f - slack_ub[entry_name, uuid, t] <= lims.max + ) + cons_lb[entry_name, uuid, t] = JuMP.@constraint( + jump_model, + f + slack_lb[entry_name, uuid, t] >= lims.min + ) + else + cons_ub[entry_name, uuid, t] = + JuMP.@constraint(jump_model, f <= lims.max) + cons_lb[entry_name, uuid, t] = + JuMP.@constraint(jump_model, f >= lims.min) + end + end + end + end + end + return +end + +_add_post_contingency_flow_constraints!( + ::OptimizationContainer, + ::PSY.System, + ::_OUTAGE_MAP, + ::Dict{Int, Bool}, + ::NetworkModel, +) = nothing + +function _add_post_contingency_slack_costs!( + container::OptimizationContainer, + monitored_components::_OUTAGE_MAP, +) + component_types = Set{DataType}() + for per_type in values(monitored_components) + union!(component_types, keys(per_type)) + end + for component_type in component_types, + slack_type in ( + PostGeneratorContingencyFlowSlackUpperBound, + PostGeneratorContingencyFlowSlackLowerBound, + ) + + has_container_key(container, slack_type, component_type) || continue + for slack in values(get_variable(container, slack_type, component_type).data) + add_to_objective_invariant_expression!( + container, + slack * CONSTRAINT_VIOLATION_SLACK_COST, + ) + end + end + return +end diff --git a/src/services_models/services_constructor.jl b/src/services_models/services_constructor.jl index 3f24f5b2..d5983389 100644 --- a/src/services_models/services_constructor.jl +++ b/src/services_models/services_constructor.jl @@ -1,9 +1,10 @@ # One `ServiceModel` per service TYPE (like `DeviceModel`). `construct_service!` runs once # per type: it gets all services of the type via `get_available_components(model, sys)`, # reads each service's contributing devices from the nested per-service map -# (`get_contributing_devices(model, service_name)`), and builds. Reserve variable and -# constraint containers are shared per `(entry type, service type)`, with each service -# filling its own slice. Group formulations are deferred to last (their members must exist). +# (`get_contributing_devices_map(model, service_name)`), and builds. Reserve award containers +# are shared per `(device type, service type)` and constraint containers per +# `(entry type, service type)`, with each service filling its own slice. Group formulations +# are deferred to last (their members must exist). # # TODO(services stability): See issue #216. @@ -47,6 +48,23 @@ device model's full available component set. Runs once per service model, before any service wires awards in, so services of the same type with contributor sets that do not nest all index an axis that holds their devices. """ +# Function barrier +function _seed_range_expression!( + container::OptimizationContainer, + sys::PSY.System, + ::Type{T}, + device_model::DeviceModel{D, W}, +) where {T <: ExpressionType, D <: PSY.Component, W <: AbstractDeviceFormulation} + has_container_key(container, T, D) && return + add_expressions!( + container, + T, + get_available_components(device_model, sys), + device_model, + ) + return +end + seed_reserve_range_expressions!( ::OptimizationContainer, ::PSY.System, @@ -81,24 +99,6 @@ function seed_reserve_range_expressions!( return end -# Function barrier: the caller's loop is uninferable, so the container work happens here -# where `T`, `D` and `W` are concrete. -function _seed_range_expression!( - container::OptimizationContainer, - sys::PSY.System, - ::Type{T}, - device_model::DeviceModel{D, W}, -) where {T <: ExpressionType, D <: PSY.Component, W <: AbstractDeviceFormulation} - has_container_key(container, T, D) && return - add_expressions!( - container, - T, - get_available_components(device_model, sys), - device_model, - ) - return -end - function construct_services!( container::OptimizationContainer, sys::PSY.System, @@ -139,6 +139,8 @@ function construct_services!( network_model, ) end + + _construct_post_contingency!(container, sys, stage, services_template, network_model) return end @@ -181,6 +183,8 @@ function construct_services!( network_model, ) end + + _construct_post_contingency!(container, sys, stage, services_template, network_model) return end @@ -204,15 +208,14 @@ function construct_service!( ts_services = [s for s in demand_services if _has_ts_requirement(model, s)] isempty(ts_services) || add_parameters!(container, RequirementTimeSeriesParameter, ts_services, model) + add_service_variables!( + container, + ActivePowerReserveVariable, + services, + model, + RangeReserve, + ) for service in services - contributing_devices = get_contributing_devices(model, PSY.get_name(service)) - add_service_variables!( - container, - ActivePowerReserveVariable, - service, - contributing_devices, - RangeReserve, - ) add_to_expression!( container, ActivePowerReserveVariable, @@ -254,7 +257,7 @@ function construct_service!( get_use_slacks(model) && add_reserve_slacks!(container, SR, demand_names) end for service in services - contributing_devices = get_contributing_devices(model, PSY.get_name(service)) + contributing_devices = get_contributing_devices_map(model, PSY.get_name(service)) if _has_reserve_demand(model, service) add_constraints!( container, @@ -263,13 +266,15 @@ function construct_service!( contributing_devices, model, ) - add_constraints!( - container, - ParticipationFractionConstraint, - service, - contributing_devices, - model, - ) + for devices in values(contributing_devices) + add_constraints!( + container, + ParticipationFractionConstraint, + service, + devices, + model, + ) + end add_to_objective_function!(container, service, model) else # Supply-only: no requirement of its own (it may serve a GroupReserve). Price any @@ -311,15 +316,14 @@ function construct_service!( # Slope/breakpoint PWL cost params for the time-series-backed ORDCs (no-op otherwise). process_stepwise_cost_reserve_parameters!(container, model, demand_services) end + add_service_variables!( + container, + ActivePowerReserveVariable, + services, + model, + StepwiseCostReserve, + ) for service in services - contributing_devices = get_contributing_devices(model, PSY.get_name(service)) - add_service_variables!( - container, - ActivePowerReserveVariable, - service, - contributing_devices, - StepwiseCostReserve, - ) add_to_expression!( container, ActivePowerReserveVariable, @@ -354,7 +358,7 @@ function construct_service!( ) end for service in services - contributing_devices = get_contributing_devices(model, PSY.get_name(service)) + contributing_devices = get_contributing_devices_map(model, PSY.get_name(service)) if _has_reserve_demand(model, service) add_constraints!( container, @@ -602,15 +606,14 @@ function construct_service!( ts_services = [s for s in services if _has_ts_requirement(model, s)] isempty(ts_services) || add_parameters!(container, RequirementTimeSeriesParameter, ts_services, model) + add_service_variables!( + container, + ActivePowerReserveVariable, + services, + model, + RampReserve, + ) for service in services - contributing_devices = get_contributing_devices(model, PSY.get_name(service)) - add_service_variables!( - container, - ActivePowerReserveVariable, - service, - contributing_devices, - RampReserve, - ) add_to_expression!( container, ActivePowerReserveVariable, @@ -645,7 +648,7 @@ function construct_service!( ) get_use_slacks(model) && add_reserve_slacks!(container, SR, service_names) for service in services - contributing_devices = get_contributing_devices(model, PSY.get_name(service)) + contributing_devices = get_contributing_devices_map(model, PSY.get_name(service)) add_constraints!( container, RequirementConstraint, @@ -653,14 +656,16 @@ function construct_service!( contributing_devices, model, ) - add_constraints!(container, RampConstraint, service, contributing_devices, model) - add_constraints!( - container, - ParticipationFractionConstraint, - service, - contributing_devices, - model, - ) + for devices in values(contributing_devices) + add_constraints!(container, RampConstraint, service, devices, model) + add_constraints!( + container, + ParticipationFractionConstraint, + service, + devices, + model, + ) + end add_to_objective_function!(container, service, model) add_feedforward_constraints!(container, model, service) end @@ -684,15 +689,14 @@ function construct_service!( ts_services = [s for s in services if _has_ts_requirement(model, s)] isempty(ts_services) || add_parameters!(container, RequirementTimeSeriesParameter, ts_services, model) + add_service_variables!( + container, + ActivePowerReserveVariable, + services, + model, + NonSpinningReserve, + ) for service in services - contributing_devices = get_contributing_devices(model, PSY.get_name(service)) - add_service_variables!( - container, - ActivePowerReserveVariable, - service, - contributing_devices, - NonSpinningReserve, - ) add_feedforward_arguments!(container, model, service) end return @@ -720,7 +724,7 @@ function construct_service!( ) get_use_slacks(model) && add_reserve_slacks!(container, SR, service_names) for service in services - contributing_devices = get_contributing_devices(model, PSY.get_name(service)) + contributing_devices = get_contributing_devices_map(model, PSY.get_name(service)) add_constraints!( container, RequirementConstraint, @@ -728,20 +732,16 @@ function construct_service!( contributing_devices, model, ) - add_constraints!( - container, - ReservePowerConstraint, - service, - contributing_devices, - model, - ) - add_constraints!( - container, - ParticipationFractionConstraint, - service, - contributing_devices, - model, - ) + for devices in values(contributing_devices) + add_constraints!(container, ReservePowerConstraint, service, devices, model) + add_constraints!( + container, + ParticipationFractionConstraint, + service, + devices, + model, + ) + end add_to_objective_function!(container, service, model) add_feedforward_constraints!(container, model, service) end @@ -1093,3 +1093,132 @@ function construct_service!( add_constraint_dual!(container, sys, model) return end + +const _RESERVE_UP = Union{PSY.OnlineReserve{PSY.ReserveUp}, PSY.OfflineReserve} + +_is_ramp_formulation(::Type{SecurityConstrainedContingencyReserve}) = false +_is_ramp_formulation(::Type{SecurityConstrainedRampReserve}) = true + +# Whether `service` is procured pre-contingency (reserve variable, requirement, ramp, +# participation, objective); otherwise it only deploys post-contingency. +_is_procured( + model::ServiceModel{<:PSY.AbstractReserve, F}, + service::PSY.AbstractReserve, +) where {F <: AbstractSecurityConstrainedReservesFormulation} = + _is_ramp_formulation(F) || _has_ts_requirement(model, service) + +function construct_service!( + container::OptimizationContainer, + sys::PSY.System, + ::ArgumentConstructStage, + model::ServiceModel{R, <:AbstractSecurityConstrainedReservesFormulation}, + devices_template::Dict{Symbol, DeviceModel}, + ::Set{<:DataType}, + ::NetworkModel{<:AbstractNetworkModel}, +) where {R <: _RESERVE_UP} + services = _services_with_contributors(model, sys) + isempty(services) && return + ts_services = [s for s in services if _has_ts_requirement(model, s)] + isempty(ts_services) || + add_parameters!(container, RequirementTimeSeriesParameter, ts_services, model) + procured = [s for s in services if _is_procured(model, s)] + isempty(procured) || add_service_variables!( + container, + ActivePowerReserveVariable, + procured, + model, + RampReserve, + ) + for service in procured + add_to_expression!( + container, + ActivePowerReserveVariable, + service, + model, + devices_template, + ) + end + for service in services + for devices in values(get_contributing_devices_map(model, PSY.get_name(service))) + _add_post_contingency_deployment!(container, sys, service, devices) + end + add_feedforward_arguments!(container, model, service) + end + return +end + +function construct_service!( + container::OptimizationContainer, + sys::PSY.System, + ::ModelConstructStage, + model::ServiceModel{R, F}, + ::Dict{Symbol, DeviceModel}, + ::Set{<:DataType}, + ::NetworkModel{<:AbstractNetworkModel}, +) where {R <: _RESERVE_UP, F <: AbstractSecurityConstrainedReservesFormulation} + services = _services_with_contributors(model, sys) + isempty(services) && return + for service in services + add_feedforward_constraints!(container, model, service) + end + procured = [s for s in services if _is_procured(model, s)] + if isempty(procured) + add_constraint_dual!(container, sys, model) + return + end + service_names = PSY.get_name.(procured) + add_constraints_container!( + container, + RequirementConstraint, + R, + service_names, + get_time_steps(container), + ) + get_use_slacks(model) && add_reserve_slacks!(container, R, service_names) + for service in procured + contributing_devices = get_contributing_devices_map(model, PSY.get_name(service)) + add_constraints!( + container, + RequirementConstraint, + service, + contributing_devices, + model, + ) + for devices in values(contributing_devices) + # Ramp limits bind spinning reserves only. + _is_ramp_formulation(F) && R <: PSY.Reserve && + add_constraints!(container, RampConstraint, service, devices, model) + add_constraints!( + container, + ParticipationFractionConstraint, + service, + devices, + model, + ) + add_constraints!( + container, + PostContingencyDeploymentConstraint, + service, + devices, + model, + ) + end + add_to_objective_function!(container, service, model) + end + add_constraint_dual!(container, sys, model) + return +end + +construct_service!( + ::OptimizationContainer, + ::PSY.System, + ::Union{ArgumentConstructStage, ModelConstructStage}, + ::ServiceModel{<:PSY.AbstractReserve, <:AbstractSecurityConstrainedReservesFormulation}, + ::Dict{Symbol, DeviceModel}, + ::Set{<:DataType}, + ::NetworkModel{<:AbstractNetworkModel}, +) = throw( + IS.ConflictingInputsError( + "Security-constrained formulations currently only support reserve-up services.", + ), +) diff --git a/src/static_injector_models/hydro_generation.jl b/src/static_injector_models/hydro_generation.jl index f785bd1c..5724a5bc 100644 --- a/src/static_injector_models/hydro_generation.jl +++ b/src/static_injector_models/hydro_generation.jl @@ -2723,7 +2723,8 @@ function add_to_expression!( isa(service, PSY.Reserve{PSY.ReserveUp}) || continue service_name = PSY.get_name(service) deployed_fraction = PSY.get_deployed_fraction(service) - variable = get_variable(container, U, typeof(service)) + variable = + get_variable(container, U, IOM.ComponentPairKey{V, typeof(service)}) for t in get_time_steps(container) add_proportional_to_jump_expression!( expression[name, t], @@ -2762,7 +2763,8 @@ function add_to_expression!( isa(service, PSY.Reserve{PSY.ReserveDown}) || continue service_name = PSY.get_name(service) deployed_fraction = PSY.get_deployed_fraction(service) - variable = get_variable(container, U, typeof(service)) + variable = + get_variable(container, U, IOM.ComponentPairKey{V, typeof(service)}) for t in get_time_steps(container) add_proportional_to_jump_expression!( expression[name, t], @@ -2831,7 +2833,7 @@ function add_to_expression!( W <: AbstractReservesFormulation, } service_name = PSY.get_name(service) - variable = get_variable(container, U, X) + variable = get_variable(container, U, IOM.ComponentPairKey{V, X}) if !has_container_key(container, T, V) add_expressions!(container, T, devices, model) end @@ -2862,7 +2864,7 @@ function add_to_expression!( W <: AbstractReservesFormulation, } service_name = PSY.get_name(service) - variable = get_variable(container, U, X) + variable = get_variable(container, U, IOM.ComponentPairKey{V, X}) if !has_container_key(container, T, V) add_expressions!(container, T, devices, model) end diff --git a/src/static_injector_models/thermal_generation.jl b/src/static_injector_models/thermal_generation.jl index 0406951a..49ea249a 100644 --- a/src/static_injector_models/thermal_generation.jl +++ b/src/static_injector_models/thermal_generation.jl @@ -1783,12 +1783,12 @@ function add_constraints!( # to the device model; the ServiceModel's contributing map carries both. offline = Tuple{String, IOM.JuMPArray, Set{String}}[] for sm in get_services(model) - _is_offline_reserve(get_component_type(sm)) || continue - variable = - get_variable(container, ActivePowerReserveVariable, get_component_type(sm)) + S = get_component_type(sm) + _is_offline_reserve(S) || continue for (service_name, dev_map) in get_contributing_devices_map(sm) members = get(dev_map, V, nothing) isnothing(members) && continue + variable = _reserve_variable(container, V, S) push!(offline, (service_name, variable, Set(PSY.get_name.(members)))) end end diff --git a/test/Project.toml b/test/Project.toml index 8fa666a7..2b934631 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -49,7 +49,8 @@ UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [sources] PowerOperationsModels = {path = ".."} InfrastructureSystems = {url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl.git", rev = "IS4"} -InfrastructureOptimizationModels = {rev = "main", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} +# TEMPORARY: repin to `rev = "main"` once IOM #171 (ac/service-containers) merges. +InfrastructureOptimizationModels = {rev = "ac/service-containers", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} PowerSystems = {url = "https://github.com/Sienna-Platform/PowerSystems.jl.git", rev = "psy6"} PowerSystemCaseBuilder = {url = "https://github.com/Sienna-Platform/PowerSystemCaseBuilder.jl.git", rev = "psy6"} # Branch pins until the OpenAPI packages are registered and the psy6-line branches merge. diff --git a/test/test_device_reserve_offers.jl b/test/test_device_reserve_offers.jl index 59d53f8f..d3e34650 100644 --- a/test/test_device_reserve_offers.jl +++ b/test/test_device_reserve_offers.jl @@ -11,6 +11,18 @@ # (instead of the flat `DEFAULT_RESERVE_COST`). These tests pin the DATA MODEL and assert the # consumer builds the 4D block variable, the award-linking constraint, and the offer-slope cost. +# Awards are stored per (device type, service type): join every device type's WIDE frame for +# `service_type` (the encoded service type, e.g. "OnlineReserve__ReserveUp"). +function _read_awards(res, service_type::String) + frames = [ + read_variable(res, key; table_format = TableFormat.WIDE) for + key in list_variable_names(res) if + startswith(key, "ActivePowerReserveVariable__") && + endswith(key, "__" * service_type) + ] + return reduce((a, b) -> innerjoin(a, b; on = :DateTime), frames) +end + # Give every contributing thermal device of `reserve` a MarketBidCost with an energy offer and a # per-device reserve OFFER curve (PiecewiseStepData, NaturalUnit) named after the service. function add_device_reserve_offers!( @@ -86,9 +98,12 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED container = get_optimization_container(model) - # The reserve award exists, keyed by service type. + # The reserve award exists, keyed by (device type, service type). @test IOM.has_container_key( - container, ActivePowerReserveVariable, OnlineReserve{ReserveUp}) + container, + ActivePowerReserveVariable, + IOM.ComponentPairKey{ThermalStandard, OnlineReserve{ReserveUp}}, + ) # The per-device reserve OFFER is now consumed: a 4D block variable keyed # (service, device, segment, time) exists for the contributing device type, plus the @@ -99,8 +114,11 @@ end container, POM.ReserveOfferLinkingConstraint, ThermalStandard) blk = IOM.get_variable(container, POM.PiecewiseLinearBlockReserveOffer, ThermalStandard) cons = IOM.get_constraint(container, POM.ReserveOfferLinkingConstraint, ThermalStandard) - award = - IOM.get_variable(container, ActivePowerReserveVariable, OnlineReserve{ReserveUp}) + award = IOM.get_variable( + container, + ActivePowerReserveVariable, + IOM.ComponentPairKey{ThermalStandard, OnlineReserve{ReserveUp}}, + ) @test !isempty(blk) sname = PSY.get_name(reserve) @@ -221,9 +239,7 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED res = IOM.OptimizationProblemOutputs(model) - awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; - table_format = TableFormat.WIDE) + awards = _read_awards(res, "OnlineReserve__ReserveUp") # WIDE columns are "__"; values are per-hour reserve awards in MW. col = "$(PSY.get_name(ordc))__$(PSY.get_name(g1))" # The linking constraint caps g1's award at the offered MW each hour; in a dummy hour that cap @@ -277,9 +293,7 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED res = IOM.OptimizationProblemOutputs(model) - awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; - table_format = TableFormat.WIDE) + awards = _read_awards(res, "OnlineReserve__ReserveUp") sname = PSY.get_name(ordc) order = sort(collect(keys(base_slope)); by = n -> base_slope[n]) cheapest, priciest = first(order), last(order) @@ -466,10 +480,7 @@ end res, "ServiceRequirementVariable__GroupReserve__ReserveUp"; table_format = TableFormat.WIDE, ) - awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; - table_format = TableFormat.WIDE, - ) + awards = _read_awards(res, "OnlineReserve__ReserveUp") sub_cols = [c for c in names(awards) if startswith(c, "GROUP_SUB_")] load_col = "GROUP_SUB_A__$(_MKT_LOAD)" @test load_col in names(awards) @@ -526,10 +537,7 @@ end res, "ServiceRequirementVariable__OfflineReserve"; table_format = TableFormat.WIDE, ) - awards = read_variable( - res, "ActivePowerReserveVariable__OfflineReserve"; - table_format = TableFormat.WIDE, - ) + awards = _read_awards(res, "OfflineReserve") load_col = "NSPIN__$(_MKT_LOAD)" @test load_col in names(awards) for t in 1:24 @@ -626,7 +634,11 @@ function _check_offline_band( ) con = IOM.get_constraint(container, POM.OfflineReserveBandConstraint, device_type) varbin = IOM.get_variable(container, POM.OnVariable, device_type) - awards = IOM.get_variable(container, POM.ActivePowerReserveVariable, OfflineReserve) + awards = IOM.get_variable( + container, + POM.ActivePowerReserveVariable, + IOM.ComponentPairKey{device_type, OfflineReserve}, + ) checked = 0 for (idx, c) in con.data name, t = idx @@ -671,14 +683,8 @@ end res = IOM.OptimizationProblemOutputs(model) on = read_variable(res, OnVariable, ThermalStandard; table_format = TableFormat.WIDE) - nspin_awards = read_variable( - res, "ActivePowerReserveVariable__OfflineReserve"; - table_format = TableFormat.WIDE, - ) - spin_awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; - table_format = TableFormat.WIDE, - ) + nspin_awards = _read_awards(res, "OfflineReserve") + spin_awards = _read_awards(res, "OnlineReserve__ReserveUp") # Standard UC: the band row is `p + online + offline <= pmax` for every t, so the # coefficient on `u` is exactly zero. This gates the row's existence and its RHS - # the surviving solve assertions below do not, since each award is separately capped @@ -769,14 +775,8 @@ end res, PowerAboveMinimumVariable, ThermalStandard; table_format = TableFormat.WIDE, ) - nspin_awards = read_variable( - res, "ActivePowerReserveVariable__OfflineReserve"; - table_format = TableFormat.WIDE, - ) - spin_awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; - table_format = TableFormat.WIDE, - ) + nspin_awards = _read_awards(res, "OfflineReserve") + spin_awards = _read_awards(res, "OnlineReserve__ReserveUp") off_name = PSY.get_name(offunit) total_off_award = 0.0 @@ -884,10 +884,7 @@ end res, PowerAboveMinimumVariable, PSY.ThermalMultiStart; table_format = TableFormat.WIDE, ) - awards = read_variable( - res, "ActivePowerReserveVariable__OfflineReserve"; - table_format = TableFormat.WIDE, - ) + awards = _read_awards(res, "OfflineReserve") # Committed multistart units honour the compact band in the solved solution. committed = 0 for d in multistarts @@ -965,10 +962,7 @@ end res, "ActivePowerVariable__InterruptiblePowerLoad"; table_format = TableFormat.WIDE, ) - awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; - table_format = TableFormat.WIDE, - ) + awards = _read_awards(res, "OnlineReserve__ReserveUp") total_award = 0.0 for t in 1:24 awarded = sum(awards[t, c] for c in _il_cols(awards)) @@ -989,10 +983,7 @@ end res, "ActivePowerVariable__InterruptiblePowerLoad"; table_format = TableFormat.WIDE, ) - awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveDown"; - table_format = TableFormat.WIDE, - ) + awards = _read_awards(res, "OnlineReserve__ReserveDown") hsl = read_parameter( res, "ActivePowerTimeSeriesParameter__InterruptiblePowerLoad"; table_format = TableFormat.WIDE, @@ -1026,14 +1017,8 @@ end res, "ActivePowerVariable__InterruptiblePowerLoad"; table_format = TableFormat.WIDE, ) - up = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; - table_format = TableFormat.WIDE, - ) - dn = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveDown"; - table_format = TableFormat.WIDE, - ) + up = _read_awards(res, "OnlineReserve__ReserveUp") + dn = _read_awards(res, "OnlineReserve__ReserveDown") hsl = read_parameter( res, "ActivePowerTimeSeriesParameter__InterruptiblePowerLoad"; table_format = TableFormat.WIDE, @@ -1101,10 +1086,7 @@ end res, "ActivePowerVariable__InterruptiblePowerLoad"; table_format = TableFormat.WIDE, ) - awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; - table_format = TableFormat.WIDE, - ) + awards = _read_awards(res, "OnlineReserve__ReserveUp") combined_total = 0.0 for t in 1:24 # One shared LB expression: a per-service headroom bug would allow up to 2*P. @@ -1154,10 +1136,7 @@ end container, POM.PiecewiseLinearBlockReserveOffer, PSY.InterruptiblePowerLoad, ) res = IOM.OptimizationProblemOutputs(model) - awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; - table_format = TableFormat.WIDE, - ) + awards = _read_awards(res, "OnlineReserve__ReserveUp") col = "$(PSY.get_name(ordc))__$(_IL_NAME)" total = 0.0 for t in 1:24 diff --git a/test/test_model_decision.jl b/test/test_model_decision.jl index 335ae2a0..f8c0865f 100644 --- a/test/test_model_decision.jl +++ b/test/test_model_decision.jl @@ -412,21 +412,29 @@ end # This test needs to be reviewed # @test isapprox(get_objective_value(res), 256937.0; atol = 10000.0) vars = res.variable_values - # Reserve variables of a type share one container keyed + # Reserve variables share one container per (device type, service type), keyed # `(service_name, device_name, time)`. - service_key = IOM.VariableKey( - ActivePowerReserveVariable, - PSY.OfflineReserve, - ) - @test service_key in keys(vars) - # That container flattens to `"service_name__device_name"` result columns + S = IOM.get_component_type( + IOM.VariableKey(ActivePowerReserveVariable, PSY.OfflineReserve), + ) + service_keys = [ + k for k in keys(vars) if + IOM.get_entry_type(k) === ActivePowerReserveVariable && + IOM.get_component_type(k) <: IOM.ComponentPairKey{<:PSY.Component, S} + ] + @test !isempty(service_keys) + # Each container flattens to `"service_name__device_name"` result columns # (WIDE format one column per flattened pair). - result = read_variable( - res, - "ActivePowerReserveVariable__OfflineReserve"; - table_format = TableFormat.WIDE, - ) - @test any(startswith(string(n), "NonSpinningReserve__") for n in names(result)) + result_columns = [ + n for k in service_keys for n in names( + read_variable( + res, + IOM.encode_key_as_string(k); + table_format = TableFormat.WIDE, + ), + ) + ] + @test any(startswith(string(n), "NonSpinningReserve__") for n in result_columns) end @testset "Test serialization/deserialization of DecisionModel outputs" begin diff --git a/test/test_services_constructor.jl b/test/test_services_constructor.jl index d40b0b4c..98225f24 100644 --- a/test/test_services_constructor.jl +++ b/test/test_services_constructor.jl @@ -112,7 +112,7 @@ end end @testset "Per-type reserve container isolates services of the same type" begin - # Two OnlineReserve{ReserveUp} services share one + # Two OnlineReserve{ReserveUp} services share one per-device-type # `(service, device, time)` ActivePowerReserveVariable container. Verify (a) each # service's requirement constraint sums only its own device variables (no # cross-service leakage) and (b) the proportional reserve cost prices each variable @@ -129,7 +129,11 @@ end IOM.ModelBuildStatus.BUILT container = get_optimization_container(model) - rv = IOM.get_variable(container, ActivePowerReserveVariable, OnlineReserve{ReserveUp}) + rv = IOM.get_variable( + container, + ActivePowerReserveVariable, + IOM.ComponentPairKey{ThermalStandard, OnlineReserve{ReserveUp}}, + ) con = IOM.get_constraint( container, RequirementConstraint, @@ -373,7 +377,11 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED container = get_optimization_container(model) - rv = IOM.get_variable(container, ActivePowerReserveVariable, OnlineReserve{ReserveUp}) + rv = IOM.get_variable( + container, + ActivePowerReserveVariable, + IOM.ComponentPairKey{ThermalStandard, OnlineReserve{ReserveUp}}, + ) con = IOM.get_constraint( container, RequirementConstraint, @@ -1317,7 +1325,11 @@ end IOM.ModelBuildStatus.BUILT container = get_optimization_container(model) - rv = IOM.get_variable(container, ActivePowerReserveVariable, OnlineReserve{ReserveUp}) + rv = IOM.get_variable( + container, + ActivePowerReserveVariable, + IOM.ComponentPairKey{ThermalStandard, OnlineReserve{ReserveUp}}, + ) con = IOM.get_constraint( container, RequirementConstraint, @@ -1368,7 +1380,11 @@ end IOM.ModelBuildStatus.BUILT container = get_optimization_container(model) - rv = IOM.get_variable(container, ActivePowerReserveVariable, OnlineReserve{ReserveUp}) + rv = IOM.get_variable( + container, + ActivePowerReserveVariable, + IOM.ComponentPairKey{ThermalStandard, OnlineReserve{ReserveUp}}, + ) con = IOM.get_constraint( container, RequirementConstraint, @@ -1446,7 +1462,11 @@ end @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED container = get_optimization_container(model) - rv = IOM.get_variable(container, ActivePowerReserveVariable, OnlineReserve{ReserveUp}) + rv = IOM.get_variable( + container, + ActivePowerReserveVariable, + IOM.ComponentPairKey{ThermalStandard, OnlineReserve{ReserveUp}}, + ) for t in IOM.get_time_steps(container) provided = sum( @@ -1747,7 +1767,7 @@ end table_format = TableFormat.WIDE, ) awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; + res, "ActivePowerReserveVariable__ThermalStandard__OnlineReserve__ReserveUp"; table_format = TableFormat.WIDE, ) sub_cols = _sub_cols(awards, "GROUP_SUB_") @@ -1762,7 +1782,7 @@ end model = _solve_group_model(sys) res = IOM.OptimizationProblemOutputs(model) awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; + res, "ActivePowerReserveVariable__ThermalStandard__OnlineReserve__ReserveUp"; table_format = TableFormat.WIDE, ) sub_a_total = sum(awards[1, c] for c in _sub_cols(awards, "GROUP_SUB_A")) @@ -1808,7 +1828,7 @@ end model = _solve_group_model(sys; include_group = false) res = IOM.OptimizationProblemOutputs(model) awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; + res, "ActivePowerReserveVariable__ThermalStandard__OnlineReserve__ReserveUp"; table_format = TableFormat.WIDE, ) for t in 1:24, c in _sub_cols(awards, "GROUP_SUB_") @@ -1864,7 +1884,7 @@ end table_format = TableFormat.WIDE, ) awards = read_variable( - res, "ActivePowerReserveVariable__OnlineReserve__ReserveUp"; + res, "ActivePowerReserveVariable__ThermalStandard__OnlineReserve__ReserveUp"; table_format = TableFormat.WIDE, ) sub_cols = _sub_cols(awards, "GROUP_SUB_") @@ -2157,7 +2177,10 @@ end reserve = IOM.get_variable( container, ActivePowerReserveVariable, - PSY.OnlineReserve{PSY.ReserveUp}, + IOM.ComponentPairKey{ + PSY.InterruptiblePowerLoad, + PSY.OnlineReserve{PSY.ReserveUp}, + }, ) for service_name in ("R1", "R2") devices = if nested @@ -2232,16 +2255,18 @@ end container = IOM.get_optimization_container(model) jump_model = IOM.get_jump_model(container) on = IOM.get_variable(container, OnVariable, PSY.InterruptiblePowerLoad) - r_up = - IOM.get_variable( - container, - ActivePowerReserveVariable, - PSY.OnlineReserve{PSY.ReserveUp}, - ) + r_up = IOM.get_variable( + container, + ActivePowerReserveVariable, + IOM.ComponentPairKey{PSY.InterruptiblePowerLoad, PSY.OnlineReserve{PSY.ReserveUp}}, + ) r_dn = IOM.get_variable( container, ActivePowerReserveVariable, - PSY.OnlineReserve{PSY.ReserveDown}, + IOM.ComponentPairKey{ + PSY.InterruptiblePowerLoad, + PSY.OnlineReserve{PSY.ReserveDown}, + }, ) for t in IOM.get_time_steps(container) 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 00000000..8711acfa --- /dev/null +++ b/test/test_static_injection_security_constrained_models.jl @@ -0,0 +1,530 @@ +# G-1 security-constrained reserves (`SecurityConstrainedContingencyReserve`). +# +# One fixture on `two_area_pjm_DA` covers the cases the formulation must keep apart: +# - an OnlineReserve and an OfflineReserve both named "Reserve1_2"; the online one is procured +# (requirement series), the offline one only deploys; +# - a RenewableDispatch named "Brighton_2", like a thermal, contributing to the online reserve; +# - three outages overlapping on "Alta_1": A (online only), B (both reserves), C (offline only); +# - a parallel line and a parallel interchange, monitored or modeled-but-unmonitored. +# The main test builds and solves it under copper plate, PTDF, and area balance, and checks +# the solution against quantities recomputed from the system data. + +const _G1_RESERVE = "Reserve1_2" +const _G1_TOL = 1e-4 + +function _attach_outage!(sys, generators, services, monitored) + outage = PSY.FixedForcedOutage(; + outage_status = 1.0, + monitored_components = monitored, + ) + for generator in generators + add_supplemental_attribute!(sys, generator, outage) + end + for service in services + add_supplemental_attribute!(sys, service, outage) + end + return outage +end + +function _add_parallel_interchange!(sys::PSY.System) + interchange = AreaInterchange(; + name = "1_2_b", + available = true, + active_power_flow = 0.0, + from_area = get_component(Area, sys, "Area1"), + to_area = get_component(Area, sys, "Area2"), + flow_limits = (from_to = 1.5, to_from = 1.5), + ) + add_component!(sys, interchange) + return interchange +end + +_default_g1_monitored(sys) = [ + get_component(Line, sys, "4_1"), + get_component(Line, sys, "2_2"), + get_component(AreaInterchange, sys, "1_2"), +] + +function g1_system(; monitored::Function = _default_g1_monitored) + sys = PSB.build_system(PSISystems, "two_area_pjm_DA"; add_reserves = true) + online = get_component(OnlineReserve{ReserveUp}, sys, _G1_RESERVE) + offline = OfflineReserve(; name = _G1_RESERVE, available = true, time_frame = 30.0) + add_service!(sys, offline, collect(get_components(ThermalStandard, sys))) + + wind = get_component(RenewableDispatch, sys, "WindBus1") + PSY.set_name!(sys, wind, "Brighton_2") + add_service!(wind, online, sys) + + add_equivalent_ac_transmission_with_parallel_circuits!( + sys, + get_component(Line, sys, "4_1"), + Line, + ) + _add_parallel_interchange!(sys) + + thermal(name) = get_component(ThermalStandard, sys, name) + monitored_components = monitored(sys) + outages = ( + a = _attach_outage!( + sys, + [thermal("Alta_1"), thermal("Park City_1")], + [online], + monitored_components, + ), + b = _attach_outage!( + sys, + [thermal("Alta_1"), thermal("Sundance_2")], + [online, offline], + monitored_components, + ), + c = _attach_outage!( + sys, + [thermal("Alta_1"), thermal("Park City_2")], + [offline], + monitored_components, + ), + ) + transform_single_time_series!(sys, Hour(24), Hour(1)) + return sys, outages +end + +function g1_template( + network::Type{<:AbstractNetworkModel}; + online_slacks::Bool = true, + offline_slacks::Bool = true, + model_interchanges::Bool = network === AreaBalanceNetworkModel, + interchange_filter = nothing, +) + template = get_thermal_dispatch_template_network(NetworkModel(network)) + set_device_model!(template, RenewableDispatch, RenewableFullDispatch) + if model_interchanges + attributes = Dict{String, Any}() + isnothing(interchange_filter) || + (attributes["filter_function"] = interchange_filter) + set_device_model!( + template, + DeviceModel(AreaInterchange, StaticBranch; attributes = attributes), + ) + end + set_service_model!( + template, + ServiceModel( + OnlineReserve{ReserveUp}, + SecurityConstrainedContingencyReserve; + use_slacks = online_slacks, + ), + ) + set_service_model!( + template, + ServiceModel( + OfflineReserve, + SecurityConstrainedContingencyReserve; + use_slacks = offline_slacks, + ), + ) + return template +end + +g1_model(template, sys) = + DecisionModel(template, sys; resolution = Hour(1), optimizer = HiGHS_optimizer) + +function aff_exprs_approx_equal(actual::JuMP.AffExpr, expected::JuMP.AffExpr) + isapprox(actual.constant, expected.constant; atol = IOM.PTDF_ZERO_TOL) || return false + return all( + isapprox( + get(actual.terms, v, 0.0), + get(expected.terms, v, 0.0); + atol = IOM.PTDF_ZERO_TOL, + ) for v in union(keys(actual.terms), keys(expected.terms)) + ) +end + +_pair_types(::Type{IOM.ComponentPairKey{D, S}}) where {D, S} = (D, S) + +# Every deployment entry as (device type, service type, service, device, outage, t) => variable. +function _deployments(container) + entries = Dict{Tuple{DataType, Type, String, String, Int, Int}, JuMP.VariableRef}() + for (key, variable) in IOM.get_variables(container) + IOM.get_entry_type(key) === PostContingencyDeploymentVariable || continue + D, S = _pair_types(IOM.get_component_type(key)) + for ((s, d, uuid, t), var) in variable.data + entries[(D, S, s, d, uuid, t)] = var + end + end + return entries +end + +_outaged(sys, outage) = + collect(PSY.get_associated_components(sys, outage; component_type = PSY.Generator)) + +function _power_value(container, generator, t) + power = IOM.get_variable(container, ActivePowerVariable, typeof(generator)) + return JuMP.value(power[PSY.get_name(generator), t]) +end + +# Σ deployment − Σ outaged power per `(outage, t, location(component))`, from the solution. +function _net_deployment(sys, container, deployments, outages, location) + net = Dict{Tuple{Int, Int, Any}, Float64}() + add!(key, value) = (net[key] = get(net, key, 0.0) + value) + for ((D, _, _, d, uuid, t), var) in deployments + add!((uuid, t, location(get_component(D, sys, d))), JuMP.value(var)) + end + for outage in outages, generator in _outaged(sys, outage), + t in IOM.get_time_steps(container) + + add!( + (IS.get_id(outage), t, location(generator)), + -_power_value(container, generator, t), + ) + end + return net +end + +function check_g1_deployment(sys, model, outages) + container = IOM.get_optimization_container(model) + time_steps = IOM.get_time_steps(container) + deployments = _deployments(container) + uuids = Dict(IS.get_id(o) => o for o in values(outages)) + + # Outaged devices do not deploy under their own outage. + outaged = Set( + (typeof(g), PSY.get_name(g), uuid) for (uuid, o) in uuids for g in _outaged(sys, o) + ) + @test !any(((D, _, _, d, uuid, _),) -> (D, d, uuid) in outaged, keys(deployments)) + + # Deployment replaces the outaged generation. + net = _net_deployment(sys, container, deployments, values(outages), _ -> :system) + @test all( + abs(get(net, (uuid, t, :system), 0.0)) < _G1_TOL for uuid in keys(uuids), + t in time_steps + ) + + # Deployment on the procured reserve stays within the award. + award(D, S) = + IOM.get_variable(container, ActivePowerReserveVariable, IOM.ComponentPairKey{D, S}) + @test all( + JuMP.value(var) <= JuMP.value(award(D, S)[s, d, t]) + _G1_TOL for + ((D, S, s, d, _, t), var) in deployments if S <: OnlineReserve + ) + + # Output plus total deployment stays within the device maximum. + totals = Dict{Tuple{DataType, String, Int, Int}, Float64}() + for ((D, _, _, d, uuid, t), var) in deployments + totals[(D, d, uuid, t)] = get(totals, (D, d, uuid, t), 0.0) + JuMP.value(var) + end + @test all( + _power_value(container, get_component(D, sys, d), t) + deployed <= + PSY.get_max_active_power(get_component(D, sys, d), PSY.SU) + _G1_TOL for + ((D, d, _, t), deployed) in totals + ) + return +end + +function check_g1_containers(model, outages) + container = IOM.get_optimization_container(model) + online_key(D) = IOM.ComponentPairKey{D, OnlineReserve{ReserveUp}} + offline_key(D) = IOM.ComponentPairKey{D, OfflineReserve} + + # Same-named contributors of different types stay apart. + for D in (ThermalStandard, RenewableDispatch) + award = IOM.get_variable(container, ActivePowerReserveVariable, online_key(D)) + @test haskey(award.data, (_G1_RESERVE, "Brighton_2", 1)) + end + + # Same-named services of different types stay apart; only the online one responds to A + # and only the offline one to C. + online = IOM.get_variable( + container, + PostContingencyDeploymentVariable, + online_key(ThermalStandard), + ) + offline = IOM.get_variable( + container, + PostContingencyDeploymentVariable, + offline_key(ThermalStandard), + ) + responds(variable, outage) = any(k -> k[3] == IS.get_id(outage), keys(variable.data)) + @test responds(online, outages.a) && !responds(offline, outages.a) + @test responds(online, outages.b) && responds(offline, outages.b) + @test !responds(online, outages.c) && responds(offline, outages.c) + + # The offline reserve has no requirement series, so it is not procured. + @test !IOM.has_container_key(container, RequirementConstraint, OfflineReserve) + @test !IOM.has_container_key( + container, + ActivePowerReserveVariable, + offline_key(ThermalStandard), + ) + @test IOM.has_container_key( + container, + PostContingencyDeploymentConstraint, + online_key(ThermalStandard), + ) + @test !IOM.has_container_key( + container, + PostContingencyDeploymentConstraint, + offline_key(ThermalStandard), + ) + return +end + +# Post-contingency flow change on each monitored entry equals the PTDF of a fresh matrix times +# the net deployment at each bus, recomputed from the solution and the system data. +function check_g1_ptdf_flows(sys, model, outages) + container = IOM.get_optimization_container(model) + network_model = IOM.get_network_model(IOM.get_template(model)) + reduction = POM.get_network_reduction(network_model) + catalog = POM.get_branch_catalog(network_model) + ptdf = PNM.VirtualPTDF(sys) + bus_axis = PNM.get_bus_axis(ptdf) + deployments = _deployments(container) + mapped_bus(component) = PNM.get_mapped_bus_number(reduction, PSY.get_bus(component)) + + net = _net_deployment(sys, container, deployments, values(outages), mapped_bus) + + post_flow = IOM.get_expression(container, PostContingencyBranchFlow, Line, "G1") + pre_flow = IOM.get_expression(container, PTDFBranchFlow, Line) + arc_map = PNM.get_name_to_arc_map(catalog, Line) + rows = Dict(name => ptdf[arc, :] for (name, arc) in arc_map) + flow_change(entry_name, uuid, t) = + sum( + r * get(net, (uuid, t, bus), 0.0) for + (r, bus) in zip(rows[entry_name], bus_axis) + ) + @test !isempty(post_flow.data) + @test all( + isapprox( + JuMP.value(flow) - JuMP.value(pre_flow[entry_name, t]), + flow_change(entry_name, uuid, t); + atol = _G1_TOL, + ) for ((entry_name, uuid, t), flow) in post_flow.data + ) + + # The parallel circuits share one reduced entry, limited once. + entries = PNM.get_component_to_reduction_name_map(catalog, Line) + @test entries["4_1"] == entries["4_1_copy"] + return +end + +function check_g1_area_flows(sys, model, outages) + container = IOM.get_optimization_container(model) + time_steps = IOM.get_time_steps(container) + deployments = _deployments(container) + deviation = IOM.get_variable( + container, + PostContingencyDeviationVariable, + AreaInterchange, + ) + area(c) = PSY.get_name(PSY.get_area(PSY.get_bus(c))) + net = _net_deployment(sys, container, deployments, values(outages), area) + + # Each area balances its net deployment against the deviations on the Area1 -> Area2 + # interchanges. + exported(uuid, t) = sum(JuMP.value(deviation[i, uuid, t]) for i in ("1_2", "1_2_b")) + uuids = [IS.get_id(o) for o in values(outages)] + @test all( + abs(get(net, (uuid, t, "Area1"), 0.0) - exported(uuid, t)) < _G1_TOL && + abs(get(net, (uuid, t, "Area2"), 0.0) + exported(uuid, t)) < _G1_TOL for + uuid in uuids, t in time_steps + ) + + # Both interchanges are modeled, but only the monitored one is limited. + limits = IOM.get_constraint( + container, + PostContingencyFlowRateConstraint, + AreaInterchange, + "G1_ub", + ) + limited = Set(k[1] for k in keys(limits.data)) + @test "1_2" in limited + @test !("1_2_b" in limited) + @test Set(axes(deviation, 1)) == Set(["1_2", "1_2_b"]) + return +end + +@testset "G-1 reserves: $(network)" for network in ( + CopperPlateNetworkModel, + PTDFNetworkModel, + AreaBalanceNetworkModel, +) + sys, outages = g1_system() + model = g1_model(g1_template(network), sys) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + check_g1_containers(model, outages) + @test solve!(model) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + check_g1_deployment(sys, model, outages) + network === PTDFNetworkModel && check_g1_ptdf_flows(sys, model, outages) + network === AreaBalanceNetworkModel && check_g1_area_flows(sys, model, outages) +end + +@testset "G-1 reserves on a degree-two reduced network" begin + # `RADIAL1-RADIAL2-i_1` is merged into a degree-two series chain. Series segments keep + # their own names but map to the chain's reduced arc. + sys = PSB.build_system(PSITestSystems, "case10_radial_series_reductions") + load = first(collect(get_components(StandardLoad, sys))) + add_time_series!( + sys, + load, + Deterministic( + "max_active_power", + Dict(DateTime("2020-01-01T00:00:00") => ones(24)), + Hour(1), + ), + ) + generators = collect(get_components(ThermalStandard, sys)) + reserve = OnlineReserve{ReserveUp}(; + name = _G1_RESERVE, + available = true, + time_frame = 0.0, + requirement = 0.0, + sustained_time = 3600, + max_output_fraction = 1.0, + max_participation_factor = 1.0, + deployed_fraction = 0.0, + ) + add_service!(sys, reserve, generators) + monitored_name = "RADIAL1-RADIAL2-i_1" + outage = _attach_outage!( + sys, + [first(generators)], + [reserve], + [get_component(Line, sys, monitored_name)], + ) + + reductions = PNM.NetworkReduction[DegreeTwoReduction()] + template = PowerOperationsProblemTemplate( + NetworkModel( + PTDFNetworkModel; + network_source = SystemNetworkSource(reductions...), + ), + ) + set_device_model!(template, ThermalStandard, ThermalBasicDispatch) + set_device_model!(template, StandardLoad, StaticPowerLoad) + set_device_model!(template, Line, StaticBranch) + set_service_model!( + template, + ServiceModel(OnlineReserve{ReserveUp}, SecurityConstrainedContingencyReserve), + ) + model = DecisionModel(template, sys; optimizer = HiGHS_optimizer) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + + # Coefficient-level: the post-contingency flow is the pre-contingency flow plus a fresh + # reduced PTDF row times the deployment and outaged power, mapped to reduced buses. + container = IOM.get_optimization_container(model) + network_model = IOM.get_network_model(IOM.get_template(model)) + catalog = POM.get_branch_catalog(network_model) + reduction = POM.get_network_reduction(network_model) + entry_name = PNM.get_component_to_reduction_name_map(catalog, Line)[monitored_name] + reduced_arc = PNM.get_name_to_arc_map(catalog, Line)[entry_name] + line_arc = PSY.get_arc(get_component(Line, sys, monitored_name)) + @test reduced_arc != + (PSY.get_number(PSY.get_from(line_arc)), PSY.get_number(PSY.get_to(line_arc))) + + ptdf = PNM.VirtualPTDF(sys; network_reductions = reductions) + bus_axis = PNM.get_bus_axis(ptdf) + row = ptdf[reduced_arc, :] + factor(component) = + row[findfirst( + ==(PNM.get_mapped_bus_number(reduction, PSY.get_bus(component))), + bus_axis, + )] + + uuid = IS.get_id(outage) + deployments = _deployments(container) + power = IOM.get_variable(container, ActivePowerVariable, ThermalStandard) + post_flow = IOM.get_expression(container, PostContingencyBranchFlow, Line, "G1") + pre_flow = IOM.get_expression(container, PTDFBranchFlow, Line) + for t in IOM.get_time_steps(container) + expected = copy(pre_flow[entry_name, t]) + for ((D, _, _, d, u, tt), var) in deployments + (u == uuid && tt == t) || continue + JuMP.add_to_expression!(expected, factor(get_component(D, sys, d)), var) + end + outaged = first(generators) + JuMP.add_to_expression!( + expected, + -factor(outaged), + power[PSY.get_name(outaged), t], + ) + @test aff_exprs_approx_equal(post_flow[entry_name, uuid, t], expected) + end +end + +@testset "G-1 reserves validation" begin + @testset "monitored interchange excluded by filter_function is rejected" begin + sys, _ = g1_system() + template = g1_template( + AreaBalanceNetworkModel; + interchange_filter = x -> PSY.get_name(x) != "1_2", + ) + @test build!(g1_model(template, sys); output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.FAILED + end + + @testset "monitored interchange without a device model is rejected" begin + sys, _ = g1_system() + template = g1_template(AreaBalanceNetworkModel; model_interchanges = false) + @test build!(g1_model(template, sys); output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.FAILED + end + + @testset "no AreaInterchange model warns and balances each area alone" begin + sys, _ = g1_system(; monitored = s -> [get_component(Line, s, "2_2")]) + model = g1_model( + g1_template(AreaBalanceNetworkModel; model_interchanges = false), + sys, + ) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + container = IOM.get_optimization_container(model) + @test !IOM.has_container_key( + container, + PostContingencyDeviationVariable, + AreaInterchange, + ) + # `build!` logs to the model's log file, not the caller's logger. + log_text = read(IOM.get_log_file(model), String) + @test length( + collect(eachmatch(r"each area must cover its own outages", log_text)), + ) == + 1 + end + + @testset "unmodeled outaged generator warns and is skipped" begin + sys, outages = g1_system() + online = get_component(OnlineReserve{ReserveUp}, sys, _G1_RESERVE) + _attach_outage!( + sys, + [get_component(RenewableDispatch, sys, "PVBus5")], + [online], + [get_component(Line, sys, "2_2")], + ) + template = g1_template(CopperPlateNetworkModel) + delete!(template.devices, :RenewableDispatch) + model = g1_model(template, sys) + @test build!(model; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + log_text = read(IOM.get_log_file(model), String) + @test occursin(r"PVBus5.*is not modeled", log_text) + end + + @testset "service models sharing an outage must agree on use_slacks" begin + sys, _ = g1_system() + template = g1_template(PTDFNetworkModel; offline_slacks = false) + @test build!(g1_model(template, sys); output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.FAILED + end + + @testset "down reserves are rejected" begin + sys, _ = g1_system() + template = g1_template(CopperPlateNetworkModel) + set_service_model!( + template, + ServiceModel(OnlineReserve{ReserveDown}, SecurityConstrainedContingencyReserve), + ) + @test build!(g1_model(template, sys); output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.FAILED + end +end