diff --git a/Project.toml b/Project.toml index c815cb8..9c38d32 100644 --- a/Project.toml +++ b/Project.toml @@ -26,7 +26,8 @@ PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6" [sources] InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl"} -PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} +# TODO: Move to main once this branch is merged +PowerSystems = {rev = "ac/hybridsystem-strip-units", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} InfrastructureOptimizationModels = {rev = "main", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} PowerNetworkMatrices = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"} diff --git a/src/PowerOperationsModels.jl b/src/PowerOperationsModels.jl index b2d47fe..8f4fc70 100644 --- a/src/PowerOperationsModels.jl +++ b/src/PowerOperationsModels.jl @@ -225,6 +225,7 @@ include("core/problem_types.jl") include("core/interfaces.jl") include("core/default_interface_methods.jl") include("core/physical_constant_definitions.jl") +include("core/reserve_traits.jl") include("core/variables.jl") include("core/expressions.jl") include("core/constraints.jl") @@ -303,6 +304,10 @@ include("services_models/reserve_group.jl") include("services_models/transmission_interface.jl") include("services_models/services_constructor.jl") +# Hybrid System Models (after services_models since they share reserve infrastructure) +include("hybrid_system_models/hybrid_systems.jl") +include("hybrid_system_models/hybridsystem_constructor.jl") + # Two-Terminal HVDC Models # NOTE: AC_branches.jl and branch_constructor.jl in twoterminal_hvdc_models/ are # identical copies of the files in ac_transmission_models/ — do NOT include them. @@ -653,6 +658,7 @@ export StorageEnergyOutput export EnergyBalanceConstraint export StateofChargeLimitsConstraint export StateofChargeTargetConstraint +export HybridEnergyTargetConstraint export StorageCyclingCharge export StorageCyclingDischarge export StorageRegularizationConstraintCharge @@ -667,19 +673,44 @@ export ReserveChargeConstraint # expressions export TotalReserveOffering -export ReserveAssignmentBalanceUpDischarge -export ReserveAssignmentBalanceUpCharge -export ReserveAssignmentBalanceDownDischarge -export ReserveAssignmentBalanceDownCharge -export ReserveDeploymentBalanceUpDischarge -export ReserveDeploymentBalanceUpCharge -export ReserveDeploymentBalanceDownDischarge -export ReserveDeploymentBalanceDownCharge # parameters export EnergyLimitParameter export EnergyTargetParameter +######## Hybrid System Formulations ######## +export AbstractHybridFormulation +export AbstractHybridFormulationWithReserves +export HybridDispatchWithReserves + +# Reserve / constraint marker traits used to parametrize hybrid + storage families. +export ReserveScale, UnscaledReserve, DeployedReserve +export ReserveSide, DischargeSide, ChargeSide + +# variables +export HybridEnergyShortageVariable +export HybridEnergySurplusVariable +export HybridRenewableActivePower +export HybridRenewableReserveVariable +export HybridStorageReservation +export HybridThermalActivePower +export HybridThermalReserveVariable + +# expressions + +# constraints +export HybridEnergyAssetBalanceConstraint +export HybridRenewableActivePowerLimitConstraint +export HybridRenewableReserveLimitConstraint +export HybridReserveAssignmentConstraint +export HybridReserveBalanceConstraint +export HybridStorageBalanceConstraint +export HybridThermalReserveLimitConstraint + +# parameters +export HybridElectricLoadTimeSeriesParameter +export HybridRenewableActivePowerTimeSeriesParameter + ################################################################################# # Exports - Constraint Types (defined in core/constraints.jl) ################################################################################# diff --git a/src/core/constraints.jl b/src/core/constraints.jl index 4608f95..4a262c2 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -1064,3 +1064,84 @@ The specified constraints are formulated as: ``` """ struct ShiftDownActivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end + +################################################################################# +# Hybrid System Constraints +################################################################################# + +""" +Couples the hybrid PCC reserve variables (out + in) to the system-level +`ActivePowerReserveVariable` of each service the hybrid participates in. +""" +struct HybridReserveAssignmentConstraint <: ConstraintType end + +""" +Couples the hybrid PCC reserve variables (out + in) to the sum of per-subcomponent reserve +allocations (thermal + renewable + charging + discharging). +""" +struct HybridReserveBalanceConstraint <: ConstraintType end + +""" +Equates the hybrid's PCC active-power injection to the sum of internal subcomponent +flows (thermal + renewable + storage discharge - storage charge - load) net of served +reserves. +""" +struct HybridEnergyAssetBalanceConstraint <: ConstraintType end + +""" +Status link between a hybrid PCC active-power variable and the reservation variable. +Parametric on [`ReserveSide`](@ref). +""" +struct HybridStatusOnConstraint{Sd <: ReserveSide} <: ConstraintType end + +""" +Bound between thermal subcomponent power and its commitment status (no-reserves case). +Parametric on `BoundDirection` (from IOM). +""" +struct HybridThermalOnVariableConstraint{B <: IOM.BoundDirection} <: ConstraintType end + +"Range constraint on thermal subcomponent power including up/down reserves." +struct HybridThermalReserveLimitConstraint <: ConstraintType end + +"Upper bound on renewable subcomponent power from the time-series forecast." +struct HybridRenewableActivePowerLimitConstraint <: ConstraintType end + +"Range constraint on renewable subcomponent power including up/down reserves." +struct HybridRenewableReserveLimitConstraint <: ConstraintType end + +"Energy balance for the storage subcomponent of a hybrid system, including reserve deployment." +struct HybridStorageBalanceConstraint <: ConstraintType end + +""" +Mutually-exclusive charge/discharge limit for the hybrid storage subcomponent +(no-reserves case). Parametric on [`ReserveSide`](@ref). +""" +struct HybridStorageStatusOnConstraint{Sd <: ReserveSide} <: ConstraintType end + +""" +Charge- or discharge-side power limit for the hybrid storage subcomponent including +reserve carve-outs. Parametric on [`ReserveSide`](@ref). +""" +struct HybridStorageReservePowerLimitConstraint{Sd <: ReserveSide} <: ConstraintType end + +""" +Bounds the absolute charge- or discharge-power step change between consecutive time +steps, penalizing oscillation. Active only when the hybrid `\"regularization\"` +attribute is set. Parametric on [`ReserveSide`](@ref). +""" +struct RegularizationConstraint{Sd <: ReserveSide} <: ConstraintType end + +""" +End-of-period energy target for the storage subcomponent of a hybrid system. +Used when the attribute `energy_target = true`. + +Like the storage [`StateofChargeTargetConstraint`](@ref), this is a soft equality with +non-negative surplus ([`HybridEnergySurplusVariable`](@ref), energy above target) and +shortage ([`HybridEnergyShortageVariable`](@ref), energy below target) slacks penalized +in the objective: + +```math +e^{st}_{T} - e^{st+} + e^{st-} = E^{st}_{T}. +``` +""" +struct HybridEnergyTargetConstraint <: ConstraintType end diff --git a/src/core/expressions.jl b/src/core/expressions.jl index 7b3284a..40695b3 100644 --- a/src/core/expressions.jl +++ b/src/core/expressions.jl @@ -78,25 +78,46 @@ of the energy balance for the system in medium term planning struct EnergyBalanceExpression <: ExpressionType end ################################################################################# -# Energy Storage Expressions +# Energy Storage / Hybrid Reserve Aggregation Expressions +# +# A single parametric family covers both the hybrid PCC boundary aggregation +# (HybridPCCReserveExpression) and the storage-subcomponent balance aggregation +# (StorageReserveBalanceExpression). The three axes are: +# D <: ReserveDirection : Up | Down +# S <: ReserveScale : UnscaledReserve (multiplier 1.0) +# | DeployedReserve (multiplier = get_deployed_fraction(s)) +# Sd <: ReserveSide : DischargeSide (PCC "Out" / storage "Discharge") +# | ChargeSide (PCC "In" / storage "Charge") ################################################################################# +""" +Per-device, per-service aggregation of the reserve quantity offered by a storage device +(or the storage subcomponent of a hybrid system). One container is created per service +participated in, and the per-component reserve variables (charging + discharging) are +summed into it. Consumed by [`HybridReserveBalanceConstraint`](@ref) and used as the +right-hand side of the system-level reserve balance. +""" struct TotalReserveOffering <: ExpressionType end -abstract type StorageReserveDischargeExpression <: ExpressionType end -abstract type StorageReserveChargeExpression <: ExpressionType end +abstract type ReserveAggregationExpression{ + D <: PSY.ReserveDirection, + S <: ReserveScale, + Sd <: ReserveSide, +} <: ExpressionType end -# Used for the Power Limits constraints -struct ReserveAssignmentBalanceUpDischarge <: StorageReserveDischargeExpression end -struct ReserveAssignmentBalanceUpCharge <: StorageReserveChargeExpression end -struct ReserveAssignmentBalanceDownDischarge <: StorageReserveDischargeExpression end -struct ReserveAssignmentBalanceDownCharge <: StorageReserveChargeExpression end +""" +Hybrid-boundary aggregation of reserve quantities offered through the discharge (out) and +charge (in) sides of a `PSY.HybridSystem`. +""" +struct HybridPCCReserveExpression{D, S, Sd} <: + ReserveAggregationExpression{D, S, Sd} end -# Used for the SoC estimates -struct ReserveDeploymentBalanceUpDischarge <: StorageReserveDischargeExpression end -struct ReserveDeploymentBalanceUpCharge <: StorageReserveChargeExpression end -struct ReserveDeploymentBalanceDownDischarge <: StorageReserveDischargeExpression end -struct ReserveDeploymentBalanceDownCharge <: StorageReserveChargeExpression end +""" +Aggregation of reserve variables allocated to the storage subcomponent of a hybrid system +(or a standalone storage device). +""" +struct StorageReserveBalanceExpression{D, S, Sd} <: + ReserveAggregationExpression{D, S, Sd} end # Method extensions for output writing should_write_resulting_value(::Type{InterfaceTotalFlow}) = true @@ -108,8 +129,10 @@ should_write_resulting_value(::Type{HydroServedReserveDownExpression}) = true should_write_resulting_value(::Type{TotalHydroFlowRateReservoirOutgoing}) = true should_write_resulting_value(::Type{TotalHydroFlowRateTurbineOutgoing}) = true -should_write_resulting_value(::Type{StorageReserveDischargeExpression}) = true -should_write_resulting_value(::Type{StorageReserveChargeExpression}) = true +should_write_resulting_value(::Type{<:StorageReserveBalanceExpression}) = true +should_write_resulting_value( + ::Type{HybridPCCReserveExpression{D, DeployedReserve, Sd}}, +) where {D <: PSY.ReserveDirection, Sd <: ReserveSide} = true # Method extensions for unit conversion convert_output_to_natural_units(::Type{InterfaceTotalFlow}) = true diff --git a/src/core/formulations.jl b/src/core/formulations.jl index 99ca898..f6b2095 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -450,3 +450,246 @@ The formulation supports the following attributes when used in a [`PowerSimulati See the [`StorageDispatchWithReserves` Mathematical Model](@ref) for the full mathematical description. """ struct StorageDispatchWithReserves <: AbstractStorageFormulation end + +############################ Hybrid System Formulations ################################### +abstract type AbstractHybridFormulation <: IOM.AbstractDeviceFormulation end +abstract type AbstractHybridFormulationWithReserves <: AbstractHybridFormulation end + +""" + HybridDispatchWithReserves + +Device formulation for a hybrid system (single point of common coupling (PCC) with +renewable, thermal, and storage subcomponents) that participates in both energy and +ancillary services markets. Implements a centralized production cost model where the +hybrid plant's net power at the PCC is constrained by ``P_{\\max,\\text{pcc}}`` and +ancillary service allocations (``sb^{\\text{out}}_{p,t}``, ``sb^{\\text{in}}_{p,t}``) are +assigned to internal assets (thermal, renewable, charge, discharge) per the +four-quadrant ancillary service model. Reserve participation is enabled by attaching a +service model to the hybrid (`set_service_model!` + `add_service!`); when no service is +attached the formulation collapses to an energy-only hybrid dispatch. + +Use with a hybrid system in a [`DeviceModel`](@ref) for unit commitment or economic +dispatch. + +**Variables:** + + - [`ActivePowerOutVariable`](@ref): + + + Domain: [0.0, ``P_{\\max,\\text{pcc}}``] + + Symbol: ``p^{\\text{out}}_t`` + + - [`ActivePowerInVariable`](@ref): + + + Domain: [0.0, ``P_{\\max,\\text{pcc}}``] + + Symbol: ``p^{\\text{in}}_t`` + + - [`ReservationVariable`](@ref) (only when `"reservation" => true`): + + + Domain: {0, 1} + + Symbol: ``u^{\\text{st}}_t`` (1 = discharge mode, 0 = charge mode) + + - [`HybridThermalActivePower`](@ref): + + + Domain: [0.0, ``P_{\\max,\\text{th}}``] when on + + Symbol: ``p^{\\text{th}}_t`` + + - [`OnVariable`](@ref): + + + Domain: {0, 1} + + Symbol: ``u^{\\text{th}}_t`` + + - [`HybridRenewableActivePower`](@ref): + + + Domain: [0.0, ``P^{*,\\text{re}}_t``] + + Symbol: ``p^{\\text{re}}_t`` + + - [`HybridStorageSubcomponentPower{ChargeSide}`](@ref): + + + Domain: [0.0, ``P_{\\max,\\text{ch}}``] + + Symbol: ``p^{\\text{ch}}_t`` + + - [`HybridStorageSubcomponentPower{DischargeSide}`](@ref): + + + Domain: [0.0, ``P_{\\max,\\text{ds}}``] + + Symbol: ``p^{\\text{ds}}_t`` + + - [`EnergyVariable`](@ref): + + + Domain: [0.0, ``E_{\\max,\\text{st}}``] + + Symbol: ``e^{\\text{st}}_t`` + + - [`HybridStorageReservation`](@ref) (only when `"storage_reservation" => true`): + + + Domain: {0, 1} + + Symbol: ``ss^{\\text{st}}_t`` (0 = charge, 1 = discharge) + + - [`HybridPCCReserveVariable{DischargeSide}`](@ref) (only when services are attached): + + + Domain: [0.0, ] + + Symbol: ``sb^{\\text{out}}_t`` + + - [`HybridPCCReserveVariable{ChargeSide}`](@ref) (only when services are attached): + + + Domain: [0.0, ] + + Symbol: ``sb^{\\text{in}}_t`` + + - [`RegularizationVariable{ChargeSide}`](@ref), [`RegularizationVariable{DischargeSide}`](@ref) + (only when `"regularization" => true`): non-negative slacks bounding step changes in + charge/discharge between consecutive time steps. + +**Time Series Parameters:** + +| Parameter | Default Time Series Name | +| :--- | :--- | +| `HybridRenewableActivePowerTimeSeriesParameter` | `"RenewableDispatch__max_active_power"` | +| `HybridElectricLoadTimeSeriesParameter` | `"PowerLoad__max_active_power"` | + +**Data requirements:** + + - **Device:** A `PSY.HybridSystem` with at least one of: thermal unit + (`PSY.get_thermal_unit`), renewable unit (`PSY.get_renewable_unit`), storage + (`PSY.get_storage`), and optionally electric load (`PSY.get_electric_load`). + - **Time series:** Forecast time series must be attached to the `PSY.HybridSystem` + itself (not its subcomponents) under the default names above (or custom names passed + when adding parameters). The subcomponent-namespaced default names + (`"RenewableDispatch__max_active_power"`, `"PowerLoad__max_active_power"`) reflect + which subcomponent each forecast describes; the subcomponent is consulted only for the + rating used to scale the parameter. + +**Static Parameters:** + + - ``P_{\\max,\\text{pcc}}`` = `PSY.get_output_active_power_limits(device).max` + - ``P_{\\max,\\text{th}}`` = `PSY.get_active_power_limits(thermal_unit).max` + - ``P_{\\min,\\text{th}}`` = `PSY.get_active_power_limits(thermal_unit).min` + - ``P_{\\max,\\text{ch}}`` = `PSY.get_input_active_power_limits(storage).max` + - ``P_{\\max,\\text{ds}}`` = `PSY.get_output_active_power_limits(storage).max` + - ``\\eta_{\\text{ch}}`` = `PSY.get_efficiency(storage).in` + - ``\\eta_{\\text{ds}}`` = `PSY.get_efficiency(storage).out` + - ``E_{\\max,\\text{st}}`` = `PSY.get_storage_level_limits(storage).max * PSY.get_storage_capacity(storage, PSY.SU) * PSY.get_conversion_factor(storage)`` + - ``E^{\\text{st}}_0`` = initial storage energy + - ``R^{*}_{p,t}`` = ancillary service deployment forecast for service ``p`` at time ``t`` + - ``F_p`` = fraction of ``P_{\\max,\\text{pcc}}`` allowed for service ``p`` + - ``N_p`` = number of periods of compliance for service ``p`` + +**Expressions:** + +Adds ``p^{\\text{out}}_t`` and ``p^{\\text{in}}_t`` to `ActivePowerBalance` for use in +network balance constraints. When services are attached, also accumulates reserve +expressions (`HybridPCCReserveExpression`) with unscaled and deployed-reserve scalings +across all four combinations of direction (up/down) and side (in/out). + +**Constraints:** + +Let ``\\mathcal{T} = \\{1, \\dots, T\\}`` denote the set of time steps. + +PCC and status. When `"reservation" => true`: +[`HybridStatusOnConstraint{DischargeSide}`](@ref), [`HybridStatusOnConstraint{ChargeSide}`](@ref). When +`"reservation" => false`: [`OutputActivePowerVariableLimitsConstraint`](@ref) and +[`InputActivePowerVariableLimitsConstraint`](@ref) (no mutual-exclusion binary). + +```math +\\begin{align*} +& 0 \\leq p^{\\text{in}}_t \\leq P_{\\max,\\text{pcc}}, \\quad 0 \\leq p^{\\text{out}}_t \\leq P_{\\max,\\text{pcc}}, \\quad \\forall t \\in \\mathcal{T} \\\\ +& u^{\\text{st}}_t \\in \\{0,1\\} \\quad \\text{(when reservation is enabled)} +\\end{align*} +``` + +Energy asset balance ([`HybridEnergyAssetBalanceConstraint`](@ref)). When services are +present, served-reserve expressions enter the balance with sign pattern +``+\\bar{r}^{\\text{out,up}} - \\bar{r}^{\\text{in,up}} - \\bar{r}^{\\text{out,dn}} + \\bar{r}^{\\text{in,dn}}``. + +```math +p^{\\text{th}}_t + p^{\\text{re}}_t + p^{\\text{ds}}_t - p^{\\text{ch}}_t - P^{\\text{ld}}_t = p^{\\text{out}}_t - p^{\\text{in}}_t, \\quad \\forall t \\in \\mathcal{T} +``` + +Thermal limits when no services are attached +([`HybridThermalOnVariableConstraint{UpperBound}`](@ref), +[`HybridThermalOnVariableConstraint{LowerBound}`](@ref)): + +```math +u^{\\text{th}}_t P_{\\min,\\text{th}} \\leq p^{\\text{th}}_t \\leq u^{\\text{th}}_t P_{\\max,\\text{th}}, \\quad u^{\\text{th}}_t \\in \\{0,1\\}, \\quad \\forall t \\in \\mathcal{T} +``` + +Renewable limit ([`HybridRenewableActivePowerLimitConstraint`](@ref)): + +```math +0 \\leq p^{\\text{re}}_t \\leq P^{*,\\text{re}}_t, \\quad \\forall t \\in \\mathcal{T} +``` + +Storage charge/discharge mutual exclusion when `"storage_reservation" => true` +([`HybridStorageStatusOnConstraint{ChargeSide}`](@ref), +[`HybridStorageStatusOnConstraint{DischargeSide}`](@ref)): + +```math +\\begin{align*} +& p^{\\text{ch}}_t \\leq (1 - ss^{\\text{st}}_t) P_{\\max,\\text{ch}}, \\quad p^{\\text{ds}}_t \\leq ss^{\\text{st}}_t P_{\\max,\\text{ds}}, \\quad \\forall t \\in \\mathcal{T} \\\\ +& ss^{\\text{st}}_t \\in \\{0,1\\} +\\end{align*} +``` + +Storage energy balance ([`HybridStorageBalanceConstraint`](@ref)): + +```math +e^{\\text{st}}_t = e^{\\text{st}}_{t-1} + \\Delta t \\left( \\eta_{\\text{ch}} p^{\\text{ch}}_t - \\frac{p^{\\text{ds}}_t}{\\eta_{\\text{ds}}} \\right), \\quad \\forall t \\in \\mathcal{T}, \\quad e^{\\text{st}}_0 = E^{\\text{st}}_0 +``` + +When ancillary services are attached: [`HybridThermalReserveLimitConstraint`](@ref), +[`HybridRenewableReserveLimitConstraint`](@ref), +[`HybridStorageReservePowerLimitConstraint{ChargeSide}`](@ref), +[`HybridStorageReservePowerLimitConstraint{DischargeSide}`](@ref), +[`ReserveCoverageConstraint`](@ref), [`ReserveCoverageConstraintEndOfPeriod`](@ref), +[`HybridReserveAssignmentConstraint`](@ref), [`HybridReserveBalanceConstraint`](@ref). + +End-of-horizon energy target (if `"energy_target" => true`), +[`HybridEnergyTargetConstraint`](@ref): + +```math +e^{\\text{st}}_T - e^{\\text{st+}} + e^{\\text{st-}} = E^{\\text{st}}_T +``` + +Charge/discharge regularization (if `"regularization" => true`), +[`RegularizationConstraint{ChargeSide}`](@ref), +[`RegularizationConstraint{DischargeSide}`](@ref): bound ``|p^{\\text{ch}}_t - +p^{\\text{ch}}_{t-1}|`` and ``|p^{\\text{ds}}_t - p^{\\text{ds}}_{t-1}|`` by a +non-negative slack carried into the objective. + +# Example + +```julia +DeviceModel( + PSY.HybridSystem, + HybridDispatchWithReserves; + attributes = Dict( + "reservation" => true, + "storage_reservation" => true, + "energy_target" => false, + "regularization" => false, + ), +) +``` + +# Attributes + + - `"reservation"` (default `true`): if `true`, adds `ReservationVariable` and uses + `HybridStatus{Out,In}OnConstraint` to mutually exclude PCC charge and discharge. + If `false`, both PCC variables are bounded by simple range constraints. + - `"storage_reservation"` (default `true`): if `true`, adds `HybridStorageReservation` + and uses the `ss`-multiplied form of the storage power-limit constraints. If + `false`, charge and discharge variables are bounded independently. + - `"energy_target"` (default `false`): adds `HybridEnergyTargetConstraint` at the + storage subcomponent. + - `"regularization"` (default `false`): adds `RegularizationVariable{ChargeSide}` and + `RegularizationVariable{DischargeSide}` plus the matching constraints, and a small + objective penalty on each, to suppress charge/discharge oscillation. + +**Objective:** + +Adds variable cost on `HybridThermalActivePower`, `HybridRenewableActivePower`, +`HybridStorageSubcomponentPower{ChargeSide}`, and `HybridStorageSubcomponentPower{DischargeSide}` from each subcomponent's +`PSY.get_operation_cost`, plus the proportional `OnVariable` cost (delegated to POM's +standard `proportional_cost` for `ThermalGenerationCost`, so a hybrid-embedded thermal +unit and a standalone copy produce identical objective coefficients). When +`"regularization" => true`, also adds a small per-time-step penalty on the +regularization slacks. +""" +struct HybridDispatchWithReserves <: AbstractHybridFormulationWithReserves end diff --git a/src/core/parameters.jl b/src/core/parameters.jl index 83f6df9..a81c09e 100644 --- a/src/core/parameters.jl +++ b/src/core/parameters.jl @@ -222,6 +222,18 @@ Parameter to record that the component changed in the availability status """ struct AvailableStatusChangeCountdownParameter <: EventParameter end +################################################################################# +# Hybrid System Parameters +################################################################################# + +"Time-series parameter for the maximum active power available from a hybrid system's +renewable subcomponent, normalized by the renewable unit's `max_active_power` rating." +struct HybridRenewableActivePowerTimeSeriesParameter <: TimeSeriesParameter end + +"Time-series parameter for a hybrid system's electric-load subcomponent demand, +normalized by the load's `max_active_power`." +struct HybridElectricLoadTimeSeriesParameter <: TimeSeriesParameter end + ################################################################################# # Method extensions for should_write_resulting_value ################################################################################# @@ -256,3 +268,6 @@ convert_output_to_natural_units(::Type{EnergyTargetTimeSeriesParameter}) = true convert_output_to_natural_units(::Type{EnergyBudgetTimeSeriesParameter}) = true convert_output_to_natural_units(::Type{InflowTimeSeriesParameter}) = false convert_output_to_natural_units(::Type{OutflowTimeSeriesParameter}) = false +convert_output_to_natural_units(::Type{HybridRenewableActivePowerTimeSeriesParameter}) = + true +convert_output_to_natural_units(::Type{HybridElectricLoadTimeSeriesParameter}) = true diff --git a/src/core/reserve_traits.jl b/src/core/reserve_traits.jl new file mode 100644 index 0000000..5ce88b1 --- /dev/null +++ b/src/core/reserve_traits.jl @@ -0,0 +1,16 @@ +# Marker singleton trait types used to parametrize hybrid/storage reserve variable, +# expression, and constraint families. These eliminate the need for paired sibling +# singletons across the codebase: a single parametric struct is used instead of +# every (Charge/Discharge) and (Unscaled/Deployed) sibling pair. + +abstract type ReserveScale end +"Reserve aggregation that uses the raw multiplier (1.0). Was Total / Assignment." +struct UnscaledReserve <: ReserveScale end +"Reserve aggregation that scales the multiplier by deployed_fraction. Was Served / Deployment." +struct DeployedReserve <: ReserveScale end + +abstract type ReserveSide end +"Discharge / outflow side of a storage or hybrid PCC. Was Out (PCC) / Discharge (storage)." +struct DischargeSide <: ReserveSide end +"Charge / inflow side of a storage or hybrid PCC. Was In (PCC) / Charge (storage)." +struct ChargeSide <: ReserveSide end diff --git a/src/core/variables.jl b/src/core/variables.jl index 7ee8fe0..4f5b72d 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -579,6 +579,90 @@ Auxiliary Variable for Storage Models that solve for total energy output """ struct StorageEnergyOutput <: AuxVariableType end +################################################################################# +# Hybrid System Variables +# +# Paired sibling variable types are parametric on `ReserveSide` (Discharge / Charge). +################################################################################# + +""" +Abstract type for variables representing flows internal to a `PSY.HybridSystem`. +""" +abstract type AbstractHybridSubcomponentVariableType <: VariableType end + +"Active power dispatched by the thermal subcomponent of a hybrid system." +struct HybridThermalActivePower <: AbstractHybridSubcomponentVariableType end + +"Active power dispatched by the renewable subcomponent of a hybrid system." +struct HybridRenewableActivePower <: AbstractHybridSubcomponentVariableType end + +""" +Active power on the storage subcomponent of a hybrid system. Parametric on +[`ReserveSide`](@ref). +""" +struct HybridStorageSubcomponentPower{Sd <: ReserveSide} <: + AbstractHybridSubcomponentVariableType end + +"Binary reservation variable for the storage subcomponent of a hybrid system." +struct HybridStorageReservation <: AbstractHybridSubcomponentVariableType end + +""" +Non-negative slack variable bounding the absolute step change in charge or discharge +power between consecutive time steps. Carried into the objective with a small fixed +penalty when the hybrid `\"regularization\"` attribute is set. +""" +struct RegularizationVariable{Sd <: ReserveSide} <: AbstractHybridSubcomponentVariableType end + +""" +Slack variable for the storage energy of a hybrid system being below its end-of-period +target. Added when the hybrid `\"energy_target\"` attribute is set and penalized in the +objective by the storage subcomponent's `energy_shortage_cost`. + +Docs abbreviation: ``e^{st-}`` +""" +struct HybridEnergyShortageVariable <: VariableType end + +""" +Slack variable for the storage energy of a hybrid system being above its end-of-period +target. Added when the hybrid `\"energy_target\"` attribute is set and penalized in the +objective by the storage subcomponent's `energy_surplus_cost`. + +Docs abbreviation: ``e^{st+}`` +""" +struct HybridEnergySurplusVariable <: VariableType end + +""" +Abstract type for hybrid reserve variables (both PCC-boundary and subcomponent). +""" +abstract type AbstractHybridReserveVariableType <: VariableType end + +""" +Reserve quantity offered to the grid through one side of a hybrid PCC. Parametric on +[`ReserveSide`](@ref). +""" +struct HybridPCCReserveVariable{Sd <: ReserveSide} <: AbstractHybridReserveVariableType end + +""" +Abstract type for per-subcomponent reserve allocations inside a hybrid system +that do not have a Discharge/Charge axis (thermal and renewable subcomponents). +""" +abstract type AbstractHybridSubcomponentInjectorReserveVariableType <: + AbstractHybridReserveVariableType end + +"Reserve allocated to the thermal subcomponent of a hybrid system." +struct HybridThermalReserveVariable <: AbstractHybridSubcomponentInjectorReserveVariableType end + +"Reserve allocated to the renewable subcomponent of a hybrid system." +struct HybridRenewableReserveVariable <: + AbstractHybridSubcomponentInjectorReserveVariableType end + +""" +Reserve allocated to one side of a hybrid system's storage subcomponent. Parametric on +[`ReserveSide`](@ref). +""" +struct HybridStorageSubcomponentReserveVariable{Sd <: ReserveSide} <: + AbstractHybridReserveVariableType end + const MULTI_START_VARIABLES = (HotStartVariable, WarmStartVariable, ColdStartVariable) should_write_resulting_value(::Type{PiecewiseLinearCostVariable}) = false diff --git a/src/energy_storage_models/storage_constructor.jl b/src/energy_storage_models/storage_constructor.jl index 3011e04..ce2c010 100644 --- a/src/energy_storage_models/storage_constructor.jl +++ b/src/energy_storage_models/storage_constructor.jl @@ -9,14 +9,14 @@ function _add_ancillary_services!( add_variables!(container, AncillaryServiceVariableCharge, devices, U) time_steps = get_time_steps(container) for exp in [ - ReserveAssignmentBalanceUpDischarge, - ReserveAssignmentBalanceUpCharge, - ReserveAssignmentBalanceDownDischarge, - ReserveAssignmentBalanceDownCharge, - ReserveDeploymentBalanceUpDischarge, - ReserveDeploymentBalanceUpCharge, - ReserveDeploymentBalanceDownDischarge, - ReserveDeploymentBalanceDownCharge, + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide}, + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, DischargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide}, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, DischargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}, ] lazy_container_addition!( container, @@ -27,10 +27,10 @@ function _add_ancillary_services!( ) end for exp in [ - ReserveAssignmentBalanceUpDischarge, - ReserveAssignmentBalanceDownDischarge, - ReserveDeploymentBalanceUpDischarge, - ReserveDeploymentBalanceDownDischarge, + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, DischargeSide}, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, DischargeSide}, ] add_to_expression!( container, @@ -41,10 +41,10 @@ function _add_ancillary_services!( ) end for exp in [ - ReserveAssignmentBalanceUpCharge, - ReserveAssignmentBalanceDownCharge, - ReserveDeploymentBalanceUpCharge, - ReserveDeploymentBalanceDownCharge, + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide}, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}, ] add_to_expression!(container, exp, AncillaryServiceVariableCharge, devices, model) end diff --git a/src/energy_storage_models/storage_models.jl b/src/energy_storage_models/storage_models.jl index bfb21ad..be2cf1b 100644 --- a/src/energy_storage_models/storage_models.jl +++ b/src/energy_storage_models/storage_models.jl @@ -209,13 +209,13 @@ end # For InputActivePower (charge), it's `P_in + down - up` — reserves swap roles because # a charging battery's net power is increased by downward reserves. _deployment_increasing_expr(::Type{<:OutputActivePowerVariableLimitsConstraint}) = - ReserveDeploymentBalanceUpDischarge + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide} _deployment_decreasing_expr(::Type{<:OutputActivePowerVariableLimitsConstraint}) = - ReserveDeploymentBalanceDownDischarge + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, DischargeSide} _deployment_increasing_expr(::Type{<:InputActivePowerVariableLimitsConstraint}) = - ReserveDeploymentBalanceDownCharge + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide} _deployment_decreasing_expr(::Type{<:InputActivePowerVariableLimitsConstraint}) = - ReserveDeploymentBalanceUpCharge + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, ChargeSide} # Reservation-binary handling: discharge active when ss=1, charge active when ss=0. _reservation_factor(::Type{<:OutputActivePowerVariableLimitsConstraint}, ss, name, t) = @@ -433,7 +433,7 @@ end ############################# Expression Logic for Ancillary Services ###################### get_variable_multiplier( ::Type{AncillaryServiceVariableCharge}, - ::Type{ReserveAssignmentBalanceDownCharge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveUp}, @@ -441,7 +441,7 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableCharge}, - ::Type{ReserveAssignmentBalanceDownCharge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveDown}, @@ -449,7 +449,7 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableCharge}, - ::Type{ReserveAssignmentBalanceUpCharge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveUp}, @@ -457,7 +457,7 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableCharge}, - ::Type{ReserveAssignmentBalanceUpCharge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveDown}, @@ -465,7 +465,9 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableDischarge}, - ::Type{ReserveAssignmentBalanceDownDischarge}, + ::Type{ + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, DischargeSide}, + }, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveUp}, @@ -473,7 +475,9 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableDischarge}, - ::Type{ReserveAssignmentBalanceDownDischarge}, + ::Type{ + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, DischargeSide}, + }, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveDown}, @@ -481,7 +485,7 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableDischarge}, - ::Type{ReserveAssignmentBalanceUpDischarge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveUp}, @@ -489,7 +493,7 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableDischarge}, - ::Type{ReserveAssignmentBalanceUpDischarge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveDown}, @@ -498,7 +502,7 @@ get_variable_multiplier( ### Deployment ### get_variable_multiplier( ::Type{AncillaryServiceVariableCharge}, - ::Type{ReserveDeploymentBalanceDownCharge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveUp}, @@ -506,7 +510,7 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableCharge}, - ::Type{ReserveDeploymentBalanceDownCharge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveDown}, @@ -514,7 +518,7 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableCharge}, - ::Type{ReserveDeploymentBalanceUpCharge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveUp}, @@ -522,7 +526,7 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableCharge}, - ::Type{ReserveDeploymentBalanceUpCharge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveDown}, @@ -530,7 +534,9 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableDischarge}, - ::Type{ReserveDeploymentBalanceDownDischarge}, + ::Type{ + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, DischargeSide}, + }, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveUp}, @@ -538,7 +544,9 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableDischarge}, - ::Type{ReserveDeploymentBalanceDownDischarge}, + ::Type{ + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, DischargeSide}, + }, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveDown}, @@ -546,7 +554,7 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableDischarge}, - ::Type{ReserveDeploymentBalanceUpDischarge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveUp}, @@ -554,7 +562,7 @@ get_variable_multiplier( get_variable_multiplier( ::Type{AncillaryServiceVariableDischarge}, - ::Type{ReserveDeploymentBalanceUpDischarge}, + ::Type{StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}}, d::PSY.Storage, ::Type{StorageDispatchWithReserves}, ::PSY.Reserve{PSY.ReserveDown}, @@ -562,16 +570,16 @@ get_variable_multiplier( #! format: off # Use 1.0 because this is to allow to reuse the code below on add_to_expression -get_fraction(::Type{ReserveAssignmentBalanceUpDischarge}, d::PSY.Reserve) = 1.0 -get_fraction(::Type{ReserveAssignmentBalanceUpCharge}, d::PSY.Reserve) = 1.0 -get_fraction(::Type{ReserveAssignmentBalanceDownDischarge}, d::PSY.Reserve) = 1.0 -get_fraction(::Type{ReserveAssignmentBalanceDownCharge}, d::PSY.Reserve) = 1.0 +get_fraction(::Type{StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide}}, d::PSY.Reserve) = 1.0 +get_fraction(::Type{StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide}}, d::PSY.Reserve) = 1.0 +get_fraction(::Type{StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, DischargeSide}}, d::PSY.Reserve) = 1.0 +get_fraction(::Type{StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide}}, d::PSY.Reserve) = 1.0 # Needs to implement served fraction in PSY -get_fraction(::Type{ReserveDeploymentBalanceUpDischarge}, d::PSY.Reserve) = PSY.get_deployed_fraction(d) -get_fraction(::Type{ReserveDeploymentBalanceUpCharge}, d::PSY.Reserve) = PSY.get_deployed_fraction(d) -get_fraction(::Type{ReserveDeploymentBalanceDownDischarge}, d::PSY.Reserve) = PSY.get_deployed_fraction(d) -get_fraction(::Type{ReserveDeploymentBalanceDownCharge}, d::PSY.Reserve) = PSY.get_deployed_fraction(d) +get_fraction(::Type{StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}}, d::PSY.Reserve) = PSY.get_deployed_fraction(d) +get_fraction(::Type{StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}}, d::PSY.Reserve) = PSY.get_deployed_fraction(d) +get_fraction(::Type{StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, DischargeSide}}, d::PSY.Reserve) = PSY.get_deployed_fraction(d) +get_fraction(::Type{StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}}, d::PSY.Reserve) = PSY.get_deployed_fraction(d) #! format: on function add_to_expression!( @@ -619,7 +627,8 @@ function add_to_expression!( devices::IS.FlattenIteratorWrapper{V}, model::DeviceModel{V, W}, ) where { - T <: StorageReserveChargeExpression, + T <: + StorageReserveBalanceExpression{<:PSY.ReserveDirection, <:ReserveScale, ChargeSide}, U <: AncillaryServiceVariableCharge, V <: PSY.Storage, W <: StorageDispatchWithReserves, @@ -651,7 +660,8 @@ function add_to_expression!( devices::IS.FlattenIteratorWrapper{V}, model::DeviceModel{V, W}, ) where { - T <: StorageReserveDischargeExpression, + T <: + StorageReserveBalanceExpression{<:PSY.ReserveDirection, <:ReserveScale, DischargeSide}, U <: AncillaryServiceVariableDischarge, V <: PSY.Storage, W <: StorageDispatchWithReserves, @@ -769,10 +779,26 @@ function add_energybalance_with_reserves!( powerin_var = get_variable(container, ActivePowerInVariable, V) powerout_var = get_variable(container, ActivePowerOutVariable, V) - r_up_ds = get_expression(container, ReserveDeploymentBalanceUpDischarge, V) - r_up_ch = get_expression(container, ReserveDeploymentBalanceUpCharge, V) - r_dn_ds = get_expression(container, ReserveDeploymentBalanceDownDischarge, V) - r_dn_ch = get_expression(container, ReserveDeploymentBalanceDownCharge, V) + r_up_ds = get_expression( + container, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}, + V, + ) + r_up_ch = get_expression( + container, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}, + V, + ) + r_dn_ds = get_expression( + container, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, DischargeSide}, + V, + ) + r_dn_ch = get_expression( + container, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}, + V, + ) constraint = add_constraints_container!(container, EnergyBalanceConstraint, V, @@ -878,13 +904,13 @@ end _reserve_assignment_power_var(::Type{ReserveDischargeConstraint}) = ActivePowerOutVariable _reserve_assignment_power_var(::Type{ReserveChargeConstraint}) = ActivePowerInVariable _reserve_assignment_up_expr(::Type{ReserveDischargeConstraint}) = - ReserveAssignmentBalanceUpDischarge + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide} _reserve_assignment_down_expr(::Type{ReserveDischargeConstraint}) = - ReserveAssignmentBalanceDownDischarge + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, DischargeSide} _reserve_assignment_up_expr(::Type{ReserveChargeConstraint}) = - ReserveAssignmentBalanceUpCharge + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide} _reserve_assignment_down_expr(::Type{ReserveChargeConstraint}) = - ReserveAssignmentBalanceDownCharge + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide} _reserve_assignment_limits(::Type{ReserveDischargeConstraint}, d) = PSY.get_output_active_power_limits(d, PSY.SU) _reserve_assignment_limits(::Type{ReserveChargeConstraint}, d) = @@ -1346,7 +1372,11 @@ function add_cycling_charge_with_reserves!( powerin_var = get_variable(container, ActivePowerInVariable, V) slack_var = get_variable(container, StorageChargeCyclingSlackVariable, V) - r_dn_ch = get_expression(container, ReserveDeploymentBalanceDownCharge, V) + r_dn_ch = get_expression( + container, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}, + V, + ) constraint = add_constraints_container!(container, StorageCyclingCharge, V, names) @@ -1435,7 +1465,11 @@ function add_cycling_discharge_with_reserves!( names = [PSY.get_name(x) for x in devices] powerout_var = get_variable(container, ActivePowerOutVariable, V) slack_var = get_variable(container, StorageDischargeCyclingSlackVariable, V) - r_up_ds = get_expression(container, ReserveDeploymentBalanceUpDischarge, V) + r_up_ds = get_expression( + container, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}, + V, + ) constraint = add_constraints_container!(container, StorageCyclingDischarge, V, names) @@ -1475,148 +1509,77 @@ function add_constraints!( return end -function add_constraints!( - container::OptimizationContainer, - ::Type{StorageRegularizationConstraintCharge}, - devices::IS.FlattenIteratorWrapper{V}, - model::DeviceModel{V, StorageDispatchWithReserves}, - network_model::NetworkModel{X}, -) where {V <: PSY.Storage, X <: AbstractPowerModel} - names = [PSY.get_name(x) for x in devices] - time_steps = get_time_steps(container) - reg_var = get_variable(container, StorageRegularizationVariableCharge, V) - powerin_var = get_variable(container, ActivePowerInVariable, V) - has_services = has_service_model(model) - - if has_services - r_up_ch = get_expression(container, ReserveDeploymentBalanceUpCharge, V) - r_dn_ch = get_expression(container, ReserveDeploymentBalanceDownCharge, V) - end - - constraint_ub = - add_constraints_container!(container, StorageRegularizationConstraintCharge, - V, - names, - time_steps; - meta = "ub", - ) - - constraint_lb = - add_constraints_container!(container, StorageRegularizationConstraintCharge, - V, - names, - time_steps; - meta = "lb", - ) - - for d in devices - name = PSY.get_name(d) - constraint_ub[name, 1] = - JuMP.@constraint(get_jump_model(container), reg_var[name, 1] == 0) - constraint_lb[name, 1] = - JuMP.@constraint(get_jump_model(container), reg_var[name, 1] == 0) - - for t in time_steps[2:end] - if has_services - constraint_ub[name, t] = JuMP.@constraint( - get_jump_model(container), - ( - powerin_var[name, t - 1] + r_dn_ch[name, t - 1] - - r_up_ch[name, t - 1] - ) - (powerin_var[name, t] + r_dn_ch[name, t] - r_up_ch[name, t]) <= - reg_var[name, t] - ) - constraint_lb[name, t] = JuMP.@constraint( - get_jump_model(container), - ( - powerin_var[name, t - 1] + r_dn_ch[name, t - 1] - - r_up_ch[name, t - 1] - ) - (powerin_var[name, t] + r_dn_ch[name, t] - r_up_ch[name, t]) >= - -reg_var[name, t] - ) - else - constraint_ub[name, t] = JuMP.@constraint( - get_jump_model(container), - powerin_var[name, t - 1] - powerin_var[name, t] <= reg_var[name, t] - ) - constraint_lb[name, t] = JuMP.@constraint( - get_jump_model(container), - powerin_var[name, t - 1] - powerin_var[name, t] >= -reg_var[name, t] - ) - end - end - end - - return -end +# Trait stubs for the unified Charge/Discharge regularization body. Sign +# convention for net injection: charge nets to (p − r_up + r_dn); discharge +# nets to (p + r_up − r_dn). Mirrors the hybrid pair in +# src/hybrid_system_models/hybrid_systems.jl. +_storage_reg_slack_var(::Type{StorageRegularizationConstraintCharge}) = + StorageRegularizationVariableCharge +_storage_reg_slack_var(::Type{StorageRegularizationConstraintDischarge}) = + StorageRegularizationVariableDischarge +_storage_reg_power_var(::Type{StorageRegularizationConstraintCharge}) = + ActivePowerInVariable +_storage_reg_power_var(::Type{StorageRegularizationConstraintDischarge}) = + ActivePowerOutVariable +_storage_reg_reserve_exprs(::Type{StorageRegularizationConstraintCharge}) = + ( + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}, + ) +_storage_reg_reserve_exprs(::Type{StorageRegularizationConstraintDischarge}) = + ( + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, DischargeSide}, + ) +_storage_reg_reserve_signs(::Type{StorageRegularizationConstraintCharge}) = (-1, +1) +_storage_reg_reserve_signs(::Type{StorageRegularizationConstraintDischarge}) = (+1, -1) function add_constraints!( container::OptimizationContainer, - ::Type{StorageRegularizationConstraintDischarge}, + ::Type{T}, devices::IS.FlattenIteratorWrapper{V}, model::DeviceModel{V, StorageDispatchWithReserves}, network_model::NetworkModel{X}, -) where {V <: PSY.Storage, X <: AbstractPowerModel} +) where { + T <: Union{ + StorageRegularizationConstraintCharge, + StorageRegularizationConstraintDischarge, + }, + V <: PSY.Storage, + X <: AbstractPowerModel, +} names = [PSY.get_name(x) for x in devices] time_steps = get_time_steps(container) - reg_var = get_variable(container, StorageRegularizationVariableDischarge, V) - powerout_var = get_variable(container, ActivePowerOutVariable, V) + reg_var = get_variable(container, _storage_reg_slack_var(T), V) + p_var = get_variable(container, _storage_reg_power_var(T), V) has_services = has_service_model(model) - if has_services - r_up_ds = get_expression(container, ReserveDeploymentBalanceUpDischarge, V) - r_dn_ds = get_expression(container, ReserveDeploymentBalanceDownDischarge, V) - end + s_up, s_dn = _storage_reg_reserve_signs(T) + UpExpr, DnExpr = _storage_reg_reserve_exprs(T) + r_up = has_services ? get_expression(container, UpExpr, V) : nothing + r_dn = has_services ? get_expression(container, DnExpr, V) : nothing constraint_ub = - add_constraints_container!(container, StorageRegularizationConstraintDischarge, - V, - names, - time_steps; - meta = "ub", - ) - + add_constraints_container!(container, T, V, names, time_steps; meta = "ub") constraint_lb = - add_constraints_container!(container, StorageRegularizationConstraintDischarge, - V, - names, - time_steps; - meta = "lb", - ) + add_constraints_container!(container, T, V, names, time_steps; meta = "lb") + jm = get_jump_model(container) for d in devices name = PSY.get_name(d) - constraint_ub[name, 1] = - JuMP.@constraint(get_jump_model(container), reg_var[name, 1] == 0) - constraint_lb[name, 1] = - JuMP.@constraint(get_jump_model(container), reg_var[name, 1] == 0) + constraint_ub[name, 1] = JuMP.@constraint(jm, reg_var[name, 1] == 0) + constraint_lb[name, 1] = JuMP.@constraint(jm, reg_var[name, 1] == 0) for t in time_steps[2:end] - if has_services - constraint_ub[name, t] = JuMP.@constraint( - get_jump_model(container), - ( - powerout_var[name, t - 1] + r_up_ds[name, t - 1] - - r_dn_ds[name, t - 1] - ) - (powerout_var[name, t] + r_up_ds[name, t] - r_dn_ds[name, t]) <= - reg_var[name, t] - ) - constraint_lb[name, t] = JuMP.@constraint( - get_jump_model(container), - ( - powerout_var[name, t - 1] + r_up_ds[name, t - 1] - - r_dn_ds[name, t - 1] - ) - (powerout_var[name, t] + r_up_ds[name, t] - r_dn_ds[name, t]) >= - -reg_var[name, t] - ) + lhs = if has_services + ( + p_var[name, t - 1] + + s_up * r_up[name, t - 1] + + s_dn * r_dn[name, t - 1] + ) - (p_var[name, t] + s_up * r_up[name, t] + s_dn * r_dn[name, t]) else - constraint_ub[name, t] = JuMP.@constraint( - get_jump_model(container), - powerout_var[name, t - 1] - powerout_var[name, t] <= reg_var[name, t] - ) - constraint_lb[name, t] = JuMP.@constraint( - get_jump_model(container), - powerout_var[name, t - 1] - powerout_var[name, t] >= -reg_var[name, t] - ) + p_var[name, t - 1] - p_var[name, t] end + constraint_ub[name, t] = JuMP.@constraint(jm, lhs <= reg_var[name, t]) + constraint_lb[name, t] = JuMP.@constraint(jm, lhs >= -reg_var[name, t]) end end return diff --git a/src/hybrid_system_models/hybrid_systems.jl b/src/hybrid_system_models/hybrid_systems.jl new file mode 100644 index 0000000..c61fa56 --- /dev/null +++ b/src/hybrid_system_models/hybrid_systems.jl @@ -0,0 +1,2396 @@ +requires_initialization(::AbstractHybridFormulation) = false + +################################################################################# +# Default time-series and attributes +################################################################################# + +function get_default_time_series_names( + ::Type{PSY.HybridSystem}, + ::Type{<:Union{FixedOutput, AbstractHybridFormulation}}, +) + return Dict{Type{<:TimeSeriesParameter}, String}( + HybridRenewableActivePowerTimeSeriesParameter => "RenewableDispatch__max_active_power", + HybridElectricLoadTimeSeriesParameter => "PowerLoad__max_active_power", + ) +end + +function get_default_attributes( + ::Type{PSY.HybridSystem}, + ::Type{<:Union{FixedOutput, AbstractHybridFormulation}}, +) + return Dict{String, Any}( + "reservation" => true, + "storage_reservation" => true, + "energy_target" => false, + "regularization" => false, + ) +end + +# Small fixed cost rate on regularization slacks. Mirrors HSS REG_COST. +const HYBRID_REGULARIZATION_COST = 1e-3 + +################################################################################# +# PCC variables — ActivePowerInVariable / ActivePowerOutVariable +################################################################################# + +get_variable_binary( + ::Type{ActivePowerInVariable}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{ActivePowerInVariable}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_input_active_power_limits(d, PSY.SU).min +get_variable_upper_bound( + ::Type{ActivePowerInVariable}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_input_active_power_limits(d, PSY.SU).max +get_variable_multiplier( + ::Type{ActivePowerInVariable}, + ::Type{<:PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = -1.0 + +get_variable_binary( + ::Type{ActivePowerOutVariable}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{ActivePowerOutVariable}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_output_active_power_limits(d, PSY.SU).min +get_variable_upper_bound( + ::Type{ActivePowerOutVariable}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_output_active_power_limits(d, PSY.SU).max +get_variable_multiplier( + ::Type{ActivePowerOutVariable}, + ::Type{<:PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = 1.0 + +get_variable_binary( + ::Type{ReactivePowerVariable}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +function get_variable_lower_bound( + ::Type{ReactivePowerVariable}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) + limits = PSY.get_reactive_power_limits(d, PSY.SU) + return limits === nothing ? nothing : limits.min +end +function get_variable_upper_bound( + ::Type{ReactivePowerVariable}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) + limits = PSY.get_reactive_power_limits(d, PSY.SU) + return limits === nothing ? nothing : limits.max +end +get_variable_multiplier( + ::Type{ReactivePowerVariable}, + ::Type{<:PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = 1.0 + +get_variable_binary( + ::Type{ReservationVariable}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = true + +get_min_max_limits( + d::PSY.HybridSystem, + ::Type{InputActivePowerVariableLimitsConstraint}, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_input_active_power_limits(d, PSY.SU) +get_min_max_limits( + d::PSY.HybridSystem, + ::Type{OutputActivePowerVariableLimitsConstraint}, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_output_active_power_limits(d, PSY.SU) +get_min_max_limits( + d::PSY.HybridSystem, + ::Type{ReactivePowerVariableLimitsConstraint}, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_reactive_power_limits(d, PSY.SU) + +################################################################################# +# Subcomponent power variables +################################################################################# + +get_variable_binary( + ::Type{HybridThermalActivePower}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{HybridThermalActivePower}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = 0.0 +get_variable_upper_bound( + ::Type{HybridThermalActivePower}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_active_power_limits(PSY.get_thermal_unit(d), PSY.SU).max + +get_variable_binary( + ::Type{HybridRenewableActivePower}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{HybridRenewableActivePower}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = 0.0 +get_variable_upper_bound( + ::Type{HybridRenewableActivePower}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_max_active_power(PSY.get_renewable_unit(d), PSY.SU) + +get_variable_binary( + ::Type{HybridStorageSubcomponentPower{ChargeSide}}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{HybridStorageSubcomponentPower{ChargeSide}}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = 0.0 +get_variable_upper_bound( + ::Type{HybridStorageSubcomponentPower{ChargeSide}}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_input_active_power_limits(PSY.get_storage(d), PSY.SU).max + +get_variable_binary( + ::Type{HybridStorageSubcomponentPower{DischargeSide}}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{HybridStorageSubcomponentPower{DischargeSide}}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = 0.0 +get_variable_upper_bound( + ::Type{HybridStorageSubcomponentPower{DischargeSide}}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_output_active_power_limits(PSY.get_storage(d), PSY.SU).max + +get_variable_binary( + ::Type{HybridStorageReservation}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = true + +get_variable_binary( + ::Type{RegularizationVariable{ChargeSide}}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{RegularizationVariable{ChargeSide}}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = 0.0 + +get_variable_binary( + ::Type{RegularizationVariable{DischargeSide}}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{RegularizationVariable{DischargeSide}}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = 0.0 + +# Storage energy state on the hybrid (uses POM's standard EnergyVariable, keyed by HybridSystem) +get_variable_binary( + ::Type{EnergyVariable}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{EnergyVariable}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = + PSY.get_storage_level_limits(PSY.get_storage(d)).min * + PSY.get_storage_capacity(PSY.get_storage(d), PSY.SU) * + PSY.get_conversion_factor(PSY.get_storage(d)) +get_variable_upper_bound( + ::Type{EnergyVariable}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = + PSY.get_storage_level_limits(PSY.get_storage(d)).max * + PSY.get_storage_capacity(PSY.get_storage(d), PSY.SU) * + PSY.get_conversion_factor(PSY.get_storage(d)) +get_variable_warm_start_value( + ::Type{EnergyVariable}, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = + PSY.get_initial_storage_capacity_level(PSY.get_storage(d)) * + PSY.get_storage_capacity(PSY.get_storage(d), PSY.SU) * + PSY.get_conversion_factor(PSY.get_storage(d)) + +# End-of-period energy-target slacks (added when `energy_target = true`). Non-negative, +# defined only at the final time step. Mirrors POM storage's slack variables at +# energy_storage_models/storage_models.jl:376-402, keyed by HybridSystem. +function add_variables!( + container::OptimizationContainer, + ::Type{T}, + devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}}, + ::Type{<:AbstractHybridFormulation}, +) where { + T <: Union{HybridEnergyShortageVariable, HybridEnergySurplusVariable}, + U <: PSY.HybridSystem, +} + @assert !isempty(devices) + time_steps = get_time_steps(container) + last_time_range = time_steps[end]:time_steps[end] + variable = add_variable_container!( + container, + T, + U, + PSY.get_name.(devices), + last_time_range, + ) + for d in devices + PSY.get_storage(d) === nothing && continue + name = PSY.get_name(d) + variable[name, time_steps[end]] = JuMP.@variable( + get_jump_model(container), + base_name = "$(T)_{$(PSY.get_name(d))}", + lower_bound = 0.0 + ) + end + return +end + +# Thermal commitment OnVariable on a hybrid (binary) +get_variable_binary( + ::Type{OnVariable}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = true +get_variable_lower_bound( + ::Type{OnVariable}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = nothing +get_variable_upper_bound( + ::Type{OnVariable}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = nothing + +################################################################################# +# Reserve variables — bounds and binary flags +################################################################################# + +get_variable_binary( + ::Type{ + <:Union{ + AbstractHybridSubcomponentInjectorReserveVariableType, + HybridStorageSubcomponentReserveVariable, + }, + }, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{ + <:Union{ + AbstractHybridSubcomponentInjectorReserveVariableType, + HybridStorageSubcomponentReserveVariable, + }, + }, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = 0.0 + +# Per-subcomponent reserve upper bounds: limited by the subcomponent's headroom × the service's max output fraction +function get_variable_upper_bound( + ::Type{HybridThermalReserveVariable}, + r::PSY.Reserve, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) + return PSY.get_max_output_fraction(r) * + PSY.get_active_power_limits(PSY.get_thermal_unit(d), PSY.SU).max +end +function get_variable_upper_bound( + ::Type{HybridRenewableReserveVariable}, + r::PSY.Reserve, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) + return PSY.get_max_output_fraction(r) * + PSY.get_max_active_power(PSY.get_renewable_unit(d), PSY.SU) +end +function get_variable_upper_bound( + ::Type{HybridStorageSubcomponentReserveVariable{ChargeSide}}, + r::PSY.Reserve, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) + return PSY.get_max_output_fraction(r) * + PSY.get_input_active_power_limits(PSY.get_storage(d), PSY.SU).max +end +function get_variable_upper_bound( + ::Type{HybridStorageSubcomponentReserveVariable{DischargeSide}}, + r::PSY.Reserve, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) + return PSY.get_max_output_fraction(r) * + PSY.get_output_active_power_limits(PSY.get_storage(d), PSY.SU).max +end + +# Hybrid PCC reserve variables — limited by the hybrid's PCC limits × max_output_fraction +get_variable_binary( + ::Type{HybridPCCReserveVariable{DischargeSide}}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{HybridPCCReserveVariable{DischargeSide}}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = 0.0 +function get_variable_upper_bound( + ::Type{HybridPCCReserveVariable{DischargeSide}}, + r::PSY.Reserve, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) + return PSY.get_max_output_fraction(r) * + PSY.get_output_active_power_limits(d, PSY.SU).max +end + +get_variable_binary( + ::Type{HybridPCCReserveVariable{ChargeSide}}, + ::Type{PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = false +get_variable_lower_bound( + ::Type{HybridPCCReserveVariable{ChargeSide}}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = 0.0 +function get_variable_upper_bound( + ::Type{HybridPCCReserveVariable{ChargeSide}}, + r::PSY.Reserve, + d::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) + return PSY.get_max_output_fraction(r) * PSY.get_input_active_power_limits(d, PSY.SU).max +end + +# Multipliers used by reserve aggregations (Out side gets +1; In side handled via separate dispatch in add_to_expression) +get_variable_multiplier( + ::Type{ + <:Union{ + AbstractHybridSubcomponentInjectorReserveVariableType, + HybridStorageSubcomponentReserveVariable, + }, + }, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulationWithReserves}, + ::PSY.Reserve, +) = 1.0 +get_variable_multiplier( + ::Type{HybridPCCReserveVariable{DischargeSide}}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulationWithReserves}, + ::PSY.Reserve, +) = 1.0 +get_variable_multiplier( + ::Type{HybridPCCReserveVariable{ChargeSide}}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulationWithReserves}, + ::PSY.Reserve, +) = 1.0 + +# When the system-side ActivePowerReserveVariable is added by the service constructor for a HybridSystem, +# direct it into the TotalReserveOffering channel keyed by HybridSystem (mirrors POM storage line 59). +get_expression_type_for_reserve( + ::Type{ActivePowerReserveVariable}, + ::Type{<:PSY.HybridSystem}, + ::Type{<:PSY.Reserve}, +) = TotalReserveOffering + +function get_variable_upper_bound( + ::Type{ActivePowerReserveVariable}, + r::PSY.Reserve, + d::PSY.HybridSystem, + ::Type{<:AbstractReservesFormulation}, +) + return PSY.get_max_output_fraction(r) * ( + PSY.get_output_active_power_limits(d, PSY.SU).max + + PSY.get_input_active_power_limits(d, PSY.SU).max + ) +end + +# Disambiguate against the generic ReserveDemandCurve method in services_models/reserves.jl. +function get_variable_upper_bound( + ::Type{ActivePowerReserveVariable}, + r::PSY.ReserveDemandCurve, + d::PSY.HybridSystem, + ::Type{<:AbstractReservesFormulation}, +) + return PSY.get_output_active_power_limits(d, PSY.SU).max + + PSY.get_input_active_power_limits(d, PSY.SU).max +end + +################################################################################# +# Time-series parameter multipliers +################################################################################# + +get_multiplier_value( + ::HybridRenewableActivePowerTimeSeriesParameter, + d::PSY.HybridSystem, + ::AbstractHybridFormulation, +) = PSY.get_max_active_power(PSY.get_renewable_unit(d), PSY.SU) + +get_multiplier_value( + ::HybridElectricLoadTimeSeriesParameter, + d::PSY.HybridSystem, + ::AbstractHybridFormulation, +) = PSY.get_max_active_power(PSY.get_electric_load(d), PSY.SU) + +get_parameter_multiplier( + ::HybridRenewableActivePowerTimeSeriesParameter, + ::PSY.HybridSystem, + ::AbstractHybridFormulation, +) = 1.0 +get_parameter_multiplier( + ::HybridElectricLoadTimeSeriesParameter, + ::PSY.HybridSystem, + ::AbstractHybridFormulation, +) = 1.0 + +################################################################################# +# Initial conditions +################################################################################# + +get_initial_conditions_device_model( + ::IOM.AbstractOptimizationModel, + model::DeviceModel{T, <:AbstractHybridFormulation}, +) where {T <: PSY.HybridSystem} = model + +initial_condition_default( + ::InitialEnergyLevel, + d::PSY.HybridSystem, + ::AbstractHybridFormulation, +) = + PSY.get_initial_storage_capacity_level(PSY.get_storage(d)) * + PSY.get_storage_capacity(PSY.get_storage(d), PSY.SU) * + PSY.get_conversion_factor(PSY.get_storage(d)) + +initial_condition_variable( + ::InitialEnergyLevel, + d::PSY.HybridSystem, + ::AbstractHybridFormulation, +) = EnergyVariable() + +function initial_conditions!( + container::OptimizationContainer, + devices::IS.FlattenIteratorWrapper{T}, + formulation::AbstractHybridFormulation, +) where {T <: PSY.HybridSystem} + storage_devices = [d for d in devices if PSY.get_storage(d) !== nothing] + if !isempty(storage_devices) + add_initial_condition!( + container, + storage_devices, + formulation, + InitialEnergyLevel(), + ) + end + return +end + +################################################################################# +# Specialized add_variables! for the per-service reserve variables. +# +# These variables are indexed by (service_type, service_name) in addition to +# (component_name, time). Mirrors POM storage's pattern at +# energy_storage_models/storage_models.jl:409–445. +################################################################################# + +function add_variables!( + container::OptimizationContainer, + ::Type{T}, + devices::U, + ::Type{F}, +) where { + T <: AbstractHybridReserveVariableType, + U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, + F <: AbstractHybridFormulation, +} where {D <: PSY.HybridSystem} + @assert !isempty(devices) + time_steps = get_time_steps(container) + services = Set{PSY.Service}() + for d in devices + union!(services, PSY.get_services(d)) + end + isempty(services) && return + for service in services + # Restrict to devices that participate in this service + participating = [d for d in devices if service in PSY.get_services(d)] + isempty(participating) && continue + variable = add_variable_container!(container, T, + D, + PSY.get_name.(participating), + time_steps; + meta = "$(typeof(service))_$(PSY.get_name(service))", + ) + for d in participating, t in time_steps + name = PSY.get_name(d) + variable[name, t] = JuMP.@variable( + get_jump_model(container), + base_name = "$(T)_$(PSY.get_name(service))_{$(name), $(t)}", + lower_bound = 0.0, + upper_bound = get_variable_upper_bound(T, service, d, F), + ) + end + end + return +end + +################################################################################# +# Objective-function multipliers (positive — we minimize cost) +################################################################################# + +objective_function_multiplier(::Type{<:VariableType}, ::Type{<:AbstractHybridFormulation}) = + OBJECTIVE_FUNCTION_POSITIVE + +################################################################################# +# PCC active-power balance: ActivePowerInVariable / ActivePowerOutVariable into +# the network's ActivePowerBalance expression. +# +# These delegate to POM's existing common_models/add_to_expression.jl methods, +# which dispatch on (ExpressionType, VariableType, AbstractDeviceFormulation). +# Because AbstractHybridFormulation <: IOM.AbstractDeviceFormulation, the +# generic methods work for HybridSystem out of the box. The methods here are +# left documented but not redefined to avoid ambiguity. +################################################################################# + +################################################################################# +# Reserve term accumulation — unified across hybrid PCC and storage subcomponent. +# +# The parametric ReserveAggregationExpression{Direction, Scale, Side} family lets +# one helper handle both the hybrid-boundary aggregation (HybridPCCReserveVariable{DischargeSide}/In +# into HybridPCCReserveExpression{...}) and the storage-subcomponent aggregation +# (HybridStorageSubcomponentReserveVariable{ChargeSide}/Discharging... into StorageReserveBalanceExpression{...}). +# Mismatched-direction services are filtered out by dispatch on the Direction parameter +# of the expression type vs the Reserve direction (ReserveUp / ReserveDown). +# The Scale parameter (UnscaledReserve / DeployedReserve) drives the multiplier scale. +################################################################################# + +# Multiplier scale: UnscaledReserve → 1.0; DeployedReserve → deployed_fraction(service). +_reserve_scale( + ::Type{<:ReserveAggregationExpression{<:PSY.ReserveDirection, UnscaledReserve}}, + ::PSY.Service, +) = 1.0 +_reserve_scale( + ::Type{<:ReserveAggregationExpression{<:PSY.ReserveDirection, DeployedReserve}}, + s::PSY.Service, +) = PSY.get_deployed_fraction(s) + +# Up-direction expressions: ReserveDown services are a no-op (skipped via dispatch). +_add_reserve_term!( + ::Type{<:ReserveAggregationExpression{PSY.ReserveUp}}, + ::OptimizationContainer, + _expression, + ::Type{<:AbstractHybridReserveVariableType}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulationWithReserves}, + ::Int, + ::PSY.Reserve{PSY.ReserveDown}, +) = nothing + +# Down-direction expressions: ReserveUp services are a no-op (skipped via dispatch). +_add_reserve_term!( + ::Type{<:ReserveAggregationExpression{PSY.ReserveDown}}, + ::OptimizationContainer, + _expression, + ::Type{<:AbstractHybridReserveVariableType}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulationWithReserves}, + ::Int, + ::PSY.Reserve{PSY.ReserveUp}, +) = nothing + +# Fallback: actually accumulate the (correct-direction) reserve term. +function _add_reserve_term!( + ::Type{T}, + container::OptimizationContainer, + expression, + ::Type{U}, + d::V, + ::Type{W}, + t::Int, + service::PSY.Service, +) where { + T <: ReserveAggregationExpression, + U <: AbstractHybridReserveVariableType, + V <: PSY.HybridSystem, + W <: AbstractHybridFormulationWithReserves, +} + name = PSY.get_name(d) + variable = + get_variable(container, U, V, "$(typeof(service))_$(PSY.get_name(service))") + mult = get_variable_multiplier(U, d, W, service) * _reserve_scale(T, service) + add_proportional_to_jump_expression!(expression[name, t], variable[name, t], mult) + return +end + +# Single add_to_expression! method covering both PCC boundary and storage subcomponent. +function add_to_expression!( + container::OptimizationContainer, + ::Type{T}, + ::Type{U}, + devices::Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + ::DeviceModel{V, W}, +) where { + T <: ReserveAggregationExpression, + U <: AbstractHybridReserveVariableType, + V <: PSY.HybridSystem, + W <: AbstractHybridFormulationWithReserves, +} + expression = get_expression(container, T, V) + for d in devices, service in PSY.get_services(d), t in get_time_steps(container) + _add_reserve_term!(T, container, expression, U, d, W, t, service) + end + return +end + +# Variant signature retained for callers that also pass a NetworkModel (hybrid PCC path). +function add_to_expression!( + container::OptimizationContainer, + ::Type{T}, + ::Type{U}, + devices::Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + T <: ReserveAggregationExpression, + U <: AbstractHybridReserveVariableType, + V <: PSY.HybridSystem, + W <: AbstractHybridFormulationWithReserves, + X <: AbstractPowerModel, +} + add_to_expression!(container, T, U, devices, model) + return +end + +# Hybrid storage subcomponent feeds TotalReserveOffering keyed by HybridSystem. +function add_to_expression!( + container::OptimizationContainer, + ::Type{TotalReserveOffering}, + ::Type{U}, + devices::Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + model::DeviceModel{V, W}, +) where { + U <: HybridStorageSubcomponentReserveVariable, + V <: PSY.HybridSystem, + W <: AbstractHybridFormulationWithReserves, +} + time_steps = get_time_steps(container) + for d in devices + name = PSY.get_name(d) + for service in PSY.get_services(d) + expression = get_expression(container, TotalReserveOffering, V, + "$(typeof(service))_$(PSY.get_name(service))") + variable = + get_variable(container, U, V, "$(typeof(service))_$(PSY.get_name(service))") + mult = get_variable_multiplier(U, d, W, service) + for t in time_steps + add_proportional_to_jump_expression!( + expression[name, t], + variable[name, t], + mult, + ) + end + end + end + return +end + +# Service-side: ActivePowerReserveVariable subtracted from per-hybrid TotalReserveOffering +# (mirrors storage_models.jl:781–808 pattern). +function add_to_expression!( + container::OptimizationContainer, + ::Type{T}, + ::Type{U}, + devices::Vector{UV}, + service_model::ServiceModel{V, W}, +) where { + T <: TotalReserveOffering, + U <: ActivePowerReserveVariable, + UV <: PSY.HybridSystem, + V <: PSY.Reserve, + W <: AbstractReservesFormulation, +} + for d in devices + name = PSY.get_name(d) + s_name = get_service_name(service_model) + expression = get_expression(container, T, UV, "$(V)_$(s_name)") + variable = get_variable(container, U, V, s_name) + for t in get_time_steps(container) + add_proportional_to_jump_expression!( + expression[name, t], + variable[name, t], + -1.0, + ) + end + end + return +end +################################################################################# +# Subcomponent (thermal / renewable) reserve accumulators — unified family. +# +# A single helper covers both thermal and renewable subcomponent reserve +# accumulation into a JuMP.AffExpr. The reserve variable type +# (U <: AbstractHybridSubcomponentInjectorReserveVariableType) selects which +# subcomponent we're aggregating, and the direction marker (Up / Down) +# filters out mismatched-direction services via dispatch. +# +# Callers in HybridThermalReserveLimitConstraint and HybridRenewableReserveLimit- +# Constraint invoke _subcomponent_reserve_expr(Up | ReserveDown, container, +# HybridThermalReserveVariable | HybridRenewableReserveVariable, d, t, services). +################################################################################# + +# Up direction: ReserveDown service is a no-op. +_subcomponent_reserve_term!( + ::Type{PSY.ReserveUp}, + ::JuMP.AffExpr, + ::OptimizationContainer, + ::Type{<:AbstractHybridSubcomponentInjectorReserveVariableType}, + ::PSY.HybridSystem, + ::Int, + ::PSY.Reserve{PSY.ReserveDown}, +) = nothing + +# Down direction: ReserveUp service is a no-op. +_subcomponent_reserve_term!( + ::Type{PSY.ReserveDown}, + ::JuMP.AffExpr, + ::OptimizationContainer, + ::Type{<:AbstractHybridSubcomponentInjectorReserveVariableType}, + ::PSY.HybridSystem, + ::Int, + ::PSY.Reserve{PSY.ReserveUp}, +) = nothing + +# Fallback: accumulate the term for the correct-direction service. +function _subcomponent_reserve_term!( + ::Type{<:PSY.ReserveDirection}, + expr::JuMP.AffExpr, + container::OptimizationContainer, + ::Type{U}, + d::V, + t::Int, + service::PSY.Service, +) where { + U <: AbstractHybridSubcomponentInjectorReserveVariableType, + V <: PSY.HybridSystem, +} + s_name = PSY.get_name(service) + s_type = typeof(service) + key = VariableKey(U, V, "$(s_type)_$s_name") + haskey(IOM.get_variables(container), key) || return + var = get_variable(container, U, V, "$(s_type)_$s_name") + add_proportional_to_jump_expression!(expr, var[PSY.get_name(d), t], 1.0) + return +end + +function _subcomponent_reserve_expr( + ::Type{Dir}, + container::OptimizationContainer, + ::Type{U}, + d::V, + t::Int, + services, +) where { + Dir <: PSY.ReserveDirection, + U <: AbstractHybridSubcomponentInjectorReserveVariableType, + V <: PSY.HybridSystem, +} + expr = JuMP.AffExpr(0.0) + for service in services + _subcomponent_reserve_term!(Dir, expr, container, U, d, t, service) + end + return expr +end + +################################################################################# +# Thermal subcomponent constraints for HybridSystem. +# +# Mirrors HSS add_constraints.jl _add_thermallimit_withreserves! (lines 1477–1506) +# for the with-reserves case, and _add_thermal_on_variable_constraints! for the +# no-reserves case. Walks PSY.get_thermal_unit(d) for the thermal unit's limits. +################################################################################# + +""" +Range constraint on the thermal subcomponent's active power, accounting for +up/down reserve allocations. Mirrors HSS `ThermalReserveLimit` (HSS +add_constraints.jl:1495–1506). +""" +function add_constraints!( + container::OptimizationContainer, + ::Type{HybridThermalReserveLimitConstraint}, + devices::U, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + W <: AbstractHybridFormulationWithReserves, + X <: AbstractPowerModel, +} where {V <: PSY.HybridSystem} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + p_th = get_variable(container, HybridThermalActivePower, V) + on_var = get_variable(container, OnVariable, V) + + con_ub = add_constraints_container!( + container, + HybridThermalReserveLimitConstraint, + V, + names, + time_steps; + meta = "ub", + ) + con_lb = add_constraints_container!( + container, + HybridThermalReserveLimitConstraint, + V, + names, + time_steps; + meta = "lb", + ) + + for d in devices, t in time_steps + name = PSY.get_name(d) + thermal_unit = PSY.get_thermal_unit(d) + thermal_unit === nothing && continue + limits = PSY.get_active_power_limits(thermal_unit, PSY.SU) + services = PSY.get_services(d) + r_up = _subcomponent_reserve_expr( + PSY.ReserveUp, + container, + HybridThermalReserveVariable, + d, + t, + services, + ) + r_dn = _subcomponent_reserve_expr( + PSY.ReserveDown, + container, + HybridThermalReserveVariable, + d, + t, + services, + ) + con_ub[name, t] = JuMP.@constraint( + get_jump_model(container), + p_th[name, t] + r_up <= limits.max * on_var[name, t] + ) + con_lb[name, t] = JuMP.@constraint( + get_jump_model(container), + p_th[name, t] - r_dn >= limits.min * on_var[name, t] + ) + end + return +end + +# Per-bound traits: which side of `lim` to use, and which JuMP relation to emit. +_thermal_on_limit(::Type{HybridThermalOnVariableConstraint{UpperBound}}, lim) = lim.max +_thermal_on_limit(::Type{HybridThermalOnVariableConstraint{LowerBound}}, lim) = lim.min +_thermal_on_relation(::Type{HybridThermalOnVariableConstraint{UpperBound}}, jm, lhs, rhs) = + JuMP.@constraint(jm, lhs <= rhs) +_thermal_on_relation(::Type{HybridThermalOnVariableConstraint{LowerBound}}, jm, lhs, rhs) = + JuMP.@constraint(jm, lhs >= rhs) + +""" +Bound link between thermal subcomponent power and its commitment status +(no-reserves case). Parametric on `BoundDirection` (from IOM): `{UpperBound}` enforces +`p_th ≤ max · on_var`, `{LowerBound}` enforces `p_th ≥ min · on_var`. +""" +function add_constraints!( + container::OptimizationContainer, + ::Type{T}, + devices::U, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + T <: HybridThermalOnVariableConstraint, + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + W <: AbstractHybridFormulation, + X <: AbstractPowerModel, +} where {V <: PSY.HybridSystem} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + p_th = get_variable(container, HybridThermalActivePower, V) + on_var = get_variable(container, OnVariable, V) + constraint = add_constraints_container!(container, T, V, names, time_steps) + jm = get_jump_model(container) + for d in devices, t in time_steps + name = PSY.get_name(d) + thermal_unit = PSY.get_thermal_unit(d) + thermal_unit === nothing && continue + bound = _thermal_on_limit(T, PSY.get_active_power_limits(thermal_unit, PSY.SU)) + constraint[name, t] = + _thermal_on_relation(T, jm, p_th[name, t], bound * on_var[name, t]) + end + return +end +################################################################################# +# Renewable subcomponent constraints for HybridSystem. +# +# - HybridRenewableActivePowerLimitConstraint: cap renewable subcomponent power +# at the time-series-derived available output (no-reserves and reserves cases +# share this; the reserve-aware variant carves out reserves in the with-reserves +# constraint). +# - HybridRenewableReserveLimitConstraint: range constraint on renewable power +# accounting for up/down reserves. The accumulator is shared with the thermal +# case via `_subcomponent_reserve_expr`, dispatching on the variable type. +################################################################################# + +""" +Cap renewable subcomponent power at the time-series-derived available output +(0 ≤ p_renewable[t] ≤ multiplier · ts[t]). +""" +function add_constraints!( + container::OptimizationContainer, + ::Type{HybridRenewableActivePowerLimitConstraint}, + devices::U, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + W <: AbstractHybridFormulation, + X <: AbstractPowerModel, +} where {V <: PSY.HybridSystem} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + p_re = get_variable(container, HybridRenewableActivePower, V) + + re_param_key = ParameterKey(HybridRenewableActivePowerTimeSeriesParameter, V) + re_param_container = if haskey(IOM.get_parameters(container), re_param_key) + get_parameter(container, HybridRenewableActivePowerTimeSeriesParameter, V) + else + nothing + end + re_multiplier = + re_param_container === nothing ? nothing : + get_multiplier_array(re_param_container) + + constraint = add_constraints_container!( + container, + HybridRenewableActivePowerLimitConstraint, + V, + names, + time_steps, + ) + + for d in devices, t in time_steps + name = PSY.get_name(d) + renewable_unit = PSY.get_renewable_unit(d) + renewable_unit === nothing && continue + if re_param_container !== nothing + re_ref = get_parameter_column_refs(re_param_container, name)[t] + constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + p_re[name, t] <= re_multiplier[name, t] * re_ref + ) + else + max_p = PSY.get_max_active_power(renewable_unit, PSY.SU) + constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + p_re[name, t] <= max_p + ) + end + end + return +end + +""" +Range constraint on renewable subcomponent power accounting for reserves. +Mirrors HSS `RenewableReserveLimit`. +""" +function add_constraints!( + container::OptimizationContainer, + ::Type{HybridRenewableReserveLimitConstraint}, + devices::U, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + W <: AbstractHybridFormulationWithReserves, + X <: AbstractPowerModel, +} where {V <: PSY.HybridSystem} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + p_re = get_variable(container, HybridRenewableActivePower, V) + + re_param_key = ParameterKey(HybridRenewableActivePowerTimeSeriesParameter, V) + re_param_container = if haskey(IOM.get_parameters(container), re_param_key) + get_parameter(container, HybridRenewableActivePowerTimeSeriesParameter, V) + else + nothing + end + re_multiplier = + re_param_container === nothing ? nothing : + get_multiplier_array(re_param_container) + + con_ub = add_constraints_container!( + container, + HybridRenewableReserveLimitConstraint, + V, + names, + time_steps; + meta = "ub", + ) + con_lb = add_constraints_container!( + container, + HybridRenewableReserveLimitConstraint, + V, + names, + time_steps; + meta = "lb", + ) + + for d in devices, t in time_steps + name = PSY.get_name(d) + renewable_unit = PSY.get_renewable_unit(d) + renewable_unit === nothing && continue + services = PSY.get_services(d) + r_up = _subcomponent_reserve_expr( + PSY.ReserveUp, + container, + HybridRenewableReserveVariable, + d, + t, + services, + ) + r_dn = _subcomponent_reserve_expr( + PSY.ReserveDown, + container, + HybridRenewableReserveVariable, + d, + t, + services, + ) + if re_param_container !== nothing + re_ref = get_parameter_column_refs(re_param_container, name)[t] + con_ub[name, t] = JuMP.@constraint( + get_jump_model(container), + p_re[name, t] + r_up <= re_multiplier[name, t] * re_ref + ) + else + max_p = PSY.get_max_active_power(renewable_unit, PSY.SU) + con_ub[name, t] = JuMP.@constraint( + get_jump_model(container), + p_re[name, t] + r_up <= max_p + ) + end + con_lb[name, t] = JuMP.@constraint( + get_jump_model(container), + p_re[name, t] - r_dn >= 0.0 + ) + end + return +end +################################################################################# +# Storage subcomponent constraints for HybridSystem (Option D core). +# +# Most methods re-emit POM's storage reserve constraint TYPES with new dispatches +# on V <: PSY.HybridSystem, substituting PSY.get_storage(hybrid) for hybrid at +# every PSY accessor. The constraint TYPES are reused (same names, same purpose, +# same shape); only the dispatch context changes. Hybrid-specific constraint +# types are introduced for the inner-storage status (charge/discharge mode) and +# the charge/discharge reserve power limits, since their math is subtly different +# from POM's storage versions. +################################################################################# + +################################################################################# +# HybridStorageBalanceConstraint — energy balance with optional reserve deployment +################################################################################# + +# With-reserves formulation: pick the body based on whether services are wired up. +function add_constraints!( + container::OptimizationContainer, + ::Type{HybridStorageBalanceConstraint}, + devices::U, + model::DeviceModel{V, W}, + network_model::NetworkModel{X}, +) where { + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + W <: AbstractHybridFormulationWithReserves, + X <: AbstractPowerModel, +} where {V <: PSY.HybridSystem} + if has_service_model(model) + _hybrid_storage_balance_with_reserves!(container, devices, model, network_model) + else + _hybrid_storage_balance_no_reserves!(container, devices, model, network_model) + end + return +end + +# Plain hybrid formulation (no reserves): always the no-reserves body. +function add_constraints!( + container::OptimizationContainer, + ::Type{HybridStorageBalanceConstraint}, + devices::U, + model::DeviceModel{V, W}, + network_model::NetworkModel{X}, +) where { + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + W <: AbstractHybridFormulation, + X <: AbstractPowerModel, +} where {V <: PSY.HybridSystem} + _hybrid_storage_balance_no_reserves!(container, devices, model, network_model) + return +end + +function _hybrid_storage_balance_no_reserves!( + container::OptimizationContainer, + devices, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where {V <: PSY.HybridSystem, W <: AbstractHybridFormulation, X <: AbstractPowerModel} + time_steps = get_time_steps(container) + resolution = get_resolution(container) + fraction_of_hour = Dates.value(Dates.Minute(resolution)) / MINUTES_IN_HOUR + names = [PSY.get_name(d) for d in devices] + initial_conditions = get_initial_condition(container, InitialEnergyLevel(), V) + energy_var = get_variable(container, EnergyVariable, V) + p_ch = get_variable(container, HybridStorageSubcomponentPower{ChargeSide}, V) + p_ds = get_variable(container, HybridStorageSubcomponentPower{DischargeSide}, V) + constraint = add_constraints_container!( + container, + HybridStorageBalanceConstraint, + V, + names, + time_steps, + ) + + for ic in initial_conditions + d = IOM.get_component(ic) + storage = PSY.get_storage(d) + storage === nothing && continue + eff = PSY.get_efficiency(storage) + name = PSY.get_name(d) + constraint[name, 1] = JuMP.@constraint( + get_jump_model(container), + energy_var[name, 1] == + get_value(ic) + + (p_ch[name, 1] * eff.in - p_ds[name, 1] / eff.out) * fraction_of_hour + ) + for t in time_steps[2:end] + constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + energy_var[name, t] == + energy_var[name, t - 1] + + (p_ch[name, t] * eff.in - p_ds[name, t] / eff.out) * fraction_of_hour + ) + end + end + return +end + +function _hybrid_storage_balance_with_reserves!( + container::OptimizationContainer, + devices, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + V <: PSY.HybridSystem, + W <: AbstractHybridFormulationWithReserves, + X <: AbstractPowerModel, +} + time_steps = get_time_steps(container) + resolution = get_resolution(container) + fraction_of_hour = Dates.value(Dates.Minute(resolution)) / MINUTES_IN_HOUR + names = [PSY.get_name(d) for d in devices] + initial_conditions = get_initial_condition(container, InitialEnergyLevel(), V) + energy_var = get_variable(container, EnergyVariable, V) + p_ch = get_variable(container, HybridStorageSubcomponentPower{ChargeSide}, V) + p_ds = get_variable(container, HybridStorageSubcomponentPower{DischargeSide}, V) + r_up_ds = get_expression( + container, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}, + V, + ) + r_up_ch = get_expression( + container, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}, + V, + ) + r_dn_ds = get_expression( + container, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, DischargeSide}, + V, + ) + r_dn_ch = get_expression( + container, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}, + V, + ) + constraint = add_constraints_container!( + container, + HybridStorageBalanceConstraint, + V, + names, + time_steps, + ) + + for ic in initial_conditions + d = IOM.get_component(ic) + storage = PSY.get_storage(d) + storage === nothing && continue + eff = PSY.get_efficiency(storage) + name = PSY.get_name(d) + constraint[name, 1] = JuMP.@constraint( + get_jump_model(container), + energy_var[name, 1] == + get_value(ic) + + ( + ((p_ch[name, 1] + r_dn_ch[name, 1] - r_up_ch[name, 1]) * eff.in) - + ((p_ds[name, 1] + r_up_ds[name, 1] - r_dn_ds[name, 1]) / eff.out) + ) * fraction_of_hour + ) + for t in time_steps[2:end] + constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + energy_var[name, t] == + energy_var[name, t - 1] + + ( + ((p_ch[name, t] + r_dn_ch[name, t] - r_up_ch[name, t]) * eff.in) - + ((p_ds[name, t] + r_up_ds[name, t] - r_dn_ds[name, t]) / eff.out) + ) * fraction_of_hour + ) + end + end + return +end + +################################################################################# +# HybridStorageStatusOnConstraint{ChargeSide} / HybridStorageStatusOnConstraint{DischargeSide} +# (no-reserves case — mutually exclusive charge/discharge via the inner storage +# reservation variable) +################################################################################# + +# Side-keyed traits shared by HybridStorageStatusOnConstraint{Sd} and +# HybridStorageReservePowerLimitConstraint{Sd} below. +# - ChargeSide : input limits, reservation factor (1 - ss). +# - DischargeSide: output limits, reservation factor ss. +const _StorageSideConstraint{Sd} = Union{ + HybridStorageStatusOnConstraint{Sd}, + HybridStorageReservePowerLimitConstraint{Sd}, +} + +_storage_side_power_var(::Type{<:_StorageSideConstraint{Sd}}) where {Sd <: ReserveSide} = + HybridStorageSubcomponentPower{Sd} +_storage_side_max(::Type{<:_StorageSideConstraint{ChargeSide}}, s) = + PSY.get_input_active_power_limits(s, PSY.SU).max +_storage_side_max(::Type{<:_StorageSideConstraint{DischargeSide}}, s) = + PSY.get_output_active_power_limits(s, PSY.SU).max +# Reservation-binary factor applied to the side limit. Charge side flips ss → (1-ss). +_storage_side_ss_factor(::Type{<:_StorageSideConstraint{ChargeSide}}, ss_val) = 1 - ss_val +_storage_side_ss_factor(::Type{<:_StorageSideConstraint{DischargeSide}}, ss_val) = ss_val + +function add_constraints!( + container::OptimizationContainer, + ::Type{T}, + devices::U, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + T <: HybridStorageStatusOnConstraint, + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + W <: AbstractHybridFormulation, + X <: AbstractPowerModel, +} where {V <: PSY.HybridSystem} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + p_var = get_variable(container, _storage_side_power_var(T), V) + ss = get_variable(container, HybridStorageReservation, V) + constraint = add_constraints_container!(container, T, V, names, time_steps) + for d in devices, t in time_steps + storage = PSY.get_storage(d) + storage === nothing && continue + name = PSY.get_name(d) + max_p = _storage_side_max(T, storage) + ss_factor = _storage_side_ss_factor(T, ss[name, t]) + constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + p_var[name, t] <= max_p * ss_factor + ) + end + return +end + +################################################################################# +# HybridStorageReservePowerLimitConstraint{ChargeSide} +# HybridStorageReservePowerLimitConstraint{DischargeSide} +# (with-reserves case — charge/discharge headroom under reservation + +# reserve-aware bounds, mirroring HSS's ChargingReservePowerLimit/ +# DischargingReservePowerLimit) +################################################################################# + +# Reserve-assignment expressions enter the bounds of the with-reserves storage +# power limits asymmetrically: charge UB picks up the down reserve (loading +# margin), charge LB subtracts the up reserve (headroom); discharge UB picks up +# the up reserve, discharge LB subtracts the down reserve. +_storage_side_ub_reserve_expr( + ::Type{HybridStorageReservePowerLimitConstraint{ChargeSide}}, +) = + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide} +_storage_side_ub_reserve_expr( + ::Type{HybridStorageReservePowerLimitConstraint{DischargeSide}}, +) = + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide} +_storage_side_lb_reserve_expr( + ::Type{HybridStorageReservePowerLimitConstraint{ChargeSide}}, +) = + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide} +_storage_side_lb_reserve_expr( + ::Type{HybridStorageReservePowerLimitConstraint{DischargeSide}}, +) = + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, DischargeSide} + +function add_constraints!( + container::OptimizationContainer, + ::Type{T}, + devices::U, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + T <: HybridStorageReservePowerLimitConstraint, + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + W <: AbstractHybridFormulationWithReserves, + X <: AbstractPowerModel, +} where {V <: PSY.HybridSystem} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + p_var = get_variable(container, _storage_side_power_var(T), V) + has_ss = haskey(IOM.get_variables(container), VariableKey(HybridStorageReservation, V)) + ss = has_ss ? get_variable(container, HybridStorageReservation, V) : nothing + r_ub = get_expression(container, _storage_side_ub_reserve_expr(T), V) + r_lb = get_expression(container, _storage_side_lb_reserve_expr(T), V) + con_ub = add_constraints_container!( + container, T, V, names, time_steps; meta = "ub") + con_lb = add_constraints_container!( + container, T, V, names, time_steps; meta = "lb") + for d in devices, t in time_steps + storage = PSY.get_storage(d) + storage === nothing && continue + name = PSY.get_name(d) + max_p = _storage_side_max(T, storage) + ub_rhs = has_ss ? max_p * _storage_side_ss_factor(T, ss[name, t]) : max_p + con_ub[name, t] = JuMP.@constraint( + get_jump_model(container), + p_var[name, t] + r_ub[name, t] <= ub_rhs + ) + con_lb[name, t] = JuMP.@constraint( + get_jump_model(container), + p_var[name, t] - r_lb[name, t] >= 0.0 + ) + end + return +end + +################################################################################# +# Charge/Discharge regularization constraints — penalize step changes in the +# charge/discharge profile via a non-negative slack. Mirrors HSS +# add_constraints.jl:1255–1424. When reserves are present, the served reserve +# expressions enter the step-change quantity so the regularization smooths the +# *net* injection profile, not the bare charge/discharge variable. +################################################################################# + +# Trait stubs for the unified Charge/Discharge regularization body. Sign +# convention for net injection: charge nets to (p − r_up + r_dn); discharge +# nets to (p + r_up − r_dn). +_reg_slack_var(::Type{RegularizationConstraint{Sd}}) where {Sd <: ReserveSide} = + RegularizationVariable{Sd} +_reg_power_var(::Type{RegularizationConstraint{Sd}}) where {Sd <: ReserveSide} = + HybridStorageSubcomponentPower{Sd} +_reg_reserve_exprs(::Type{RegularizationConstraint{Sd}}) where {Sd <: ReserveSide} = ( + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, Sd}, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, Sd}, +) +_reg_reserve_signs(::Type{RegularizationConstraint{ChargeSide}}) = (-1, +1) +_reg_reserve_signs(::Type{RegularizationConstraint{DischargeSide}}) = (+1, -1) + +function _hybrid_served_reserve_pair(container, ::Type{T}, V, name, t) where {T} + UpExpr, DnExpr = _reg_reserve_exprs(T) + if has_container_key(container, UpExpr, V) && + has_container_key(container, DnExpr, V) + up = get_expression(container, UpExpr, V)[name, t] + dn = get_expression(container, DnExpr, V)[name, t] + return up, dn + end + return 0.0, 0.0 +end + +# Per-formulation: does the regularization body include served-reserve terms? +# With-reserves formulation: include them only if the device model has a service model. +# Plain hybrid formulation: never include them (no axes to integrate against). +_regularization_has_services(::Type{<:AbstractHybridFormulationWithReserves}, model) = + has_service_model(model) +_regularization_has_services(::Type{<:AbstractHybridFormulation}, _model) = false + +function add_constraints!( + container::OptimizationContainer, + ::Type{T}, + devices::U, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + T <: RegularizationConstraint, + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + W <: AbstractHybridFormulation, + X <: AbstractPowerModel, +} where {V <: PSY.HybridSystem} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + reg_var = get_variable(container, _reg_slack_var(T), V) + p_var = get_variable(container, _reg_power_var(T), V) + has_services = _regularization_has_services(W, model) + s_up, s_dn = _reg_reserve_signs(T) + con_ub = add_constraints_container!( + container, T, V, names, time_steps; meta = "ub") + con_lb = add_constraints_container!( + container, T, V, names, time_steps; meta = "lb") + jm = get_jump_model(container) + t1 = first(time_steps) + for d in devices + PSY.get_storage(d) === nothing && continue + name = PSY.get_name(d) + # First time step: pin slack to zero (no previous step to compare against). + con_ub[name, t1] = JuMP.@constraint(jm, reg_var[name, t1] == 0) + con_lb[name, t1] = JuMP.@constraint(jm, reg_var[name, t1] == 0) + for t in time_steps[2:end] + if has_services + up_prev, dn_prev = + _hybrid_served_reserve_pair(container, T, V, name, t - 1) + up_t, dn_t = _hybrid_served_reserve_pair(container, T, V, name, t) + lhs = + (p_var[name, t - 1] + s_up * up_prev + s_dn * dn_prev) - + (p_var[name, t] + s_up * up_t + s_dn * dn_t) + else + lhs = p_var[name, t - 1] - p_var[name, t] + end + con_ub[name, t] = JuMP.@constraint(jm, lhs <= reg_var[name, t]) + con_lb[name, t] = JuMP.@constraint(jm, lhs >= -reg_var[name, t]) + end + end + return +end + +################################################################################# +# Reuse POM's ReserveCoverageConstraint{,EndOfPeriod} types with V <: HybridSystem +# dispatches. Bodies mirror storage_models.jl:1038–1108, substituting +# PSY.get_storage(d) for d at every PSY accessor. +################################################################################# + +const _ReserveCoverageT = + Union{ReserveCoverageConstraint, ReserveCoverageConstraintEndOfPeriod} + +# Container setup: dispatch on service type. Up → "_discharge" suffix (storage discharges +# to deliver up-reserve); Down → "_charge" suffix. `lazy_container_addition!` is idempotent, +# so calling these methods more than once per (T, V, service) is safe. +function _init_coverage_container!( + container::OptimizationContainer, + ::Type{T}, + ::Type{V}, + names::Vector{String}, + time_steps::UnitRange{Int}, + service::PSY.Reserve{PSY.ReserveUp}, +) where {T <: _ReserveCoverageT, V <: PSY.HybridSystem} + return lazy_container_addition!( + container, T, V, names, time_steps; + meta = "$(typeof(service))_$(PSY.get_name(service))_discharge", + ) +end + +function _init_coverage_container!( + container::OptimizationContainer, + ::Type{T}, + ::Type{V}, + names::Vector{String}, + time_steps::UnitRange{Int}, + service::PSY.Reserve{PSY.ReserveDown}, +) where {T <: _ReserveCoverageT, V <: PSY.HybridSystem} + return lazy_container_addition!( + container, T, V, names, time_steps; + meta = "$(typeof(service))_$(PSY.get_name(service))_charge", + ) +end + +_init_coverage_container!( + ::OptimizationContainer, + ::Type{<:_ReserveCoverageT}, + ::Type{<:PSY.HybridSystem}, + ::Vector{String}, + ::UnitRange{Int}, + ::PSY.Service, +) = nothing # subsumes the `(service isa PSY.Reserve) || continue` guard + +# Constraint emission: dispatch on service type. Up uses HybridStorageSubcomponentReserveVariable{DischargeSide} +# bounded by SoC; Down uses HybridStorageSubcomponentReserveVariable{ChargeSide} bounded by (soc_max − SoC). +# Sustained-time accessors exist only on PSY.Reserve, so the param computation lives +# inside the per-direction helpers — the PSY.Service fallback never touches them. +function _emit_coverage_constraint!( + container::OptimizationContainer, + ::Type{T}, + ::Type{V}, + ic::InitialCondition, + energy_var, + ci_name::String, + eff_in::Float64, + inv_eff_out::Float64, + fraction_of_hour::Float64, + resolution::Dates.Period, + storage::PSY.Storage, + time_steps::UnitRange{Int}, + service::PSY.Reserve{PSY.ReserveUp}, +) where {T <: _ReserveCoverageT, V <: PSY.HybridSystem} + s_type = typeof(service) + s_name = PSY.get_name(service) + num_periods = PSY.get_sustained_time(service) / Dates.value(Dates.Second(resolution)) + sustained_param_discharge = inv_eff_out * fraction_of_hour * num_periods + reserve_var = + get_variable( + container, + HybridStorageSubcomponentReserveVariable{DischargeSide}, + V, + "$(s_type)_$s_name", + ) + con = get_constraint(container, T, V, "$(s_type)_$(s_name)_discharge") + jm = get_jump_model(container) + soc_min = + PSY.get_storage_level_limits(storage).min * + PSY.get_storage_capacity(storage, PSY.SU) * + PSY.get_conversion_factor(storage) + if time_offset(T) == -1 + con[ci_name, 1] = JuMP.@constraint( + jm, + sustained_param_discharge * reserve_var[ci_name, 1] <= get_value(ic) - soc_min + ) + for t in time_steps[2:end] + con[ci_name, t] = JuMP.@constraint( + jm, + sustained_param_discharge * reserve_var[ci_name, t] <= + energy_var[ci_name, t - 1] - soc_min + ) + end + else # EndOfPeriod + for t in time_steps + con[ci_name, t] = JuMP.@constraint( + jm, + sustained_param_discharge * reserve_var[ci_name, t] <= + energy_var[ci_name, t] - soc_min + ) + end + end + end + return +end + +function _emit_coverage_constraint!( + container::OptimizationContainer, + ::Type{T}, + ::Type{V}, + ic::InitialCondition, + energy_var, + ci_name::String, + eff_in::Float64, + inv_eff_out::Float64, + fraction_of_hour::Float64, + resolution::Dates.Period, + storage::PSY.Storage, + time_steps::UnitRange{Int}, + service::PSY.Reserve{PSY.ReserveDown}, +) where {T <: _ReserveCoverageT, V <: PSY.HybridSystem} + s_type = typeof(service) + s_name = PSY.get_name(service) + num_periods = PSY.get_sustained_time(service) / Dates.value(Dates.Second(resolution)) + sustained_param_charge = eff_in * fraction_of_hour * num_periods + reserve_var = + get_variable( + container, + HybridStorageSubcomponentReserveVariable{ChargeSide}, + V, + "$(s_type)_$s_name", + ) + con = get_constraint(container, T, V, "$(s_type)_$(s_name)_charge") + soc_max = + PSY.get_storage_level_limits(storage).max * + PSY.get_storage_capacity(storage, PSY.SU) * + PSY.get_conversion_factor(storage) + jm = get_jump_model(container) + if time_offset(T) == -1 + con[ci_name, 1] = JuMP.@constraint( + jm, + sustained_param_charge * reserve_var[ci_name, 1] <= soc_max - get_value(ic) + ) + for t in time_steps[2:end] + con[ci_name, t] = JuMP.@constraint( + jm, + sustained_param_charge * reserve_var[ci_name, t] <= + soc_max - energy_var[ci_name, t - 1] + ) + end + else # EndOfPeriod + for t in time_steps + con[ci_name, t] = JuMP.@constraint( + jm, + sustained_param_charge * reserve_var[ci_name, t] <= + soc_max - energy_var[ci_name, t] + ) + end + end + return +end + +_emit_coverage_constraint!( + ::OptimizationContainer, + ::Type{<:_ReserveCoverageT}, + ::Type{<:PSY.HybridSystem}, + ::InitialCondition, + _energy_var, + ::String, + ::Float64, + ::Float64, + ::Float64, + ::Dates.Period, + ::PSY.Storage, + ::UnitRange{Int}, + ::PSY.Service, +) = nothing + +function add_constraints!( + container::OptimizationContainer, + ::Type{T}, + devices::Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + model::DeviceModel{V, HybridDispatchWithReserves}, + network_model::NetworkModel{X}, +) where { + T <: _ReserveCoverageT, + V <: PSY.HybridSystem, + X <: AbstractPowerModel, +} + time_steps = get_time_steps(container) + resolution = get_resolution(container) + fraction_of_hour = Dates.value(Dates.Minute(resolution)) / MINUTES_IN_HOUR + names = [PSY.get_name(d) for d in devices] + initial_conditions = get_initial_condition(container, InitialEnergyLevel(), V) + energy_var = get_variable(container, EnergyVariable, V) + + services_set = Set{PSY.Service}() + for ic in initial_conditions + d = IOM.get_component(ic) + union!(services_set, PSY.get_services(d)) + end + + for service in services_set + _init_coverage_container!(container, T, V, names, time_steps, service) + end + + for ic in initial_conditions + d = IOM.get_component(ic) + storage = PSY.get_storage(d) + storage === nothing && continue + ci_name = PSY.get_name(d) + eff_in = PSY.get_efficiency(storage).in + inv_eff_out = 1.0 / PSY.get_efficiency(storage).out + for service in PSY.get_services(d) + _emit_coverage_constraint!( + container, T, V, ic, energy_var, ci_name, + eff_in, inv_eff_out, fraction_of_hour, resolution, + storage, time_steps, service, + ) + end + end + return +end + +################################################################################# +################################################################################# +# HybridEnergyTargetConstraint on hybrids with energy_target=true. A soft equality +# (e_T - e^+ + e^- = E_T) with non-negative surplus/shortage slacks penalized in the +# objective. Mirrors the storage StateofChargeTargetConstraint; the target RHS is +# scaled to absolute energy units to match the hybrid EnergyVariable. +################################################################################# +################################################################################# + +function add_constraints!( + container::OptimizationContainer, + ::Type{HybridEnergyTargetConstraint}, + devices::Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, + model::DeviceModel{V, HybridDispatchWithReserves}, + ::NetworkModel{X}, +) where {V <: PSY.HybridSystem, X <: AbstractPowerModel} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + energy_var = get_variable(container, EnergyVariable, V) + surplus_var = get_variable(container, HybridEnergySurplusVariable, V) + shortage_var = get_variable(container, HybridEnergyShortageVariable, V) + constraint = add_constraints_container!( + container, + HybridEnergyTargetConstraint, + V, + names, + [last(time_steps)], + ) + for d in devices + storage = PSY.get_storage(d) + storage === nothing && continue + name = PSY.get_name(d) + target = + PSY.get_storage_target(storage) * + PSY.get_storage_capacity(storage, PSY.SU) * + PSY.get_conversion_factor(storage) + t_end = last(time_steps) + constraint[name, t_end] = JuMP.@constraint( + get_jump_model(container), + energy_var[name, t_end] - surplus_var[name, t_end] + + shortage_var[name, t_end] == target + ) + end + return +end + +################################################################################# +# Hybrid PCC ↔ subcomponent balance and reserve plumbing. +# +# These constraints are genuinely new — they have no analogue in POM's storage, +# thermal, or renewable code. They tie: +# - PCC active-power variables (ActivePowerOutVariable / ActivePowerInVariable) +# to the reservation variable (mutually exclusive charge/discharge at the +# hybrid boundary) +# - Internal subcomponent flows (thermal + renewable + storage discharge - +# storage charge - load) to the PCC injection +# - Per-subcomponent reserve allocations to the hybrid-boundary reserve +# variables, and the hybrid-boundary reserve variables to the system-level +# ActivePowerReserveVariable +################################################################################# + +# Plain range constraints on the PCC variables, used when `reservation = false`. +# When `reservation = true` the PCC mutual-exclusion is enforced by +# `HybridStatusOnConstraint{DischargeSide}` / `HybridStatusOnConstraint{ChargeSide}` instead. +function add_constraints!( + container::OptimizationContainer, + ::Type{T}, + ::Type{U}, + devices::IS.FlattenIteratorWrapper{V}, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + T <: Union{ + OutputActivePowerVariableLimitsConstraint, + InputActivePowerVariableLimitsConstraint, + }, + U <: Union{ActivePowerOutVariable, ActivePowerInVariable}, + V <: PSY.HybridSystem, + W <: AbstractHybridFormulation, + X <: AbstractPowerModel, +} + add_range_constraints!(container, T, U, devices, model, X) + return +end + +""" +Couple the hybrid PCC active-power variable to the reservation binary so that +only one direction is active at a time. `HybridStatusOnConstraint{DischargeSide}` enforces +`p_out ≤ reservation·max_out` (out-mode when reservation=1); `HybridStatusOnConstraint{ChargeSide}` +enforces `p_in ≤ (1 − reservation)·max_in` (in-mode when reservation=0). With +ancillary services attached, the asymmetric reserve expressions enter both +bounds — Out side picks up Out{ReserveUp,Down}; In side picks up In{ReserveDown,Up} — mirroring +HSS `_add_constraints_status{out,in}_withreserves!`. +""" +# Side-keyed traits for HybridStatusOnConstraint{Sd}. The reserve-expression mapping is +# asymmetric: DischargeSide UB picks up Out-ReserveUp, In side UB picks up In-Down (and vice-versa +# for LB). The reservation-binary factor is `reservation` for DischargeSide, `(1-reservation)` +# for ChargeSide (mirrors the storage Charge/Discharge ss_factor trait). +_pcc_power_var(::Type{HybridStatusOnConstraint{DischargeSide}}) = ActivePowerOutVariable +_pcc_power_var(::Type{HybridStatusOnConstraint{ChargeSide}}) = ActivePowerInVariable +_pcc_max_limit(::Type{HybridStatusOnConstraint{DischargeSide}}, d) = + PSY.get_output_active_power_limits(d, PSY.SU).max +_pcc_max_limit(::Type{HybridStatusOnConstraint{ChargeSide}}, d) = + PSY.get_input_active_power_limits(d, PSY.SU).max +_pcc_reserve_ub_expr(::Type{HybridStatusOnConstraint{DischargeSide}}) = + HybridPCCReserveExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide} +_pcc_reserve_ub_expr(::Type{HybridStatusOnConstraint{ChargeSide}}) = + HybridPCCReserveExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide} +_pcc_reserve_lb_expr(::Type{HybridStatusOnConstraint{DischargeSide}}) = + HybridPCCReserveExpression{PSY.ReserveDown, UnscaledReserve, DischargeSide} +_pcc_reserve_lb_expr(::Type{HybridStatusOnConstraint{ChargeSide}}) = + HybridPCCReserveExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide} +_pcc_reservation_factor(::Type{HybridStatusOnConstraint{DischargeSide}}, r_val) = r_val +_pcc_reservation_factor(::Type{HybridStatusOnConstraint{ChargeSide}}, r_val) = 1 - r_val + +# Helper: do PCC status constraints carry reserve terms? Type-dispatched, no body-level <:. +_pcc_has_reserves(::Type{<:AbstractHybridFormulationWithReserves}, model) = + has_service_model(model) +_pcc_has_reserves(::Type{<:AbstractHybridFormulation}, _model) = false + +function add_constraints!( + container::OptimizationContainer, + ::Type{T}, + devices::IS.FlattenIteratorWrapper{V}, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + T <: HybridStatusOnConstraint, + V <: PSY.HybridSystem, + W <: AbstractHybridFormulation, + X <: AbstractPowerModel, +} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + p_var = get_variable(container, _pcc_power_var(T), V) + reservation = get_variable(container, ReservationVariable, V) + constraint = add_constraints_container!(container, T, V, names, time_steps) + + has_reserves = _pcc_has_reserves(W, model) + r_ub, r_lb, con_lb = if has_reserves + ( + get_expression(container, _pcc_reserve_ub_expr(T), V), + get_expression(container, _pcc_reserve_lb_expr(T), V), + add_constraints_container!(container, T, V, names, time_steps; meta = "lb"), + ) + else + (nothing, nothing, nothing) + end + + for d in devices, t in time_steps + name = PSY.get_name(d) + max_p = _pcc_max_limit(T, d) + rhs_factor = _pcc_reservation_factor(T, reservation[name, t]) + if has_reserves + constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + p_var[name, t] + r_ub[name, t] <= rhs_factor * max_p + ) + con_lb[name, t] = JuMP.@constraint( + get_jump_model(container), + p_var[name, t] - r_lb[name, t] >= 0.0 + ) + else + constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + p_var[name, t] <= rhs_factor * max_p + ) + end + end + return +end + +""" +Energy asset balance: the hybrid's PCC injection equals the sum of subcomponent +injections (thermal + renewable + storage discharge - storage charge - load). +When ancillary services are attached, served (deployed-fraction) reserve expressions +also enter the balance with sign pattern `+out_up - in_up - out_down + in_down`, +mirroring HSS `_add_constraints_energyassetbalance_with_reserves!`. +""" +function add_constraints!( + container::OptimizationContainer, + ::Type{HybridEnergyAssetBalanceConstraint}, + devices::IS.FlattenIteratorWrapper{V}, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where {V <: PSY.HybridSystem, W <: AbstractHybridFormulation, X <: AbstractPowerModel} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + p_out = get_variable(container, ActivePowerOutVariable, V) + p_in = get_variable(container, ActivePowerInVariable, V) + constraint = add_constraints_container!( + container, + HybridEnergyAssetBalanceConstraint, + V, + names, + time_steps, + ) + + # Optional subcomponent variables — only present when the hybrid has them + p_th = if haskey(IOM.get_variables(container), VariableKey(HybridThermalActivePower, V)) + get_variable(container, HybridThermalActivePower, V) + else + nothing + end + p_re = + if haskey(IOM.get_variables(container), VariableKey(HybridRenewableActivePower, V)) + get_variable(container, HybridRenewableActivePower, V) + else + nothing + end + p_ch = + if haskey( + IOM.get_variables(container), + VariableKey(HybridStorageSubcomponentPower{ChargeSide}, V), + ) + get_variable(container, HybridStorageSubcomponentPower{ChargeSide}, V) + else + nothing + end + p_ds = + if haskey( + IOM.get_variables(container), + VariableKey(HybridStorageSubcomponentPower{DischargeSide}, V), + ) + get_variable(container, HybridStorageSubcomponentPower{DischargeSide}, V) + else + nothing + end + + load_param_container = + if haskey( + IOM.get_parameters(container), + ParameterKey(HybridElectricLoadTimeSeriesParameter, V), + ) + get_parameter(container, HybridElectricLoadTimeSeriesParameter, V) + else + nothing + end + load_multiplier = if load_param_container === nothing + nothing + else + get_multiplier_array(load_param_container) + end + + has_reserves = W <: AbstractHybridFormulationWithReserves && has_service_model(model) + serv_out_up, serv_out_dn, serv_in_up, serv_in_dn = if has_reserves + ( + get_expression( + container, + HybridPCCReserveExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}, + V, + ), + get_expression( + container, + HybridPCCReserveExpression{PSY.ReserveDown, DeployedReserve, DischargeSide}, + V, + ), + get_expression( + container, + HybridPCCReserveExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}, + V, + ), + get_expression( + container, + HybridPCCReserveExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}, + V, + ), + ) + else + (nothing, nothing, nothing, nothing) + end + + for d in devices, t in time_steps + name = PSY.get_name(d) + rhs = JuMP.AffExpr(0.0) + if p_th !== nothing && PSY.get_thermal_unit(d) !== nothing + add_proportional_to_jump_expression!(rhs, p_th[name, t], 1.0) + end + if p_re !== nothing && PSY.get_renewable_unit(d) !== nothing + add_proportional_to_jump_expression!(rhs, p_re[name, t], 1.0) + end + if p_ds !== nothing && PSY.get_storage(d) !== nothing + add_proportional_to_jump_expression!(rhs, p_ds[name, t], 1.0) + end + if p_ch !== nothing && PSY.get_storage(d) !== nothing + add_proportional_to_jump_expression!(rhs, p_ch[name, t], -1.0) + end + if load_param_container !== nothing && PSY.get_electric_load(d) !== nothing + load_ref = get_parameter_column_refs(load_param_container, name)[t] + add_proportional_to_jump_expression!(rhs, load_ref, -load_multiplier[name, t]) + end + if has_reserves + add_proportional_to_jump_expression!(rhs, serv_out_up[name, t], 1.0) + add_proportional_to_jump_expression!(rhs, serv_in_dn[name, t], 1.0) + add_proportional_to_jump_expression!(rhs, serv_out_dn[name, t], -1.0) + add_proportional_to_jump_expression!(rhs, serv_in_up[name, t], -1.0) + end + constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + p_out[name, t] - p_in[name, t] == rhs + ) + end + return +end + +""" +Couple the hybrid PCC reserve variables (Out + In, summed across subcomponents) +to the system-level `ActivePowerReserveVariable` for each service the hybrid +participates in. +""" +function add_constraints!( + container::OptimizationContainer, + ::Type{HybridReserveAssignmentConstraint}, + devices::IS.FlattenIteratorWrapper{V}, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + V <: PSY.HybridSystem, + W <: AbstractHybridFormulationWithReserves, + X <: AbstractPowerModel, +} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + + services = Set{PSY.Service}() + for d in devices + union!(services, PSY.get_services(d)) + end + + for service in services + s_name = PSY.get_name(service) + s_type = typeof(service) + constraint = + add_constraints_container!(container, HybridReserveAssignmentConstraint, V, + names, time_steps; + meta = "$(s_type)_$s_name") + # System-level reserve variable for this service + sys_reserve = get_variable(container, ActivePowerReserveVariable, s_type, s_name) + # Per-hybrid reserve variables for this service + r_out = get_variable( + container, + HybridPCCReserveVariable{DischargeSide}, + V, + "$(s_type)_$s_name", + ) + r_in = get_variable( + container, + HybridPCCReserveVariable{ChargeSide}, + V, + "$(s_type)_$s_name", + ) + for d in devices, t in time_steps + name = PSY.get_name(d) + (service in PSY.get_services(d)) || continue + constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + r_out[name, t] + r_in[name, t] == sys_reserve[name, t] + ) + end + end + return +end + +""" +Couple the hybrid PCC reserve variables (Out + In) to the sum of per-subcomponent +reserve allocations (thermal + renewable + charging + discharging). +""" +function add_constraints!( + container::OptimizationContainer, + ::Type{HybridReserveBalanceConstraint}, + devices::IS.FlattenIteratorWrapper{V}, + model::DeviceModel{V, W}, + ::NetworkModel{X}, +) where { + V <: PSY.HybridSystem, + W <: AbstractHybridFormulationWithReserves, + X <: AbstractPowerModel, +} + time_steps = get_time_steps(container) + names = [PSY.get_name(d) for d in devices] + + services = Set{PSY.Service}() + for d in devices + union!(services, PSY.get_services(d)) + end + + for service in services + s_name = PSY.get_name(service) + s_type = typeof(service) + constraint = + add_constraints_container!(container, HybridReserveBalanceConstraint, V, names, + time_steps; + meta = "$(s_type)_$s_name") + r_out = get_variable( + container, + HybridPCCReserveVariable{DischargeSide}, + V, + "$(s_type)_$s_name", + ) + r_in = get_variable( + container, + HybridPCCReserveVariable{ChargeSide}, + V, + "$(s_type)_$s_name", + ) + for d in devices, t in time_steps + name = PSY.get_name(d) + (service in PSY.get_services(d)) || continue + rhs = JuMP.AffExpr(0.0) + for var_t in (HybridThermalReserveVariable, HybridRenewableReserveVariable, + HybridStorageSubcomponentReserveVariable{ChargeSide}, + HybridStorageSubcomponentReserveVariable{DischargeSide}) + key = VariableKey(var_t, V, "$(s_type)_$s_name") + if haskey(IOM.get_variables(container), key) + var = get_variable(container, key) + add_proportional_to_jump_expression!(rhs, var[name, t], 1.0) + end + end + constraint[name, t] = JuMP.@constraint( + get_jump_model(container), + r_out[name, t] + r_in[name, t] == rhs + ) + end + end + return +end +################################################################################# +# Objective function for HybridSystem. +# +# The hybrid envelope itself carries `MarketBidCost(nothing)`. The variable +# costs come from the subcomponents: +# - Thermal subcomponent → ThermalGenerationCost (PSY 5.x) +# - Renewable subcomponent → RenewableGenerationCost +# - Storage subcomponent → StorageCost +# +# We walk into each subcomponent's operation_cost, then delegate to IOM's +# add_variable_cost_to_objective! with the subcomponent's cost data and the +# hybrid-specific variable type. The variable is keyed by the hybrid's name +# (not the subcomponent's), but the cost data drives the linear/piecewise terms. +################################################################################# + +function _add_hybrid_subcomponent_variable_cost!( + container::OptimizationContainer, + ::Type{V}, + devices, + accessor::Function, + ::Type{W}, +) where {V <: VariableType, W <: AbstractHybridFormulation} + for d in devices + sub = accessor(d) + sub === nothing && continue + op_cost = PSY.get_operation_cost(sub) + # Use IOM's add_variable_cost_to_objective! with the hybrid device + # but the subcomponent's cost data. The dispatch on (V, HybridSystem, W) + # is what variable_cost(op_cost, V, HybridSystem, W) needs to see. + add_variable_cost_to_objective!(container, V, d, op_cost, W) + end + return +end + +# Hybrid `OnVariable` proportional cost — delegate to the standalone thermal +# `proportional_cost` so a hybrid-embedded thermal unit and a standalone copy with the +# same `ThermalGenerationCost` produce identical objective coefficients +# (`onvar_cost + vom_constant + fixed`). Implemented via the same IOM cost-term +# helpers (`add_cost_term_invariant!` / `add_cost_term_variant!`) that +# `add_proportional_cost_maybe_time_variant!` uses, so the time-variant fuel-cost +# branch lights up when the embedded thermal cost is backed by a time series. +function add_proportional_cost!( + container::OptimizationContainer, + ::Type{OnVariable}, + devices::Vector{D}, + ::Type{W}, +) where {D <: PSY.HybridSystem, W <: AbstractHybridFormulation} + multiplier = objective_function_multiplier(OnVariable, W) + on_var = get_variable(container, OnVariable, D) + for d in devices + thermal = PSY.get_thermal_unit(d) + thermal === nothing && continue + thermal_cost = PSY.get_operation_cost(thermal) + thermal_cost === nothing && continue + # Select the variant- vs invariant-aware cost-term writer once per device, then + # call it inside the time loop. Mirrors IOM.add_proportional_cost_maybe_time_variant!. + add_cost_term! = if IOM.is_time_variant_proportional(thermal_cost) + add_cost_term_variant! + else + add_cost_term_invariant! + end + name = PSY.get_name(d) + for t in get_time_steps(container) + cost_term = proportional_cost( + container, thermal_cost, OnVariable, thermal, + ThermalBasicUnitCommitment, t, + ) + iszero(cost_term) && continue + rate = cost_term * multiplier + add_cost_term!( + container, on_var[name, t], rate, ProductionCostExpression, D, name, t, + ) + end + end + return +end + +function objective_function!( + container::OptimizationContainer, + devices::U, + model::DeviceModel{D, W}, + ::Type{<:AbstractPowerModel}, +) where { + U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}, + W <: AbstractHybridFormulation, +} where {D <: PSY.HybridSystem} + devices_vec = collect(devices) + hybrids_with_thermal = [d for d in devices_vec if PSY.get_thermal_unit(d) !== nothing] + hybrids_with_renewable = + [d for d in devices_vec if PSY.get_renewable_unit(d) !== nothing] + hybrids_with_storage = [d for d in devices_vec if PSY.get_storage(d) !== nothing] + + # Thermal: variable cost on HybridThermalActivePower; OnVariable proportional cost + # routed through POM's standard add_proportional_cost! pathway so hybrids match + # standalone thermal exactly (onvar_cost + vom_constant + fixed). + if !isempty(hybrids_with_thermal) + _add_hybrid_subcomponent_variable_cost!(container, HybridThermalActivePower, + hybrids_with_thermal, PSY.get_thermal_unit, W) + add_proportional_cost!(container, OnVariable, hybrids_with_thermal, W) + end + + # Renewable: variable cost on HybridRenewableActivePower (typically a curtailment cost) + if !isempty(hybrids_with_renewable) + _add_hybrid_subcomponent_variable_cost!(container, HybridRenewableActivePower, + hybrids_with_renewable, PSY.get_renewable_unit, W) + end + + # Storage: variable costs on charge/discharge, plus optional regularization penalty. + if !isempty(hybrids_with_storage) + _add_hybrid_subcomponent_variable_cost!(container, + HybridStorageSubcomponentPower{ChargeSide}, + hybrids_with_storage, PSY.get_storage, W) + _add_hybrid_subcomponent_variable_cost!(container, + HybridStorageSubcomponentPower{DischargeSide}, + hybrids_with_storage, PSY.get_storage, W) + if get_attribute(model, "regularization") + _add_hybrid_regularization_cost!( + container, RegularizationVariable{ChargeSide}, hybrids_with_storage, W) + _add_hybrid_regularization_cost!( + container, RegularizationVariable{DischargeSide}, hybrids_with_storage, W) + end + if get_attribute(model, "energy_target") + _add_hybrid_energy_target_cost!( + container, HybridEnergySurplusVariable, hybrids_with_storage, W) + _add_hybrid_energy_target_cost!( + container, HybridEnergyShortageVariable, hybrids_with_storage, W) + end + end + return +end + +# Routes regularization slacks through IOM's add_cost_term_invariant! so the penalty +# lands in both the objective and (when present) the production-cost expression. +function _add_hybrid_regularization_cost!( + container::OptimizationContainer, + ::Type{V}, + devices::Vector{D}, + ::Type{W}, +) where {V <: VariableType, D <: PSY.HybridSystem, W <: AbstractHybridFormulation} + multiplier = objective_function_multiplier(V, W) + rate = HYBRID_REGULARIZATION_COST * multiplier + var = get_variable(container, V, D) + for d in devices + PSY.get_storage(d) === nothing && continue + name = PSY.get_name(d) + for t in get_time_steps(container) + add_cost_term_invariant!( + container, var[name, t], rate, ProductionCostExpression, D, name, t, + ) + end + end + return +end + +# Penalizes the end-of-period energy-target slacks. The per-unit cost comes from the +# storage subcomponent's `StorageCost` (energy_surplus_cost / energy_shortage_cost), +# mirroring POM storage's StateofChargeTargetConstraint objective handling. Slacks live +# only at the final time step. +function _add_hybrid_energy_target_cost!( + container::OptimizationContainer, + ::Type{V}, + devices::Vector{D}, + ::Type{W}, +) where {V <: VariableType, D <: PSY.HybridSystem, W <: AbstractHybridFormulation} + multiplier = objective_function_multiplier(V, W) + var = get_variable(container, V, D) + t_end = last(get_time_steps(container)) + for d in devices + storage = PSY.get_storage(d) + storage === nothing && continue + name = PSY.get_name(d) + op_cost = PSY.get_operation_cost(storage) + cost_term = proportional_cost(op_cost, V, d, W) * multiplier + add_cost_term_invariant!( + container, var[name, t_end], cost_term, ProductionCostExpression, D, name, t_end, + ) + end + return +end + +################################################################################# +# IOM.variable_cost dispatches — reach into subcomponent cost types +################################################################################# + +# Thermal subcomponent variable cost +IOM.variable_cost( + cost::PSY.ThermalGenerationCost, + ::Type{HybridThermalActivePower}, + ::Type{<:PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_variable(cost) + +# Renewable subcomponent variable cost (typically a curtailment penalty) +IOM.variable_cost( + cost::PSY.RenewableGenerationCost, + ::Type{HybridRenewableActivePower}, + ::Type{<:PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_curtailment_cost(cost) + +# Storage subcomponent variable costs +IOM.variable_cost( + cost::PSY.StorageCost, + ::Type{HybridStorageSubcomponentPower{ChargeSide}}, + ::Type{<:PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_charge_variable_cost(cost) + +IOM.variable_cost( + cost::PSY.StorageCost, + ::Type{HybridStorageSubcomponentPower{DischargeSide}}, + ::Type{<:PSY.HybridSystem}, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_discharge_variable_cost(cost) + +# End-of-period energy-target slack penalties, pulled from the storage subcomponent's +# StorageCost. Mirrors POM storage (energy_storage_models/storage_models.jl:74-75). +proportional_cost( + cost::PSY.StorageCost, + ::Type{HybridEnergySurplusVariable}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_energy_surplus_cost(cost) + +proportional_cost( + cost::PSY.StorageCost, + ::Type{HybridEnergyShortageVariable}, + ::PSY.HybridSystem, + ::Type{<:AbstractHybridFormulation}, +) = PSY.get_energy_shortage_cost(cost) diff --git a/src/hybrid_system_models/hybridsystem_constructor.jl b/src/hybrid_system_models/hybridsystem_constructor.jl new file mode 100644 index 0000000..2b280d3 --- /dev/null +++ b/src/hybrid_system_models/hybridsystem_constructor.jl @@ -0,0 +1,601 @@ +################################################################################# +# Two-stage construct_device! for HybridSystem with HybridDispatchWithReserves. +# +# Argument stage: variables, parameters, expression containers, expression +# accumulation, feedforward arguments. +# Model stage: constraints, feedforward constraints, objective function, +# dual recording. +# +# Mirrors energy_storage_models/storage_constructor.jl in structure. We provide +# both AbstractPowerModel (with reactive power) and AbstractActivePowerModel +# (without) variants. +################################################################################# + +function _filter_hybrids(devices) + devices_vec = collect(devices) + return ( + all = devices_vec, + with_thermal = [d for d in devices_vec if PSY.get_thermal_unit(d) !== nothing], + with_renewable = [d for d in devices_vec if PSY.get_renewable_unit(d) !== nothing], + with_storage = [d for d in devices_vec if PSY.get_storage(d) !== nothing], + with_load = [d for d in devices_vec if PSY.get_electric_load(d) !== nothing], + ) +end + +function _add_hybrid_reserve_arguments!( + container::OptimizationContainer, + devices, + hybrids_with_storage, + hybrids_with_thermal, + hybrids_with_renewable, + model::DeviceModel{T, D}, + network_model::NetworkModel{S}, +) where {T <: PSY.HybridSystem, D <: HybridDispatchWithReserves, S <: AbstractPowerModel} + time_steps = get_time_steps(container) + + # Hybrid PCC reserve variables + add_variables!(container, HybridPCCReserveVariable{DischargeSide}, devices, D) + add_variables!(container, HybridPCCReserveVariable{ChargeSide}, devices, D) + + # Allocate hybrid-boundary aggregation expression containers + for E in ( + HybridPCCReserveExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide}, + HybridPCCReserveExpression{PSY.ReserveDown, UnscaledReserve, DischargeSide}, + HybridPCCReserveExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide}, + HybridPCCReserveExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide}, + HybridPCCReserveExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}, + HybridPCCReserveExpression{PSY.ReserveDown, DeployedReserve, DischargeSide}, + HybridPCCReserveExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}, + HybridPCCReserveExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}, + ) + lazy_container_addition!(container, E, T, PSY.get_name.(devices), time_steps) + end + + # Accumulate Out/In reserve variables into Total* and Served* expressions + for E in ( + HybridPCCReserveExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide}, + HybridPCCReserveExpression{PSY.ReserveDown, UnscaledReserve, DischargeSide}, + HybridPCCReserveExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}, + HybridPCCReserveExpression{PSY.ReserveDown, DeployedReserve, DischargeSide}, + ) + add_to_expression!( + container, + E, + HybridPCCReserveVariable{DischargeSide}, + devices, + model, + network_model, + ) + end + for E in ( + HybridPCCReserveExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide}, + HybridPCCReserveExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide}, + HybridPCCReserveExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}, + HybridPCCReserveExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}, + ) + add_to_expression!( + container, + E, + HybridPCCReserveVariable{ChargeSide}, + devices, + model, + network_model, + ) + end + + # Per-subcomponent reserve variables + if !isempty(hybrids_with_thermal) + add_variables!(container, HybridThermalReserveVariable, hybrids_with_thermal, D) + end + if !isempty(hybrids_with_renewable) + add_variables!(container, HybridRenewableReserveVariable, hybrids_with_renewable, D) + end + if !isempty(hybrids_with_storage) + add_variables!( + container, + HybridStorageSubcomponentReserveVariable{ChargeSide}, + hybrids_with_storage, + D, + ) + add_variables!( + container, + HybridStorageSubcomponentReserveVariable{DischargeSide}, + hybrids_with_storage, + D, + ) + + # Storage-side reserve expression containers, keyed by HybridSystem + for E in ( + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide}, + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide}, + StorageReserveBalanceExpression{ + PSY.ReserveDown, + UnscaledReserve, + DischargeSide, + }, + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide}, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}, + StorageReserveBalanceExpression{ + PSY.ReserveDown, + DeployedReserve, + DischargeSide, + }, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}, + ) + lazy_container_addition!( + container, + E, + T, + PSY.get_name.(hybrids_with_storage), + time_steps, + ) + end + + # Wire HybridStorageSubcomponentReserveVariable{DischargeSide} into Discharge expressions + for E in ( + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, DischargeSide}, + StorageReserveBalanceExpression{ + PSY.ReserveDown, + UnscaledReserve, + DischargeSide, + }, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, DischargeSide}, + StorageReserveBalanceExpression{ + PSY.ReserveDown, + DeployedReserve, + DischargeSide, + }, + ) + add_to_expression!( + container, + E, + HybridStorageSubcomponentReserveVariable{DischargeSide}, + hybrids_with_storage, + model, + ) + end + for E in ( + StorageReserveBalanceExpression{PSY.ReserveUp, UnscaledReserve, ChargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, UnscaledReserve, ChargeSide}, + StorageReserveBalanceExpression{PSY.ReserveUp, DeployedReserve, ChargeSide}, + StorageReserveBalanceExpression{PSY.ReserveDown, DeployedReserve, ChargeSide}, + ) + add_to_expression!( + container, + E, + HybridStorageSubcomponentReserveVariable{ChargeSide}, + hybrids_with_storage, + model, + ) + end + end + + # TotalReserveOffering aggregation per service, keyed by HybridSystem. Created for + # EVERY hybrid that participates in a reserve service, storage-less hybrids included: + # get_expression_type_for_reserve routes all hybrids' ActivePowerReserveVariable into + # TotalReserveOffering, so the container must exist regardless of storage. + services = Set{PSY.Service}() + for d in devices + union!(services, PSY.get_services(d)) + end + for s in services + lazy_container_addition!(container, TotalReserveOffering, T, + PSY.get_name.(devices), time_steps; + meta = "$(typeof(s))_$(PSY.get_name(s))") + end + # Only storage hybrids have subcomponent reserve variables to aggregate into the + # offering; storage-less hybrids feed it via the PCC reserve path instead. + if !isempty(hybrids_with_storage) + for v in ( + HybridStorageSubcomponentReserveVariable{ChargeSide}, + HybridStorageSubcomponentReserveVariable{DischargeSide}, + ) + add_to_expression!( + container, + TotalReserveOffering, + v, + hybrids_with_storage, + model, + ) + end + end + return +end + +_maybe_add_reactive_power_variable!( + container::OptimizationContainer, + devices::U, + ::Type{D}, + ::Type{<:AbstractPowerModel}, +) where { + D <: AbstractHybridFormulation, + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, +} where {V <: PSY.HybridSystem} = + add_variables!(container, ReactivePowerVariable, devices, D) + +_maybe_add_reactive_power_balance!( + container::OptimizationContainer, + devices::U, + model::DeviceModel{V, D}, + network_model::NetworkModel{<:AbstractPowerModel}, +) where { + D <: AbstractHybridFormulation, + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, +} where {V <: PSY.HybridSystem} = + add_to_expression!( + container, + ReactivePowerBalance, + ReactivePowerVariable, + devices, + model, + network_model, + ) + +_maybe_add_reactive_power_limits!( + container::OptimizationContainer, + devices::U, + model::DeviceModel{V, D}, + network_model::NetworkModel{<:AbstractPowerModel}, +) where { + D <: AbstractHybridFormulation, + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, +} where {V <: PSY.HybridSystem} = + add_constraints!( + container, + ReactivePowerVariableLimitsConstraint, + ReactivePowerVariable, + devices, + model, + network_model, + ) + +_maybe_add_reactive_power_variable!( + ::OptimizationContainer, + ::U, + ::Type{D}, + ::Type{<:AbstractActivePowerModel}, +) where { + D <: AbstractHybridFormulation, + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, +} where {V <: PSY.HybridSystem} = nothing + +_maybe_add_reactive_power_balance!( + ::OptimizationContainer, + ::U, + ::DeviceModel{V, D}, + ::NetworkModel{<:AbstractActivePowerModel}, +) where { + D <: AbstractHybridFormulation, + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, +} where {V <: PSY.HybridSystem} = nothing + +_maybe_add_reactive_power_limits!( + ::OptimizationContainer, + ::U, + ::DeviceModel{V, D}, + ::NetworkModel{<:AbstractActivePowerModel}, +) where { + D <: AbstractHybridFormulation, + U <: Union{Vector{V}, IS.FlattenIteratorWrapper{V}}, +} where {V <: PSY.HybridSystem} = nothing + +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ArgumentConstructStage, + model::DeviceModel{T, D}, + network_model::NetworkModel{S}, +) where {T <: PSY.HybridSystem, D <: HybridDispatchWithReserves, S <: AbstractPowerModel} + devices = get_available_components(model, sys) + grouped = _filter_hybrids(devices) + + # PCC variables + add_variables!(container, ActivePowerOutVariable, devices, D) + add_variables!(container, ActivePowerInVariable, devices, D) + _maybe_add_reactive_power_variable!(container, devices, D, S) + if get_attribute(model, "reservation") + add_variables!(container, ReservationVariable, devices, D) + end + + add_to_expression!( + container, + ActivePowerBalance, + ActivePowerInVariable, + devices, + model, + network_model, + ) + add_to_expression!( + container, + ActivePowerBalance, + ActivePowerOutVariable, + devices, + model, + network_model, + ) + _maybe_add_reactive_power_balance!(container, devices, model, network_model) + + # Subcomponent variables + if !isempty(grouped.with_thermal) + add_variables!(container, HybridThermalActivePower, grouped.with_thermal, D) + add_variables!(container, OnVariable, grouped.with_thermal, D) + end + if !isempty(grouped.with_renewable) + add_variables!(container, HybridRenewableActivePower, grouped.with_renewable, D) + add_parameters!( + container, + HybridRenewableActivePowerTimeSeriesParameter, + grouped.with_renewable, + model, + ) + end + if !isempty(grouped.with_storage) + add_variables!( + container, + HybridStorageSubcomponentPower{ChargeSide}, + grouped.with_storage, + D, + ) + add_variables!( + container, + HybridStorageSubcomponentPower{DischargeSide}, + grouped.with_storage, + D, + ) + add_variables!(container, EnergyVariable, grouped.with_storage, D) + if get_attribute(model, "storage_reservation") + add_variables!(container, HybridStorageReservation, grouped.with_storage, D) + end + if get_attribute(model, "regularization") + add_variables!( + container, + RegularizationVariable{ChargeSide}, + grouped.with_storage, + D, + ) + add_variables!( + container, + RegularizationVariable{DischargeSide}, + grouped.with_storage, + D, + ) + end + if get_attribute(model, "energy_target") + add_variables!(container, HybridEnergySurplusVariable, grouped.with_storage, D) + add_variables!(container, HybridEnergyShortageVariable, grouped.with_storage, D) + end + initial_conditions!(container, devices, D()) + end + if !isempty(grouped.with_load) + add_parameters!( + container, + HybridElectricLoadTimeSeriesParameter, + grouped.with_load, + model, + ) + end + + if has_service_model(model) + _add_hybrid_reserve_arguments!(container, devices, + grouped.with_storage, grouped.with_thermal, grouped.with_renewable, + model, network_model) + end + + add_feedforward_arguments!(container, model, devices) + add_event_arguments!(container, devices, model, network_model) + return +end + +function construct_device!( + container::OptimizationContainer, + sys::PSY.System, + ::ModelConstructStage, + model::DeviceModel{T, D}, + network_model::NetworkModel{S}, +) where {T <: PSY.HybridSystem, D <: HybridDispatchWithReserves, S <: AbstractPowerModel} + devices = get_available_components(model, sys) + grouped = _filter_hybrids(devices) + + # PCC reactive-power limits (active-power limits handled via the asset balance + status constraints) + _maybe_add_reactive_power_limits!(container, devices, model, network_model) + + # PCC ↔ subcomponent plumbing. With reservation, mutual-exclusion via ReservationVariable; + # without, simple range constraints on the PCC variables (no binary). + if get_attribute(model, "reservation") + add_constraints!( + container, + HybridStatusOnConstraint{DischargeSide}, + devices, + model, + network_model, + ) + add_constraints!( + container, + HybridStatusOnConstraint{ChargeSide}, + devices, + model, + network_model, + ) + else + add_constraints!( + container, + OutputActivePowerVariableLimitsConstraint, + ActivePowerOutVariable, + devices, + model, + network_model, + ) + add_constraints!( + container, + InputActivePowerVariableLimitsConstraint, + ActivePowerInVariable, + devices, + model, + network_model, + ) + end + add_constraints!( + container, + HybridEnergyAssetBalanceConstraint, + devices, + model, + network_model, + ) + + # Thermal subcomponent + if !isempty(grouped.with_thermal) + if has_service_model(model) + add_constraints!( + container, + HybridThermalReserveLimitConstraint, + grouped.with_thermal, + model, + network_model, + ) + else + add_constraints!( + container, + HybridThermalOnVariableConstraint{UpperBound}, + grouped.with_thermal, + model, + network_model, + ) + add_constraints!( + container, + HybridThermalOnVariableConstraint{LowerBound}, + grouped.with_thermal, + model, + network_model, + ) + end + end + + # Renewable subcomponent + if !isempty(grouped.with_renewable) + add_constraints!( + container, + HybridRenewableActivePowerLimitConstraint, + grouped.with_renewable, + model, + network_model, + ) + if has_service_model(model) + add_constraints!( + container, + HybridRenewableReserveLimitConstraint, + grouped.with_renewable, + model, + network_model, + ) + end + end + + # Storage subcomponent + if !isempty(grouped.with_storage) + add_constraints!( + container, + HybridStorageBalanceConstraint, + grouped.with_storage, + model, + network_model, + ) + if get_attribute(model, "energy_target") + add_constraints!( + container, + HybridEnergyTargetConstraint, + grouped.with_storage, + model, + network_model, + ) + end + if has_service_model(model) + add_constraints!( + container, + ReserveCoverageConstraint, + grouped.with_storage, + model, + network_model, + ) + add_constraints!( + container, + ReserveCoverageConstraintEndOfPeriod, + grouped.with_storage, + model, + network_model, + ) + add_constraints!( + container, + HybridStorageReservePowerLimitConstraint{ChargeSide}, + grouped.with_storage, + model, + network_model, + ) + add_constraints!( + container, + HybridStorageReservePowerLimitConstraint{DischargeSide}, + grouped.with_storage, + model, + network_model, + ) + elseif get_attribute(model, "storage_reservation") + # No reserves attached: enforce mutual exclusion via HybridStorageReservation. + # When `storage_reservation` is false, charge/discharge are bounded + # independently by their variable upper bounds — no extra constraint needed. + add_constraints!( + container, + HybridStorageStatusOnConstraint{ChargeSide}, + grouped.with_storage, + model, + network_model, + ) + add_constraints!( + container, + HybridStorageStatusOnConstraint{DischargeSide}, + grouped.with_storage, + model, + network_model, + ) + end + if get_attribute(model, "regularization") + add_constraints!( + container, + RegularizationConstraint{ChargeSide}, + grouped.with_storage, + model, + network_model, + ) + add_constraints!( + container, + RegularizationConstraint{DischargeSide}, + grouped.with_storage, + model, + network_model, + ) + end + end + + # Hybrid-boundary reserve coupling + if has_service_model(model) + add_constraints!( + container, + HybridReserveAssignmentConstraint, + devices, + model, + network_model, + ) + add_constraints!( + container, + HybridReserveBalanceConstraint, + devices, + model, + network_model, + ) + end + + add_feedforward_constraints!(container, model, devices) + add_event_constraints!(container, devices, model, network_model) + objective_function!(container, devices, model, S) + add_constraint_dual!(container, sys, model) + return +end diff --git a/test/Project.toml b/test/Project.toml index a1d00b2..5e1cf87 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -21,6 +21,7 @@ PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd" PowerOperationsModels = "bed98974-b02a-5e2f-9ee0-a103f5c450dd" PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e" PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd" +PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" @@ -33,7 +34,8 @@ UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [sources] InfrastructureSystems = {rev = "IS4", url = "https://github.com/Sienna-Platform/InfrastructureSystems.jl"} -PowerSystems = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} +# TODO: Move to main once this branch is merged +PowerSystems = {rev = "ac/hybridsystem-strip-units", url = "https://github.com/Sienna-Platform/PowerSystems.jl"} InfrastructureOptimizationModels = {rev = "main", url = "https://github.com/Sienna-Platform/InfrastructureOptimizationModels.jl"} PowerSystemCaseBuilder = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerSystemCaseBuilder.jl"} PowerNetworkMatrices = {rev = "psy6", url = "https://github.com/Sienna-Platform/PowerNetworkMatrices.jl"} diff --git a/test/includes.jl b/test/includes.jl index abfc6c8..5ea2677 100644 --- a/test/includes.jl +++ b/test/includes.jl @@ -56,6 +56,7 @@ include("test_utils/mbc_system_utils.jl") include("test_utils/mbc_math_helpers.jl") include("test_utils/iec_test_systems.jl") include("test_utils/hydro_testing_utils.jl") +include("test_utils/hybrid_test_utils.jl") ENV["RUNNING_SIENNA_TESTS"] = "true" ENV["SIENNA_RANDOM_SEED"] = 1234 # Set a fixed seed for reproducibility in tests diff --git a/test/runtests.jl b/test/runtests.jl index 8b924be..1e8242b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,6 +14,7 @@ const LOG_FILE = "power-simulations-test.log" const DISABLED_TEST_FILES = [ # Can generate with ls -1 test | grep "test_.*.jl" # "test_device_branch_constructors.jl", # "test_device_hvdc.jl", +# "test_device_hybrid_constructors.jl", # "test_device_hydro_constructors.jl", # "test_device_lcc.jl", # "test_device_load_constructors.jl", diff --git a/test/test_device_hybrid_constructors.jl b/test/test_device_hybrid_constructors.jl new file mode 100644 index 0000000..65323c8 --- /dev/null +++ b/test/test_device_hybrid_constructors.jl @@ -0,0 +1,378 @@ +# Tests for HybridSystem device formulations. + +# Helpers shared across testsets ---------------------------------------------- + +const _NON_HYBRID_RESERVES = ("Spin_Up_R1", "Spin_Up_R2") + +# These tests assert structural facts (variable presence/absence) and directional +# objective invariants, neither of which needs the full 24-hour RTS-GMLC horizon. +# A short horizon keeps every model feasible while cutting each MILP solve from +# ~50-80s to well under a second. +const _HYBRID_HORIZON = Hour(3) + +function _build_hybrid_test_system(; + with_reserves::Bool = true, + with_thermal::Bool = true, + with_renewable::Bool = true, + with_storage::Bool = true, + with_load::Bool = true, + energy_target::Bool = false, +) + sys = PSB.build_system(PSB.PSITestSystems, "test_RTS_GMLC_sys") + modify_ren_curtailment_cost!(sys) + hybrid = add_hybrid_to_chuhsi_bus!(sys; + with_thermal = with_thermal, + with_renewable = with_renewable, + with_storage = with_storage, + with_load = with_load, + energy_target = energy_target, + ) + if with_reserves + for s in PSY.get_components(PSY.VariableReserve, sys) + s_name = PSY.get_name(s) + any(occursin(prefix, s_name) for prefix in _NON_HYBRID_RESERVES) && continue + PSY.add_service!(hybrid, s, sys) + end + end + return sys, hybrid +end + +function _build_hybrid_template( + sys; + attributes::Dict{String, Any} = Dict{String, Any}(), + with_reserves::Bool = true, +) + template = PowerOperationsProblemTemplate(POM.CopperPlatePowerModel) + POM.set_device_model!(template, PSY.ThermalStandard, POM.ThermalStandardUnitCommitment) + POM.set_device_model!(template, PSY.RenewableDispatch, POM.RenewableFullDispatch) + POM.set_device_model!(template, PSY.PowerLoad, POM.StaticPowerLoad) + POM.set_device_model!(template, + POM.DeviceModel(PSY.HybridSystem, POM.HybridDispatchWithReserves; + attributes = attributes), + ) + if with_reserves + for service in PSY.get_components(PSY.VariableReserve, sys) + s_name = PSY.get_name(service) + any(occursin(prefix, s_name) for prefix in _NON_HYBRID_RESERVES) && continue + POM.set_service_model!(template, + POM.ServiceModel(typeof(service), POM.RangeReserve, s_name), + ) + end + end + return template +end + +function _build_and_solve(template, sys) + m = POM.DecisionModel(template, sys; + optimizer = HiGHS_optimizer, horizon = _HYBRID_HORIZON) + @test POM.build!(m; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + @test POM.solve!(m) == IOM.RunStatus.SUCCESSFULLY_FINALIZED + return m +end + +_var_keys(model) = keys(IOM.get_variables(IOM.get_optimization_container(model))) +_obj(model) = IOM.get_optimization_container(model).optimizer_stats.objective_value + +# --------------------------------------------------------------------------- +# 9a. Build+solve smoke tests for each attribute / structural perturbation. +# --------------------------------------------------------------------------- + +@testset "Test HybridSystem DispatchWithReserves DeviceModel" begin + sys, _ = _build_hybrid_test_system() + template = _build_hybrid_template(sys) + _build_and_solve(template, sys) +end + +@testset "HybridDispatchWithReserves: reservation = false" begin + sys, _ = _build_hybrid_test_system() + template = _build_hybrid_template(sys; + attributes = Dict{String, Any}("reservation" => false)) + m = _build_and_solve(template, sys) + # ReservationVariable must NOT have been created. + @test !any( + k -> + IOM.get_entry_type(k) === POM.ReservationVariable && + IOM.get_component_type(k) === PSY.HybridSystem, _var_keys(m)) +end + +@testset "HybridDispatchWithReserves: storage_reservation = false" begin + sys, _ = _build_hybrid_test_system() + template = _build_hybrid_template(sys; + attributes = Dict{String, Any}("storage_reservation" => false)) + m = _build_and_solve(template, sys) + @test !any( + k -> + IOM.get_entry_type(k) === POM.HybridStorageReservation && + IOM.get_component_type(k) === PSY.HybridSystem, _var_keys(m)) +end + +@testset "HybridDispatchWithReserves: regularization = true" begin + sys, _ = _build_hybrid_test_system() + template = _build_hybrid_template(sys; + attributes = Dict{String, Any}("regularization" => true)) + m = _build_and_solve(template, sys) + @test any( + k -> IOM.get_entry_type(k) === POM.RegularizationVariable{ChargeSide}, + _var_keys(m), + ) + @test any( + k -> IOM.get_entry_type(k) === POM.RegularizationVariable{DischargeSide}, + _var_keys(m), + ) +end + +@testset "HybridDispatchWithReserves: energy_target = true (soft equality + slacks)" begin + sys, _ = _build_hybrid_test_system(; energy_target = true) + template = _build_hybrid_template(sys; + attributes = Dict{String, Any}("energy_target" => true)) + m = _build_and_solve(template, sys) + + # Both end-of-period slack variables must be created, keyed by HybridSystem. This is + # the check that would have caught the original port dropping the slacks. + @test any( + k -> + IOM.get_entry_type(k) === POM.HybridEnergySurplusVariable && + IOM.get_component_type(k) === PSY.HybridSystem, _var_keys(m)) + @test any( + k -> + IOM.get_entry_type(k) === POM.HybridEnergyShortageVariable && + IOM.get_component_type(k) === PSY.HybridSystem, _var_keys(m)) + + # The target is a soft EQUALITY (e_T - e^+ + e^- = E_T), not a one-sided floor (>=). + container = IOM.get_optimization_container(m) + con_key = IOM.ConstraintKey(POM.HybridEnergyTargetConstraint, PSY.HybridSystem) + target_cons = IOM.get_constraints(container)[con_key] + @test !isempty(target_cons) + @test all(JuMP.constraint_object(c).set isa MOI.EqualTo for c in target_cons) +end + +@testset "HybridDispatchWithReserves: energy_target = false omits slacks" begin + sys, _ = _build_hybrid_test_system() + template = _build_hybrid_template(sys) + m = _build_and_solve(template, sys) + @test !any( + k -> IOM.get_entry_type(k) === POM.HybridEnergySurplusVariable, _var_keys(m)) + @test !any( + k -> IOM.get_entry_type(k) === POM.HybridEnergyShortageVariable, _var_keys(m)) +end + +@testset "HybridDispatchWithReserves: no reserves attached" begin + sys, _ = _build_hybrid_test_system(; with_reserves = false) + template = _build_hybrid_template(sys; with_reserves = false) + m = _build_and_solve(template, sys) + # When no service model is attached, hybrid reserve variables should not exist. + @test !any( + k -> IOM.get_entry_type(k) === POM.HybridPCCReserveVariable{DischargeSide}, + _var_keys(m), + ) + @test !any( + k -> IOM.get_entry_type(k) === POM.HybridPCCReserveVariable{ChargeSide}, + _var_keys(m), + ) +end + +@testset "HybridDispatchWithReserves: hybrid with no subcomponents (build only)" begin + sys = PSB.build_system(PSB.PSITestSystems, "test_RTS_GMLC_sys") + bus = PSY.get_component(PSY.ACBus, sys, "Chuhsi") + # A bare hybrid envelope: no thermal, renewable, storage, or load attached. + hybrid = PSY.HybridSystem(; + name = string(PSY.get_number(bus)) * "_BareHybrid", + available = true, + status = true, + bus = bus, + active_power = 0.0, + reactive_power = 0.0, + base_power = 100.0, + operation_cost = PSY.MarketBidCost(nothing), + thermal_unit = nothing, + electric_load = nothing, + storage = nothing, + renewable_unit = nothing, + interconnection_impedance = 0.0 + 0.0im, + interconnection_rating = nothing, + input_active_power_limits = (min = 0.0, max = 1.0), + output_active_power_limits = (min = 0.0, max = 1.0), + reactive_power_limits = nothing, + ) + PSY.add_component!(sys, hybrid) + template = _build_hybrid_template(sys; with_reserves = false) + m = POM.DecisionModel(template, sys; + optimizer = HiGHS_optimizer, horizon = _HYBRID_HORIZON) + @test POM.build!(m; output_dir = mktempdir(; cleanup = true)) == + IOM.ModelBuildStatus.BUILT + # Subcomponent variables must be absent for a bare envelope. + @test !any(k -> IOM.get_entry_type(k) === POM.HybridThermalActivePower, _var_keys(m)) + @test !any( + k -> IOM.get_entry_type(k) === POM.HybridStorageSubcomponentPower{ChargeSide}, + _var_keys(m), + ) +end + +# --------------------------------------------------------------------------- +# 9b. Comparison tests — verify outputs respond sensibly to structural changes. +# These tests use directional inequalities (not exact magnitudes) so they are +# robust to test-system-specific cost calibration. +# --------------------------------------------------------------------------- + +@testset "Comparison: reserves vs. no reserves" begin + # Same hybrid structure; one model has reserves attached, the other doesn't. + sys_r, _ = _build_hybrid_test_system(; with_reserves = true) + sys_n, _ = _build_hybrid_test_system(; with_reserves = false) + m_r = _build_and_solve(_build_hybrid_template(sys_r), sys_r) + m_n = _build_and_solve(_build_hybrid_template(sys_n; with_reserves = false), sys_n) + + obj_r = _obj(m_r) + obj_n = _obj(m_n) + # Reserves are extra constraints (and may carry slack penalties); the system + # under reserves cannot solve cheaper than the unconstrained one *at the true + # optimum*. Both models are solved only to HiGHS's default MIP relative gap + # (~1e-4), so the reported objectives can flip by that much; compare with a + # gap-aware relative tolerance rather than an absolute one. + @test obj_r >= obj_n * (1 - 1e-3) + + # Some reserve provision must occur in the reserves case if any service has a + # positive requirement. Sum non-zero values across all hybrid reserve variables. + container_r = IOM.get_optimization_container(m_r) + total_reserve_provision = 0.0 + for (key, var_arr) in IOM.get_variables(container_r) + IOM.get_component_type(key) === PSY.HybridSystem || continue + if IOM.get_entry_type(key) in + ( + POM.HybridPCCReserveVariable{DischargeSide}, + POM.HybridPCCReserveVariable{ChargeSide}, + ) + total_reserve_provision += sum(JuMP.value, var_arr) + end + end + @test total_reserve_provision >= 0.0 # always non-negative by variable bounds +end + +@testset "Comparison: with-thermal vs. without-thermal" begin + # The standalone "318_CC_1" stays in the system whether or not the hybrid + # references it as a subcomponent, so the objective-ordering invariant doesn't + # cleanly hold here. We assert structural changes instead: variable presence, + # finite/positive objective. + sys_t, _ = _build_hybrid_test_system() + sys_nt, _ = _build_hybrid_test_system(; with_thermal = false) + + m_t = _build_and_solve(_build_hybrid_template(sys_t), sys_t) + m_nt = _build_and_solve(_build_hybrid_template(sys_nt), sys_nt) + + @test isfinite(_obj(m_t)) && _obj(m_t) > 0 + @test isfinite(_obj(m_nt)) && _obj(m_nt) > 0 + + @test any(k -> IOM.get_entry_type(k) === POM.HybridThermalActivePower, _var_keys(m_t)) + @test !any(k -> IOM.get_entry_type(k) === POM.HybridThermalActivePower, _var_keys(m_nt)) +end + +@testset "Comparison: with-storage vs. without-storage" begin + # Without storage there are no per-storage reserve variables to aggregate, so + # disable reserves on this comparison. We assert structural changes instead + # of objective ordering (the system already has alternate sources of + # arbitrage, so the directional invariant is system-dependent). + sys_s, _ = _build_hybrid_test_system(; with_reserves = false) + sys_ns, _ = _build_hybrid_test_system(; + with_reserves = false, with_storage = false) + + m_s = _build_and_solve(_build_hybrid_template(sys_s; with_reserves = false), sys_s) + m_ns = _build_and_solve(_build_hybrid_template(sys_ns; with_reserves = false), sys_ns) + + @test isfinite(_obj(m_s)) && _obj(m_s) > 0 + @test isfinite(_obj(m_ns)) && _obj(m_ns) > 0 + + @test any( + k -> IOM.get_entry_type(k) === POM.HybridStorageSubcomponentPower{ChargeSide}, + _var_keys(m_s), + ) + @test !any( + k -> IOM.get_entry_type(k) === POM.HybridStorageSubcomponentPower{ChargeSide}, + _var_keys(m_ns), + ) +end + +@testset "Storage-less hybrid with reserves builds (C4 regression)" begin + # Regression: TotalReserveOffering containers used to be created only for storage + # hybrids, but get_expression_type_for_reserve routes *every* hybrid's + # ActivePowerReserveVariable into TotalReserveOffering. A storage-less hybrid with + # reserves attached previously hit a missing-container error during service + # construction; it must now build and solve. + sys, _ = _build_hybrid_test_system(; with_reserves = true, with_storage = false) + template = _build_hybrid_template(sys; with_reserves = true) + m = _build_and_solve(template, sys) + @test isfinite(_obj(m)) && _obj(m) > 0 + + # No storage-subcomponent reserve variables exist for a storage-less hybrid... + @test !any( + k -> + IOM.get_entry_type(k) === + POM.HybridStorageSubcomponentReserveVariable{ChargeSide}, + _var_keys(m), + ) + # ...but the PCC reserve variables that feed TotalReserveOffering are still present. + @test any( + k -> IOM.get_entry_type(k) === POM.HybridPCCReserveVariable{DischargeSide}, + _var_keys(m), + ) +end + +@testset "Comparison: energy_target on vs. off" begin + sys_on, _ = _build_hybrid_test_system(; energy_target = true) + sys_off, _ = _build_hybrid_test_system(; energy_target = true) + + m_on = _build_and_solve( + _build_hybrid_template(sys_on; + attributes = Dict{String, Any}("energy_target" => true)), + sys_on, + ) + m_off = _build_and_solve( + _build_hybrid_template(sys_off; + attributes = Dict{String, Any}("energy_target" => false)), + sys_off, + ) + + # Enabling the energy target adds the soft-equality constraint and penalizes the + # surplus/shortage slacks, which can only weakly raise the true objective. Confirms + # the penalty is wired into the objective (the original port had no penalty at all). + obj_on, obj_off = _obj(m_on), _obj(m_off) + @test isfinite(obj_on) && obj_on > 0 + @test isfinite(obj_off) && obj_off > 0 + @test obj_on >= obj_off * (1 - 1e-3) +end + +@testset "Comparison: regularization on vs. off" begin + sys_on, _ = _build_hybrid_test_system() + sys_off, _ = _build_hybrid_test_system() + + m_on = _build_and_solve( + _build_hybrid_template(sys_on; + attributes = Dict{String, Any}("regularization" => true)), + sys_on, + ) + m_off = _build_and_solve( + _build_hybrid_template(sys_off; + attributes = Dict{String, Any}("regularization" => false)), + sys_off, + ) + + # Adding the regularization penalty can only weakly raise the optimal *true* + # objective, but HiGHS's default MIP relative gap (~1e-4) lets the solver + # accept any feasible solution within that tolerance, so the directional + # inequality can flip on either side. Use a generous relative tolerance to + # account for that — we just want to confirm the two objectives sit in the + # same ballpark (both finite, both positive, within MIP gap of each other). + obj_on, obj_off = _obj(m_on), _obj(m_off) + @test isfinite(obj_on) && obj_on > 0 + @test isfinite(obj_off) && obj_off > 0 + @test isapprox(obj_on, obj_off; rtol = 1e-3) + + # The slack variables exist only in the on case. + @test any( + k -> IOM.get_entry_type(k) === POM.RegularizationVariable{ChargeSide}, + _var_keys(m_on), + ) + @test !any( + k -> IOM.get_entry_type(k) === POM.RegularizationVariable{ChargeSide}, + _var_keys(m_off), + ) +end diff --git a/test/test_utils/hybrid_test_utils.jl b/test/test_utils/hybrid_test_utils.jl new file mode 100644 index 0000000..89bde36 --- /dev/null +++ b/test/test_utils/hybrid_test_utils.jl @@ -0,0 +1,137 @@ +# Test fixtures for HybridSystem device tests. +# Ports HybridSystemsSimulations.jl/test/test_utils/function_utils.jl +# (modify_ren_curtailment_cost! and add_hybrid_to_chuhsi_bus!) to PSY 5.3. + +""" +Set a flat curtailment penalty on every RenewableDispatch in `sys`. PSY 5.x +RenewableGenerationCost wraps a CostCurve(LinearCurve(...)). +""" +function modify_ren_curtailment_cost!(sys::PSY.System; cost = 15.0) + for ren in PSY.get_components(PSY.RenewableDispatch, sys) + PSY.set_operation_cost!( + ren, + PSY.RenewableGenerationCost(PSY.CostCurve(PSY.LinearCurve(cost))), + ) + end + return +end + +""" +Build an EnergyReservoirStorage device sized for hybrid testing. PSY 5.x +constructor; default StorageCost(nothing) is fine for our test (no storage costs +contribute to the objective). Pass `storage_target` (a ratio of capacity) and a +`storage_cost` with non-zero surplus/shortage costs to exercise the energy target. +""" +function _build_hybrid_storage( + bus::PSY.ACBus, + energy_capacity, + rating, + eff_in, + eff_out; + storage_target = 0.0, + storage_cost = PSY.StorageCost(nothing), +) + name = string(PSY.get_number(bus)) * "_BATTERY" + return PSY.EnergyReservoirStorage(; + name = name, + available = true, + bus = bus, + prime_mover_type = PSY.PrimeMovers.BA, + storage_technology_type = PSY.StorageTech.OTHER_CHEM, + storage_capacity = energy_capacity, + storage_level_limits = (min = 0.05, max = 1.0), + initial_storage_capacity_level = 0.5, + rating = rating, + active_power = 0.0, + input_active_power_limits = (min = 0.0, max = rating), + output_active_power_limits = (min = 0.0, max = rating), + efficiency = (in = eff_in, out = eff_out), + reactive_power = 0.0, + reactive_power_limits = nothing, + base_power = 100.0, + operation_cost = storage_cost, + storage_target = storage_target, + ) +end + +""" +Add a HybridSystem to bus "Chuhsi" of an RTS-GMLC system, composed of: + - thermal subcomponent: existing "318_CC_1" generator + - renewable subcomponent: existing "317_WIND_1" generator + - electric load subcomponent: existing "Clark" load + - storage subcomponent: a fresh EnergyReservoirStorage built on bus Chuhsi + +Mirrors HSS test_utils/function_utils.jl:add_hybrid_to_chuhsi_bus!. +""" +function add_hybrid_to_chuhsi_bus!( + sys::PSY.System; + with_thermal::Bool = true, + with_renewable::Bool = true, + with_storage::Bool = true, + with_load::Bool = true, + energy_target::Bool = false, +) + bus = PSY.get_component(PSY.ACBus, sys, "Chuhsi") + bus === nothing && error("add_hybrid_to_chuhsi_bus!: bus 'Chuhsi' not found in system") + # With energy_target on, set an end-of-period target (ratio of capacity) and + # non-zero surplus/shortage penalties on the storage subcomponent so the slacks + # carry a real cost in the objective. + storage_target = energy_target ? 0.8 : 0.0 + storage_cost = + if energy_target + PSY.StorageCost(; energy_shortage_cost = 1000.0, energy_surplus_cost = 1000.0) + else + PSY.StorageCost(nothing) + end + bat = + if with_storage + _build_hybrid_storage( + bus, 4.0, 2.0, 0.93, 0.93; + storage_target = storage_target, + storage_cost = storage_cost, + ) + else + nothing + end + + # Subcomponents borrowed from adjacent existing components in RTS-GMLC. + renewable = + with_renewable ? + PSY.get_component(PSY.StaticInjection, sys, "317_WIND_1") : nothing + thermal = + with_thermal ? + PSY.get_component(PSY.StaticInjection, sys, "318_CC_1") : nothing + load = with_load ? PSY.get_component(PSY.PowerLoad, sys, "Clark") : nothing + for (flag, name, cmp) in ( + (with_renewable, "317_WIND_1", renewable), + (with_thermal, "318_CC_1", thermal), + (with_load, "Clark", load), + ) + if flag && cmp === nothing + error("add_hybrid_to_chuhsi_bus!: component '$name' not found") + end + end + + hybrid_name = string(PSY.get_number(bus)) * "_Hybrid" + hybrid = PSY.HybridSystem(; + name = hybrid_name, + available = true, + status = true, + bus = bus, + active_power = 1.0, + reactive_power = 0.0, + base_power = 100.0, + operation_cost = PSY.MarketBidCost(nothing), + thermal_unit = thermal, + electric_load = load, + storage = bat, + renewable_unit = renewable, + interconnection_impedance = 0.0 + 0.0im, + interconnection_rating = nothing, + input_active_power_limits = (min = 0.0, max = 10.0), + output_active_power_limits = (min = 0.0, max = 10.0), + reactive_power_limits = nothing, + ) + PSY.add_component!(sys, hybrid) + return hybrid +end