diff --git a/Project.toml b/Project.toml index b5ae2462fb..065acbf52e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PowerSimulations" uuid = "e690365d-45e2-57bb-ac84-44ba829e73c4" authors = ["Jose Daniel Lara", "Clayton Barrows", "Daniel Thom", "Dheepak Krishnamurthy", "Sourabh Dalvi"] -version = "0.27.7" +version = "0.27.8" [deps] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" diff --git a/docs/make.jl b/docs/make.jl index dda574be59..ce9640a385 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -15,9 +15,9 @@ pages = OrderedDict( "modeler_guide/psi_structure.md", "modeler_guide/problem_templates.md", "modeler_guide/running_a_simulation.md", + "modeler_guide/read_results.md", "modeler_guide/simulation_recorder.md", "modeler_guide/logging.md", - "modeler_guide/tips_and_tricks.md", "modeler_guide/debugging_infeasible_models.md", "modeler_guide/parallel_simulations.md", "modeler_guide/modeling_faq.md", @@ -31,12 +31,16 @@ pages = OrderedDict( "Troubleshooting" => "code_base_developer_guide/troubleshooting.md", ], "Formulation Library" => Any[ + "Introduction" => "formulation_library/Introduction.md", "General" => "formulation_library/General.md", + "Network" => "formulation_library/Network.md", "Thermal Generation" => "formulation_library/ThermalGen.md", "Renewable Generation" => "formulation_library/RenewableGen.md", "Load" => "formulation_library/Load.md", - "Network" => "formulation_library/Network.md", "Branch" => "formulation_library/Branch.md", + "Services" => "formulation_library/Service.md", + "Feedforwards" => "formulation_library/Feedforward.md", + "Piecewise Linear Cost" => "formulation_library/Piecewise.md", ], "API Reference" => "api/PowerSimulations.md", ) diff --git a/docs/src/api/PowerSimulations.md b/docs/src/api/PowerSimulations.md index c038a7d6c1..18f6bbb941 100644 --- a/docs/src/api/PowerSimulations.md +++ b/docs/src/api/PowerSimulations.md @@ -9,14 +9,38 @@ end ### Table of Contents -1. [Device Models](#device-models) -2. [Decision Models](#decision-models) -3. [Emulation Models](#emulation-models) -4. [Service Models](#service-models) -5. [Simulation Models](#simulation-models) -6. [Variables](#variables) -7. [Constraints](#constraints) -8. [Parameters](#parameters) +* [Device Models](#Device-Models) + * [Formulations](#Formulations) + * [Problem Templates](#Problem-Templates) +* [Decision Models](#Decision-Models) +* [Emulation Models](#Emulation-Models) +* [Service Models](#Service-Models) +* [Simulation Models](#Simulation-Models) +* [Variables](#Variables) + * [Common Variables](#Common-Variables) + * [Thermal Unit Variables](#Thermal-Unit-Variables) + * [Storage Unit Variables](#Storage-Unit-Variables) + * [Branches and Network Variables](#Branches-and-Network-Variables) + * [Services Variables](#Services-Variables) + * [Feedforward Variables](#Feedforward-Variables) +* [Constraints](#Constraints) + * [Common Constraints](#Common-Constraints) + * [Network Constraints](#Network-Constraints) + * [Power Variable Limit Constraints](#Power-Variable-Limit-Constraints) + * [Services Constraints](#Services-Constraints) + * [Thermal Unit Constraints](#Thermal-Unit-Constraints) + * [Renewable Unit Constraints](#Renewable-Unit-Constraints) + * [Branches Constraints](#Branches-Constraints) + * [Feedforward Constraints](#Feedforward-Constraints) +* [Parameters](#Parameters) + * [Time Series Parameters](#Time-Series-Parameters) + * [Variable Value Parameters](#Variable-Value-Parameters) + * [Objective Function Parameters](#Objective-Function-Parameters) + +```@raw html +  +  +``` # Device Models @@ -34,22 +58,13 @@ Refer to the [Formulations Page](@ref formulation_library) for each Abstract Dev Refer to the [Problem Templates Page](@ref op_problem_template) for available `ProblemTemplate`s. -### Problem Templates - -Refer to the [Problem Templates Page](https://nrel-siip.github.io/PowerSimulations.jl/latest/modeler_guide/problem_templates/) for available `ProblemTemplate`s. ```@raw html     ``` -# Service Models - -List of structures and methods for Service models - -```@docs -ServiceModel -``` +--- # Decision Models @@ -66,6 +81,8 @@ solve!(::DecisionModel)   ``` +--- + # Emulation Models ```@docs @@ -81,6 +98,23 @@ run!(::EmulationModel)   ``` +--- + +# Service Models + +List of structures and methods for Service models + +```@docs +ServiceModel +``` + +```@raw html +  +  +``` + +--- + # Simulation Models Refer to the [Simulations Page](@ref running_a_simulation) to explanations on how to setup a Simulation, with Sequencing and Feedforwards. @@ -99,6 +133,8 @@ execute!(::Simulation)   ``` +--- + # Variables For a list of variables for each device refer to its Formulations page. @@ -122,6 +158,7 @@ HotStartVariable WarmStartVariable ColdStartVariable PowerAboveMinimumVariable +PowerOutput ``` ### Storage Unit Variables @@ -134,6 +171,8 @@ ReservationVariable ```@docs FlowActivePowerVariable +FlowActivePowerSlackUpperBound +FlowActivePowerSlackLowerBound FlowActivePowerFromToVariable FlowActivePowerToFromVariable FlowReactivePowerFromToVariable @@ -145,21 +184,23 @@ VoltageMagnitude VoltageAngle ``` -### Regulation and Services Variables +### Services Variables ```@docs ActivePowerReserveVariable ServiceRequirementVariable -DeltaActivePowerUpVariable -DeltaActivePowerDownVariable -AdditionalDeltaActivePowerUpVariable -AdditionalDeltaActivePowerDownVariable -AreaMismatchVariable -SteadyStateFrequencyDeviation -SmoothACE SystemBalanceSlackUp SystemBalanceSlackDown ReserveRequirementSlack +InterfaceFlowSlackUp +InterfaceFlowSlackDown +``` + +### Feedforward Variables + +```@docs +UpperBoundFeedForwardSlack +LowerBoundFeedForwardSlack ``` ```@raw html @@ -167,6 +208,8 @@ ReserveRequirementSlack   ``` +--- + # Constraints ### Common Constraints @@ -180,10 +223,7 @@ PieceWiseLinearCostConstraint ```@docs AreaDispatchBalanceConstraint -AreaParticipationAssignmentConstraint -BalanceAuxConstraint CopperPlateBalanceConstraint -FrequencyResponseConstraint NodalBalanceActiveConstraint NodalBalanceReactiveConstraint ``` @@ -198,13 +238,11 @@ InputActivePowerVariableLimitsConstraint OutputActivePowerVariableLimitsConstraint ``` -### Regulation and Services Constraints +### Services Constraints ```@docs -ParticipationAssignmentConstraint -RegulationLimitsConstraint RequirementConstraint -ReserveEnergyCoverageConstraint +ParticipationFractionConstraint ReservePowerConstraint ``` @@ -215,7 +253,6 @@ ActiveRangeICConstraint CommitmentConstraint DurationConstraint RampConstraint -RampLimitConstraint StartupInitialConditionConstraint StartupTimeLimitTemperatureConstraint ``` @@ -224,41 +261,40 @@ StartupTimeLimitTemperatureConstraint ```@docs EqualityConstraint - ``` ### Branches Constraints ```@docs -AbsoluteValueConstraint -FlowLimitFromToConstraint -FlowLimitToFromConstraint +FlowLimitConstraint FlowRateConstraint FlowRateConstraintFromTo FlowRateConstraintToFrom -HVDCDirection HVDCLossesAbsoluteValue HVDCPowerBalance NetworkFlowConstraint RateLimitConstraint -RateLimitConstraintFromTo -RateLimitConstraintToFrom PhaseAngleControlLimit ``` ### Feedforward Constraints ```@docs -FeedforwardSemiContinousConstraint -FeedforwardIntegralLimitConstraint +FeedforwardSemiContinuousConstraint FeedforwardUpperBoundConstraint FeedforwardLowerBoundConstraint -FeedforwardEnergyTargetConstraint ``` +```@raw html +  +  +``` + +--- + # Parameters -## Time Series Parameters +### Time Series Parameters ```@docs ActivePowerTimeSeriesParameter @@ -266,15 +302,13 @@ ReactivePowerTimeSeriesParameter RequirementTimeSeriesParameter ``` -## Variable Value Parameters +### Variable Value Parameters ```@docs UpperBoundValueParameter LowerBoundValueParameter OnStatusParameter -EnergyLimitParameter FixValueParameter -EnergyTargetParameter ``` ### Objective Function Parameters diff --git a/docs/src/formulation_library/Branch.md b/docs/src/formulation_library/Branch.md index f55d6c37e7..2a70cae817 100644 --- a/docs/src/formulation_library/Branch.md +++ b/docs/src/formulation_library/Branch.md @@ -1,67 +1,356 @@ # `PowerSystems.Branch` Formulations +!!! note + The use of reactive power variables and constraints will depend on the network model used, i.e., whether it uses (or does not use) reactive power. If the network model is purely active power-based, reactive power variables and related constraints are not created. -Valid `DeviceModel`s for subtypes of `Branch` include the following: - -```@eval -using PowerSimulations -using PowerSystems -using DataFrames -using Latexify -combos = PowerSimulations.generate_device_formulation_combinations() -filter!(x -> x["device_type"] <: Branch, combos) -combo_table = DataFrame( - "Valid DeviceModel" => ["`DeviceModel($(c["device_type"]), $(c["formulation"]))`" for c in combos], - "Device Type" => ["[$(c["device_type"])](https://nrel-Sienna.github.io/PowerSystems.jl/stable/model_library/generated_$(c["device_type"])/)" for c in combos], - "Formulation" => ["[$(c["formulation"])](@ref)" for c in combos], - ) -mdtable(combo_table, latex = false) -``` +### Table of contents ---- +1. [`StaticBranch`](#StaticBranch) +2. [`StaticBranchBounds`](#StaticBranchBounds) +3. [`StaticBranchUnbounded`](#StaticBranchUnbounded) +4. [`HVDCTwoTerminalUnbounded`](#HVDCTwoTerminalUnbounded) +5. [`HVDCTwoTerminalLossless`](#HVDCTwoTerminalLossless) +6. [`HVDCTwoTerminalDispatch`](#HVDCTwoTerminalDispatch) +7. [`PhaseAngleControl`](#PhaseAngleControl) +8. [Valid configurations](#Valid-configurations) ## `StaticBranch` +Formulation valid for `PTDFPowerModel` Network model + ```@docs StaticBranch ``` +**Variables:** + +- [`FlowActivePowerVariable`](@ref): + - Bounds: ``(-\infty,\infty)`` + - Symbol: ``f`` +If Slack variables are enabled: +- [`FlowActivePowerSlackUpperBound`](@ref): + - Bounds: [0.0, ] + - Default proportional cost: 2e5 + - Symbol: ``f^\text{sl,up}`` +- [`FlowActivePowerSlackLowerBound`](@ref): + - Bounds: [0.0, ] + - Default proportional cost: 2e5 + - Symbol: ``f^\text{sl,lo}`` + +**Static Parameters** + +- ``R^\text{max}`` = `PowerSystems.get_rate(branch)` + +**Objective:** + +Add a large proportional cost to the objective function if rate constraint slack variables are used ``+ (f^\text{sl,up} + f^\text{sl,lo}) \cdot 2 \cdot 10^5`` + +**Expressions:** + +No expressions are used. + +**Constraints:** + +For each branch ``b \in \{1,\dots, B\}`` (in a system with ``N`` buses) the constraints are given by: + +```math +\begin{aligned} +& f_t = \sum_{i=1}^N \text{PTDF}_{i,b} \cdot \text{Bal}_{i,t}, \quad \forall t \in \{1,\dots, T\}\\ +& f_t - f_t^\text{sl,up} \le R^\text{max},\quad \forall t \in \{1,\dots, T\} \\ +& f_t + f_t^\text{sl,lo} \ge -R^\text{max},\quad \forall t \in \{1,\dots, T\} +\end{aligned} +``` +on which ``\text{PTDF}`` is the ``N \times B`` system Power Transfer Distribution Factors (PTDF) matrix, and ``\text{Bal}_{i,t}`` is the active power bus balance expression (i.e. ``\text{Generation}_{i,t} - \text{Demand}_{i,t}``) at bus ``i`` at time-step ``t``. + --- ## `StaticBranchBounds` +Formulation valid for `PTDFPowerModel` Network model + ```@docs StaticBranchBounds ``` +**Variables:** + +- [`FlowActivePowerVariable`](@ref): + - Bounds: ``\left[-R^\text{max},R^\text{max}\right]`` + - Symbol: ``f`` + +**Static Parameters** + +- ``R^\text{max}`` = `PowerSystems.get_rate(branch)` + +**Objective:** + +No cost is added to the objective function. + +**Expressions:** + +No expressions are used. + +**Constraints:** + +For each branch ``b \in \{1,\dots, B\}`` (in a system with ``N`` buses) the constraints are given by: + +```math +\begin{aligned} +& f_t = \sum_{i=1}^N \text{PTDF}_{i,b} \cdot \text{Bal}_{i,t}, \quad \forall t \in \{1,\dots, T\} +\end{aligned} +``` +on which ``\text{PTDF}`` is the ``N \times B`` system Power Transfer Distribution Factors (PTDF) matrix, and ``\text{Bal}_{i,t}`` is the active power bus balance expression (i.e. ``\text{Generation}_{i,t} - \text{Demand}_{i,t}``) at bus ``i`` at time-step ``t``. + --- ## `StaticBranchUnbounded` +Formulation valid for `PTDFPowerModel` Network model + ```@docs StaticBranchUnbounded ``` +- [`FlowActivePowerVariable`](@ref): + - Bounds: ``(-\infty,\infty)`` + - Symbol: ``f`` + + +**Objective:** + +No cost is added to the objective function. + +**Expressions:** + +No expressions are used. + +**Constraints:** + +For each branch ``b \in \{1,\dots, B\}`` (in a system with ``N`` buses) the constraints are given by: + +```math +\begin{aligned} +& f_t = \sum_{i=1}^N \text{PTDF}_{i,b} \cdot \text{Bal}_{i,t}, \quad \forall t \in \{1,\dots, T\} +\end{aligned} +``` +on which ``\text{PTDF}`` is the ``N \times B`` system Power Transfer Distribution Factors (PTDF) matrix, and ``\text{Bal}_{i,t}`` is the active power bus balance expression (i.e. ``\text{Generation}_{i,t} - \text{Demand}_{i,t}``) at bus ``i`` at time-step ``t``. + +--- + +## `HVDCTwoTerminalUnbounded` + +Formulation valid for `PTDFPowerModel` Network model + +```@docs +HVDCTwoTerminalUnbounded +``` + +This model assumes that it can transfer power from two AC buses without losses and no limits. + +**Variables:** + +- [`FlowActivePowerVariable`](@ref): + - Bounds: ``\left(-\infty,\infty\right)`` + - Symbol: ``f`` + + +**Objective:** + +No cost is added to the objective function. + +**Expressions:** + +The variable `FlowActivePowerVariable` ``f`` is added to the nodal balance expression `ActivePowerBalance`, by adding the flow ``f`` in the receiving bus and subtracting it from the sending bus. This is used then to compute the AC flows using the PTDF equation. + +**Constraints:** + +No constraints are added. + + --- ## `HVDCTwoTerminalLossless` +Formulation valid for `PTDFPowerModel` Network model + ```@docs HVDCTwoTerminalLossless ``` +This model assumes that it can transfer power from two AC buses without losses. + +**Variables:** + +- [`FlowActivePowerVariable`](@ref): + - Bounds: ``\left(-\infty,\infty\right)`` + - Symbol: ``f`` + + +**Static Parameters** + +- ``R^\text{from,min}`` = `PowerSystems.get_active_power_limits_from(branch).min` +- ``R^\text{from,max}`` = `PowerSystems.get_active_power_limits_from(branch).max` +- ``R^\text{to,min}`` = `PowerSystems.get_active_power_limits_to(branch).min` +- ``R^\text{to,max}`` = `PowerSystems.get_active_power_limits_to(branch).max` + +**Objective:** + +No cost is added to the objective function. + +**Expressions:** + +The variable `FlowActivePowerVariable` ``f`` is added to the nodal balance expression `ActivePowerBalance`, by adding the flow ``f`` in the receiving bus and subtracting it from the sending bus. This is used then to compute the AC flows using the PTDF equation. + +**Constraints:** + +```math +\begin{align*} +& R^\text{min} \le f_t \le R^\text{max},\quad \forall t \in \{1,\dots, T\} \\ +\end{align*} +``` +where: +```math +\begin{align*} +& R^\text{min} = \begin{cases} + \min\left(R^\text{from,min}, R^\text{to,min}\right), & \text{if } R^\text{from,min} \ge 0 \text{ and } R^\text{to,min} \ge 0 \\ + \max\left(R^\text{from,min}, R^\text{to,min}\right), & \text{if } R^\text{from,min} \le 0 \text{ and } R^\text{to,min} \le 0 \\ + R^\text{from,min},& \text{if } R^\text{from,min} \le 0 \text{ and } R^\text{to,min} \ge 0 \\ + R^\text{to,min},& \text{if } R^\text{from,min} \ge 0 \text{ and } R^\text{to,min} \le 0 + \end{cases} +\end{align*} +``` +and +```math +\begin{align*} +& R^\text{max} = \begin{cases} + \min\left(R^\text{from,max}, R^\text{to,max}\right), & \text{if } R^\text{from,max} \ge 0 \text{ and } R^\text{to,max} \ge 0 \\ + \max\left(R^\text{from,max}, R^\text{to,max}\right), & \text{if } R^\text{from,max} \le 0 \text{ and } R^\text{to,max} \le 0 \\ + R^\text{from,max},& \text{if } R^\text{from,max} \le 0 \text{ and } R^\text{to,max} \ge 0 \\ + R^\text{to,max},& \text{if } R^\text{from,max} \ge 0 \text{ and } R^\text{to,max} \le 0 + \end{cases} +\end{align*} +``` + --- -## `HVDCTwoTerminalDispatch` + +## `HVDCTwoTerminalDispatch` + +Formulation valid for `PTDFPowerModel` Network model ```@docs HVDCTwoTerminalDispatch ``` +**Variables** + +- [`FlowActivePowerToFromVariable`](@ref): + - Symbol: ``f^\text{to-from}`` +- [`FlowActivePowerFromToVariable`](@ref): + - Symbol: ``f^\text{from-to}`` +- [`HVDCLosses`](@ref): + - Symbol: ``\ell`` +- [`HVDCFlowDirectionVariable`](@ref) + - Bounds: ``\{0,1\}`` + - Symbol: ``u^\text{dir}`` + +**Static Parameters** + +- ``R^\text{from,min}`` = `PowerSystems.get_active_power_limits_from(branch).min` +- ``R^\text{from,max}`` = `PowerSystems.get_active_power_limits_from(branch).max` +- ``R^\text{to,min}`` = `PowerSystems.get_active_power_limits_to(branch).min` +- ``R^\text{to,max}`` = `PowerSystems.get_active_power_limits_to(branch).max` +- ``L_0`` = `PowerSystems.get_loss(branch).l0` +- ``L_1`` = `PowerSystems.get_loss(branch).l1` + +**Objective:** + +No cost is added to the objective function. + +**Expressions:** + +Each `FlowActivePowerToFromVariable` ``f^\text{to-from}`` and `FlowActivePowerFromToVariable` ``f^\text{from-to}`` is added to the nodal balance expression `ActivePowerBalance`, by adding the respective flow in the receiving bus and subtracting it from the sending bus. That is, ``f^\text{to-from}`` adds the flow to the `from` bus, and subtracts the flow from the `to` bus, while ``f^\text{from-to}`` adds the flow to the `to` bus, and subtracts the flow from the `from` bus This is used then to compute the AC flows using the PTDF equation. + +In addition, the `HVDCLosses` are subtracted to the `from` bus in the `ActivePowerBalance` expression. + +**Constraints:** + +```math +\begin{align*} +& R^\text{from,min} \le f_t^\text{from-to} \le R^\text{from,max}, \forall t \in \{1,\dots, T\} \\ +& R^\text{to,min} \le f_t^\text{to-from} \le R^\text{to,max},\quad \forall t \in \{1,\dots, T\} \\ +& f_t^\text{to-from} - f_t^\text{from-to} \le L_1 \cdot f_t^\text{to-from} - L_0,\quad \forall t \in \{1,\dots, T\} \\ +& f_t^\text{from-to} - f_t^\text{to-from} \ge L_1 \cdot f_t^\text{from-to} + L_0,\quad \forall t \in \{1,\dots, T\} \\ +& f_t^\text{from-to} - f_t^\text{to-from} \ge - M^\text{big} (1 - u^\text{dir}_t),\quad \forall t \in \{1,\dots, T\} \\ +& f_t^\text{to-from} - f_t^\text{from-to} \ge - M^\text{big} u^\text{dir}_t,\quad \forall t \in \{1,\dots, T\} \\ +& f_t^\text{to-from} - f_t^\text{from-to} \le \ell_t,\quad \forall t \in \{1,\dots, T\} \\ +& f_t^\text{from-to} - f_t^\text{to-from} \le \ell_t,\quad \forall t \in \{1,\dots, T\} +\end{align*} +``` + --- -## `HVDCTwoTerminalUnbounded` +## `PhaseAngleControl` + +Formulation valid for `PTDFPowerModel` Network model ```@docs -HVDCTwoTerminalUnbounded +PhaseAngleControl ``` + +**Variables:** + +- [`FlowActivePowerVariable`](@ref): + - Bounds: ``(-\infty,\infty)`` + - Symbol: ``f`` +- [`PhaseShifterAngle`](@ref): + - Symbol: ``\theta^\text{shift}`` + +**Static Parameters** + +- ``R^\text{max}`` = `PowerSystems.get_rate(branch)` +- ``\Theta^\text{min}`` = `PowerSystems.get_phase_angle_limits(branch).min` +- ``\Theta^\text{max}`` = `PowerSystems.get_phase_angle_limits(branch).max` +- ``X`` = `PowerSystems.get_x(branch)` (series reactance) + +**Objective:** + +No changes to objective function + +**Expressions:** + +Adds to the `ActivePowerBalance` expression the term ``-\theta^\text{shift} /X`` to the `from` bus and ``+\theta^\text{shift} /X`` to the `to` bus, that the `PhaseShiftingTransformer` is connected. + +**Constraints:** + +For each branch ``b \in \{1,\dots, B\}`` (in a system with ``N`` buses) the constraints are given by: + +```math +\begin{aligned} +& f_t = \sum_{i=1}^N \text{PTDF}_{i,b} \cdot \text{Bal}_{i,t} + \frac{\theta^\text{shift}_t}{X}, \quad \forall t \in \{1,\dots, T\}\\ +& -R^\text{max} \le f_t \le R^\text{max},\quad \forall t \in \{1,\dots, T\} +\end{aligned} +``` +on which ``\text{PTDF}`` is the ``N \times B`` system Power Transfer Distribution Factors (PTDF) matrix, and ``\text{Bal}_{i,t}`` is the active power bus balance expression (i.e. ``\text{Generation}_{i,t} - \text{Demand}_{i,t}``) at bus ``i`` at time-step ``t``. + + +--- + +## Valid configurations + +Valid `DeviceModel`s for subtypes of `Branch` include the following: + +```@eval +using PowerSimulations +using PowerSystems +using DataFrames +using Latexify +combos = PowerSimulations.generate_device_formulation_combinations() +filter!(x -> (x["device_type"] <: Branch) && (x["device_type"] != TModelHVDCLine), combos) +combo_table = DataFrame( + "Valid DeviceModel" => ["`DeviceModel($(c["device_type"]), $(c["formulation"]))`" for c in combos], + "Device Type" => ["[$(c["device_type"])](https://nrel-Sienna.github.io/PowerSystems.jl/stable/model_library/generated_$(c["device_type"])/)" for c in combos], + "Formulation" => ["[$(c["formulation"])](@ref)" for c in combos], + ) +mdtable(combo_table, latex = false) +``` \ No newline at end of file diff --git a/docs/src/formulation_library/Feedforward.md b/docs/src/formulation_library/Feedforward.md new file mode 100644 index 0000000000..bdda721f36 --- /dev/null +++ b/docs/src/formulation_library/Feedforward.md @@ -0,0 +1,161 @@ +# [FeedForward Formulations](@id ff_formulations) + +**FeedForwards** are the mechanism to define how information is shared between models. Specifically, a FeedForward defines what to do with information passed with an inter-stage chronology in a Simulation. The most common FeedForward is the `SemiContinuousFeedForward` that affects the semi-continuous range constraints of thermal generators in the economic dispatch problems based on the value of the (already solved) unit-commitment variables. + +The creation of a FeedForward requires at least specifying the `component_type` on which the FeedForward will be applied. The `source` variable specifies which variable will be taken from the problem solved, for example, the commitment variable of the thermal unit in the unit commitment problem. Finally, the `affected_values` specify which variables will be affected in the problem to be solved, for example, the next economic dispatch problem. + +### Table of contents + +1. [`SemiContinuousFeedforward`](#SemiContinuousFeedForward) +2. [`FixValueFeedforward`](#FixValueFeedforward) +3. [`UpperBoundFeedforward`](#UpperBoundFeedforward) +4. [`LowerBoundFeedforward`](#LowerBoundFeedforward) + +--- + +## `SemiContinuousFeedforward` + +```@docs +SemiContinuousFeedforward +``` + +**Variables:** + +No variables are created + +**Parameters:** + +- ``\text{on}^\text{th}`` = `OnStatusParameter` obtained from the source variable, typically the commitment variable of the unit commitment problem ``u^\text{th}``. + +**Objective:** + +No changes to the objective function. + +**Expressions:** + +Adds ``-\text{on}^\text{th}P^\text{th,max}`` to the `ActivePowerRangeExpressionUB` expression and ``-\text{on}^\text{th}P^\text{th,min}`` to the `ActivePowerRangeExpressionLB` expression. + +**Constraints:** + +Limits the `ActivePowerRangeExpressionUB` and `ActivePowerRangeExpressionLB` by zero as: + +```math +\begin{align*} +& \text{ActivePowerRangeExpressionUB}_t := p_t^\text{th} - \text{on}_t^\text{th}P^\text{th,max} \le 0, \quad \forall t\in \{1, \dots, T\} \\ +& \text{ActivePowerRangeExpressionLB}_t := p_t^\text{th} - \text{on}_t^\text{th}P^\text{th,min} \ge 0, \quad \forall t\in \{1, \dots, T\} +\end{align*} +``` + +Thus, if the commitment parameter is zero, the dispatch is limited to zero, forcing to turn off the generator without introducing binary variables in the economic dispatch problem. + +--- + +## `FixValueFeedforward` + +```@docs +FixValueFeedforward +``` + +**Variables:** + +No variables are created + +**Parameters:** + +The parameter `FixValueParameter` is used to match the result obtained from the source variable (from the simulation state). + +**Objective:** + +No changes to the objective function. + +**Expressions:** + +No changes on expressions. + +**Constraints:** + +Set the `VariableType` from the `affected_values` to be equal to the source parameter store in `FixValueParameter` + +```math +\begin{align*} +& \text{AffectedVariable}_t = \text{SourceVariableParameter}_t, \quad \forall t \in \{1,\dots, T\} +\end{align*} +``` + +--- + +## `UpperBoundFeedforward` + +```@docs +UpperBoundFeedforward +``` + +**Variables:** + +If slack variables are enabled: +- [`UpperBoundFeedForwardSlack`](@ref) + - Bounds: [0.0, ] + - Default proportional cost: 1e6 + - Symbol: ``p^\text{ff,ubsl}`` + + +**Parameters:** + +The parameter `UpperBoundValueParameter` stores the result obtained from the source variable (from the simulation state) that will be used as an upper bound to the affected variable. + +**Objective:** + +The slack variable is added to the objective function using its large default cost ``+ p^\text{ff,ubsl} \cdot 10^6`` + +**Expressions:** + +No changes on expressions. + +**Constraints:** + +Set the `VariableType` from the `affected_values` to be lower than the source parameter store in `UpperBoundValueParameter`. + +```math +\begin{align*} +& \text{AffectedVariable}_t - p_t^\text{ff,ubsl} \le \text{SourceVariableParameter}_t, \quad \forall t \in \{1,\dots, T\} +\end{align*} +``` + +--- + +## `LowerBoundFeedforward` + +```@docs +LowerBoundFeedforward +``` + +**Variables:** + +If slack variables are enabled: +- [`LowerBoundFeedForwardSlack`](@ref) + - Bounds: [0.0, ] + - Default proportional cost: 1e6 + - Symbol: ``p^\text{ff,lbsl}`` + + +**Parameters:** + +The parameter `LowerBoundValueParameter` stores the result obtained from the source variable (from the simulation state) that will be used as a lower bound to the affected variable. + +**Objective:** + +The slack variable is added to the objective function using its large default cost ``+ p^\text{ff,lbsl} \cdot 10^6`` + +**Expressions:** + +No changes on expressions. + +**Constraints:** + +Set the `VariableType` from the `affected_values` to be greater than the source parameter store in `LowerBoundValueParameter`. + +```math +\begin{align*} +& \text{AffectedVariable}_t + p_t^\text{ff,lbsl} \ge \text{SourceVariableParameter}_t, \quad \forall t \in \{1,\dots, T\} +\end{align*} +``` \ No newline at end of file diff --git a/docs/src/formulation_library/General.md b/docs/src/formulation_library/General.md index edc2fb3f1b..ceb5e1b9da 100644 --- a/docs/src/formulation_library/General.md +++ b/docs/src/formulation_library/General.md @@ -15,11 +15,11 @@ No variables are created for `DeviceModel(<:DeviceType, FixedOutput)` **Static Parameters:** - ThermalGen: - - ``Pg^\text{max}`` = `PowerSystems.get_max_active_power(device)` - - ``Qg^\text{max}`` = `PowerSystems.get_max_reactive_power(device)` + - ``P^\text{th,max}`` = `PowerSystems.get_max_active_power(device)` + - ``Q^\text{th,max}`` = `PowerSystems.get_max_reactive_power(device)` - Storage: - - ``Pg^\text{max}`` = `PowerSystems.get_max_active_power(device)` - - ``Qg^\text{max}`` = `PowerSystems.get_max_reactive_power(device)` + - ``P^\text{st,max}`` = `PowerSystems.get_max_active_power(device)` + - ``Q^\text{st,max}`` = `PowerSystems.get_max_reactive_power(device)` **Time Series Parameters:** @@ -48,7 +48,7 @@ No objective terms are created for `DeviceModel(<:DeviceType, FixedOutput)` **Expressions:** -Adds the active and reactive parameters listed for specific device types above to the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations) +Adds the active and reactive parameters listed for specific device types above to the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations). **Constraints:** @@ -104,9 +104,9 @@ where - For `variable_cost::PiecewiseLinearData`, ``f(x)`` is the piecewise linear function obtained by connecting the `(x, y)` points `get_points(variable_cost)` in order. - For `variable_cost = PiecewiseLinearSlopeData([x0, x1, x2, ...], y0, [s0, s1, s2, ...])`, ``f(x)`` is the piecewise linear function obtained by starting at `(x0, y0)`, drawing a segment at slope `s0` to `x=x1`, drawing a segment at slope `s1` to `x=x2`, etc. -___ +--- -### `StorageManagementCost` +## `StorageCost` Adds an objective function cost term according to: @@ -118,7 +118,7 @@ Adds an objective function cost term according to: **Impact of different cost configurations:** -The following table describes all possible configuration of the `StorageManagementCost` with the target constraint in hydro or storage device models. Cases 1(a) & 2(a) will have no impact of the models operations and the target constraint will be rendered useless. In most cases that have no energy target and a non-zero value for ``C^{value}``, if this cost is too high (``C^{value} >> 0``) or too low (``C^{value} <<0``) can result in either the model holding on to stored energy till the end or the model not storing any energy in the device. This is caused by the fact that when energy target is zero, we have ``E_t = - E^{shortage}_t``, and ``- E^{shortage}_t * C^{value}`` in the objective function is replaced by ``E_t * C^{value}``, thus resulting in ``C^{value}`` to be seen as the cost of stored energy. +The following table describes all possible configurations of the `StorageCost` with the target constraint in hydro or storage device models. Cases 1(a) & 2(a) will not impact the model's operations, and the target constraint will be rendered useless. In most cases that have no energy target and a non-zero value for ``C^{value}``, if this cost is too high (``C^{value} >> 0``) or too low (``C^{value} <<0``) can result in either the model holding on to stored energy till the end of the model not storing any energy in the device. This is caused by the fact that when the energy target is zero, we have ``E_t = - E^{shortage}_t``, and ``- E^{shortage}_t * C^{value}`` in the objective function is replaced by ``E_t * C^{value}``, thus resulting in ``C^{value}`` to be seen as the cost of stored energy. | Case | Energy Target | Energy Shortage Cost | Energy Value / Energy Surplus cost | Effect | diff --git a/docs/src/formulation_library/Introduction.md b/docs/src/formulation_library/Introduction.md new file mode 100644 index 0000000000..47a8425d4e --- /dev/null +++ b/docs/src/formulation_library/Introduction.md @@ -0,0 +1,67 @@ +# [Formulations Introduction](@id formulation_intro) + +PowerSimulations.jl enables modularity in its formulations by assigning a `DeviceModel` to each `PowerSystems.jl` component type existing in a defined system. + +`PowerSimulations.jl` has a multiple `AbstractDeviceFormulation` subtypes that can be applied to different `PowerSystems.jl` device types, each dispatching to different methods for populating the optimization problem **variables**, **objective function**, **expressions** and **constraints**. + +## Example Formulation + +For example a typical optimization problem in a `DecisionModel` in `PowerSimulations.jl` with three `DeviceModel` has the abstract form of: + +```math +\begin{align*} + &\min_{\boldsymbol{x}}~ \text{Objective\_DeviceModelA} + \text{Objective\_DeviceModelB} + \text{Objective\_DeviceModelC} \\ + & ~~\text{s.t.} \\ + & \hspace{0.9cm} \text{Constraints\_NetworkModel} \\ + & \hspace{0.9cm} \text{Constraints\_DeviceModelA} \\ + & \hspace{0.9cm} \text{Constraints\_DeviceModelB} \\ + & \hspace{0.9cm} \text{Constraints\_DeviceModelC} +\end{align*} +``` + +Suppose this is a system with the following characteristics: +- Horizon: 48 hours +- Interval: 24 hours +- Resolution: 1 hour +- Three Buses: 1, 2 and 3 +- One `ThermalStandard` (device A) unit at bus 1 +- One `RenewableDispatch` (device B) unit at bus 2 +- One `PowerLoad` (device C) at bus 3 +- Three `Line` that connects all the buses + +Now, we assign the following `DeviceModel` to each `PowerSystems.jl` with: + +| Type | Formulation | +| ----------- | ----------- | +| Network | `CopperPlatePowerModel` | +| `ThermalStandard` | `ThermalDispatchNoMin` | +| `RenewableDispatch` | `RenewableFullDispatch` | +| `PowerLoad` | `StaticPowerLoad` | + +Note that we did not assign any `DeviceModel` to `Line` since the `CopperPlatePowerModel` used for the network assumes that everything is lumped in the same node (like a copper plate with infinite capacity), and hence there are no flows between buses that branches can limit. + +Each `DeviceModel` formulation is described in specific in their respective page, but the overall optimization problem will end-up as: + +```math +\begin{align*} + &\min_{\boldsymbol{p}^\text{th}, \boldsymbol{p}^\text{re}}~ \sum_{t=1}^{48} C^\text{th} p_t^\text{th} - C^\text{re} p_t^\text{re} \\ + & ~~\text{s.t.} \\ + & \hspace{0.9cm} p_t^\text{th} + p_t^\text{re} = P_t^\text{load}, \quad \forall t \in {1,\dots, 48} \\ + & \hspace{0.9cm} 0 \le p_t^\text{th} \le P^\text{th,max} \\ + & \hspace{0.9cm} 0 \le p_t^\text{re} \le \text{ActivePowerTimeSeriesParameter}_t +\end{align*} +``` + +Note that the `StaticPowerLoad` does not impose any cost to the objective function or constraint but adds its power demand to the supply-balance demand of the `CopperPlatePowerModel` used. Since we are using the `ThermalDispatchNoMin` formulation for the thermal generation, the lower bound for the power is 0, instead of ``P^\text{th,min}``. In addition, we are assuming a linear cost ``C^\text{th}``. Finally, the `RenewableFullDispatch` formulation allows the dispatch of the renewable unit between 0 and its maximum injection time series ``p_t^\text{re,param}``. + +# Nomenclature + +In the formulations described in the other pages, the nomenclature is as follows: +- Lowercase letters are used for variables, e.g., ``p`` for power. +- Uppercase letters are used for parameters, e.g., ``C`` for costs. +- Subscripts are used for indexing, e.g., ``(\cdot)_t`` for indexing at time ``t``. +- Superscripts are used for descriptions, e.g., ``(\cdot)^\text{th}`` to describe a thermal (th) variable/parameter. +- Bold letters are used for vectors, e.g., ``\boldsymbol{p} = \{p\}_{1,\dots,24}``. + + + diff --git a/docs/src/formulation_library/Load.md b/docs/src/formulation_library/Load.md index c3bcbabb3c..dcb7b9b8d0 100644 --- a/docs/src/formulation_library/Load.md +++ b/docs/src/formulation_library/Load.md @@ -1,21 +1,16 @@ # `PowerSystems.ElectricLoad` Formulations -Valid `DeviceModel`s for subtypes of `ElectricLoad` include the following: +Electric load formulations define the optimization models that describe load units (demand) mathematical model in different operational settings, such as economic dispatch and unit commitment. -```@eval -using PowerSimulations -using PowerSystems -using DataFrames -using Latexify -combos = PowerSimulations.generate_device_formulation_combinations() -filter!(x -> x["device_type"] <: ElectricLoad, combos) -combo_table = DataFrame( - "Valid DeviceModel" => ["`DeviceModel($(c["device_type"]), $(c["formulation"]))`" for c in combos], - "Device Type" => ["[$(c["device_type"])](https://nrel-Sienna.github.io/PowerSystems.jl/stable/model_library/generated_$(c["device_type"])/)" for c in combos], - "Formulation" => ["[$(c["formulation"])](@ref)" for c in combos], - ) -mdtable(combo_table, latex = false) -``` +!!! note + The use of reactive power variables and constraints will depend on the network model used, i.e., whether it uses (or does not use) reactive power. If the network model is purely active power-based, reactive power variables and related constraints are not created. + +### Table of contents + +1. [`StaticPowerLoad`](#StaticPowerLoad) +2. [`PowerLoadInterruption`](#PowerLoadInterruption) +3. [`PowerLoadDispatch`](#PowerLoadDispatch) +4. [Valid configurations](#Valid-configurations) --- @@ -31,6 +26,8 @@ No variables are created **Time Series Parameters:** +Uses the `max_active_power` timeseries parameter to determine the demand value at each time-step + ```@eval using PowerSimulations using PowerSystems @@ -46,7 +43,7 @@ mdtable(combo_table, latex = false) **Expressions:** -Subtracts the parameters listed above from the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations) +Subtracts the parameters listed above from the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations). **Constraints:** @@ -65,12 +62,19 @@ PowerLoadInterruption - [`ActivePowerVariable`](@ref): - Bounds: [0.0, ] - Default initial value: 0.0 + - Symbol: ``p^\text{ld}`` - [`ReactivePowerVariable`](@ref): - Bounds: [0.0, ] - Default initial value: 0.0 + - Symbol: ``q^\text{ld}`` - [`OnVariable`](@ref): - - Bounds: {0,1} + - Bounds: ``\{0,1\}`` - Default initial value: 1 + - Symbol: ``u^\text{ld}`` + +**Static Parameters:** +- ``P^\text{ld,max}`` = `PowerSystems.get_max_active_power(device)` +- ``Q^\text{ld,max}`` = `PowerSystems.get_max_reactive_power(device)` **Time Series Parameters:** @@ -89,25 +93,22 @@ mdtable(combo_table, latex = false) **Objective:** -Creates an objective function term based on the [`FunctionData` Options](@ref) where the quantity term is defined as ``Pg``. +Creates an objective function term based on the [`FunctionData` Options](@ref) where the quantity term is defined as ``p^\text{ld}``. + **Expressions:** -- Adds ``Pg`` and ``Qg`` terms and to the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations) -- Subtracts the time series parameters listed above terms from the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations) +- Subtract``p^\text{ld}`` and ``q^\text{ld}`` terms and to the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations) **Constraints:** -``Pg`` and ``Qg`` represent the "unserved" active and reactive power loads - ```math \begin{aligned} -& Pg_t \le ActivePowerTimeSeriesParameter_t\\ -& Pg_t - u_t ActivePowerTimeSeriesParameter_t \le 0 \\ -& Qg_t \le ReactivePowerTimeSeriesParameter_t\\ -& Qg_t - u_t ReactivePowerTimeSeriesParameter_t\le 0 +& p_t^\text{ld} \le u_t^\text{ld} \cdot \text{ActivePowerTimeSeriesParameter}_t, \quad \forall t \in \{1,\dots, T\} \\ +& q_t^\text{re} = \text{pf} \cdot p_t^\text{re}, \quad \forall t \in \{1,\dots, T\} \end{aligned} ``` +on which ``\text{pf} = \sin(\arctan(Q^\text{ld,max}/P^\text{ld,max}))``. --- @@ -122,9 +123,15 @@ PowerLoadDispatch - [`ActivePowerVariable`](@ref): - Bounds: [0.0, ] - Default initial value: `PowerSystems.get_active_power(device)` + - Symbol: ``p^\text{ld}`` - [`ReactivePowerVariable`](@ref): - Bounds: [0.0, ] - Default initial value: `PowerSystems.get_reactive_power(device)` + - Symbol: ``q^\text{ld}`` + +**Static Parameters:** +- ``P^\text{ld,max}`` = `PowerSystems.get_max_active_power(device)` +- ``Q^\text{ld,max}`` = `PowerSystems.get_max_reactive_power(device)` **Time Series Parameters:** @@ -143,20 +150,38 @@ mdtable(combo_table, latex = false) **Objective:** -Creates an objective function term based on the [`FunctionData` Options](@ref) where the quantity term is defined as ``Pg``. +Creates an objective function term based on the [`FunctionData` Options](@ref) where the quantity term is defined as ``p^\text{ld}``. + **Expressions:** -- Adds ``Pg`` and ``Qg`` terms and to the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations) -- Subtracts the time series parameters listed above terms from the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations) +- Subtract``p^\text{ld}`` and ``q^\text{ld}`` terms and to the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations) **Constraints:** -``Pg`` and ``Qg`` represent the "unserved" active and reactive power loads - ```math \begin{aligned} -& Pg_t \le ActivePowerTimeSeriesParameter_t\\ -& Qg_t \le ReactivePowerTimeSeriesParameter_t\\ +& p_t^\text{ld} \le \text{ActivePowerTimeSeriesParameter}_t, \quad \forall t \in \{1,\dots, T\}\\ +& q_t^\text{ld} = \text{pf} \cdot p_t^\text{ld}, \quad \forall t \in \{1,\dots, T\}\\ \end{aligned} ``` +on which ``\text{pf} = \sin(\arctan(Q^\text{ld,max}/P^\text{ld,max}))``. + +## Valid configurations + +Valid `DeviceModel`s for subtypes of `ElectricLoad` include the following: + +```@eval +using PowerSimulations +using PowerSystems +using DataFrames +using Latexify +combos = PowerSimulations.generate_device_formulation_combinations() +filter!(x -> x["device_type"] <: ElectricLoad, combos) +combo_table = DataFrame( + "Valid DeviceModel" => ["`DeviceModel($(c["device_type"]), $(c["formulation"]))`" for c in combos], + "Device Type" => ["[$(c["device_type"])](https://nrel-Sienna.github.io/PowerSystems.jl/stable/model_library/generated_$(c["device_type"])/)" for c in combos], + "Formulation" => ["[$(c["formulation"])](@ref)" for c in combos], + ) +mdtable(combo_table, latex = false) +``` diff --git a/docs/src/formulation_library/Network.md b/docs/src/formulation_library/Network.md index e5f5e742e4..1a4b01d67f 100644 --- a/docs/src/formulation_library/Network.md +++ b/docs/src/formulation_library/Network.md @@ -1,3 +1,131 @@ # [Network Formulations](@id network_formulations) -TODO +Network formulations are used to describe how the network and buses are handled when constructing constraints. The most common constraint decided by the network formulation is the supply-demand balance constraint. Available Network Models are: + +| Formulation | Description | +| ----- | ---- | +| `CopperPlatePowerModel` | Copper plate connection between all components, i.e. infinite transmission capacity | +| `AreaBalancePowerModel` | Network model approximation to represent inter-area flow with each area represented as a single node | +| `PTDFPowerModel` | Uses the PTDF factor matrix to compute the fraction of power transferred in the network across the branches | + +[`PowerModels.jl`](https://github.com/lanl-ansi/PowerModels.jl) available formulations: +- Exact non-convex models: `ACPPowerModel`, `ACRPowerModel`, `ACTPowerModel`. +- Linear approximations: `DCPPowerModel`, `NFAPowerModel`. +- Quadratic approximations: `DCPLLPowerModel`, `LPACCPowerModel` +- Quadratic relaxations: `SOCWRPowerModel`, `SOCWRConicPowerModel`, `SOCBFPowerModel`, `SOCBFConicPowerModel`, `QCRMPowerModel`, `QCLSPowerModel`. +- SDP relaxations: `SDPWRMPowerModel`, `SparseSDPWRMPowerModel`. + +All of these formulations are described in the [PowerModels.jl documentation](https://lanl-ansi.github.io/PowerModels.jl/stable/formulation-details/) and will not be described here. + +--- + +## `CopperPlatePowerModel` + +```@docs +CopperPlatePowerModel +``` + +**Variables:** + +If Slack variables are enabled: +- [`SystemBalanceSlackUp`](@ref): + - Bounds: [0.0, ] + - Default initial value: 0.0 + - Default proportional cost: 1e6 + - Symbol: ``p^\text{sl,up}`` +- [`SystemBalanceSlackDown`](@ref): + - Bounds: [0.0, ] + - Default initial value: 0.0 + - Default proportional cost: 1e6 + - Symbol: ``p^\text{sl,dn}`` + +**Objective:** + +Add a large proportional cost to the objective function if slack variables are used ``+ (p^\text{sl,up} + p^\text{sl,dn}) \cdot 10^6`` + +**Expressions:** + +Adds ``p^\text{sl,up}`` and ``p^\text{sl,dn}`` terms to the respective active power balance expressions `ActivePowerBalance` created by this `CopperPlatePowerModel` network formulation. + +**Constraints:** + +Adds the `CopperPlateBalanceConstraint` to balance the active power of all components available in the system + +```math +\begin{align} +& \sum_{c \in \text{components}} p_t^c = 0, \quad \forall t \in \{1, \dots, T\} +\end{align} +``` + +--- + +## `AreaBalancePowerModel` + +```@docs +AreaBalancePowerModel +``` + +**Variables:** + +Slack variables are not supported for `AreaBalancePowerModel` + +**Objective:** + +No changes to the objective function. + +**Expressions:** + +Creates `ActivePowerBalance` expressions for each bus that then are used to balance active power for all buses within a single area. + +**Constraints:** + +Adds the `AreaDispatchBalanceConstraint` to balance the active power of all components available in an area. + +```math +\begin{align} +& \sum_{c \in \text{components}_a} p_t^c = 0, \quad \forall a\in \{1,\dots, A\}, t \in \{1, \dots, T\} +\end{align} +``` + +--- + +## `PTDFPowerModel` + +```@docs +PTDFPowerModel +``` + +**Variables:** + +If Slack variables are enabled: +- [`SystemBalanceSlackUp`](@ref): + - Bounds: [0.0, ] + - Default initial value: 0.0 + - Default proportional cost: 1e6 + - Symbol: ``p^\text{sl,up}`` +- [`SystemBalanceSlackDown`](@ref): + - Bounds: [0.0, ] + - Default initial value: 0.0 + - Default proportional cost: 1e6 + - Symbol: ``p^\text{sl,dn}`` + +**Objective:** + +Add a large proportional cost to the objective function if slack variables are used ``+ (p^\text{sl,up} + p^\text{sl,dn}) \cdot 10^6`` + +**Expressions:** + +Adds ``p^\text{sl,up}`` and ``p^\text{sl,dn}`` terms to the respective system-wide active power balance expressions `ActivePowerBalance` created by this `CopperPlatePowerModel` network formulation. In addition, it creates `ActivePowerBalance` expressions for each bus to be used in the calculation of branch flows. + +**Constraints:** + +Adds the `CopperPlateBalanceConstraint` to balance the active power of all components available in the system + +```math +\begin{align} +& \sum_{c \in \text{components}} p_t^c = 0, \quad \forall t \in \{1, \dots, T\} +\end{align} +``` + +In addition creates `NodalBalanceActiveConstraint` for HVDC buses balance, if DC components are connected to an HVDC network. + diff --git a/docs/src/formulation_library/Piecewise.md b/docs/src/formulation_library/Piecewise.md new file mode 100644 index 0000000000..2167769162 --- /dev/null +++ b/docs/src/formulation_library/Piecewise.md @@ -0,0 +1,77 @@ +# [Piecewise linear cost functions](@id pwl_cost) + +The choice for piecewise-linear (PWL) cost representation in `PowerSimulations.jl` is equivalent to the so-called λ-model from the paper [_The Impacts of Convex Piecewise Linear Cost Formulations on AC Optimal Power Flow_](https://www.sciencedirect.com/science/article/pii/S0378779621001723). The SOS constraints in each model are only implemented if the data for PWL is not convex. + +## Special Ordered Set (SOS) Constraints + +A special ordered set (SOS) is an ordered set of variables used as an additional way to specify integrality conditions in an optimization model. + +- Special Ordered Sets of type 1 (SOS1) are a set of variables, at most one of which can take a non-zero value, all others being at 0. They most frequently applications is in a a set of variables that are actually binary variables: in other words, we have to choose at most one from a set of possibilities. +- Special Ordered Sets of type 2 (SOS2) are an ordered set of non-negative variables, of which at most two can be non-zero, and if two are non-zero these must be consecutive in their ordering. Special Ordered Sets of type 2 are typically used to model non-linear functions of a variable in a linear model, such as non-convex quadratic functions using PWL functions. + +## Standard representation of PWL costs + +Piecewise-linear costs are defined by a sequence of points representing the line segments for each generator: ``(P_k^\text{max}, C_k)`` on which we assume ``C_k`` is the cost of generating ``P_k^\text{max}`` power, and ``k \in \{1,\dots, K\}`` are the number of segments each generator cost function has. + +!!! note + `PowerSystems` has more options to specify cost functions for each thermal unit. Independent of which form of the cost data is provided, `PowerSimulations.jl` will internally transform the data to use the λ-model formulation. See **TODO: ADD PSY COST DOCS** for more information. + +### Commitment formulation + + With this the standard representation of PWL costs for a thermal unit commitment is given by: + +```math +\begin{align*} + \min_{\substack{p_{t}, \delta_{k,t}}} + & \sum_{t \in \mathcal{T}} \left(\sum_{k \in \mathcal{K}} C_{k,t} \delta_{k,t} \right) \Delta t\\ + & \sum_{k \in \mathcal{K}} P_{k}^{\text{max}} \delta_{k,t} = p_{t} & \forall t \in \mathcal{T}\\ + & \sum_{k \in \mathcal{K}} \delta_{k,t} = u_{t} & \forall t \in \mathcal{T}\\ + & P^{\text{min}} u_{t} \leq p_{t} \leq P^{\text{max}} u_{t} & \forall t \in \mathcal{T}\\ + &\left \{\delta_{1,t}, \dots, \delta_{K,t} \right \} \in \text{SOS}_{2} & \forall t \in \mathcal{T} +\end{align*} +``` +on which ``\delta_{k,t} \in [0,1]`` is the interpolation variable, ``p`` is the active power of the generator and ``u \in \{0,1\}`` is the commitment variable of the generator. In the case of a PWL convex costs, i.e. increasing slopes, the SOS constraint is omitted. + +### Dispatch formulation + +```math +\begin{align*} + \min_{\substack{p_{t}, \delta_{k,t}}} + & \sum_{t \in \mathcal{T}} \left(\sum_{k \in \mathcal{K}} C_{k,t} \delta_{k,t} \right) \Delta t\\ + & \sum_{k \in \mathcal{K}} P_{k}^{\text{max}} \delta_{k,t} = p_{t} & \forall t \in \mathcal{T}\\ + & \sum_{k \in \mathcal{K}} \delta_{k,t} = \text{on}_{t} & \forall t \in \mathcal{T}\\ + & P^{\text{min}} \text{on}_{t} \leq p_{t} \leq P^{\text{max}} \text{on}_{t} & \forall t \in \mathcal{T}\\ + &\left \{\delta_{i,t}, \dots, \delta_{k,t} \right \} \in \text{SOS}_{2} & \forall t \in \mathcal{T} +\end{align*} +``` +on which ``\delta_{k,t} \in [0,1]`` is the interpolation variable, ``p`` is the active power of the generator and ``\text{on} \in \{0,1\}`` is the parameter that decides if the generator is available or not. In the case of a PWL convex costs, i.e. increasing slopes, the SOS constraint is omitted. + +## Compact representation of PWL costs + +### Commitment Formulation + +```math +\begin{align*} + \min_{\substack{p_{t}, \delta_{k,t}}} + & \sum_{t \in \mathcal{T}} \left(\sum_{k \in \mathcal{K}} C_{k,t} \delta_{k,t} \right) \Delta t\\ + & \sum_{k \in \mathcal{K}} P_{k}^{\text{max}} \delta_{k,t} = P^{\text{min}} u_{t} + \Delta p_{t} & \forall t \in \mathcal{T}\\ + & \sum_{k \in \mathcal{K}} \delta_{k,t} = u_{t} & \forall t \in \mathcal{T}\\ + & 0 \leq \Delta p_{t} \leq \left( P^{\text{max}} - P^{\text{min}} \right)u_{t} & \forall t \in \mathcal{T}\\ + &\left \{\delta_{i,t} \dots \delta_{k,t} \right \} \in \text{SOS}_{2} & \forall t \in \mathcal{T} +\end{align*} +``` +on which ``\delta_{k,t} \in [0,1]`` is the interpolation variable, ``\Delta p`` is the active power of the generator above the minimum power and ``u \in \{0,1\}`` is the commitment variable of the generator. In the case of a PWL convex costs, i.e. increasing slopes, the SOS constraint is omitted. + +### Dispatch formulation + +```math +\begin{align*} + \min_{\substack{p_{t}, \delta_{k,t}}} + & \sum_{t \in \mathcal{T}} \left(\sum_{k \in \mathcal{K}} C_{k,t} \delta_{k,t} \right) \Delta t\\ + & \sum_{k \in \mathcal{K}} P_{k}^{\text{max}} \delta_{k,t} = P^{\text{min}} \text{on}_{t} + \Delta p_{t} & \forall t \in \mathcal{T}\\ + & \sum_{k \in \mathcal{K}} \delta_{k,t} = \text{on}_{t} & \forall t \in \mathcal{T}\\ + & 0 \leq \Delta p_{t} \leq \left( P^{\text{max}} - P^{\text{min}} \right)\text{on}_{t} & \forall t \in \mathcal{T}\\ + &\left \{\delta_{i,t} \dots \delta_{k,t} \right \} \in \text{SOS}_{2} & \forall t \in \mathcal{T} +\end{align*} +``` +on which ``\delta_{k,t} \in [0,1]`` is the interpolation variable, ``\Delta p`` is the active power of the generator above the minimum power and ``u \in \{0,1\}`` is the commitment variable of the generator. In the case of a PWL convex costs, i.e. increasing slopes, the SOS constraint is omitted. \ No newline at end of file diff --git a/docs/src/formulation_library/RenewableGen.md b/docs/src/formulation_library/RenewableGen.md index dd2de3122d..5dfca92c3e 100644 --- a/docs/src/formulation_library/RenewableGen.md +++ b/docs/src/formulation_library/RenewableGen.md @@ -1,21 +1,18 @@ # `PowerSystems.RenewableGen` Formulations -Valid `DeviceModel`s for subtypes of `RenewableGen` include the following: +Renewable generation formulations define the optimization models that describe renewable units mathematical model in different operational settings, such as economic dispatch and unit commitment. -```@eval -using PowerSimulations -using PowerSystems -using DataFrames -using Latexify -combos = PowerSimulations.generate_device_formulation_combinations() -filter!(x -> x["device_type"] <: RenewableGen, combos) -combo_table = DataFrame( - "Valid DeviceModel" => ["`DeviceModel($(c["device_type"]), $(c["formulation"]))`" for c in combos], - "Device Type" => ["[$(c["device_type"])](https://nrel-Sienna.github.io/PowerSystems.jl/stable/model_library/generated_$(c["device_type"])/)" for c in combos], - "Formulation" => ["[$(c["formulation"])](@ref)" for c in combos], - ) -mdtable(combo_table, latex = false) -``` +!!! note + The use of reactive power variables and constraints will depend on the network model used, i.e., whether it uses (or does not use) reactive power. If the network model is purely active power-based, reactive power variables and related constraints are not created. + +!!! note + Reserve variables for services are not included in the formulation, albeit their inclusion change the variables, expressions, constraints and objective functions created. A detailed description of the implications in the optimization models is described in the [Service formulation](@ref service_formulations) section. + +### Table of contents + +1. [`RenewableFullDispatch`](#RenewableFullDispatch) +2. [`RenewableConstantPowerFactor`](#RenewableConstantPowerFactor) +3. [Valid configurations](#Valid-configurations) --- @@ -29,19 +26,21 @@ RenewableFullDispatch - [`ActivePowerVariable`](@ref): - Bounds: [0.0, ] - - Default initial value: `PowerSystems.get_active_power(device)` + - Symbol: ``p^\text{re}`` - [`ReactivePowerVariable`](@ref): - Bounds: [0.0, ] - - Default initial value: `PowerSystems.get_reactive_power(device)` + - Symbol: ``q^\text{re}`` **Static Parameters:** -- ``Pg^\text{min}`` = `PowerSystems.get_active_power_limits(device).min` -- ``Qg^\text{min}`` = `PowerSystems.get_reactive_power_limits(device).min` -- ``Qg^\text{max}`` = `PowerSystems.get_reactive_power_limits(device).max` +- ``P^\text{re,min}`` = `PowerSystems.get_active_power_limits(device).min` +- ``Q^\text{re,min}`` = `PowerSystems.get_reactive_power_limits(device).min` +- ``Q^\text{re,max}`` = `PowerSystems.get_reactive_power_limits(device).max` **Time Series Parameters:** +Uses the `max_active_power` timeseries parameter to limit the available active power at each time-step. + ```@eval using PowerSimulations using PowerSystems @@ -57,18 +56,19 @@ mdtable(combo_table, latex = false) **Objective:** -Creates an objective function term based on the [`FunctionData` Options](@ref) where the quantity term is defined as ``- Pg_t`` to incentivize generation from `RenewableGen` devices. +Creates an objective function term based on the [`FunctionData` Options](@ref) where the quantity term is defined as ``- p^\text{re}`` to incentivize generation from `RenewableGen` devices. + **Expressions:** -Adds ``Pg`` and ``Qg`` terms to the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations) +Adds ``p^\text{re}`` and ``q^\text{re}`` terms to the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations). **Constraints:** ```math \begin{aligned} -& Pg^\text{min} \le Pg_t \le ActivePowerTimeSeriesParameter_t \\ -& Qg^\text{min} \le Qg_t \le Qg^\text{max} +& P^\text{re,min} \le p_t^\text{re} \le \text{ActivePowerTimeSeriesParameter}_t, \quad \forall t \in \{1,\dots, T\} \\ +& Q^\text{re,min} \le q_t^\text{re} \le Q^\text{re,max}, \quad \forall t \in \{1,\dots, T\} \end{aligned} ``` @@ -85,16 +85,18 @@ RenewableConstantPowerFactor - [`ActivePowerVariable`](@ref): - Bounds: [0.0, ] - Default initial value: `PowerSystems.get_active_power(device)` + - Symbol: ``p^\text{re}`` - [`ReactivePowerVariable`](@ref): - Bounds: [0.0, ] - Default initial value: `PowerSystems.get_reactive_power(device)` + - Symbol: ``q^\text{re}`` **Static Parameters:** -- ``Pg^\text{min}`` = `PowerSystems.get_active_power_limits(device).min` -- ``Qg^\text{min}`` = `PowerSystems.get_reactive_power_limits(device).min` -- ``Qg^\text{max}`` = `PowerSystems.get_reactive_power_limits(device).max` -- ``pf`` = `PowerSystems.get_power_factor(device)` +- ``P^\text{re,min}`` = `PowerSystems.get_active_power_limits(device).min` +- ``Q^\text{re,min}`` = `PowerSystems.get_reactive_power_limits(device).min` +- ``Q^\text{re,max}`` = `PowerSystems.get_reactive_power_limits(device).max` +- ``\text{pf}`` = `PowerSystems.get_power_factor(device)` **Time Series Parameters:** @@ -113,18 +115,39 @@ mdtable(combo_table, latex = false) **Objective:** -Creates an objective function term based on the [`FunctionData` Options](@ref) where the quantity term is defined as ``- Pg_t`` to incentivize generation from `RenewableGen` devices. +Creates an objective function term based on the [`FunctionData` Options](@ref) where the quantity term is defined as ``- p_t^\text{re}`` to incentivize generation from `RenewableGen` devices. **Expressions:** -Adds ``Pg`` and ``Qg`` terms to the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations) +Adds ``p^\text{re}`` and ``q^\text{re}`` terms to the respective active and reactive power balance expressions created by the selected [Network Formulations](@ref network_formulations) **Constraints:** ```math \begin{aligned} -& Pg^\text{min} \le Pg_t \le ActivePowerTimeSeriesParameter_t \\ -& Qg^\text{min} \le Qg_t \le Qg^\text{max} \\ -& Qg_t = pf * Pg_t +& P^\text{re,min} \le p_t^\text{re} \le \text{ActivePowerTimeSeriesParameter}_t, \quad \forall t \in \{1,\dots, T\} \\ +& q_t^\text{re} = \text{pf} \cdot p_t^\text{re}, \quad \forall t \in \{1,\dots, T\} \end{aligned} ``` + +--- + +## Valid configurations + +Valid `DeviceModel`s for subtypes of `RenewableGen` include the following: + +```@eval +using PowerSimulations +using PowerSystems +using DataFrames +using Latexify +combos = PowerSimulations.generate_device_formulation_combinations() +filter!(x -> x["device_type"] <: RenewableGen, combos) +combo_table = DataFrame( + "Valid DeviceModel" => ["`DeviceModel($(c["device_type"]), $(c["formulation"]))`" for c in combos], + "Device Type" => ["[$(c["device_type"])](https://nrel-Sienna.github.io/PowerSystems.jl/stable/model_library/generated_$(c["device_type"])/)" for c in combos], + "Formulation" => ["[$(c["formulation"])](@ref)" for c in combos], + ) +mdtable(combo_table, latex = false) +``` + diff --git a/docs/src/formulation_library/Service.md b/docs/src/formulation_library/Service.md index f4331eba63..61aaf1d335 100644 --- a/docs/src/formulation_library/Service.md +++ b/docs/src/formulation_library/Service.md @@ -1,3 +1,495 @@ -# `PowerSystems.Service` Formulations +# [`PowerSystems.Service` Formulations](@id service_formulations) -TODO +`Services` (or ancillary services) are models used to ensure that there is necessary support to the power grid from generators to consumers, in order to ensure reliable operation of the system. + +The most common application for ancillary services are reserves, i.e., generation (or load) that is not currently being used, but can be quickly made available in case of unexpected changes of grid conditions, for example a sudden loss of load or generation. + +A key challenge of adding services to a system, from a mathematical perspective, is specifying which units contribute to the specified requirement of a service, that implies the creation of new variables (such as reserve variables) and modification of constraints. + +In this documentation, we first specify the available `Services` in the grid, and what requirements impose in the system, and later we discuss the implication on device formulations for specific units. + +### Table of contents + +1. [`RangeReserve`](#RangeReserve) +2. [`StepwiseCostReserve`](#StepwiseCostReserve) +3. [`GroupReserve`](#GroupReserve) +4. [`RampReserve`](#RampReserve) +5. [`NonSpinningReserve`](#NonSpinningReserve) +6. [`ConstantMaxInterfaceFlow`](#ConstantMaxInterfaceFlow) +7. [Changes on Expressions](#Changes-on-Expressions-due-to-Service-models) + +--- + +## `RangeReserve` + +```@docs +RangeReserve +``` + +For each service ``s`` of the model type `RangeReserve` the following variables are created: + +**Variables**: + +- [`ActivePowerReserveVariable`](@ref): + - Bounds: [0.0, ] + - Default proportional cost: ``1.0 / \text{SystemBasePower}`` + - Symbol: ``r_{d}`` for ``d`` in contributing devices to the service ``s`` +If slacks are enabled: +- [`ReserveRequirementSlack`](@ref): + - Bounds: [0.0, ] + - Default proportional cost: 1e5 + - Symbol: ``r^\text{sl}`` + +Depending on the `PowerSystems.jl` type associated to the `RangeReserve` formulation model, the parameters are: + +**Static Parameters** + +- ``\text{PF}`` = `PowerSystems.get_max_participation_factor(service)` + +For a `StaticReserve` `PowerSystems` type: +- ``\text{Req}`` = `PowerSystems.get_requirement(service)` + +**Time Series Parameters** + +For a `VariableReserve` `PowerSystems` type: +```@eval +using PowerSimulations +using PowerSystems +using DataFrames +using Latexify +combos = PowerSimulations.get_default_time_series_names(VariableReserve, RangeReserve) +combo_table = DataFrame( + "Parameter" => map(x -> "[`$x`](@ref)", collect(keys(combos))), + "Default Time Series Name" => map(x -> "`$x`", collect(values(combos))), + ) +mdtable(combo_table, latex = false) +``` + +**Relevant Methods:** + +- ``\mathcal{D}_s`` = `PowerSystems.get_contributing_devices(system, service)`: Set (vector) of all contributing devices to the service ``s`` in the system. + +**Objective:** + +Add a large proportional cost to the objective function if slack variables are used ``+ r^\text{sl} \cdot 10^5``. In addition adds the default cost for `ActivePowerReserveVariables` as a proportional cost. + +**Expressions:** + +Adds the `ActivePowerReserveVariable` for upper/lower bound expressions of contributing devices. + +For `ReserveUp` types, the variable is added to `ActivePowerRangeExpressionUB`, such that this expression considers both the `ActivePowerVariable` and its reserve variable. Similarly, For `ReserveDown` types, the variable is added to `ActivePowerRangeExpressionLB`, such that this expression considers both the `ActivePowerVariable` and its reserve variable + + +*Example*: for a thermal unit ``d`` contributing to two different `ReserveUp` ``s_1, s_2`` services (e.g. Reg-Up and Spin): +```math +\text{ActivePowerRangeExpressionUB}_{t} = p_t^\text{th} + r_{s_1,t} + r_{s_2, t} \le P^\text{th,max} +``` +similarly if ``s_3`` is a `ReserveDown` service (e.g. Reg-Down): +```math +\text{ActivePowerRangeExpressionLB}_{t} = p_t^\text{th} - r_{s_3,t} \ge P^\text{th,min} +``` + +**Constraints:** + +A RangeReserve implements two fundamental constraints. The first is that the sum of all reserves of contributing devices must be larger than the `RangeReserve` requirement. Thus, for a service ``s``: + +```math +\sum_{d\in\mathcal{D}_s} r_{d,t} + r_t^\text{sl} \ge \text{Req},\quad \forall t\in \{1,\dots, T\} \quad \text{(for a StaticReserve)} \\ +\sum_{d\in\mathcal{D}_s} r_{d,t} + r_t^\text{sl} \ge \text{RequirementTimeSeriesParameter}_{t},\quad \forall t\in \{1,\dots, T\} \quad \text{(for a VariableReserve)} +``` + +In addition, there is a restriction on how much each contributing device ``d`` can contribute to the requirement, based on the max participation factor allowed. + +```math +r_{d,t} \le \text{Req} \cdot \text{PF} ,\quad \forall d\in \mathcal{D}_s, \forall t\in \{1,\dots, T\} \quad \text{(for a StaticReserve)} \\ +r_{d,t} \le \text{RequirementTimeSeriesParameter}_{t} \cdot \text{PF}\quad \forall d\in \mathcal{D}_s, \forall t\in \{1,\dots, T\}, \quad \text{(for a VariableReserve)} +``` + +--- + +## `StepwiseCostReserve` + +Service must be used with `ReserveDemandCurve` `PowerSystems.jl` type. This service model is used to model ORDC (Operating Reserve Demand Curve) in ERCOT. + +```@docs +StepwiseCostReserve +``` + +For each service ``s`` of the model type `ReserveDemandCurve` the following variables are created: + +**Variables**: + +- [`ActivePowerReserveVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``r_{d}`` for ``d`` in contributing devices to the service ``s`` +- [`ServiceRequirementVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``\text{req}`` + +**Time Series Parameters** + +For a `ReserveDemandCurve` `PowerSystems` type: +```@eval +using PowerSimulations +using PowerSystems +using DataFrames +using Latexify +combos = PowerSimulations.get_default_time_series_names(ReserveDemandCurve, StepwiseCostReserve) +combo_table = DataFrame( + "Parameter" => map(x -> "[`$x`](@ref)", collect(keys(combos))), + "Default Time Series Name" => map(x -> "`$x`", collect(values(combos))), + ) +mdtable(combo_table, latex = false) +``` + +**Relevant Methods:** + +- ``\mathcal{D}_s`` = `PowerSystems.get_contributing_devices(system, service)`: Set (vector) of all contributing devices to the service ``s`` in the system. + +**Objective:** + +The `ServiceRequirementVariable` is added as a piecewise linear cost based on the decreasing offers listed in the `variable_cost` time series. These decreasing cost represent the scarcity prices of not having sufficient reserves. For example, if the variable ``\text{req} = 0``, then a really high cost is paid for not having enough reserves, and if ``\text{req}`` is larger, then a lower cost (or even zero) is paid. + +**Expressions:** + +Adds the `ActivePowerReserveVariable` for upper/lower bound expressions of contributing devices. + +For `ReserveUp` types, the variable is added to `ActivePowerRangeExpressionUB`, such that this expression considers both the `ActivePowerVariable` and its reserve variable. Similarly, For `ReserveDown` types, the variable is added to `ActivePowerRangeExpressionLB`, such that this expression considers both the `ActivePowerVariable` and its reserve variable + + +*Example*: for a thermal unit ``d`` contributing to two different `ReserveUp` ``s_1, s_2`` services (e.g. Reg-Up and Spin): +```math +\text{ActivePowerRangeExpressionUB}_{t} = p_t^\text{th} + r_{s_1,t} + r_{s_2, t} \le P^\text{th,max} +``` +similarly if ``s_3`` is a `ReserveDown` service (e.g. Reg-Down): +```math +\text{ActivePowerRangeExpressionLB}_{t} = p_t^\text{th} - r_{s_3,t} \ge P^\text{th,min} +``` + +**Constraints:** + +A `StepwiseCostReserve` implements a single constraint, such that the sum of all reserves of contributing devices must be larger than the `ServiceRequirementVariable` variable. Thus, for a service ``s``: + +```math +\sum_{d\in\mathcal{D}_s} r_{d,t} \ge \text{req}_t,\quad \forall t\in \{1,\dots, T\} +``` + +## `GroupReserve` + +Service must be used with `StaticReserveGroup` `PowerSystems.jl` type. This service model is used to model an aggregation of services. + +```@docs +GroupReserve +``` + +For each service ``s`` of the model type `GroupReserve` the following variables are created: + +**Variables**: + +No variables are created, but the services associated with the `GroupReserve` must have created variables. + +**Static Parameters** + +- ``\text{Req}`` = `PowerSystems.get_requirement(service)` + +**Relevant Methods:** + +- ``\mathcal{S}_s`` = `PowerSystems.get_contributing_services(system, service)`: Set (vector) of all contributing services to the group service ``s`` in the system. +- ``\mathcal{D}_{s_i}`` = `PowerSystems.get_contributing_devices(system, service_aux)`: Set (vector) of all contributing devices to the service ``s_i`` in the system. + +**Objective:** + +Does not modify the objective function, besides the changes to the objective function due to the other services associated to the group service. + +**Expressions:** + +No changes, besides the changes to the expressions due to the other services associated to the group service. + +**Constraints:** + +A GroupReserve implements that the sum of all reserves of contributing devices, of all contributing services, must be larger than the `GroupReserve` requirement. Thus, for a `GroupReserve` service ``s``: + +```math +\sum_{d\in\mathcal{D}_{s_i}} \sum_{i \in \mathcal{S}_s} r_{d,t} \ge \text{Req},\quad \forall t\in \{1,\dots, T\} +``` + +--- + +## `RampReserve` + +```@docs +RampReserve +``` + +For each service ``s`` of the model type `RampReserve` the following variables are created: + +**Variables**: + +- [`ActivePowerReserveVariable`](@ref): + - Bounds: [0.0, ] + - Default proportional cost: ``1.0 / \text{SystemBasePower}`` + - Symbol: ``r_{d}`` for ``d`` in contributing devices to the service ``s`` +If slacks are enabled: +- [`ReserveRequirementSlack`](@ref): + - Bounds: [0.0, ] + - Default proportional cost: 1e5 + - Symbol: ``r^\text{sl}`` + +`RampReserve` only accepts `VariableReserve` `PowerSystems.jl` type. With that, the parameters are: + +**Static Parameters** + +- ``\text{TF}`` = `PowerSystems.get_time_frame(service)` +- ``R^\text{th,up}`` = `PowerSystems.get_ramp_limits(device).up` for thermal contributing devices +- ``R^\text{th,dn}`` = `PowerSystems.get_ramp_limits(device).down` for thermal contributing devices + + +**Time Series Parameters** + +For a `VariableReserve` `PowerSystems` type: +```@eval +using PowerSimulations +using PowerSystems +using DataFrames +using Latexify +combos = PowerSimulations.get_default_time_series_names(VariableReserve, RampReserve) +combo_table = DataFrame( + "Parameter" => map(x -> "[`$x`](@ref)", collect(keys(combos))), + "Default Time Series Name" => map(x -> "`$x`", collect(values(combos))), + ) +mdtable(combo_table, latex = false) +``` + +**Relevant Methods:** + +- ``\mathcal{D}_s`` = `PowerSystems.get_contributing_devices(system, service)`: Set (vector) of all contributing devices to the service ``s`` in the system. + +**Objective:** + +Add a large proportional cost to the objective function if slack variables are used ``+ r^\text{sl} \cdot 10^5``. In addition adds the default cost for `ActivePowerReserveVariables` as a proportional cost. + +**Expressions:** + +Adds the `ActivePowerReserveVariable` for upper/lower bound expressions of contributing devices. + +For `ReserveUp` types, the variable is added to `ActivePowerRangeExpressionUB`, such that this expression considers both the `ActivePowerVariable` and its reserve variable. Similarly, For `ReserveDown` types, the variable is added to `ActivePowerRangeExpressionLB`, such that this expression considers both the `ActivePowerVariable` and its reserve variable + + +*Example*: for a thermal unit ``d`` contributing to two different `ReserveUp` ``s_1, s_2`` services (e.g. Reg-Up and Spin): +```math +\text{ActivePowerRangeExpressionUB}_{t} = p_t^\text{th} + r_{s_1,t} + r_{s_2, t} \le P^\text{th,max} +``` +similarly if ``s_3`` is a `ReserveDown` service (e.g. Reg-Down): +```math +\text{ActivePowerRangeExpressionLB}_{t} = p_t^\text{th} - r_{s_3,t} \ge P^\text{th,min} +``` + +**Constraints:** + +A RampReserve implements three fundamental constraints. The first is that the sum of all reserves of contributing devices must be larger than the `RampReserve` requirement. Thus, for a service ``s``: + +```math +\sum_{d\in\mathcal{D}_s} r_{d,t} + r_t^\text{sl} \ge \text{RequirementTimeSeriesParameter}_{t},\quad \forall t\in \{1,\dots, T\} +``` + +Finally, there is a restriction based on the ramp limits of the contributing devices: + +```math +r_{d,t} \le R^\text{th,up} \cdot \text{TF}\quad \forall d\in \mathcal{D}_s, \forall t\in \{1,\dots, T\}, \quad \text{(for ReserveUp)} \\ +r_{d,t} \le R^\text{th,dn} \cdot \text{TF}\quad \forall d\in \mathcal{D}_s, \forall t\in \{1,\dots, T\}, \quad \text{(for ReserveDown)} +``` + +--- + +## `NonSpinningReserve` + +```@docs +NonSpinningReserve +``` + +For each service ``s`` of the model type `NonSpinningReserve`, the following variables are created: + +**Variables**: + +- [`ActivePowerReserveVariable`](@ref): + - Bounds: [0.0, ] + - Default proportional cost: ``1.0 / \text{SystemBasePower}`` + - Symbol: ``r_{d}`` for ``d`` in contributing devices to the service ``s`` +If slacks are enabled: +- [`ReserveRequirementSlack`](@ref): + - Bounds: [0.0, ] + - Default proportional cost: 1e5 + - Symbol: ``r^\text{sl}`` + +`NonSpinningReserve` only accepts `VariableReserve` `PowerSystems.jl` type. With that, the parameters are: + +**Static Parameters** + +- ``\text{PF}`` = `PowerSystems.get_max_participation_factor(service)` +- ``\text{TF}`` = `PowerSystems.get_time_frame(service)` +- ``P^\text{th,min}`` = `PowerSystems.get_active_power_limits(device).min` for thermal contributing devices +- ``T^\text{st,up}`` = `PowerSystems.get_time_limits(d).up` for thermal contributing devices +- ``R^\text{th,up}`` = `PowerSystems.get_ramp_limits(device).down` for thermal contributing devices + +Other parameters: + +- ``\Delta T``: Resolution of the problem in minutes. + +**Time Series Parameters** + +For a `VariableReserve` `PowerSystems` type: +```@eval +using PowerSimulations +using PowerSystems +using DataFrames +using Latexify +combos = PowerSimulations.get_default_time_series_names(VariableReserve, NonSpinningReserve) +combo_table = DataFrame( + "Parameter" => map(x -> "[`$x`](@ref)", collect(keys(combos))), + "Default Time Series Name" => map(x -> "`$x`", collect(values(combos))), + ) +mdtable(combo_table, latex = false) +``` + +**Relevant Methods:** + +- ``\mathcal{D}_s`` = `PowerSystems.get_contributing_devices(system, service)`: Set (vector) of all contributing devices to the service ``s`` in the system. + +**Objective:** + +Add a large proportional cost to the objective function if slack variables are used ``+ r^\text{sl} \cdot 10^5``. In addition adds the default cost for `ActivePowerReserveVariables` as a proportional cost. + +**Expressions:** + +Adds the `ActivePowerReserveVariable` for upper/lower bound expressions of contributing devices. + +For `ReserveUp` types, the variable is added to `ActivePowerRangeExpressionUB`, such that this expression considers both the `ActivePowerVariable` and its reserve variable. Similarly, For `ReserveDown` types, the variable is added to `ActivePowerRangeExpressionLB`, such that this expression considers both the `ActivePowerVariable` and its reserve variable + + +*Example*: for a thermal unit ``d`` contributing to two different `ReserveUp` ``s_1, s_2`` services (e.g. Reg-Up and Spin): +```math +\text{ActivePowerRangeExpressionUB}_{t} = p_t^\text{th} + r_{s_1,t} + r_{s_2, t} \le P^\text{th,max} +``` +similarly if ``s_3`` is a `ReserveDown` service (e.g. Reg-Down): +```math +\text{ActivePowerRangeExpressionLB}_{t} = p_t^\text{th} - r_{s_3,t} \ge P^\text{th,min} +``` + +**Constraints:** + +A NonSpinningReserve implements three fundamental constraints. The first is that the sum of all reserves of contributing devices must be larger than the `NonSpinningReserve` requirement. Thus, for a service ``s``: + +```math +\sum_{d\in\mathcal{D}_s} r_{d,t} + r_t^\text{sl} \ge \text{RequirementTimeSeriesParameter}_{t},\quad \forall t\in \{1,\dots, T\} +``` + +In addition, there is a restriction on how much each contributing device ``d`` can contribute to the requirement, based on the max participation factor allowed. + +```math +r_{d,t} \le \text{RequirementTimeSeriesParameter}_{t} \cdot \text{PF}\quad \forall d\in \mathcal{D}_s, \forall t\in \{1,\dots, T\}, +``` + +Finally, there is a restriction based on the reserve response time for the non-spinning reserve if the unit is off. To do so, compute ``R^\text{limit}_d`` as the reserve response limit as: +```math +R^\text{limit}_d = \begin{cases} +0 & \text{ if TF } \le T^\text{st,up}_d \\ +P^\text{th,min}_d + (\text{TF}_s - T^\text{st,up}_d) \cdot R^\text{th,up}_d \Delta T \cdot R^\text{th,up}_d & \text{ if TF } > T^\text{st,up}_d +\end{cases}, \quad \forall d\in \mathcal{D}_s +``` + +Then, the constraint depends on the commitment variable ``u_t^\text{th}`` as: + +```math +r_{d,t} \le (1 - u_{d,t}^\text{th}) \cdot R^\text{limit}_d, \quad \forall d \in \mathcal{D}_s, \forall t \in \{1,\dots, T\} +``` + +--- + +## `ConstantMaxInterfaceFlow` + +This Service model only accepts the `PowerSystems.jl` `TransmissionInterface` type to properly function. It is used to model a collection of branches that make up an interface or corridor with a maximum transfer of power. + +```@docs +ConstantMaxInterfaceFlow +``` + +**Variables** + +If slacks are used: +- [`InterfaceFlowSlackUp`](@ref): + - Bounds: [0.0, ] + - Symbol: ``f^\text{sl,up}`` +- [`InterfaceFlowSlackDown`](@ref): + - Bounds: [0.0, ] + - Symbol: ``f^\text{sl,dn}`` + +**Static Parameters** + +- ``F^\text{max}`` = `PowerSystems.get_active_power_flow_limits(service).max` +- ``F^\text{min}`` = `PowerSystems.get_active_power_flow_limits(service).min` +- ``C^\text{flow}`` = `PowerSystems.get_violation_penalty(service)` +- ``\mathcal{M}_s`` = `PowerSystems.get_direction_mapping(service)`. Dictionary of contributing branches with its specified direction (``\text{Dir}_d = 1`` or ``\text{Dir}_d = -1``) with respect to the interface. + +**Relevant Methods** + +- ``\mathcal{D}_s`` = `PowerSystems.get_contributing_devices(system, service)`: Set (vector) of all contributing branches to the service ``s`` in the system. + +**Objective:** + +Add the violation penalty proportional cost to the objective function if slack variables are used ``+ (f^\text{sl,up} + f^\text{sl,dn}) \cdot C^\text{flow}``. + +**Expressions:** + +Creates the expression `InterfaceTotalFlow` to keep track of all `FlowActivePowerVariable` of contributing branches to the transmission interface. + +**Constraints:** + +It adds the constraint to limit the `InterfaceTotalFlow` by the specified bounds of the service ``s``: + +```math +F^\text{min} \le f^\text{sl,up}_t - f^\text{sl,dn}_t + \sum_{d\in\mathcal{D}_s} \text{Dir}_d f_{d,t} \le F^\text{max}, \quad \forall t \in \{1,\dots,T\} +``` + +## Changes on Expressions due to Service models + +It is important to note that by adding a service to a Optimization Problem, variables for each contributing device must be created. For example, for every contributing generator ``d \in \mathcal{D}`` that is participating in services ``s_1,s_2,s_3``, it is required to create three set of `ActivePowerReserveVariable` variables: + +```math +r_{s_1,d,t},~ r_{s_2,d,t},~ r_{s_3,d,t},\quad \forall d \in \mathcal{D}, \forall t \in \{1,\dots, T\} +``` + +### Changes on UpperBound (UB) and LowerBound (LB) limits + +Each contributing generator ``d`` has active power limits that the reserve variables affect. In simple terms, the limits are implemented using expressions `ActivePowerRangeExpressionUB` and `ActivePowerRangeExpressionLB` as: + +```math +\text{ActivePowerRangeExpressionUB}_t \le P^\text{max} \\ +\text{ActivePowerRangeExpressionLB}_t \ge P^\text{min} +``` +`ReserveUp` type variables contribute to the upper bound expression, while `ReserveDown` variables contribute to the lower bound expressions. So if ``s_1,s_2`` are `ReserveUp` services, and ``s_3`` is a `ReserveDown` service, then for a thermal generator ``d`` using a `ThermalStandardDispatch`: + +```math +\begin{align*} +& p_{d,t}^\text{th} + r_{s_1,d,t} + r_{s_2,d,t} \le P^\text{th,max},\quad \forall d\in \mathcal{D}^\text{th}, \forall t \in \{1,\dots,T\} \\ +& p_{d,t}^\text{th} - r_{s_3,d,t} \ge P^\text{th,min},\quad \forall d\in \mathcal{D}^\text{th}, \forall t \in \{1,\dots,T\} +\end{align*} +``` + +while for a renewable generator ``d`` using a `RenewableFullDispatch`: + +```math +\begin{align*} +& p_{d,t}^\text{re} + r_{s_1,d,t} + r_{s_2,d,t} \le \text{ActivePowerTimeSeriesParameter}_t,\quad \forall d\in \mathcal{D}^\text{re}, \forall t \in \{1,\dots,T\}\\ +& p_{d,t}^\text{re} - r_{s_3,d,t} \ge 0,\quad \forall d\in \mathcal{D}^\text{re}, \forall t \in \{1,\dots,T\} +\end{align*} +``` + +### Changes in Ramp limits + +For the case of Ramp Limits (of formulation that model these limits), the reserve variables only affect the current time, and not the previous time. Then, for the same example as before: +```math +\begin{align*} +& p_{d,t}^\text{th} + r_{s_1,d,t} + r_{s_2,d,t} - p_{d,t-1}^\text{th}\le R^\text{th,up},\quad \forall d\in \mathcal{D}^\text{th}, \forall t \in \{1,\dots,T\}\\ +& p_{d,t}^\text{th} - r_{s_3,d,t} - p_{d,t-1}^\text{th} \ge -R^\text{th,dn},\quad \forall d\in \mathcal{D}^\text{th}, \forall t \in \{1,\dots,T\} +\end{align*} +``` diff --git a/docs/src/formulation_library/ThermalGen.md b/docs/src/formulation_library/ThermalGen.md index d80072ff2b..e179c8c8e1 100644 --- a/docs/src/formulation_library/ThermalGen.md +++ b/docs/src/formulation_library/ThermalGen.md @@ -1,21 +1,29 @@ # `ThermalGen` Formulations -Valid `DeviceModel`s for subtypes of `ThermalGen` include the following: +Thermal generation formulations define the optimization models that describe thermal units mathematical model in different operational settings, such as economic dispatch and unit commitment. -```@eval -using PowerSimulations -using PowerSystems -using DataFrames -using Latexify -combos = PowerSimulations.generate_device_formulation_combinations() -filter!(x -> x["device_type"] <: ThermalGen, combos) -combo_table = DataFrame( - "Valid DeviceModel" => ["`DeviceModel($(c["device_type"]), $(c["formulation"]))`" for c in combos], - "Device Type" => ["[$(c["device_type"])](https://nrel-Sienna.github.io/PowerSystems.jl/stable/model_library/generated_$(c["device_type"])/)" for c in combos], - "Formulation" => ["[$(c["formulation"])](@ref)" for c in combos], - ) -mdtable(combo_table, latex = false) -``` + +!!! note + Thermal units can include multiple terms added to the objective function, such as no-load cost, turn-on/off cost, fixed cost and variable cost. In addition, variable costs can be linear, quadratic or piecewise-linear formulations. These methods are properly described in the [cost function page](@ref pwl_cost). + + +!!! note + The use of reactive power variables and constraints will depend on the network model used, i.e., whether it uses (or does not use) reactive power. If the network model is purely active power-based, reactive power variables and related constraints are not created. + +!!! note + Reserve variables for services are not included in the formulation, albeit their inclusion change the variables, expressions, constraints and objective functions created. A detailed description of the implications in the optimization models is described in the [Service formulation](@ref service_formulations) section. + +### Table of Contents + +1. [`ThermalBasicDispatch`](#ThermalBasicDispatch) +2. [`ThermalDispatchNoMin`](#ThermalDispatchNoMin) +3. [`ThermalCompactDispatch`](#ThermalCompactDispatch) +4. [`ThermalStandardDispatch`](#ThermalStandardDispatch) +5. [`ThermalBasicUnitCommitment`](#ThermalBasicUnitCommitment) +6. [`ThermalBasicCompactUnitCommitment`](#ThermalBasicCompactUnitCommitment) +7. [`ThermalStandardUnitCommitment`](#ThermalStandardUnitCommitment) +8. [`ThermalMultiStartUnitCommitment`](#ThermalMultiStartUnitCommitment) +9. [Valid configurations](#Valid-configurations) --- @@ -24,38 +32,244 @@ mdtable(combo_table, latex = false) ```@docs ThermalBasicDispatch ``` +**Variables:** -TODO +- [`ActivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``p^\text{th}`` +- [`ReactivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``q^\text{th}`` + +**Static Parameters:** + +- ``P^\text{th,min}`` = `PowerSystems.get_active_power_limits(device).min` +- ``P^\text{th,max}`` = `PowerSystems.get_active_power_limits(device).max` +- ``Q^\text{th,min}`` = `PowerSystems.get_reactive_power_limits(device).min` +- ``Q^\text{th,max}`` = `PowerSystems.get_reactive_power_limits(device).max` + +**Objective:** + +Add a cost to the objective function depending on the defined cost structure of the thermal unit by adding it to its `ProductionCostExpression`. + +**Expressions:** + +Adds ``p^\text{th}`` to the `ActivePowerBalance` expression and ``q^\text{th}`` to the `ReactivePowerBalance`, to be used in the supply-balance constraint depending on the network model used. + +**Constraints:** + +For each thermal unit creates the range constraints for its active and reactive power depending on its static parameters. + +```math +\begin{align*} +& P^\text{th,min} \le p^\text{th}_t \le P^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& Q^\text{th,min} \le q^\text{th}_t \le Q^\text{th,max}, \quad \forall t\in \{1, \dots, T\} +\end{align*} +``` --- -## `ThermalCompactDispatch` + +## `ThermalDispatchNoMin` ```@docs -ThermalCompactDispatch +ThermalDispatchNoMin ``` -TODO +**Variables:** + +- [`ActivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``p^\text{th}`` +- [`ReactivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``q^\text{th}`` + +**Static Parameters:** + +- ``P^\text{th,max}`` = `PowerSystems.get_active_power_limits(device).max` +- ``Q^\text{th,min}`` = `PowerSystems.get_reactive_power_limits(device).min` +- ``Q^\text{th,max}`` = `PowerSystems.get_reactive_power_limits(device).max` + +**Objective:** + +Add a cost to the objective function depending on the defined cost structure of the thermal unit by adding it to its `ProductionCostExpression`. + +**Expressions:** + +Adds ``p^\text{th}`` to the `ActivePowerBalance` expression and ``q^\text{th}`` to the `ReactivePowerBalance`, to be used in the supply-balance constraint depending on the network model used. + +**Constraints:** + +For each thermal unit creates the range constraints for its active and reactive power depending on its static parameters. + +```math +\begin{align} +& 0 \le p^\text{th}_t \le P^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& Q^\text{th,min} \le q^\text{th}_t \le Q^\text{th,max}, \quad \forall t\in \{1, \dots, T\} +\end{align} +``` --- -## `ThermalDispatchNoMin` +## `ThermalCompactDispatch` ```@docs -ThermalDispatchNoMin +ThermalCompactDispatch ``` -TODO +**Variables:** + +- [`PowerAboveMinimumVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``\Delta p^\text{th}`` +- [`ReactivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``q^\text{th}`` + +**Auxiliary Variables:** +- [`PowerOutput`](@ref): + - Symbol: ``P^\text{th}`` + - Definition: ``P^\text{th} = \text{on}^\text{th}P^\text{min} + \Delta p^\text{th}`` + +**Static Parameters:** + +- ``P^\text{th,min}`` = `PowerSystems.get_active_power_limits(device).min` +- ``P^\text{th,max}`` = `PowerSystems.get_active_power_limits(device).max` +- ``Q^\text{th,min}`` = `PowerSystems.get_reactive_power_limits(device).min` +- ``Q^\text{th,max}`` = `PowerSystems.get_reactive_power_limits(device).max` +- ``R^\text{th,up}`` = `PowerSystems.get_ramp_limits(device).up` +- ``R^\text{th,dn}`` = `PowerSystems.get_ramp_limits(device).down` + +**Variable Value Parameters:** + +- ``\text{on}^\text{th}``: Used in feedforwards to define if the unit is on/off at each time-step from another problem. If no feedforward is used, the parameter takes a {0,1} value if the unit is available or not. + +**Objective:** + +Add a cost to the objective function depending on the defined cost structure of the thermal unit by adding it to its `ProductionCostExpression`. + +**Expressions:** + +Adds ``\text{on}^\text{th}P^\text{th,min} + \Delta p^\text{th}`` to the `ActivePowerBalance` expression and ``q^\text{th}`` to the `ReactivePowerBalance`, to be used in the supply-balance constraint depending on the network model used. + +**Constraints:** + +For each thermal unit creates the range constraints for its active and reactive power depending on its static parameters. It also implements ramp constraints for the active power variable. + +```math +\begin{align*} +& 0 \le \Delta p^\text{th}_t \le \text{on}^\text{th}_t\left(P^\text{th,max} - P^\text{th,min}\right), \quad \forall t\in \{1, \dots, T\} \\ +& \text{on}^\text{th}_t Q^\text{th,min} \le q^\text{th}_t \le \text{on}^\text{th}_t Q^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& -R^\text{th,dn} \le \Delta p_1^\text{th} - \Delta p^\text{th, init} \le R^\text{th,up} \\ +& -R^\text{th,dn} \le \Delta p_t^\text{th} - \Delta p_{t-1}^\text{th} \le R^\text{th,up}, \quad \forall t\in \{2, \dots, T\} +\end{align*} +``` --- + ## `ThermalStandardDispatch` ```@docs ThermalStandardDispatch ``` -TODO +**Variables:** + +- [`ActivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``p^\text{th}`` +- [`ReactivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``q^\text{th}`` + +**Static Parameters:** + +- ``P^\text{th,min}`` = `PowerSystems.get_active_power_limits(device).min` +- ``P^\text{th,max}`` = `PowerSystems.get_active_power_limits(device).max` +- ``Q^\text{th,min}`` = `PowerSystems.get_reactive_power_limits(device).min` +- ``Q^\text{th,max}`` = `PowerSystems.get_reactive_power_limits(device).max` +- ``R^\text{th,up}`` = `PowerSystems.get_ramp_limits(device).up` +- ``R^\text{th,dn}`` = `PowerSystems.get_ramp_limits(device).down` + +**Objective:** + +Add a cost to the objective function depending on the defined cost structure of the thermal unit by adding it to its `ProductionCostExpression`. + +**Expressions:** + +Adds ``p^\text{th}`` to the `ActivePowerBalance` expression and ``q^\text{th}`` to the `ReactivePowerBalance`, to be used in the supply-balance constraint depending on the network model used. + +**Constraints:** + +For each thermal unit creates the range constraints for its active and reactive power depending on its static parameters. + +```math +\begin{align*} +& P^\text{th,min} \le p^\text{th}_t \le P^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& Q^\text{th,min} \le q^\text{th}_t \le Q^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& -R^\text{th,dn} \le p_1^\text{th} - p^\text{th, init} \le R^\text{th,up} \\ +& -R^\text{th,dn} \le p_t^\text{th} - p_{t-1}^\text{th} \le R^\text{th,up}, \quad \forall t\in \{2, \dots, T\} +\end{align*} +``` + +--- + +## `ThermalBasicUnitCommitment` + +```@docs +ThermalBasicUnitCommitment +``` + +**Variables:** + +- [`ActivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``p^\text{th}`` +- [`ReactivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``q^\text{th}`` +- [`OnVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``u_t^\text{th}`` +- [`StartVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``v_t^\text{th}`` +- [`StopVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``w_t^\text{th}`` + + +**Static Parameters:** + +- ``P^\text{th,min}`` = `PowerSystems.get_active_power_limits(device).min` +- ``P^\text{th,max}`` = `PowerSystems.get_active_power_limits(device).max` +- ``Q^\text{th,min}`` = `PowerSystems.get_reactive_power_limits(device).min` +- ``Q^\text{th,max}`` = `PowerSystems.get_reactive_power_limits(device).max` + + +**Objective:** + +Add a cost to the objective function depending on the defined cost structure of the thermal unit by adding it to its `ProductionCostExpression`. + +**Expressions:** + +Adds ``p^\text{th}`` to the `ActivePowerBalance` expression and ``q^\text{th}`` to the `ReactivePowerBalance`, to be used in the supply-balance constraint depending on the network model used. + +**Constraints:** + +For each thermal unit creates the range constraints for its active and reactive power depending on its static parameters. In addition, it creates the commitment constraint to turn on/off the device. + +```math +\begin{align*} +& u_t^\text{th} P^\text{th,min} \le p^\text{th}_t \le u_t^\text{th} P^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& u_t^\text{th} Q^\text{th,min} \le q^\text{th}_t \le u_t^\text{th} Q^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& u_1^\text{th} = u^\text{th,init} + v_1^\text{th} - w_1^\text{th} \\ +& u_t^\text{th} = u_{t-1}^\text{th} + v_t^\text{th} - w_t^\text{th}, \quad \forall t \in \{2,\dots,T\} \\ +& v_t^\text{th} + w_t^\text{th} \le 1, \quad \forall t \in \{1,\dots,T\} +\end{align*} +``` --- @@ -65,7 +279,60 @@ TODO ThermalBasicCompactUnitCommitment ``` -TODO + +**Variables:** + +- [`PowerAboveMinimumVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``\Delta p^\text{th}`` +- [`ReactivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``q^\text{th}`` +- [`OnVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``u_t^\text{th}`` +- [`StartVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``v_t^\text{th}`` +- [`StopVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``w_t^\text{th}`` + +**Auxiliary Variables:** +- [`PowerOutput`](@ref): + - Symbol: ``P^\text{th}`` + - Definition: ``P^\text{th} = u^\text{th}P^\text{min} + \Delta p^\text{th}`` + + +**Static Parameters:** + +- ``P^\text{th,min}`` = `PowerSystems.get_active_power_limits(device).min` +- ``P^\text{th,max}`` = `PowerSystems.get_active_power_limits(device).max` +- ``Q^\text{th,min}`` = `PowerSystems.get_reactive_power_limits(device).min` +- ``Q^\text{th,max}`` = `PowerSystems.get_reactive_power_limits(device).max` + + +**Objective:** + +Add a cost to the objective function depending on the defined cost structure of the thermal unit by adding it to its `ProductionCostExpression`. + +**Expressions:** + +Adds ``u^\text{th}P^\text{th,min} + \Delta p^\text{th}`` to the `ActivePowerBalance` expression and ``q^\text{th}`` to the `ReactivePowerBalance`, to be used in the supply-balance constraint depending on the network model used. + +**Constraints:** + +For each thermal unit creates the range constraints for its active and reactive power depending on its static parameters. In addition, it creates the commitment constraint to turn on/off the device. + +```math +\begin{align*} +& 0 \le \Delta p^\text{th}_t \le u^\text{th}_t\left(P^\text{th,max} - P^\text{th,min}\right), \quad \forall t\in \{1, \dots, T\} \\ +& u_t^\text{th} Q^\text{th,min} \le q^\text{th}_t \le u_t^\text{th} Q^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& u_1^\text{th} = u^\text{th,init} + v_1^\text{th} - w_1^\text{th} \\ +& u_t^\text{th} = u_{t-1}^\text{th} + v_t^\text{th} - w_t^\text{th}, \quad \forall t \in \{2,\dots,T\} \\ +& v_t^\text{th} + w_t^\text{th} \le 1, \quad \forall t \in \{1,\dots,T\} +\end{align*} +``` --- @@ -75,36 +342,335 @@ TODO ThermalCompactUnitCommitment ``` -TODO +**Variables:** + +- [`PowerAboveMinimumVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``\Delta p^\text{th}`` +- [`ReactivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``q^\text{th}`` +- [`OnVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``u_t^\text{th}`` +- [`StartVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``v_t^\text{th}`` +- [`StopVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``w_t^\text{th}`` + +**Auxiliary Variables:** +- [`PowerOutput`](@ref): + - Symbol: ``P^\text{th}`` + - Definition: ``P^\text{th} = u^\text{th}P^\text{min} + \Delta p^\text{th}`` +- [`TimeDurationOn`](@ref): + - Symbol: ``V_t^\text{th}`` + - Definition: Computed post optimization by adding consecutive turned on variable ``u_t^\text{th}`` +- [`TimeDurationOff`](@ref): + - Symbol: ``W_t^\text{th}`` + - Definition: Computed post optimization by adding consecutive turned off variable ``1 - u_t^\text{th}`` + +**Static Parameters:** + +- ``P^\text{th,min}`` = `PowerSystems.get_active_power_limits(device).min` +- ``P^\text{th,max}`` = `PowerSystems.get_active_power_limits(device).max` +- ``Q^\text{th,min}`` = `PowerSystems.get_reactive_power_limits(device).min` +- ``Q^\text{th,max}`` = `PowerSystems.get_reactive_power_limits(device).max` +- ``R^\text{th,up}`` = `PowerSystems.get_ramp_limits(device).up` +- ``R^\text{th,dn}`` = `PowerSystems.get_ramp_limits(device).down` +- ``D^\text{min,up}`` = `PowerSystems.get_time_limits(device).up` +- ``D^\text{min,dn}`` = `PowerSystems.get_time_limits(device).down` + + +**Objective:** + +Add a cost to the objective function depending on the defined cost structure of the thermal unit by adding it to its `ProductionCostExpression`. + +**Expressions:** + +Adds ``u^\text{th}P^\text{th,min} + \Delta p^\text{th}`` to the `ActivePowerBalance` expression and ``q^\text{th}`` to the `ReactivePowerBalance`, to be used in the supply-balance constraint depending on the network model used. + +**Constraints:** + +For each thermal unit creates the range constraints for its active and reactive power depending on its static parameters. It also creates the commitment constraint to turn on/off the device. + +```math +\begin{align*} +& 0 \le \Delta p^\text{th}_t \le u^\text{th}_t\left(P^\text{th,max} - P^\text{th,min}\right), \quad \forall t\in \{1, \dots, T\} \\ +& u_t^\text{th} Q^\text{th,min} \le q^\text{th}_t \le u_t^\text{th} Q^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& -R^\text{th,dn} \le \Delta p_1^\text{th} - \Delta p^\text{th, init} \le R^\text{th,up} \\ +& -R^\text{th,dn} \le \Delta p_t^\text{th} - \Delta p_{t-1}^\text{th} \le R^\text{th,up}, \quad \forall t\in \{2, \dots, T\} \\ +& u_1^\text{th} = u^\text{th,init} + v_1^\text{th} - w_1^\text{th} \\ +& u_t^\text{th} = u_{t-1}^\text{th} + v_t^\text{th} - w_t^\text{th}, \quad \forall t \in \{2,\dots,T\} \\ +& v_t^\text{th} + w_t^\text{th} \le 1, \quad \forall t \in \{1,\dots,T\} +\end{align*} +``` ---- +In addition, this formulation adds duration constraints, i.e. minimum-up time and minimum-down time constraints. The duration constraints are added over the start times looking backwards. -## `ThermalMultiStartUnitCommitment` +The duration times ``D^\text{min,up}`` and ``D^\text{min,dn}`` are processed to be used in multiple of the time-steps, given the resolution of the specific problem. In addition, parameters ``D^\text{init,up}`` and ``D^\text{init,dn}`` are used to identify how long the unit was on or off, respectively, before the simulation started. -```@docs -ThermalMultiStartUnitCommitment +Minimum up-time constraint for ``t \in \{1,\dots T\}``: +```math +\begin{align*} +& \text{If } t \leq D^\text{min,up} - D^\text{init,up} \text{ and } D^\text{init,up} > 0: \\ +& 1 + \sum_{i=t-D^\text{min,up} + 1}^t v_i^\text{th} \leq u_t^\text{th} \quad \text{(for } i \text{ in the set of time steps).} \\ +& \text{Otherwise:} \\ +& \sum_{i=t-D^\text{min,up} + 1}^t v_i^\text{th} \leq u_t^\text{th} +\end{align*} ``` -TODO +Minimum down-time constraint for ``t \in \{1,\dots T\}``: +```math +\begin{align*} +& \text{If } t \leq D^\text{min,dn} - D^\text{init,dn} \text{ and } D^\text{init,up} > 0: \\ +& 1 + \sum_{i=t-D^\text{min,dn} + 1}^t w_i^\text{th} \leq 1 - u_t^\text{th} \quad \text{(for } i \text{ in the set of time steps).} \\ +& \text{Otherwise:} \\ +& \sum_{i=t-D^\text{min,dn} + 1}^t w_i^\text{th} \leq 1 - u_t^\text{th} +\end{align*} +``` --- -## `ThermalBasicUnitCommitment` +## `ThermalStandardUnitCommitment` ```@docs -ThermalBasicUnitCommitment +ThermalStandardUnitCommitment +``` + +**Variables:** + +- [`ActivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``p^\text{th}`` +- [`ReactivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``q^\text{th}`` +- [`OnVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``u_t^\text{th}`` +- [`StartVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``v_t^\text{th}`` +- [`StopVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``w_t^\text{th}`` + +**Auxiliary Variables:** +- [`TimeDurationOn`](@ref): + - Symbol: ``V_t^\text{th}`` + - Definition: Computed post optimization by adding consecutive turned on variable ``u_t^\text{th}`` +- [`TimeDurationOff`](@ref): + - Symbol: ``W_t^\text{th}`` + - Definition: Computed post optimization by adding consecutive turned off variable ``1 - u_t^\text{th}`` + +**Static Parameters:** + +- ``P^\text{th,min}`` = `PowerSystems.get_active_power_limits(device).min` +- ``P^\text{th,max}`` = `PowerSystems.get_active_power_limits(device).max` +- ``Q^\text{th,min}`` = `PowerSystems.get_reactive_power_limits(device).min` +- ``Q^\text{th,max}`` = `PowerSystems.get_reactive_power_limits(device).max` +- ``R^\text{th,up}`` = `PowerSystems.get_ramp_limits(device).up` +- ``R^\text{th,dn}`` = `PowerSystems.get_ramp_limits(device).down` +- ``D^\text{min,up}`` = `PowerSystems.get_time_limits(device).up` +- ``D^\text{min,dn}`` = `PowerSystems.get_time_limits(device).down` + + +**Objective:** + +Add a cost to the objective function depending on the defined cost structure of the thermal unit by adding it to its `ProductionCostExpression`. + +**Expressions:** + +Adds ``p^\text{th}`` to the `ActivePowerBalance` expression and ``q^\text{th}`` to the `ReactivePowerBalance`, to be used in the supply-balance constraint depending on the network model used. + +**Constraints:** + +For each thermal unit creates the range constraints for its active and reactive power depending on its static parameters. It also creates the commitment constraint to turn on/off the device. + +```math +\begin{align*} +& u^\text{th}_t P^\text{th,min} \le p^\text{th}_t \le u^\text{th}_t P^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& u_t^\text{th} Q^\text{th,min} \le q^\text{th}_t \le u_t^\text{th} Q^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& -R^\text{th,dn} \le p_1^\text{th} - p^\text{th, init} \le R^\text{th,up} \\ +& -R^\text{th,dn} \le p_t^\text{th} - p_{t-1}^\text{th} \le R^\text{th,up}, \quad \forall t\in \{2, \dots, T\} \\ +& u_1^\text{th} = u^\text{th,init} + v_1^\text{th} - w_1^\text{th} \\ +& u_t^\text{th} = u_{t-1}^\text{th} + v_t^\text{th} - w_t^\text{th}, \quad \forall t \in \{2,\dots,T\} \\ +& v_t^\text{th} + w_t^\text{th} \le 1, \quad \forall t \in \{1,\dots,T\} +\end{align*} +``` + +In addition, this formulation adds duration constraints, i.e. minimum-up time and minimum-down time constraints. The duration constraints are added over the start times looking backwards. + +The duration times ``D^\text{min,up}`` and ``D^\text{min,dn}`` are processed to be used in multiple of the time-steps, given the resolution of the specific problem. In addition, parameters ``D^\text{init,up}`` and ``D^\text{init,dn}`` are used to identify how long the unit was on or off, respectively, before the simulation started. + +Minimum up-time constraint for ``t \in \{1,\dots T\}``: +```math +\begin{align*} +& \text{If } t \leq D^\text{min,up} - D^\text{init,up} \text{ and } D^\text{init,up} > 0: \\ +& 1 + \sum_{i=t-D^\text{min,up} + 1}^t v_i^\text{th} \leq u_t^\text{th} \quad \text{(for } i \text{ in the set of time steps).} \\ +& \text{Otherwise:} \\ +& \sum_{i=t-D^\text{min,up} + 1}^t v_i^\text{th} \leq u_t^\text{th} +\end{align*} +``` + +Minimum down-time constraint for ``t \in \{1,\dots T\}``: +```math +\begin{align*} +& \text{If } t \leq D^\text{min,dn} - D^\text{init,dn} \text{ and } D^\text{init,up} > 0: \\ +& 1 + \sum_{i=t-D^\text{min,dn} + 1}^t w_i^\text{th} \leq 1 - u_t^\text{th} \quad \text{(for } i \text{ in the set of time steps).} \\ +& \text{Otherwise:} \\ +& \sum_{i=t-D^\text{min,dn} + 1}^t w_i^\text{th} \leq 1 - u_t^\text{th} +\end{align*} ``` -TODO --- -## `ThermalStandardUnitCommitment` +## `ThermalMultiStartUnitCommitment` ```@docs -ThermalStandardUnitCommitment +ThermalMultiStartUnitCommitment +``` + + +**Variables:** + +- [`PowerAboveMinimumVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``\Delta p^\text{th}`` +- [`ReactivePowerVariable`](@ref): + - Bounds: [0.0, ] + - Symbol: ``q^\text{th}`` +- [`OnVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``u_t^\text{th}`` +- [`StartVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``v_t^\text{th}`` +- [`StopVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``w_t^\text{th}`` +- [`ColdStartVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``x_t^\text{th}`` +- [`WarmStartVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``y_t^\text{th}`` +- [`HotStartVariable`](@ref): + - Bounds: ``\{0,1\}`` + - Symbol: ``z_t^\text{th}`` + +**Auxiliary Variables:** +- [`PowerOutput`](@ref): + - Symbol: ``P^\text{th}`` + - Definition: ``P^\text{th} = u^\text{th}P^\text{min} + \Delta p^\text{th}`` +- [`TimeDurationOn`](@ref): + - Symbol: ``V_t^\text{th}`` + - Definition: Computed post optimization by adding consecutive turned on variable ``u_t^\text{th}`` +- [`TimeDurationOff`](@ref): + - Symbol: ``W_t^\text{th}`` + - Definition: Computed post optimization by adding consecutive turned off variable ``1 - u_t^\text{th}`` + +**Static Parameters:** + +- ``P^\text{th,min}`` = `PowerSystems.get_active_power_limits(device).min` +- ``P^\text{th,max}`` = `PowerSystems.get_active_power_limits(device).max` +- ``Q^\text{th,min}`` = `PowerSystems.get_reactive_power_limits(device).min` +- ``Q^\text{th,max}`` = `PowerSystems.get_reactive_power_limits(device).max` +- ``R^\text{th,up}`` = `PowerSystems.get_ramp_limits(device).up` +- ``R^\text{th,dn}`` = `PowerSystems.get_ramp_limits(device).down` +- ``D^\text{min,up}`` = `PowerSystems.get_time_limits(device).up` +- ``D^\text{min,dn}`` = `PowerSystems.get_time_limits(device).down` +- ``D^\text{cold}`` = `PowerSystems.get_start_time_limits(device).cold` +- ``D^\text{warm}`` = `PowerSystems.get_start_time_limits(device).warm` +- ``D^\text{hot}`` = `PowerSystems.get_start_time_limits(device).hot` +- ``P^\text{th,startup}`` = `PowerSystems.get_power_trajectory(device).startup` +- ``P^\text{th, shdown}`` = `PowerSystems.get_power_trajectory(device).shutdown` + + +**Objective:** + +Add a cost to the objective function depending on the defined cost structure of the thermal unit by adding it to its `ProductionCostExpression`. + +**Expressions:** + +Adds ``u^\text{th}P^\text{th,min} + \Delta p^\text{th}`` to the `ActivePowerBalance` expression and ``q^\text{th}`` to the `ReactivePowerBalance`, to be used in the supply-balance constraint depending on the network model used. + +**Constraints:** + +For each thermal unit creates the range constraints for its active and reactive power depending on its static parameters. It also creates the commitment constraint to turn on/off the device. + +```math +\begin{align*} +& 0 \le \Delta p^\text{th}_t \le u^\text{th}_t\left(P^\text{th,max} - P^\text{th,min}\right), \quad \forall t\in \{1, \dots, T\} \\ +& u_t^\text{th} Q^\text{th,min} \le q^\text{th}_t \le u_t^\text{th} Q^\text{th,max}, \quad \forall t\in \{1, \dots, T\} \\ +& -R^\text{th,dn} \le \Delta p_1^\text{th} - \Delta p^\text{th, init} \le R^\text{th,up} \\ +& -R^\text{th,dn} \le \Delta p_t^\text{th} - \Delta p_{t-1}^\text{th} \le R^\text{th,up}, \quad \forall t\in \{2, \dots, T\} \\ +& u_1^\text{th} = u^\text{th,init} + v_1^\text{th} - w_1^\text{th} \\ +& u_t^\text{th} = u_{t-1}^\text{th} + v_t^\text{th} - w_t^\text{th}, \quad \forall t \in \{2,\dots,T\} \\ +& v_t^\text{th} + w_t^\text{th} \le 1, \quad \forall t \in \{1,\dots,T\} \\ +& \max\{P^\text{th,max} - P^\text{th,shdown}, 0\} \cdot w_1^\text{th} \le u^\text{th,init} (P^\text{th,max} - P^\text{th,min}) - P^\text{th,init} +\end{align*} +``` + +In addition, this formulation adds duration constraints, i.e. minimum-up time and minimum-down time constraints. The duration constraints are added over the start times looking backwards. + +The duration times ``D^\text{min,up}`` and ``D^\text{min,dn}`` are processed to be used in multiple of the time-steps, given the resolution of the specific problem. In addition, parameters ``D^\text{init,up}`` and ``D^\text{init,dn}`` are used to identify how long the unit was on or off, respectively, before the simulation started. + +Minimum up-time constraint for ``t \in \{1,\dots T\}``: +```math +\begin{align*} +& \text{If } t \leq D^\text{min,up} - D^\text{init,up} \text{ and } D^\text{init,up} > 0: \\ +& 1 + \sum_{i=t-D^\text{min,up} + 1}^t v_i^\text{th} \leq u_t^\text{th} \quad \text{(for } i \text{ in the set of time steps).} \\ +& \text{Otherwise:} \\ +& \sum_{i=t-D^\text{min,up} + 1}^t v_i^\text{th} \leq u_t^\text{th} +\end{align*} +``` + +Minimum down-time constraint for ``t \in \{1,\dots T\}``: +```math +\begin{align*} +& \text{If } t \leq D^\text{min,dn} - D^\text{init,dn} \text{ and } D^\text{init,up} > 0: \\ +& 1 + \sum_{i=t-D^\text{min,dn} + 1}^t w_i^\text{th} \leq 1 - u_t^\text{th} \quad \text{(for } i \text{ in the set of time steps).} \\ +& \text{Otherwise:} \\ +& \sum_{i=t-D^\text{min,dn} + 1}^t w_i^\text{th} \leq 1 - u_t^\text{th} +\end{align*} +``` + +Finally, multi temperature start/stop constraints are implemented using the following constraints: + +```math +\begin{align*} +& v_t^\text{th} = x_t^\text{th} + y_t^\text{th} + z_t^\text{th}, \quad \forall t \in \{1, \dots, T\} \\ +& z_t^\text{th} \le \sum_{i \in [D^\text{hot}, D^\text{warm})}w_{t-i}^\text{th}, \quad \forall t \in \{D^\text{warm}, \dots, T\} \\ +& y_t^\text{th} \le \sum_{i \in [D^\text{warm}, D^\text{cold})}w_{t-i}^\text{th}, \quad \forall t \in \{D^\text{cold}, \dots, T\} \\ +& (D^\text{warm} - 1) z_t^\text{th} + (1 - z_t^\text{th}) M^\text{big} \ge \sum_{i=1}^t (1 - u_i^\text{th}) + D^\text{init,hot}, \quad \forall t \in \{1, \dots, T\} \\ +& D^\text{hot} z_t^\text{th} \le \sum_{i=1}^t (1 - u_i^\text{th}) + D^\text{init,hot}, \quad \forall t \in \{1, \dots, T\} \\ +& (D^\text{cold} - 1) y_t^\text{th} + (1 - y_t^\text{th}) M^\text{big} \ge \sum_{i=1}^t (1 - u_i^\text{th}) + D^\text{init,warm}, \quad \forall t \in \{1, \dots, T\} \\ +& D^\text{warm} y_t^\text{th} \le \sum_{i=1}^t (1 - u_i^\text{th}) + D^\text{init,warm}, \quad \forall t \in \{1, \dots, T\} \\ +\end{align*} ``` -TODO --- + +## Valid configurations + +Valid `DeviceModel`s for subtypes of `ThermalGen` include the following: + +```@eval +using PowerSimulations +using PowerSystems +using DataFrames +using Latexify +combos = PowerSimulations.generate_device_formulation_combinations() +filter!(x -> x["device_type"] <: ThermalGen, combos) +combo_table = DataFrame( + "Valid DeviceModel" => ["`DeviceModel($(c["device_type"]), $(c["formulation"]))`" for c in combos], + "Device Type" => ["[$(c["device_type"])](https://nrel-Sienna.github.io/PowerSystems.jl/stable/model_library/generated_$(c["device_type"])/)" for c in combos], + "Formulation" => ["[$(c["formulation"])](@ref)" for c in combos], + ) +mdtable(combo_table, latex = false) +``` diff --git a/docs/src/modeler_guide/debugging_infeasible_models.md b/docs/src/modeler_guide/debugging_infeasible_models.md index c96c7ff78b..cc52f7a2ec 100644 --- a/docs/src/modeler_guide/debugging_infeasible_models.md +++ b/docs/src/modeler_guide/debugging_infeasible_models.md @@ -5,4 +5,168 @@ Getting infeasible solutions to models is a common occurrence in operations simu ## Adding slacks to the model +One of the most common infeasibility issues observed is due to not enough generation to supply demand, or conversely, excessive fixed (non-curtailable) generation in a low demand scenario. + +The recommended solution for any of these cases is adding slack variables to the network model, for example: + +```julia +template_uc = ProblemTemplate( + NetworkModel( + CopperPlatePowerModel, + use_slacks=true, + ), + ) +``` +will add slack variables to the `ActivePowerBalance` expression. + +In this case, if the problem is now feasible, the user can check the solution of the variables `SystemBalanceSlackUp` and `SystemBalanceSlackDown`, and if one value is greater than zero, it represents that not enough generation (for Slack Up) or not enough demand (for Slack Down) in the optimization problem. + +### Services cases + +In many scenarios, certain units are also required to provide reserve requirements, e.g. thermal units mandated to provide up-regulation. In such scenarios, it is also possible to add slack variables, by specifying the service model (`RangeReserve`) for the specific service type (`VariableReserve{ReserveUp}`) as: +```julia +set_service_model!( + template_uc, + ServiceModel( + VariableReserve{ReserveUp}, + RangeReserve; + use_slacks=true + ), +) +``` +Again, if the problem is now feasible, check the solution of `ReserveRequirementSlack` variable, and if it is larger than zero in a specific time-step, then it is evidence that there is not enough reserve available to satisfy the requirement. + ## Getting the infeasibility conflict + +Some solvers allows to identify which constraints and variables are producing the infeasibility, by finding the irreducible infeasible set (IIS), that is the subset of constraints and variable bounds that will become feasible if any single constraint or variable bound is removed. + +To enable this feature in `PowerSimulations` the keyword argument `calculate_conflict` must be set to `true`, when creating the `DecisionModel`. Note that not all solvers allow the computation of the IIS, but most commercial solvers have this capability. It is also recommended to enable the keyword argument `store_variable_names=true` to help understanding which variables are with infeasibility issues. + +The following code creates a decision model with the `Xpress` optimizer, and enabling the `calculate_conflict=true` keyword argument. + +```julia +DecisionModel( + template_ed, + sys_rts_rt; + name="ED", + optimizer=optimizer_with_attributes(Xpress.Optimizer, "MIPRELSTOP" => 1e-2), + optimizer_solve_log_print=true, + calculate_conflict=true, + store_variable_names=true, +) +``` + +Here is an example on how the IIS will be displayed as: + +```raw +Error: Constraints participating in conflict basis (IIS) +│ +│ ┌──────────────────────────────────────┐ +│ │ CopperPlateBalanceConstraint__System │ +│ ├──────────────────────────────────────┤ +│ │ (113, 26) │ +│ └──────────────────────────────────────┘ +│ ┌──────────────────────────────────┐ +│ │ EnergyAssetBalance__HybridSystem │ +│ ├──────────────────────────────────┤ +│ │ ("317_Hybrid", 26) │ +│ └──────────────────────────────────┘ +│ ┌─────────────────────────────────────────────┐ +│ │ PieceWiseLinearCostConstraint__HybridSystem │ +│ ├─────────────────────────────────────────────┤ +│ │ ("317_Hybrid", 26) │ +│ └─────────────────────────────────────────────┘ +│ ┌────────────────────────────────────────────────┐ +│ │ PieceWiseLinearCostConstraint__ThermalStandard │ +│ ├────────────────────────────────────────────────┤ +│ │ ("202_STEAM_3", 26) │ +│ │ ("101_STEAM_3", 26) │ +│ │ ("118_CC_1", 26) │ +│ │ ("202_STEAM_4", 26) │ +│ │ ("315_CT_6", 26) │ +│ │ ("201_STEAM_3", 26) │ +│ │ ("102_STEAM_4", 26) │ +│ └────────────────────────────────────────────────┘ +│ ┌──────────────────────────────────────────────────────────────────────┐ +│ │ ActivePowerVariableTimeSeriesLimitsConstraint__RenewableDispatch__ub │ +│ ├──────────────────────────────────────────────────────────────────────┤ +│ │ ("122_WIND_1", 26) │ +│ │ ("324_PV_3", 26) │ +│ │ ("312_PV_1", 26) │ +│ │ ("102_PV_1", 26) │ +│ │ ("101_PV_1", 26) │ +│ │ ("324_PV_2", 26) │ +│ │ ("313_PV_2", 26) │ +│ │ ("104_PV_1", 26) │ +│ │ ("101_PV_2", 26) │ +│ │ ("309_WIND_1", 26) │ +│ │ ("310_PV_2", 26) │ +│ │ ("113_PV_1", 26) │ +│ │ ("314_PV_1", 26) │ +│ │ ("324_PV_1", 26) │ +│ │ ("103_PV_1", 26) │ +│ │ ("303_WIND_1", 26) │ +│ │ ("314_PV_2", 26) │ +│ │ ("102_PV_2", 26) │ +│ │ ("314_PV_3", 26) │ +│ │ ("320_PV_1", 26) │ +│ │ ("101_PV_3", 26) │ +│ │ ("319_PV_1", 26) │ +│ │ ("314_PV_4", 26) │ +│ │ ("310_PV_1", 26) │ +│ │ ("215_PV_1", 26) │ +│ │ ("313_PV_1", 26) │ +│ │ ("101_PV_4", 26) │ +│ │ ("119_PV_1", 26) │ +│ └──────────────────────────────────────────────────────────────────────┘ +│ ┌─────────────────────────────────────────────────────────────────────────────┐ +│ │ FeedforwardSemiContinuousConstraint__ThermalStandard__ActivePowerVariable_ub │ +│ ├─────────────────────────────────────────────────────────────────────────────┤ +│ │ ("322_CT_6", 26) │ +│ │ ("321_CC_1", 26) │ +│ │ ("223_CT_4", 26) │ +│ │ ("213_CT_1", 26) │ +│ │ ("223_CT_6", 26) │ +│ │ ("123_CT_1", 26) │ +│ │ ("113_CT_3", 26) │ +│ │ ("302_CT_3", 26) │ +│ │ ("215_CT_4", 26) │ +│ │ ("301_CT_4", 26) │ +│ │ ("113_CT_2", 26) │ +│ │ ("221_CC_1", 26) │ +│ │ ("223_CT_5", 26) │ +│ │ ("315_CT_7", 26) │ +│ │ ("215_CT_5", 26) │ +│ │ ("113_CT_1", 26) │ +│ │ ("307_CT_2", 26) │ +│ │ ("213_CT_2", 26) │ +│ │ ("113_CT_4", 26) │ +│ │ ("218_CC_1", 26) │ +│ │ ("213_CC_3", 26) │ +│ │ ("323_CC_2", 26) │ +│ │ ("322_CT_5", 26) │ +│ │ ("207_CT_2", 26) │ +│ │ ("123_CT_5", 26) │ +│ │ ("123_CT_4", 26) │ +│ │ ("207_CT_1", 26) │ +│ │ ("301_CT_3", 26) │ +│ │ ("302_CT_4", 26) │ +│ │ ("307_CT_1", 26) │ +│ └─────────────────────────────────────────────────────────────────────────────┘ +│ ┌───────────────────────────────────────────────────────┐ +│ │ RenewableActivePowerLimitConstraint__HybridSystem__ub │ +│ ├───────────────────────────────────────────────────────┤ +│ │ ("317_Hybrid", 26) │ +│ └───────────────────────────────────────────────────────┘ +│ ┌───────────────────────────────────────┐ +│ │ ThermalOnVariableUb__HybridSystem__ub │ +│ ├───────────────────────────────────────┤ +│ │ ("317_Hybrid", 26) │ +│ └───────────────────────────────────────┘ + + Error: Serializing Infeasible Problem at /var/folders/1v/t69qyl0n5059n6c1nn7sp8zm7g8s6z/T/jl_jNSREb/compact_sim/problems/ED/infeasible_ED_2020-10-06T15:00:00.json +``` + +Note that the IIS clearly identify that the issue is happening at time step 26, and constraints are related with the `CopperPlateBalanceConstraint__System`, with multiple upper bound constraints, for the hybrid system, renewable units and thermal units. This highlights that there may not be enough generation in the system. Indeed, by enabling system slacks, the problem become feasible. + +Finally, the infeasible model is exported in a `json` file that can be loaded directly in `JuMP` to be explored. More information about this is [available here](https://jump.dev/JuMP.jl/stable/moi/submodules/FileFormats/overview/#Read-from-file). \ No newline at end of file diff --git a/docs/src/modeler_guide/definitions.md b/docs/src/modeler_guide/definitions.md index 8a4946f85f..a8effa9ab5 100644 --- a/docs/src/modeler_guide/definitions.md +++ b/docs/src/modeler_guide/definitions.md @@ -1,15 +1,44 @@ # Definitions +## A + +* *Attributes*: Certain device formulations can be customized by specifying attributes that will include/remove certain variables, expressions and/or constraints. For example, in `StorageSystemsSimulations.jl`, the device formulation of `StorageDispatchWithReserves` can be specified with the following dictionary of attributes: +```julia +set_device_model!( + template, + DeviceModel( + GenericBattery, + StorageDispatchWithReserves; + attributes=Dict{String, Any}( + "reservation" => false, + "cycling_limits" => false, + "energy_target" => false, + "complete_coverage" => false, + "regularization" => false, + ), + ), +) +``` +Changing the attributes between `true` or `false` can enable/disable multiple aspects of the formulation. + +## C + +* *Chronologies:* In `PowerSimulations.jl`, chronologies define where information is flowing. There are two types of chronologies. 1) **inter-stage chronologies** (`InterProblemChronology`) that define how information flows between stages. e.g. day-ahead solutions are used to inform economic dispatch problems; and 2) **intra-stage chronologies** (`IntraProblemChronology`) that define how information flows between multiple executions of a single stage. e.g. the dispatch setpoints of the first period of an economic dispatch problem are constrained by the ramping limits from setpoints in the final period of the previous problem. + ## D -* *Decision Problem*: A decision problem calculates the desired system operation based on forecasts of uncertain inputs and information about the state of the system. The output of a decision problem represents the policies used to drive the set-points of the system's devices, like generators or switches, and depends on the purpose of the problem. See the [Decision Model Tutorial](op_problem_tutorial) to learn more about solving individual problems. +* *Decision Problem*: A decision problem calculates the desired system operation based on forecasts of uncertain inputs and information about the state of the system. The output of a decision problem represents the policies used to drive the set-points of the system's devices, like generators or switches, and depends on the purpose of the problem. See the [Decision Model Tutorial](@ref op_problem_tutorial) to learn more about solving individual problems. -* *Device Formulation*: The model of a device that is incorporated into a large system optimization models. For instance, the storage device model used inside of a Unit Commitment (UC) problem. A device model needs to follow some requirements to be integrated into operation problems. +* *Device Formulation*: The model of a device that is incorporated into a large system optimization models. For instance, the storage device model used inside of a Unit Commitment (UC) problem. A device model needs to follow some requirements to be integrated into operation problems. For more information about valid `DeviceModel`s and their mathematical representations, check out the [Formulation Library](@ref formulation_intro). ## E * *Emulation Problem*: An emulation problem is used to mimic the system's behavior subject to an incoming decision and the realization of a forecasted inputs. The solution of the emulator produces outputs representative of the system performance when operating subject the policies resulting from the decision models. +## F + +* *FeedForward*: The definition of exactly what information is passed using the defined chronologies is accomplished using FeedForwards. Specifically, a FeedForward is used to define what to do with information being passed with an inter-stage chronology in a Simulation. The most common FeedForward is the `SemiContinuousFeedForward` that affects the semi-continuous range constraints of thermal generators in the economic dispatch problems based on the value of the (already solved) unit-commitment variables. + ## H * *Horizon*: The number of steps in the look-ahead of a decision problem. For instance, a Day-Ahead problem usually has a 48 step horizon. Check the time [Time Series Data Section in PowerSystems.jl](https://nrel-sienna.github.io/PowerSystems.jl/stable/modeler_guide/time_series/) @@ -20,4 +49,18 @@ ## R -* *Resolution*: The amount of time between timesteps in a simulation. For instance 1-hour or 5-minutes. In Julia these are defined using the syntax `Hour(1)` and `Minute(5)`. Check the time [Time Series Data Section in PowerSystems.jl](https://nrel-sienna.github.io/PowerSystems.jl/stable/modeler_guide/time_series/) +* *Resolution*: The amount of time between time steps in a simulation. For instance 1-hour or 5-minutes. In Julia these are defined using the syntax `Hour(1)` and `Minute(5)`. Check the time [Time Series Data Section in PowerSystems.jl](https://nrel-sienna.github.io/PowerSystems.jl/stable/modeler_guide/time_series/) + +* *Results vs Realized Results*: In `PowerSimulations.jl` the term *results* is used to refer to the solution of all optimization problems in a *Simulation*. When using `read_variable(results, Variable)` in a `DecisionModel` of a simulation, the output is a dictionary with the values of such variable for every optimization problem solved, while `read_realized_variable(results, Variable)` will return the values of the specified interval and number of steps in the simulation. See the [Read Results page](@ref read_results) for more details. + +## S + +* *Service Formulation*: The model of a service that is incorporated into a large system optimization models. `Services` (or ancillary services) are models used to ensure that there is necessary support to the power grid from generators to consumers, in order to ensure reliable operation of the system. The most common application for ancillary services are reserves, i.e., generation (or load) that is not currently being used, but can be quickly made available in case of unexpected changes of grid conditions, for example a sudden loss of load or generation. A service model needs to follow some requirements to be integrated into operation problems. For more information about valid `ServiceModel`s and their mathematical representations, check out the [Formulation Library](@ref service_formulations). + +* *Simulation*: A simulation is a pre-determined sequence of decision problems in a way that solving it, resembles the solution procedures commonly used by operators. The most common simulation model is the solution of a Unit Commitment and Economic Dispatch sequence of problems. + +* *Solver*: A solver is a software package that incorporates algorithms for finding solutions to one or more classes of optimization problem. For example, FICO Xpress is a commercial optimization solver for linear programming (LP), convex quadratic programming (QP) problems, convex quadratically constrained quadratic programming (QCQP), second-order cone programming (SOCP) and their mixed integer counterparts. **A solver is required to be specified** in order to solve any computer optimization problem. + +## T + +* *Template*: A `ProblemTemplate` is just a collection of `DeviceModel`s that allows the user to specify the formulations of each set of devices (by device type) independently so that the modeler can adjust the level of detail according to the question of interest and the available data. For more information about valid `DeviceModel`s and their mathematical representations, check out the [Formulation Library](@ref formulation_intro). \ No newline at end of file diff --git a/docs/src/modeler_guide/problem_templates.md b/docs/src/modeler_guide/problem_templates.md index 52eff0a56a..aea9e98c44 100644 --- a/docs/src/modeler_guide/problem_templates.md +++ b/docs/src/modeler_guide/problem_templates.md @@ -3,7 +3,7 @@ Templates are used to specify the modeling properties of the devices and network that are going to he used to specify a problem. A `ProblemTemplate` is just a collection of `DeviceModel`s that allows the user to specify the formulations of each set of devices (by device type) independently so that the modeler can adjust the level of detail according to the question of interest and the available data. -For more information about valid `DeviceModel`s and their mathematical representations, check out the [Formulation Library](@ref formulation_library). +For more information about valid `DeviceModel`s and their mathematical representations, check out the [Formulation Library](@ref formulation_intro). ## Building a `ProblemTemplate` diff --git a/docs/src/modeler_guide/read_results.md b/docs/src/modeler_guide/read_results.md new file mode 100644 index 0000000000..6fcec568f9 --- /dev/null +++ b/docs/src/modeler_guide/read_results.md @@ -0,0 +1,201 @@ +# [Read results](@id read_results) + +Once a `DecisionModel` is solved via `solve!(model)` or a Simulation is executed (and solved) via `execute!(simulation)`, the results are stored and can be accessed directly in the REPL for result exploration and plotting. + +## Read results of a Decision Problem + +Once a `DecisionModel` is solved, results are accessed using `ProblemResults(model)` as follows: + +```julia +# The DecisionModel is already constructed +build!(model, output_dir = mktempdir()) +solve!(model) + +results = ProblemResults(model) +``` + +The output will showcase the available expressions, parameters and variables to read. For example it will look like: + +```raw +Start: 2020-01-01T00:00:00 +End: 2020-01-03T23:00:00 +Resolution: 60 minutes + +PowerSimulations Problem Auxiliary variables Results +┌──────────────────────────────────────────┐ +│ CumulativeCyclingCharge__HybridSystem │ +│ CumulativeCyclingDischarge__HybridSystem │ +└──────────────────────────────────────────┘ + +PowerSimulations Problem Expressions Results +┌─────────────────────────────────────────────┐ +│ ProductionCostExpression__RenewableDispatch │ +│ ProductionCostExpression__ThermalStandard │ +└─────────────────────────────────────────────┘ + +PowerSimulations Problem Duals Results +┌──────────────────────────────────────┐ +│ CopperPlateBalanceConstraint__System │ +└──────────────────────────────────────┘ + +PowerSimulations Problem Parameters Results +┌────────────────────────────────────────────────────────────────────────┐ +│ ActivePowerTimeSeriesParameter__RenewableFix │ +│ RenewablePowerTimeSeries__HybridSystem │ +│ RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Spin_Up_R3 │ +│ RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Reg_Up │ +│ ActivePowerTimeSeriesParameter__PowerLoad │ +│ ActivePowerTimeSeriesParameter__RenewableDispatch │ +│ RequirementTimeSeriesParameter__VariableReserve__ReserveDown__Reg_Down │ +│ ActivePowerTimeSeriesParameter__HydroDispatch │ +│ RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Spin_Up_R1 │ +│ RequirementTimeSeriesParameter__VariableReserve__ReserveUp__Spin_Up_R2 │ +└────────────────────────────────────────────────────────────────────────┘ + +PowerSimulations Problem Variables Results +┌────────────────────────────────────────────────────────────────────┐ +│ ActivePowerOutVariable__HybridSystem │ +│ ReservationVariable__HybridSystem │ +│ RenewablePower__HybridSystem │ +│ ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R1 │ +│ SystemBalanceSlackUp__System │ +│ BatteryEnergyShortageVariable__HybridSystem │ +│ ActivePowerReserveVariable__VariableReserve__ReserveUp__Reg_Up │ +│ StopVariable__ThermalStandard │ +│ BatteryStatus__HybridSystem │ +│ BatteryDischarge__HybridSystem │ +│ ActivePowerInVariable__HybridSystem │ +│ DischargeRegularizationVariable__HybridSystem │ +│ BatteryCharge__HybridSystem │ +│ ActivePowerVariable__RenewableDispatch │ +│ ActivePowerReserveVariable__VariableReserve__ReserveDown__Reg_Down │ +│ EnergyVariable__HybridSystem │ +│ OnVariable__HybridSystem │ +│ BatteryEnergySurplusVariable__HybridSystem │ +│ SystemBalanceSlackDown__System │ +│ ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R2 │ +│ ThermalPower__HybridSystem │ +│ ActivePowerVariable__ThermalStandard │ +│ StartVariable__ThermalStandard │ +│ ActivePowerReserveVariable__VariableReserve__ReserveUp__Spin_Up_R3 │ +│ OnVariable__ThermalStandard │ +│ ChargeRegularizationVariable__HybridSystem │ +└────────────────────────────────────────────────────────────────────┘ +``` + +Then the following code can be used to read results: + +```julia +# Read active power of Thermal Standard +thermal_active_power = read_variable(results, "ActivePowerVariable__ThermalStandard") + +# Read max active power parameter of RenewableDispatch +renewable_param = read_parameter(results, "ActivePowerTimeSeriesParameter__RenewableDispatch") + +# Read cost expressions of ThermalStandard units +cost_thermal = read_expression(results, "ProductionCostExpression__ThermalStandard") + +# Read dual variables +dual_balance_constraint = read_dual(results, "CopperPlateBalanceConstraint__System") + +# Read auxiliary variables +aux_var_result = read_aux_variable(results, "CumulativeCyclingCharge__HybridSystem") +``` + +Results will be in the form of DataFrames that can be easily explored. + +## Read results of a Simulation + +```julia +# The Simulation is already constructed +build!(sim) +execute!(sim; enable_progress_bar=true) + +results_sim = SimulationResults(sim) +``` + +As an example, the `SimulationResults` printing will look like: + +```raw +Decision Problem Results +┌──────────────┬─────────────────────┬──────────────┬─────────────────────────┐ +│ Problem Name │ Initial Time │ Resolution │ Last Solution Timestamp │ +├──────────────┼─────────────────────┼──────────────┼─────────────────────────┤ +│ ED │ 2020-10-02T00:00:00 │ 60 minutes │ 2020-10-09T23:00:00 │ +│ UC │ 2020-10-02T00:00:00 │ 1440 minutes │ 2020-10-09T00:00:00 │ +└──────────────┴─────────────────────┴──────────────┴─────────────────────────┘ + +Emulator Results +┌─────────────────┬───────────┐ +│ Name │ Emulator │ +│ Resolution │ 5 minutes │ +│ Number of steps │ 2304 │ +└─────────────────┴───────────┘ +``` + +With this, it is possible to obtain results of each `DecisionModel` and `EmulationModel` as follows: + +```julia +# Use the Problem Name for Decision Problems +results_uc = get_decision_problem_results(results_sim, "UC") +results_ed = get_decision_problem_results(results_sim, "ED") +results_emulator = get_emulation_problem_results(results_sim) +``` + +Once we have each decision (or emulation) problem results, we can explore directly using the approach for Decision Models, mentioned in the previous section. + +### Reading solutions for all simulation steps + +In this case, using `read_variable` (or read expression, parameter or dual), will return a dictionary of all steps (of that Decision Problem). For example, the following code: + +```julia +thermal_active_power = read_variable(results_uc, "ActivePowerVariable__ThermalStandard") +``` +will return: +``` +DataStructures.SortedDict{Any, Any, Base.Order.ForwardOrdering} with 8 entries: + DateTime("2020-10-02T00:00:00") => 72×54 DataFrame… + DateTime("2020-10-03T00:00:00") => 72×54 DataFrame… + DateTime("2020-10-04T00:00:00") => 72×54 DataFrame… + DateTime("2020-10-05T00:00:00") => 72×54 DataFrame… + DateTime("2020-10-06T00:00:00") => 72×54 DataFrame… + DateTime("2020-10-07T00:00:00") => 72×54 DataFrame… + DateTime("2020-10-08T00:00:00") => 72×54 DataFrame… + DateTime("2020-10-09T00:00:00") => 72×54 DataFrame… +``` +That is, a sorted dictionary for each simulation step, using as a key the initial timestamp for that specific simulation step. + +Note that in this case, each DataFrame, has a dimension of ``72 \times 54``, since the horizon is 72 hours (number of rows), but the interval is only 24 hours. Indeed, note the initial timestamp of each simulation step is the beginning of each day, i.e. 24 hours. Finally, there 54 columns, since this example system has 53 `ThermalStandard` units (plus 1 column for the timestamps). The user is free to explore the solution of any simulation step as needed. + +### Reading the "realized" solution (i.e. the interval) + +Using `read_realized_variable` (or read realized expression, parameter or dual), will return the DataFrame of the realized solution of any specific variable. That is, it will concatenate the corresponding simulation step with the specified interval of that step, to construct a single DataFrame with the "realized solution" of the entire simulation. + +For example, the code: +```julia +th_realized_power = read_realized_variable(results_uc, "ActivePowerVariable__ThermalStandard") +``` +will return: +```raw +92×54 DataFrame + Row │ DateTime 322_CT_6 321_CC_1 202_STEAM_3 223_CT_4 123_STEAM_2 213_CT_1 223_CT_6 313_CC_1 101_STEAM_3 123_C ⋯ + │ DateTime Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float ⋯ +─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── + 1 │ 2020-10-02T00:00:00 0.0 293.333 0.0 0.0 0.0 0.0 0.0 231.667 76.0 0.0 ⋯ + 2 │ 2020-10-02T01:00:00 0.0 267.552 0.0 0.0 0.0 0.0 0.0 231.667 76.0 0.0 + 3 │ 2020-10-02T02:00:00 0.0 234.255 0.0 0.0 -4.97544e-11 0.0 0.0 231.667 76.0 0.0 + 4 │ 2020-10-02T03:00:00 0.0 249.099 0.0 0.0 -4.97544e-11 0.0 0.0 231.667 76.0 0.0 + 5 │ 2020-10-02T04:00:00 0.0 293.333 0.0 0.0 -4.97544e-11 0.0 0.0 231.667 76.0 0.0 ⋯ + 6 │ 2020-10-02T05:00:00 0.0 293.333 1.27578e-11 0.0 -4.97544e-11 0.0 0.0 293.333 76.0 0.0 + ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ + 187 │ 2020-10-09T18:00:00 0.0 293.333 76.0 0.0 155.0 0.0 0.0 318.843 76.0 0.0 + 188 │ 2020-10-09T19:00:00 0.0 293.333 76.0 0.0 124.0 0.0 0.0 293.333 76.0 0.0 + 189 │ 2020-10-09T20:00:00 0.0 293.333 60.6667 0.0 124.0 0.0 0.0 0.0 76.0 0.0 ⋯ + 190 │ 2020-10-09T21:00:00 -7.65965e-12 293.333 60.6667 0.0 124.0 0.0 0.0 0.0 76.0 0.0 + 191 │ 2020-10-09T22:00:00 0.0 0.0 60.6667 0.0 124.0 0.0 0.0 0.0 76.0 7.156 + 192 │ 2020-10-09T23:00:00 0.0 0.0 60.6667 0.0 117.81 0.0 0.0 0.0 76.0 0.0 + 44 columns and 180 rows omitted +``` +In this case, the 8 simulation steps of 24 hours (192 hours), in a single DataFrame, to enable easy exploration of the realized results for the user. + + diff --git a/docs/src/modeler_guide/running_a_simulation.md b/docs/src/modeler_guide/running_a_simulation.md index b2174d9bad..80df408917 100644 --- a/docs/src/modeler_guide/running_a_simulation.md +++ b/docs/src/modeler_guide/running_a_simulation.md @@ -7,10 +7,112 @@ Check out the [Operations Problem Tutorial](@ref op_problem_tutorial) ## Feedforward -TODO +The definition of exactly what information is passed using the defined chronologies is accomplished using FeedForwards. + +Specifically, a FeedForward is used to define what to do with information being passed with an inter-stage chronology in a Simulation. The most common FeedForward is the `SemiContinuousFeedForward` that affects the semi-continuous range constraints of thermal generators in the economic dispatch problems based on the value of the (already solved) unit-commitment variables. + +The creation of a FeedForward requires at least to specify the `component_type` on which the FeedForward will be applied. The `source` variable specify which variable will be taken from the problem solved, for example the commitment variable of the thermal unit in the unit commitment problem. Finally, the `affected_values` specify which variables will be affected in the problem to be solved, for example the next economic dispatch problem. + +The following code specify the creation of semi-continuous range constraints on the `ActivePowerVariable` based on the solution of the commitment variable `OnVariable` for all `ThermalStandard` units. + +```julia +SemiContinuousFeedforward( + component_type=ThermalStandard, + source=OnVariable, + affected_values=[ActivePowerVariable], +) +``` + +## Chronologies + +In PowerSimulations, chronologies define where information is flowing. There are two types +of chronologies. + +- inter-stage chronologies: Define how information flows between stages. e.g. day-ahead solutions are used to inform economic dispatch problems +- intra-stage chronologies: Define how information flows between multiple executions of a single stage. e.g. the dispatch setpoints of the first period of an economic dispatch problem are constrained by the ramping limits from setpoints in the final period of the previous problem. ## Sequencing In a typical simulation pipeline, we want to connect daily (24-hours) day-ahead unit commitment problems, with multiple economic dispatch problems. Usually, our day-ahead unit commitment problem will have an hourly (1-hour) resolution, while the economic dispatch will have a 5-minute resolution. -Depending on your problem, it is common to use a 2-day look-ahead for unit commitment problems, so in this case, the Day-Ahead problem will have: resolution = Hour(1) with interval = Hour(24) and horizon = Hour(48). In the case of the economic dispatch problem, it is common to use a look-ahead of two hours. Thus, the Real-Time problem will have: resolution = Minute(5), with interval = Minute(5) (we only store the first operating point) and horizon = Hour(24) (24 time steps of 5 minutes are 120 minutes, that is 2 hours). + +Depending on your problem, it is common to use a 2-day look-ahead for unit commitment problems, so in this case, the Day-Ahead problem will have: resolution = Hour(1) with interval = Hour(24) and horizon = Hour(48). In the case of the economic dispatch problem, it is common to use a look-ahead of two hours. Thus, the Real-Time problem will have: resolution = Minute(5), with interval = Minute(5) (we only store the first operating point) and horizon = 24 (24 time steps of 5 minutes are 120 minutes, that is 2 hours). + +## Simulation Setup + +The following code creates the entire simulation pipeline: + +```julia +# We assume that the templates for UC and ED are ready +# sys_da has the resolution of 1 hour: +# with the 24 hours interval and horizon of 48 hours. +# sys_rt has the resolution of 5 minutes: +# with a 5-minute interval and horizon of 2 hours (24 time steps) + +# Create the UC Decision Model +decision_model_uc = DecisionModel( + template_uc, + sys_da; + name="UC", + optimizer=optimizer_with_attributes( + Xpress.Optimizer, + "MIPRELSTOP" => 1e-1, + ), +) + +# Create the ED Decision Model +decision_model_ed = DecisionModel( + template_ed, + sys_rt; + name="ED", + optimizer=optimizer_with_attributes(Xpress.Optimizer), +) + +# Specify the SimulationModels using a Vector of decision_models: UC, ED +sim_models = SimulationModels( + decision_models=[ + decision_model_uc, + decision_model_ed, + ], +) + +# Create the FeedForwards: +semi_ff = SemiContinuousFeedforward( + component_type=ThermalStandard, + source=OnVariable, + affected_values=[ActivePowerVariable], +) + +# Specify the sequencing: +sim_sequence = SimulationSequence( + # Specify the vector of decision models: sim_models + models=sim_models, + # Specify a Dict of feedforwards on which the FF applies + # based on the DecisionModel name, in this case "ED" + feedforwards=Dict( + "ED" => [semi_ff], + ), + # Specify the chronology, in this case inter-stage + ini_cond_chronology=InterProblemChronology(), +) + +# Construct the simulation: +sim = Simulation( + name="compact_sim", + steps=10, # 10 days + models=sim_models, + sequence=sim_sequence, + # Specify the start_time as a DateTime: e.g. DateTime("2020-10-01T00:00:00") + initial_time=start_time, + # Specify a temporary folder to avoid storing logs if not needed + simulation_folder=mktempdir(cleanup=true), +) + +# Build the decision models and simulation setup +build!(sim) + +# Execute the simulation using the Optimizer specified in each DecisionModel +execute!(sim, enable_progress_bar=true) +``` + +Check the [PCM tutorial](@ref pcm_tutorial) for a more detailed tutorial on executing a simulation in a production cost modeling (PCM) environment. diff --git a/docs/src/tutorials/basics_of_developing_models.md b/docs/src/tutorials/basics_of_developing_models.md index a027918d6a..654ec25da0 100644 --- a/docs/src/tutorials/basics_of_developing_models.md +++ b/docs/src/tutorials/basics_of_developing_models.md @@ -1,3 +1,3 @@ # Basics of Developing Operation Models -Check the page [PowerSimulations Structure](@ref) for more background on PowerSimulations.jl +Check the page PowerSimulations Structure for more background on PowerSimulations.jl diff --git a/src/PowerSimulations.jl b/src/PowerSimulations.jl index 23dacdd660..3183187939 100644 --- a/src/PowerSimulations.jl +++ b/src/PowerSimulations.jl @@ -219,6 +219,8 @@ export ReserveRequirementSlack export VoltageMagnitude export VoltageAngle export FlowActivePowerVariable +export FlowActivePowerSlackUpperBound +export FlowActivePowerSlackLowerBound export FlowActivePowerFromToVariable export FlowActivePowerToFromVariable export FlowReactivePowerFromToVariable @@ -227,6 +229,8 @@ export PowerAboveMinimumVariable export PhaseShifterAngle export UpperBoundFeedForwardSlack export LowerBoundFeedForwardSlack +export InterfaceFlowSlackUp +export InterfaceFlowSlackDown # Auxiliary variables export TimeDurationOn @@ -246,7 +250,7 @@ export CopperPlateBalanceConstraint export DurationConstraint export EnergyBalanceConstraint export EqualityConstraint -export FeedforwardSemiContinousConstraint +export FeedforwardSemiContinuousConstraint export FeedforwardUpperBoundConstraint export FeedforwardLowerBoundConstraint export FeedforwardIntegralLimitConstraint @@ -265,6 +269,7 @@ export FlowReactivePowerToFromConstraint export FrequencyResponseConstraint export HVDCPowerBalance export HVDCLosses +export HVDCFlowDirectionVariable export InputActivePowerVariableLimitsConstraint export NetworkFlowConstraint export NodalBalanceActiveConstraint @@ -527,8 +532,11 @@ include("simulation/simulation.jl") include("simulation/simulation_results_export.jl") include("simulation/simulation_results.jl") -include("devices_models/devices/common/objective_functions.jl") -include("devices_models/devices/common/market_bid_objective_function.jl") +include("devices_models/devices/common/objective_function/common.jl") +include("devices_models/devices/common/objective_function/linear_curve.jl") +include("devices_models/devices/common/objective_function/quadratic_curve.jl") +include("devices_models/devices/common/objective_function/market_bid.jl") +include("devices_models/devices/common/objective_function/piecewise_linear.jl") include("devices_models/devices/common/range_constraint.jl") include("devices_models/devices/common/add_variable.jl") include("devices_models/devices/common/add_auxiliary_variable.jl") diff --git a/src/core/constraints.jl b/src/core/constraints.jl index e695bdd824..6cce8e8ede 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -1,62 +1,416 @@ struct AbsoluteValueConstraint <: ConstraintType end +""" +Struct to create the constraint for starting up ThermalMultiStart units. +For more information check [ThermalGen Formulations](@ref ThermalGen-Formulations) for ThermalMultiStartUnitCommitment. + +The specified constraint is formulated as: + +```math +\\max\\{P^\\text{th,max} - P^\\text{th,shdown}, 0\\} \\cdot w_1^\\text{th} \\le u^\\text{th,init} (P^\\text{th,max} - P^\\text{th,min}) - P^\\text{th,init} +``` +""" struct ActiveRangeICConstraint <: ConstraintType end +""" +Struct to create the constraint to balance power across specified areas. +For more information check [Network Formulations](@ref network_formulations). + +The specified constraint is generally formulated as: + +```math +\\sum_{c \\in \\text{components}_a} p_t^c = 0, \\quad \\forall a\\in \\{1,\\dots, A\\}, t \\in \\{1, \\dots, T\\} +``` +""" struct AreaDispatchBalanceConstraint <: ConstraintType end struct AreaParticipationAssignmentConstraint <: ConstraintType end struct BalanceAuxConstraint <: ConstraintType end +""" +Struct to create the commitment constraint between the on, start, and stop variables. +For more information check [ThermalGen Formulations](@ref ThermalGen-Formulations). + +The specified constraints are formulated as: + +```math +u_1^\\text{th} = u^\\text{th,init} + v_1^\\text{th} - w_1^\\text{th} \\\\ +u_t^\\text{th} = u_{t-1}^\\text{th} + v_t^\\text{th} - w_t^\\text{th}, \\quad \\forall t \\in \\{2,\\dots,T\\} \\\\ +v_t^\\text{th} + w_t^\\text{th} \\le 1, \\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" struct CommitmentConstraint <: ConstraintType end +""" +Struct to create the constraint to balance power in the copperplate model. +For more information check [Network Formulations](@ref network_formulations). + +The specified constraint is generally formulated as: + +```math +\\sum_{c \\in \\text{components}} p_t^c = 0, \\quad \\forall t \\in \\{1, \\dots, T\\} +``` +""" struct CopperPlateBalanceConstraint <: ConstraintType end +""" +Struct to create the duration constraint for commitment formulations, i.e. min-up and min-down. + +For more information check [ThermalGen Formulations](@ref ThermalGen-Formulations). +""" struct DurationConstraint <: ConstraintType end struct EnergyBalanceConstraint <: ConstraintType end + +""" +Struct to create the constraint that sets the reactive power to the power factor +in the RenewableConstantPowerFactor formulation for renewable units. + +For more information check [RenewableGen Formulations](@ref PowerSystems.RenewableGen-Formulations). + +The specified constraint is formulated as: + +```math +q_t^\\text{re} = \\text{pf} \\cdot p_t^\\text{re}, \\quad \\forall t \\in \\{1,\\dots, T\\} +``` +""" struct EqualityConstraint <: ConstraintType end -struct FeedforwardSemiContinousConstraint <: ConstraintType end +""" +Struct to create the constraint for semicontinuous feedforward limits. + +For more information check [Feedforward Formulations](@ref ff_formulations). + +The specified constraint is formulated as: + +```math +\\begin{align*} +& \\text{ActivePowerRangeExpressionUB}_t := p_t^\\text{th} - \\text{on}_t^\\text{th}P^\\text{th,max} \\le 0, \\quad \\forall t\\in \\{1, \\dots, T\\} \\\\ +& \\text{ActivePowerRangeExpressionLB}_t := p_t^\\text{th} - \\text{on}_t^\\text{th}P^\\text{th,min} \\ge 0, \\quad \\forall t\\in \\{1, \\dots, T\\} +\\end{align*} +``` +""" +struct FeedforwardSemiContinuousConstraint <: ConstraintType end struct FeedforwardIntegralLimitConstraint <: ConstraintType end +""" +Struct to create the constraint for upper bound feedforward limits. + +For more information check [Feedforward Formulations](@ref ff_formulations). + +The specified constraint is formulated as: + +```math +\\begin{align*} +& \\text{AffectedVariable}_t - p_t^\\text{ff,ubsl} \\le \\text{SourceVariableParameter}_t, \\quad \\forall t \\in \\{1,\\dots, T\\} +\\end{align*} +``` +""" struct FeedforwardUpperBoundConstraint <: ConstraintType end +""" +Struct to create the constraint for lower bound feedforward limits. + +For more information check [Feedforward Formulations](@ref ff_formulations). + +The specified constraint is formulated as: + +```math +\\begin{align*} +& \\text{AffectedVariable}_t + p_t^\\text{ff,lbsl} \\ge \\text{SourceVariableParameter}_t, \\quad \\forall t \\in \\{1,\\dots, T\\} +\\end{align*} +``` +""" struct FeedforwardLowerBoundConstraint <: ConstraintType end struct FeedforwardEnergyTargetConstraint <: ConstraintType end struct FlowActivePowerConstraint <: ConstraintType end #not being used struct FlowActivePowerFromToConstraint <: ConstraintType end #not being used struct FlowActivePowerToFromConstraint <: ConstraintType end #not being used -struct FlowLimitConstraint <: ConstraintType end #not being used +""" +Struct to create the constraint that set the flow limits through a PhaseShiftingTransformer. + +For more information check [Branch Formulations](@ref PowerSystems.Branch-Formulations). + +The specified constraint is formulated as: + +```math +-R^\\text{max} \\le f_t \\le R^\\text{max}, \\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" +struct FlowLimitConstraint <: ConstraintType end struct FlowLimitFromToConstraint <: ConstraintType end struct FlowLimitToFromConstraint <: ConstraintType end +""" +Struct to create the constraint that set the flow limits through an HVDC two-terminal branch. + +For more information check [Branch Formulations](@ref PowerSystems.Branch-Formulations). + +The specified constraint is formulated as: + +```math +R^\\text{min} \\le f_t \\le R^\\text{max}, \\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" struct FlowRateConstraint <: ConstraintType end +""" +Struct to create the constraint that set the flow from-to limits through an HVDC two-terminal branch. + +For more information check [Branch Formulations](@ref PowerSystems.Branch-Formulations). + +The specified constraint is formulated as: + +```math +R^\\text{from,min} \\le f_t^\\text{from-to} \\le R^\\text{from,max}, \\forall t \\in \\{1,\\dots, T\\} +``` +""" struct FlowRateConstraintFromTo <: ConstraintType end +""" +Struct to create the constraint that set the flow to-from limits through an HVDC two-terminal branch. + +For more information check [Branch Formulations](@ref PowerSystems.Branch-Formulations). + +The specified constraint is formulated as: + +```math +R^\\text{to,min} \\le f_t^\\text{to-from} \\le R^\\text{to,max},\\quad \\forall t \\in \\{1,\\dots, T\\} +``` +""" struct FlowRateConstraintToFrom <: ConstraintType end struct FlowReactivePowerConstraint <: ConstraintType end #not being used struct FlowReactivePowerFromToConstraint <: ConstraintType end #not being used struct FlowReactivePowerToFromConstraint <: ConstraintType end #not being used +""" +Struct to create the constraints that set the power balance across a lossy HVDC two-terminal line. + +For more information check [Branch Formulations](@ref PowerSystems.Branch-Formulations). + +The specified constraints are formulated as: + +```math +\\begin{align*} +& f_t^\\text{to-from} - f_t^\\text{from-to} \\le L_1 \\cdot f_t^\\text{to-from} - L_0,\\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ +& f_t^\\text{from-to} - f_t^\\text{to-from} \\ge L_1 \\cdot f_t^\\text{from-to} + L_0,\\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ +& f_t^\\text{from-to} - f_t^\\text{to-from} \\ge - M^\\text{big} (1 - u^\\text{dir}_t),\\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ +& f_t^\\text{to-from} - f_t^\\text{from-to} \\ge - M^\\text{big} u^\\text{dir}_t,\\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ +\\end{align*} +``` +""" struct HVDCPowerBalance <: ConstraintType end struct FrequencyResponseConstraint <: ConstraintType end +""" +Struct to create the constraint the AC branch flows depending on the network model. +For more information check [Branch Formulations](@ref PowerSystems.Branch-Formulations). + +The specified constraint depends on the network model chosen. The most common application is the StaticBranch in a PTDF Network Model: + +```math +f_t = \\sum_{i=1}^N \\text{PTDF}_{i,b} \\cdot \\text{Bal}_{i,t}, \\quad \\forall t \\in \\{1,\\dots, T\\} +``` +""" struct NetworkFlowConstraint <: ConstraintType end +""" +Struct to create the constraint to balance active power in nodal formulation. +For more information check [Network Formulations](@ref network_formulations). + +The specified constraint depends on the network model chosen. +""" struct NodalBalanceActiveConstraint <: ConstraintType end +""" +Struct to create the constraint to balance reactive power in nodal formulation. +For more information check [Network Formulations](@ref network_formulations). + +The specified constraint depends on the network model chosen. +""" struct NodalBalanceReactiveConstraint <: ConstraintType end struct ParticipationAssignmentConstraint <: ConstraintType end +""" +Struct to create the constraint to participation assignments limits in the active power reserves. +For more information check [Service Formulations](@ref service_formulations). + +The constraint is as follows: + +```math +r_{d,t} \\le \\text{Req} \\cdot \\text{PF} ,\\quad \\forall d\\in \\mathcal{D}_s, \\forall t\\in \\{1,\\dots, T\\} \\quad \\text{(for a StaticReserve)} \\\\ +r_{d,t} \\le \\text{RequirementTimeSeriesParameter}_{t} \\cdot \\text{PF}\\quad \\forall d\\in \\mathcal{D}_s, \\forall t\\in \\{1,\\dots, T\\}, \\quad \\text{(for a VariableReserve)} +``` +""" struct ParticipationFractionConstraint <: ConstraintType end +""" +Struct to create the PieceWiseLinearCostConstraint associated with a specified variable. + +See [Piecewise linear cost functions](@ref pwl_cost) for more information. +""" struct PieceWiseLinearCostConstraint <: ConstraintType end + +""" +Struct to create the PieceWiseLinearBlockOfferConstraint associated with a specified variable. + +See [Piecewise linear cost functions](@ref pwl_cost) for more information. +""" +struct PieceWiseLinearBlockOfferConstraint <: ConstraintType end + +""" +Struct to create the PieceWiseLinearUpperBoundConstraint associated with a specified variable. + +See [Piecewise linear cost functions](@ref pwl_cost) for more information. +""" +struct PieceWiseLinearUpperBoundConstraint <: ConstraintType end + +""" +Struct to create the RampConstraint associated with a specified thermal device or reserve service. + +For thermal units, see more information in [Thermal Formulations](@ref ThermalGen-Formulations). The constraint is as follows: +```math +-R^\\text{th,dn} \\le p_t^\\text{th} - p_{t-1}^\\text{th} \\le R^\\text{th,up}, \\quad \\forall t\\in \\{1, \\dots, T\\} +``` + +For Ramp Reserve, see more information in [Service Formulations](@ref service_formulations). The constraint is as follows: + +```math +r_{d,t} \\le R^\\text{th,up} \\cdot \\text{TF}\\quad \\forall d\\in \\mathcal{D}_s, \\forall t\\in \\{1,\\dots, T\\}, \\quad \\text{(for ReserveUp)} \\\\ +r_{d,t} \\le R^\\text{th,dn} \\cdot \\text{TF}\\quad \\forall d\\in \\mathcal{D}_s, \\forall t\\in \\{1,\\dots, T\\}, \\quad \\text{(for ReserveDown)} +``` +""" struct RampConstraint <: ConstraintType end struct RampLimitConstraint <: ConstraintType end struct RangeLimitConstraint <: ConstraintType end +""" +Struct to create the constraint that set the AC flow limits through branches. + +For more information check [Branch Formulations](@ref PowerSystems.Branch-Formulations). + +The specified constraint is formulated as: + +```math +\\begin{align*} +& f_t - f_t^\\text{sl,up} \\le R^\\text{max},\\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ +& f_t + f_t^\\text{sl,lo} \\ge -R^\\text{max},\\quad \\forall t \\in \\{1,\\dots, T\\} +\\end{align*} +``` +""" struct RateLimitConstraint <: ConstraintType end struct RateLimitConstraintFromTo <: ConstraintType end struct RateLimitConstraintToFrom <: ConstraintType end struct RegulationLimitsConstraint <: ConstraintType end +""" +Struct to create the constraint for satisfying active power reserve requirements. +For more information check [Service Formulations](@ref service_formulations). + +The constraint is as follows: + +```math +\\sum_{d\\in\\mathcal{D}_s} r_{d,t} + r_t^\\text{sl} \\ge \\text{Req},\\quad \\forall t\\in \\{1,\\dots, T\\} \\quad \\text{(for a StaticReserve)} \\\\ +\\sum_{d\\in\\mathcal{D}_s} r_{d,t} + r_t^\\text{sl} \\ge \\text{RequirementTimeSeriesParameter}_{t},\\quad \\forall t\\in \\{1,\\dots, T\\} \\quad \\text{(for a VariableReserve)} +``` +""" struct RequirementConstraint <: ConstraintType end struct ReserveEnergyCoverageConstraint <: ConstraintType end +""" +Struct to create the constraint for ensuring that NonSpinning Reserve can be delivered from turn-off thermal units. + +For more information check [Service Formulations](@ref service_formulations) for NonSpinningReserve. + +The constraint is as follows: + +```math +r_{d,t} \\le (1 - u_{d,t}^\\text{th}) \\cdot R^\\text{limit}_d, \\quad \\forall d \\in \\mathcal{D}_s, \\forall t \\in \\{1,\\dots, T\\} +``` +""" struct ReservePowerConstraint <: ConstraintType end struct SACEPIDAreaConstraint <: ConstraintType end struct StartTypeConstraint <: ConstraintType end +""" +Struct to create the start-up initial condition constraints for ThermalMultiStart. + +For more information check [ThermalGen Formulations](@ref ThermalGen-Formulations) for ThermalMultiStartUnitCommitment. +""" struct StartupInitialConditionConstraint <: ConstraintType end +""" +Struct to create the start-up time limit constraints for ThermalMultiStart. + +For more information check [ThermalGen Formulations](@ref ThermalGen-Formulations) for ThermalMultiStartUnitCommitment. +""" struct StartupTimeLimitTemperatureConstraint <: ConstraintType end +""" +Struct to create the constraint that set the angle limits through a PhaseShiftingTransformer. + +For more information check [Branch Formulations](@ref PowerSystems.Branch-Formulations). + +The specified constraint is formulated as: + +```math +\\Theta^\\text{min} \\le \\theta^\\text{shift}_t \\le \\Theta^\\text{max}, \\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" struct PhaseAngleControlLimit <: ConstraintType end +""" +Struct to create the constraints that set the losses through a lossy HVDC two-terminal line. + +For more information check [Branch Formulations](@ref PowerSystems.Branch-Formulations). + +The specified constraints are formulated as: + +```math +\\begin{align*} +& f_t^\\text{to-from} - f_t^\\text{from-to} \\le \\ell_t,\\quad \\forall t \\in \\{1,\\dots, T\\} \\\\ +& f_t^\\text{from-to} - f_t^\\text{to-from} \\le \\ell_t,\\quad \\forall t \\in \\{1,\\dots, T\\} +\\end{align*} +``` +""" struct HVDCLossesAbsoluteValue <: ConstraintType end struct HVDCDirection <: ConstraintType end struct InterfaceFlowLimit <: ConstraintType end abstract type PowerVariableLimitsConstraint <: ConstraintType end +""" +Struct to create the constraint to limit active power input expressions. +For more information check [Device Formulations](@ref formulation_intro). + +The specified constraint depends on the UpperBound and LowerBound expressions, but +in its most basic formulation is of the form: + +```math +P^\\text{min} \\le p_t^\\text{in} \\le P^\\text{max}, \\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" struct InputActivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end +""" +Struct to create the constraint to limit active power output expressions. +For more information check [Device Formulations](@ref formulation_intro). + +The specified constraint depends on the UpperBound and LowerBound expressions, but +in its most basic formulation is of the form: + +```math +P^\\text{min} \\le p_t^\\text{out} \\le P^\\text{max}, \\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" struct OutputActivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end +""" +Struct to create the constraint to limit active power expressions. +For more information check [Device Formulations](@ref formulation_intro). + +The specified constraint depends on the UpperBound and LowerBound expressions, but +in its most basic formulation is of the form: + +```math +P^\\text{min} \\le p_t \\le P^\\text{max}, \\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" struct ActivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end +""" +Struct to create the constraint to limit reactive power expressions. +For more information check [Device Formulations](@ref formulation_intro). + +The specified constraint depends on the UpperBound and LowerBound expressions, but +in its most basic formulation is of the form: + +```math +Q^\\text{min} \\le q_t \\le Q^\\text{max}, \\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" struct ReactivePowerVariableLimitsConstraint <: PowerVariableLimitsConstraint end +""" +Struct to create the constraint to limit active power expressions by a time series parameter. +For more information check [Device Formulations](@ref formulation_intro). + +The specified constraint depends on the UpperBound expressions, but +in its most basic formulation is of the form: + +```math +p_t \\le \\text{ActivePowerTimeSeriesParameter}_t, \\quad \\forall t \\in \\{1,\\dots,T\\} +``` +""" struct ActivePowerVariableTimeSeriesLimitsConstraint <: PowerVariableLimitsConstraint end abstract type EventConstraint <: ConstraintType end diff --git a/src/core/definitions.jl b/src/core/definitions.jl index d7dd810df4..6d89cbe94c 100644 --- a/src/core/definitions.jl +++ b/src/core/definitions.jl @@ -45,6 +45,7 @@ const MISSING_INITIAL_CONDITIONS_TIME_COUNT = 999.0 const SECONDS_IN_MINUTE = 60.0 const MINUTES_IN_HOUR = 60.0 const SECONDS_IN_HOUR = 3600.0 +const MILLISECONDS_IN_HOUR = 3600000.0 const MAX_START_STAGES = 3 const OBJECTIVE_FUNCTION_POSITIVE = 1.0 const OBJECTIVE_FUNCTION_NEGATIVE = -1.0 diff --git a/src/core/formulations.jl b/src/core/formulations.jl index 2afdff9183..0bdffc9ad4 100644 --- a/src/core/formulations.jl +++ b/src/core/formulations.jl @@ -33,7 +33,7 @@ Formulation type to enable standard dispatch with a range and enforce intertempo """ struct ThermalStandardDispatch <: AbstractThermalDispatchFormulation end """ -Formulation type to enable basic dispatch without any intertemporal constraints and relaxed minimum generation. *may not work with PWL cost definitions* +Formulation type to enable basic dispatch without any intertemporal constraints and relaxed minimum generation. *May not work with non-convex PWL cost definitions* """ struct ThermalDispatchNoMin <: AbstractThermalDispatchFormulation end """ @@ -58,7 +58,7 @@ abstract type AbstractLoadFormulation <: AbstractDeviceFormulation end abstract type AbstractControllablePowerLoadFormulation <: AbstractLoadFormulation end """ -Formulation type to add a time series parameter for non-dispatchable `ElectricLoad` withdrawls to power balance constraints +Formulation type to add a time series parameter for non-dispatchable `ElectricLoad` withdrawals to power balance constraints """ struct StaticPowerLoad <: AbstractLoadFormulation end @@ -145,7 +145,9 @@ LossLess InterconnectingConverter Model """ struct LossLessConverter <: AbstractConverterFormulation end -# TODO: Think if this an ok abstraction for future use cases +""" +LossLess Line Abstract Model +""" struct LossLessLine <: AbstractBranchFormulation end ############################## Network Model Formulations ################################## @@ -153,11 +155,11 @@ struct LossLessLine <: AbstractBranchFormulation end abstract type AbstractPTDFModel <: PM.AbstractDCPModel end """ -Linear active power approximation using the power transfer distribution factor ((PTDF)[https://nrel-sienna.github.io/PowerNetworkMatrices.jl/stable/tutorials/tutorial_PTDF_matrix/]) matrix. +Linear active power approximation using the power transfer distribution factor [PTDF](https://nrel-sienna.github.io/PowerNetworkMatrices.jl/stable/tutorials/tutorial_PTDF_matrix/) matrix. """ struct PTDFPowerModel <: AbstractPTDFModel end """ -Infinate capacity approximation of network flow to represent entire system with a single node. +Infinite capacity approximation of network flow to represent entire system with a single node. """ struct CopperPlatePowerModel <: PM.AbstractActivePowerModel end """ @@ -227,10 +229,28 @@ abstract type AbstractAGCFormulation <: AbstractServiceFormulation end struct PIDSmoothACE <: AbstractAGCFormulation end +""" +Struct to add reserves to be larger than a specified requirement for an aggregated collection of services +""" struct GroupReserve <: AbstractReservesFormulation end + +""" +Struct for to add reserves to be larger than a specified requirement +""" struct RangeReserve <: AbstractReservesFormulation end +""" +Struct for to add reserves to be larger than a variable requirement depending of costs +""" struct StepwiseCostReserve <: AbstractReservesFormulation end +""" +Struct to add reserves to be larger than a specified requirement, with ramp constraints +""" struct RampReserve <: AbstractReservesFormulation end +""" +Struct to add non spinning reserve requirements larger than specified requirement +""" struct NonSpinningReserve <: AbstractReservesFormulation end - +""" +Struct to add a constant maximum transmission flow for specified interface +""" struct ConstantMaxInterfaceFlow <: AbstractServiceFormulation end diff --git a/src/core/optimization_container.jl b/src/core/optimization_container.jl index 5cd63b65b3..dace2ba0b9 100644 --- a/src/core/optimization_container.jl +++ b/src/core/optimization_container.jl @@ -846,7 +846,7 @@ function add_variable_container!( ::T, ::Type{U}; meta = IS.Optimization.CONTAINER_KEY_EMPTY_META, -) where {T <: PieceWiseLinearCostVariable, U <: Union{PSY.Component, PSY.System}} +) where {T <: SparseVariableType, U <: Union{PSY.Component, PSY.System}} var_key = VariableKey(T, U, meta) _assign_container!(container.variables, var_key, _get_pwl_variables_container()) return container.variables[var_key] diff --git a/src/core/settings.jl b/src/core/settings.jl index bd1a08806c..d9dd095cbc 100644 --- a/src/core/settings.jl +++ b/src/core/settings.jl @@ -44,11 +44,10 @@ function Settings( store_variable_names = false, ext = Dict{String, Any}(), ) - - # if time_series_cache_size > 0 && PSY.stores_time_series_in_memory(sys) - # @info "Overriding time_series_cache_size because time series is stored in memory" - # time_series_cache_size = 0 - # end + if time_series_cache_size > 0 && PSY.stores_time_series_in_memory(sys) + @info "Overriding time_series_cache_size because time series is stored in memory" + time_series_cache_size = 0 + end if isa(optimizer, MOI.OptimizerWithAttributes) || optimizer === nothing optimizer_ = optimizer diff --git a/src/core/variables.jl b/src/core/variables.jl index de535d5854..42d5bc0560 100644 --- a/src/core/variables.jl +++ b/src/core/variables.jl @@ -1,56 +1,56 @@ """ Struct to dispatch the creation of Active Power Variables -Docs abbreviation: ``Pg`` +Docs abbreviation: ``p`` """ struct ActivePowerVariable <: VariableType end """ Struct to dispatch the creation of Active Power Variables above minimum power for Thermal Compact formulations -Docs abbreviation: ``\\hat{Pg}`` +Docs abbreviation: ``\\Delta p`` """ struct PowerAboveMinimumVariable <: VariableType end """ Struct to dispatch the creation of Active Power Input Variables for 2-directional devices. For instance storage or pump-hydro -Docs abbreviation: ``Pg^{in}`` +Docs abbreviation: ``p^\\text{in}`` """ struct ActivePowerInVariable <: VariableType end """ Struct to dispatch the creation of Active Power Output Variables for 2-directional devices. For instance storage or pump-hydro -Docs abbreviation: ``Pg^{out}`` +Docs abbreviation: ``p^\\text{out}`` """ struct ActivePowerOutVariable <: VariableType end """ Struct to dispatch the creation of Hot Start Variable for Thermal units with temperature considerations -Docs abbreviation: TODO +Docs abbreviation: ``z^\\text{th}`` """ struct HotStartVariable <: VariableType end """ Struct to dispatch the creation of Warm Start Variable for Thermal units with temperature considerations -Docs abbreviation: TODO +Docs abbreviation: ``y^\\text{th}`` """ struct WarmStartVariable <: VariableType end """ Struct to dispatch the creation of Cold Start Variable for Thermal units with temperature considerations -Docs abbreviation: TODO +Docs abbreviation: ``x^\\text{th}`` """ struct ColdStartVariable <: VariableType end """ Struct to dispatch the creation of a variable for energy storage level (state of charge) -Docs abbreviation: ``E`` +Docs abbreviation: ``e`` """ struct EnergyVariable <: VariableType end @@ -66,37 +66,42 @@ struct OnVariable <: VariableType end """ Struct to dispatch the creation of Reactive Power Variables -Docs abbreviation: ``Qg`` +Docs abbreviation: ``q`` """ struct ReactivePowerVariable <: VariableType end """ Struct to dispatch the creation of binary storage charge reservation variable -Docs abbreviation: ``r`` +Docs abbreviation: ``u^\\text{st}`` """ struct ReservationVariable <: VariableType end """ Struct to dispatch the creation of Active Power Reserve Variables -Docs abbreviation: ``Pr`` +Docs abbreviation: ``r`` """ struct ActivePowerReserveVariable <: VariableType end +""" +Struct to dispatch the creation of Service Requirement Variables + +Docs abbreviation: ``\\text{req}`` +""" struct ServiceRequirementVariable <: VariableType end """ Struct to dispatch the creation of Binary Start Variables -Docs abbreviation: TODO +Docs abbreviation: ``v`` """ struct StartVariable <: VariableType end """ Struct to dispatch the creation of Binary Stop Variables -Docs abbreviation: TODO +Docs abbreviation: ``w`` """ struct StopVariable <: VariableType end @@ -114,34 +119,59 @@ struct AdditionalDeltaActivePowerDownVariable <: VariableType end struct SmoothACE <: VariableType end +""" +Struct to dispatch the creation of System-wide slack up variables. Used when there is not enough generation. + +Docs abbreviation: ``p^\\text{sl,up}`` +""" struct SystemBalanceSlackUp <: VariableType end +""" +Struct to dispatch the creation of System-wide slack down variables. Used when there is not enough load curtailment. + +Docs abbreviation: ``p^\\text{sl,dn}`` +""" struct SystemBalanceSlackDown <: VariableType end +""" +Struct to dispatch the creation of Reserve requirement slack variables. Used when there is not reserves in the system to satisfy the requirement. + +Docs abbreviation: ``r^\\text{sl}`` +""" struct ReserveRequirementSlack <: VariableType end +""" +Struct to dispatch the creation of active power flow upper bound slack variables. Used when there is not enough flow through the branch in the forward direction. + +Docs abbreviation: ``f^\\text{sl,up}`` +""" struct FlowActivePowerSlackUpperBound <: VariableType end +""" +Struct to dispatch the creation of active power flow lower bound slack variables. Used when there is not enough flow through the branch in the reverse direction. + +Docs abbreviation: ``f^\\text{sl,lo}`` +""" struct FlowActivePowerSlackLowerBound <: VariableType end """ Struct to dispatch the creation of Voltage Magnitude Variables for AC formulations -Docs abbreviation: TODO +Docs abbreviation: ``v`` """ struct VoltageMagnitude <: VariableType end """ Struct to dispatch the creation of Voltage Angle Variables for AC/DC formulations -Docs abbreviation: TODO +Docs abbreviation: ``\\theta`` """ struct VoltageAngle <: VariableType end """ Struct to dispatch the creation of bidirectional Active Power Flow Variables -Docs abbreviation: ``P`` +Docs abbreviation: ``f`` """ struct FlowActivePowerVariable <: VariableType end @@ -151,35 +181,35 @@ struct FlowActivePowerVariable <: VariableType end """ Struct to dispatch the creation of unidirectional Active Power Flow Variables -Docs abbreviation: ``\\overrightarrow{P}`` +Docs abbreviation: ``f^\\text{from-to}`` """ struct FlowActivePowerFromToVariable <: VariableType end """ Struct to dispatch the creation of unidirectional Active Power Flow Variables -Docs abbreviation: ``\\overleftarrow{P}`` +Docs abbreviation: ``f^\\text{to-from}`` """ struct FlowActivePowerToFromVariable <: VariableType end """ Struct to dispatch the creation of unidirectional Reactive Power Flow Variables -Docs abbreviation: ``\\overrightarrow{Q}`` +Docs abbreviation: ``f^\\text{q,from-to}`` """ struct FlowReactivePowerFromToVariable <: VariableType end """ Struct to dispatch the creation of unidirectional Reactive Power Flow Variables -Docs abbreviation: ``\\overleftarrow{Q}`` +Docs abbreviation: ``f^\\text{q,to-from}`` """ struct FlowReactivePowerToFromVariable <: VariableType end """ Struct to dispatch the creation of Phase Shifters Variables -Docs abbreviation: TODO +Docs abbreviation: ``\\theta^\\text{shift}`` """ struct PhaseShifterAngle <: VariableType end @@ -187,36 +217,63 @@ struct PhaseShifterAngle <: VariableType end """ Struct to dispatch the creation of HVDC Losses Auxiliary Variables -Docs abbreviation: TODO +Docs abbreviation: ``\\ell`` """ struct HVDCLosses <: VariableType end """ Struct to dispatch the creation of HVDC Flow Direction Auxiliary Variables -Docs abbreviation: TODO +Docs abbreviation: ``u^\\text{dir}`` """ struct HVDCFlowDirectionVariable <: VariableType end +abstract type SparseVariableType <: VariableType end + """ Struct to dispatch the creation of piecewise linear cost variables for objective function -Docs abbreviation: TODO +Docs abbreviation: ``\\delta`` +""" +struct PieceWiseLinearCostVariable <: SparseVariableType end + +""" +Struct to dispatch the creation of piecewise linear block offer variables for objective function + +Docs abbreviation: ``\\delta`` +""" +struct PieceWiseLinearBlockOffer <: SparseVariableType end + """ -struct PieceWiseLinearCostVariable <: VariableType end +Struct to dispatch the creation of Interface Flow Slack Up variables +Docs abbreviation: ``f^\\text{sl,up}`` +""" struct InterfaceFlowSlackUp <: VariableType end +""" +Struct to dispatch the creation of Interface Flow Slack Down variables +Docs abbreviation: ``f^\\text{sl,dn}`` +""" struct InterfaceFlowSlackDown <: VariableType end +""" +Struct to dispatch the creation of Slack variables for UpperBoundFeedforward + +Docs abbreviation: ``p^\\text{ff,ubsl}`` +""" struct UpperBoundFeedForwardSlack <: VariableType end +""" +Struct to dispatch the creation of Slack variables for LowerBoundFeedforward +Docs abbreviation: ``p^\\text{ff,lbsl}`` +""" struct LowerBoundFeedForwardSlack <: VariableType end const START_VARIABLES = (HotStartVariable, WarmStartVariable, ColdStartVariable) should_write_resulting_value(::Type{PieceWiseLinearCostVariable}) = false - +should_write_resulting_value(::Type{PieceWiseLinearBlockOffer}) = false convert_result_to_natural_units(::Type{ActivePowerVariable}) = true convert_result_to_natural_units(::Type{PowerAboveMinimumVariable}) = true convert_result_to_natural_units(::Type{ActivePowerInVariable}) = true diff --git a/src/devices_models/devices/common/add_to_expression.jl b/src/devices_models/devices/common/add_to_expression.jl index c3127f79ba..b0c248a1a0 100644 --- a/src/devices_models/devices/common/add_to_expression.jl +++ b/src/devices_models/devices/common/add_to_expression.jl @@ -1012,7 +1012,7 @@ function add_to_expression!( end expression = get_expression(container, T(), V) for d in devices - mult in get_expression_multiplier(U(), T(), d, W()) + mult = get_expression_multiplier(U(), T(), d, W()) for t in get_time_steps(container) name = PSY.get_name(d) _add_to_jump_expression!( @@ -1043,7 +1043,8 @@ function add_to_expression!( add_expressions!(container, T, devices, model) end expression = get_expression(container, T(), V) - for d in devices, mult in get_expression_multiplier(U(), T(), d, W()) + for d in devices + mult = get_expression_multiplier(U(), T(), d, W()) for t in get_time_steps(container) name = PSY.get_name(d) _add_to_jump_expression!(expression[name, t], parameter_array[name, t], -mult) diff --git a/src/devices_models/devices/common/market_bid_objective_function.jl b/src/devices_models/devices/common/market_bid_objective_function.jl deleted file mode 100644 index 7a22a134ca..0000000000 --- a/src/devices_models/devices/common/market_bid_objective_function.jl +++ /dev/null @@ -1,102 +0,0 @@ -function _add_variable_cost_to_objective!( - container::OptimizationContainer, - ::T, - component::PSY.Component, - op_cost::PSY.MarketBidCost, - ::U, -) where {T <: VariableType, U <: AbstractDeviceFormulation} - component_name = PSY.get_name(component) - @debug "Market Bid" _group = LOG_GROUP_COST_FUNCTIONS component_name - time_steps = get_time_steps(container) - initial_time = get_initial_time(container) - variable_cost_forecast = PSY.get_variable_cost( - component, - op_cost; - start_time = initial_time, - len = length(time_steps), - ) - variable_cost_forecast_values = TimeSeries.values(variable_cost_forecast) - parameter_container = _get_cost_function_parameter_container( - container, - CostFunctionParameter(), - component, - T(), - U(), - eltype(variable_cost_forecast_values), - ) - pwl_cost_expressions = - _add_pwl_term!(container, component, variable_cost_forecast_values, T(), U()) - jump_model = get_jump_model(container) - for t in time_steps - set_multiplier!( - parameter_container, - # Using 1.0 here since we want to reuse the existing code that adds the mulitpler - # of base power times the time delta. - 1.0, - component_name, - t, - ) - set_parameter!( - parameter_container, - jump_model, - variable_cost_forecast_values[t], - component_name, - t, - ) - add_to_expression!( - container, - ProductionCostExpression, - pwl_cost_expressions[t], - component, - t, - ) - add_to_objective_variant_expression!(container, pwl_cost_expressions[t]) - end - - # Service Cost Bid - ancillary_services = PSY.get_ancillary_service_offers(op_cost) - for service in ancillary_services - _add_service_bid_cost!(container, component, service) - end - return -end - -function _add_service_bid_cost!( - container::OptimizationContainer, - component::PSY.Component, - service::T, -) where {T <: PSY.Reserve{<:PSY.ReserveDirection}} - time_steps = get_time_steps(container) - initial_time = get_initial_time(container) - base_power = get_base_power(container) - forecast_data = PSY.get_services_bid( - component, - PSY.get_operation_cost(component), - service; - start_time = initial_time, - len = length(time_steps), - ) - forecast_data_values = PSY.get_cost.(TimeSeries.values(forecast_data)) - # Single Price Bid - if eltype(forecast_data_values) == Float64 - data_values = forecast_data_values - # Single Price/Quantity Bid - elseif eltype(forecast_data_values) == Vector{NTuple{2, Float64}} - data_values = [v[1][1] for v in forecast_data_values] - else - error("$(eltype(forecast_data_values)) not supported for MarketBidCost") - end - - reserve_variable = - get_variable(container, ActivePowerReserveVariable(), T, PSY.get_name(service)) - component_name = PSY.get_name(component) - for t in time_steps - add_to_objective_invariant_expression!( - container, - data_values[t] * base_power * reserve_variable[component_name, t], - ) - end - return -end - -function _add_service_bid_cost!(::OptimizationContainer, ::PSY.Component, ::PSY.Service) end diff --git a/src/devices_models/devices/common/objective_function/common.jl b/src/devices_models/devices/common/objective_function/common.jl new file mode 100644 index 0000000000..082cb79932 --- /dev/null +++ b/src/devices_models/devices/common/objective_function/common.jl @@ -0,0 +1,236 @@ +################################## +#### ActivePowerVariable Cost #### +################################## + +function add_variable_cost!( + container::OptimizationContainer, + ::U, + devices::IS.FlattenIteratorWrapper{T}, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} + for d in devices + op_cost_data = PSY.get_operation_cost(d) + _add_variable_cost_to_objective!(container, U(), d, op_cost_data, V()) + end + return +end + +################################## +#### Start/Stop Variable Cost #### +################################## + +function add_shut_down_cost!( + container::OptimizationContainer, + ::U, + devices::IS.FlattenIteratorWrapper{T}, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} + multiplier = objective_function_multiplier(U(), V()) + for d in devices + op_cost_data = PSY.get_operation_cost(d) + cost_term = shut_down_cost(op_cost_data, d, V()) + iszero(cost_term) && continue + for t in get_time_steps(container) + _add_proportional_term!(container, U(), d, cost_term * multiplier, t) + end + end + return +end + +################################## +####### Proportional Cost ######## +################################## + +function add_proportional_cost!( + container::OptimizationContainer, + ::U, + devices::IS.FlattenIteratorWrapper{T}, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} + multiplier = objective_function_multiplier(U(), V()) + for d in devices + op_cost_data = PSY.get_operation_cost(d) + cost_term = proportional_cost(op_cost_data, U(), d, V()) + iszero(cost_term) && continue + for t in get_time_steps(container) + _add_proportional_term!(container, U(), d, cost_term * multiplier, t) + end + end + return +end + +################################## +######## OnVariable Cost ######### +################################## + +function add_proportional_cost!( + container::OptimizationContainer, + ::U, + devices::IS.FlattenIteratorWrapper{T}, + ::V, +) where {T <: PSY.ThermalGen, U <: OnVariable, V <: AbstractCompactUnitCommitment} + multiplier = objective_function_multiplier(U(), V()) + for d in devices + op_cost_data = PSY.get_operation_cost(d) + cost_term = proportional_cost(op_cost_data, U(), d, V()) + iszero(cost_term) && continue + for t in get_time_steps(container) + exp = _add_proportional_term!(container, U(), d, cost_term * multiplier, t) + add_to_expression!(container, ProductionCostExpression, exp, d, t) + end + end + return +end + +function _add_variable_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + op_cost::PSY.OperationalCost, + ::U, +) where {T <: VariableType, U <: AbstractDeviceFormulation} + variable_cost_data = variable_cost(op_cost, T(), component, U()) + _add_variable_cost_to_objective!(container, T(), component, variable_cost_data, U()) + return +end + +function add_start_up_cost!( + container::OptimizationContainer, + ::U, + devices::IS.FlattenIteratorWrapper{T}, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} + for d in devices + op_cost_data = PSY.get_operation_cost(d) + _add_start_up_cost_to_objective!(container, U(), d, op_cost_data, V()) + end + return +end + +function _add_start_up_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.ThermalGen, + op_cost::Union{PSY.ThermalGenerationCost, PSY.MarketBidCost}, + ::U, +) where {T <: VariableType, U <: AbstractDeviceFormulation} + cost_term = start_up_cost(op_cost, component, U()) + iszero(cost_term) && return + multiplier = objective_function_multiplier(T(), U()) + for t in get_time_steps(container) + _add_proportional_term!(container, T(), component, cost_term * multiplier, t) + end + return +end + +const MULTI_START_COST_MAP = Dict{DataType, Int}( + HotStartVariable => 1, + WarmStartVariable => 2, + ColdStartVariable => 3, +) + +function _add_start_up_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.ThermalMultiStart, + op_cost::PSY.ThermalGenerationCost, + ::U, +) where {T <: VariableType, U <: ThermalMultiStartUnitCommitment} + cost_terms = start_up_cost(op_cost, component, U()) + cost_term = cost_terms[MULTI_START_COST_MAP[T]] + iszero(cost_term) && return + multiplier = objective_function_multiplier(T(), U()) + for t in get_time_steps(container) + _add_proportional_term!(container, T(), component, cost_term * multiplier, t) + end + return +end + +function _get_cost_function_parameter_container( + container::OptimizationContainer, + ::S, + component::T, + ::U, + ::V, + cost_type::Type{W}, +) where { + S <: ObjectiveFunctionParameter, + T <: PSY.Component, + U <: VariableType, + V <: Union{AbstractDeviceFormulation, AbstractServiceFormulation}, + W, +} + if has_container_key(container, S, T) + return get_parameter(container, S(), T) + else + container_axes = axes(get_variable(container, U(), T)) + if has_container_key(container, OnStatusParameter, T) + sos_val = SOSStatusVariable.PARAMETER + else + sos_val = sos_status(component, V()) + end + return add_param_container!( + container, + S(), + T, + U, + sos_val, + uses_compact_power(component, V()), + W, + container_axes..., + ) + end +end + +function _add_proportional_term!( + container::OptimizationContainer, + ::T, + component::U, + linear_term::Float64, + time_period::Int, +) where {T <: VariableType, U <: PSY.Component} + component_name = PSY.get_name(component) + @debug "Linear Variable Cost" _group = LOG_GROUP_COST_FUNCTIONS component_name + variable = get_variable(container, T(), U)[component_name, time_period] + lin_cost = variable * linear_term + add_to_objective_invariant_expression!(container, lin_cost) + return lin_cost +end + +function _add_quadratic_term!( + container::OptimizationContainer, + ::T, + component::U, + q_terms::NTuple{2, Float64}, + expression_multiplier::Float64, + time_period::Int, +) where {T <: VariableType, U <: PSY.Component} + component_name = PSY.get_name(component) + @debug "$component_name Quadratic Variable Cost" _group = LOG_GROUP_COST_FUNCTIONS component_name + var = get_variable(container, T(), U)[component_name, time_period] + q_cost_ = var .^ 2 * q_terms[1] + var * q_terms[2] + q_cost = q_cost_ * expression_multiplier + add_to_objective_invariant_expression!(container, q_cost) + return q_cost +end + +################################################## +################## Fuel Cost ##################### +################################################## + +function _get_fuel_cost_value( + ::OptimizationContainer, + fuel_cost::Float64, + ::Int, +) + return fuel_cost +end + +function _get_fuel_cost_value( + container::OptimizationContainer, + fuel_cost::IS.TimeSeriesKey, + time_period::Int, +) + error("Not implemented yet fuel cost") + return fuel_cost +end diff --git a/src/devices_models/devices/common/objective_function/linear_curve.jl b/src/devices_models/devices/common/objective_function/linear_curve.jl new file mode 100644 index 0000000000..06d47c9b29 --- /dev/null +++ b/src/devices_models/devices/common/objective_function/linear_curve.jl @@ -0,0 +1,165 @@ +# Add proportional terms to objective function and expression +function _add_linearcurve_variable_term_to_model!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + proportional_term_per_unit::Float64, + time_period::Int, +) where {T <: VariableType} + resolution = get_resolution(container) + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR + linear_cost = _add_proportional_term!( + container, + T(), + component, + proportional_term_per_unit * dt, + time_period, + ) + add_to_expression!( + container, + ProductionCostExpression, + linear_cost, + component, + time_period, + ) + return +end + +# Dispatch for vector of proportional terms +function _add_linearcurve_variable_cost!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + proportional_terms_per_unit::Vector{Float64}, +) where {T <: VariableType} + for t in get_time_steps(container) + _add_linearcurve_variable_term_to_model!( + container, + T(), + component, + proportional_terms_per_unit[t], + t, + ) + end + return +end + +# Dispatch for scalar proportional terms +function _add_linearcurve_variable_cost!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + proportional_term_per_unit::Float64, +) where {T <: VariableType} + for t in get_time_steps(container) + _add_linearcurve_variable_term_to_model!( + container, + T(), + component, + proportional_term_per_unit, + t, + ) + end + return +end + +""" +Adds to the cost function cost terms for sum of variables with common factor to be used for cost expression for optimization_container model. + +# Arguments + + - container::OptimizationContainer : the optimization_container model built in PowerSimulations + - var_key::VariableKey: The variable name + - component_name::String: The component_name of the variable container + - cost_component::PSY.CostCurve{PSY.LinearCurve} : container for cost to be associated with variable +""" +function _add_variable_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + cost_function::PSY.CostCurve{PSY.LinearCurve}, + ::U, +) where {T <: VariableType, U <: AbstractDeviceFormulation} + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + value_curve = PSY.get_value_curve(cost_function) + power_units = PSY.get_power_units(cost_function) + cost_component = PSY.get_function_data(value_curve) + proportional_term = PSY.get_proportional_term(cost_component) + proportional_term_per_unit = get_proportional_cost_per_system_unit( + proportional_term, + power_units, + base_power, + device_base_power, + ) + multiplier = objective_function_multiplier(T(), U()) + _add_linearcurve_variable_cost!( + container, + T(), + component, + multiplier * proportional_term_per_unit, + ) + return +end + +function _add_fuel_linear_variable_cost!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + fuel_curve::Float64, + fuel_cost::Float64, +) where {T <: VariableType} + _add_linearcurve_variable_cost!(container, T(), component, fuel_curve * fuel_cost) +end + +function _add_fuel_linear_variable_cost!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + fuel_curve::Float64, + fuel_cost::IS.TimeSeriesKey, +) where {T <: VariableType} + error("Not implemented yet") + _add_linearcurve_variable_cost!(container, T(), component, fuel_curve) +end + +""" +Adds to the cost function cost terms for sum of variables with common factor to be used for cost expression for optimization_container model. + +# Arguments + + - container::OptimizationContainer : the optimization_container model built in PowerSimulations + - var_key::VariableKey: The variable name + - component_name::String: The component_name of the variable container + - cost_component::PSY.FuelCurve{PSY.LinearCurve} : container for cost to be associated with variable +""" +function _add_variable_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + cost_function::PSY.FuelCurve{PSY.LinearCurve}, + ::U, +) where {T <: VariableType, U <: AbstractDeviceFormulation} + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + value_curve = PSY.get_value_curve(cost_function) + power_units = PSY.get_power_units(cost_function) + cost_component = PSY.get_function_data(value_curve) + proportional_term = PSY.get_proportional_term(cost_component) + fuel_curve_per_unit = get_proportional_cost_per_system_unit( + proportional_term, + power_units, + base_power, + device_base_power, + ) + fuel_cost = PSY.get_fuel_cost(cost_function) + # Multiplier is not necessary here. There is no negative cost for fuel curves. + _add_fuel_linear_variable_cost!( + container, + T(), + component, + fuel_curve_per_unit, + fuel_cost, + ) + return +end diff --git a/src/devices_models/devices/common/objective_function/market_bid.jl b/src/devices_models/devices/common/objective_function/market_bid.jl new file mode 100644 index 0000000000..65f68a28c6 --- /dev/null +++ b/src/devices_models/devices/common/objective_function/market_bid.jl @@ -0,0 +1,495 @@ +################################################## +################# PWL Variables ################## +################################################## + +# For Market Bid +function _add_pwl_variables!( + container::OptimizationContainer, + ::Type{T}, + component_name::String, + time_period::Int, + cost_data::PSY.PiecewiseStepData, +) where {T <: PSY.Component} + var_container = lazy_container_addition!(container, PieceWiseLinearBlockOffer(), T) + # length(PiecewiseStepData) gets number of segments, here we want number of points + break_points = PSY.get_x_coords(cost_data) + pwlvars = Array{JuMP.VariableRef}(undef, length(break_points)) + for i in 1:(length(break_points) - 1) + pwlvars[i] = + var_container[(component_name, i, time_period)] = JuMP.@variable( + get_jump_model(container), + base_name = "PieceWiseLinearBlockOffer_$(component_name)_supply_{pwl_$(i), $time_period}", + lower_bound = 0.0, + ) + end + return pwlvars +end + +################################################## +################# PWL Constraints ################ +################################################## + +""" +Implement the constraints for PWL Block Offer variables. That is: + +```math +\\sum_{k\\in\\mathcal{K}} \\delta_{k,t} = p_t \\\\ +\\sum_{k\\in\\mathcal{K}} \\delta_{k,t} <= P_{k+1,t}^{max} - P_{k,t}^{max} +``` +""" +function _add_pwl_constraint!( + container::OptimizationContainer, + component::T, + ::U, + break_points::Vector{Float64}, + period::Int, +) where {T <: PSY.Component, U <: VariableType} + variables = get_variable(container, U(), T) + const_container = lazy_container_addition!( + container, + PieceWiseLinearBlockOfferConstraint(), + T, + axes(variables)..., + ) + len_cost_data = length(break_points) - 1 + jump_model = get_jump_model(container) + pwl_vars = get_variable(container, PieceWiseLinearBlockOffer(), T) + name = PSY.get_name(component) + const_container[name, period] = JuMP.@constraint( + jump_model, + variables[name, period] == + sum(pwl_vars[name, ix, period] for ix in 1:len_cost_data) + ) + + #= + const_upperbound_container = lazy_container_addition!( + container, + PieceWiseLinearUpperBoundConstraint(), + T, + axes(pwl_vars)...; + ) + =# + + # TODO: Parameter for this + for ix in 1:len_cost_data + JuMP.@constraint( + jump_model, + pwl_vars[name, ix, period] <= break_points[ix + 1] - break_points[ix] + ) + end + return +end + +################################################## +################ PWL Expressions ################# +################################################## + +function _get_pwl_cost_expression( + container::OptimizationContainer, + component::T, + time_period::Int, + cost_data::PSY.PiecewiseStepData, + multiplier::Float64, +) where {T <: PSY.Component} + name = PSY.get_name(component) + pwl_var_container = get_variable(container, PieceWiseLinearBlockOffer(), T) + gen_cost = JuMP.AffExpr(0.0) + cost_data = PSY.get_y_coords(cost_data) + for (i, cost) in enumerate(cost_data) + JuMP.add_to_expression!( + gen_cost, + cost * multiplier * pwl_var_container[(name, i, time_period)], + ) + end + return gen_cost +end + +function _get_pwl_cost_expression( + container::OptimizationContainer, + component::T, + time_period::Int, + cost_function::PSY.MarketBidCost, + ::PSY.PiecewiseStepData, + ::U, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} + incremental_curve = PSY.get_incremental_offer_curves(cost_function) + value_curve = PSY.get_value_curve(incremental_curve) + power_units = PSY.get_power_units(incremental_curve) + cost_component = PSY.get_function_data(value_curve) + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + cost_data_normalized = get_piecewise_incrementalcurve_per_system_unit( + cost_component, + power_units, + base_power, + device_base_power, + ) + resolution = get_resolution(container) + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR + return _get_pwl_cost_expression( + container, + component, + time_period, + cost_data_normalized, + dt, + ) +end + +#= +# For Market Bid +function _add_pwl_variables!( + container::OptimizationContainer, + ::Type{T}, + component_name::String, + time_period::Int, + cost_data::PSY.PiecewiseStepData, +) where {T <: PSY.Component} + var_container = lazy_container_addition!(container, PieceWiseLinearCostVariable(), T) + # length(PiecewiseStepData) gets number of segments, here we want number of points + pwlvars = Array{JuMP.VariableRef}(undef, length(cost_data) + 1) + for i in 1:(length(cost_data) + 1) + pwlvars[i] = + var_container[(component_name, i, time_period)] = JuMP.@variable( + get_jump_model(container), + base_name = "PieceWiseLinearCostVariable_$(component_name)_{pwl_$(i), $time_period}", + ) + end + return pwlvars +end + +# For Market Bid # +function _get_pwl_cost_expression( + container::OptimizationContainer, + component::T, + time_period::Int, + cost_data::PSY.PiecewiseStepData, + multiplier::Float64, +) where {T <: PSY.Component} + # TODO: This functions needs to be reimplemented for the new model. The code is repeated + # because the internals will be different + name = PSY.get_name(component) + pwl_var_container = get_variable(container, PieceWiseLinearCostVariable(), T) + gen_cost = JuMP.AffExpr(0.0) + cost_data = PSY.get_y_coords(cost_data) + for (i, cost) in enumerate(cost_data) + JuMP.add_to_expression!( + gen_cost, + cost * multiplier * pwl_var_container[(name, i, time_period)], + ) + end + return gen_cost +end + +function _add_pwl_term!( + container::OptimizationContainer, + component::T, + cost_data::AbstractVector{PSY.LinearFunctionData}, + ::U, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} + multiplier = objective_function_multiplier(U(), V()) + resolution = get_resolution(container) + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR + base_power = get_base_power(container) + # Re-scale breakpoints by Basepower + time_steps = get_time_steps(container) + cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) + for t in time_steps + proportional_value = + PSY.get_proportional_term(cost_data[t]) * multiplier * base_power * dt + cost_expressions[t] = + _add_proportional_term!(container, U(), component, proportional_value, t) + end + return cost_expressions +end +=# + +############################################### +######## MarketBidCost: Fixed Curves ########## +############################################### + +""" +Add PWL cost terms for data coming from the MarketBidCost +with a fixed incremental offer curve +""" +function _add_pwl_term!( + container::OptimizationContainer, + component::T, + cost_function::PSY.MarketBidCost, + ::PSY.CostCurve{PSY.PiecewiseIncrementalCurve}, + ::U, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} + name = PSY.get_name(component) + incremental_offer_curve = PSY.get_incremental_offer_curves(cost_function) + value_curve = PSY.get_value_curve(incremental_offer_curve) + cost_component = PSY.get_function_data(value_curve) + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + power_units = PSY.get_power_units(incremental_offer_curve) + + data = get_piecewise_incrementalcurve_per_system_unit( + cost_component, + power_units, + base_power, + device_base_power, + ) + + compact_status = validate_compact_pwl_data(component, data, base_power) + if !uses_compact_power(component, V()) && compact_status == COMPACT_PWL_STATUS.VALID + error( + "The data provided is not compatible with formulation $V. Use a formulation compatible with Compact Cost Functions", + ) + # data = _convert_to_full_variable_cost(data, component) + elseif uses_compact_power(component, V()) && compact_status != COMPACT_PWL_STATUS.VALID + @warn( + "The cost data provided is not in compact form. Will attempt to convert. Errors may occur." + ) + data = convert_to_compact_variable_cost(data) + else + @debug uses_compact_power(component, V()) compact_status name T V + end + + cost_is_convex = PSY.is_convex(data) + if !cost_is_convex + error("MarketBidCost for component $(name) is non-convex") + end + + break_points = PSY.get_x_coords(data) + time_steps = get_time_steps(container) + pwl_cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) + for t in time_steps + _add_pwl_variables!(container, T, name, t, data) + _add_pwl_constraint!(container, component, U(), break_points, t) + pwl_cost = + _get_pwl_cost_expression(container, component, t, cost_function, data, U(), V()) + pwl_cost_expressions[t] = pwl_cost + end + return pwl_cost_expressions +end + +#= +""" +Add PWL cost terms for data coming from the MarketBidCost +with a timeseries incremental offer curve +""" +function _add_pwl_term!( + container::OptimizationContainer, + component::T, + cost_function::PSY.MarketBidCost, + ::PSY.TimeSeriesKey, + ::U, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} + name = PSY.get_name(component) + value_curve = PSY.get_value_curve(incremental_offer_curve) + cost_component = PSY.get_function_data(value_curve) + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + power_units = PSY.get_power_units(cost_function) + + data = get_piecewise_incrementalcurve_per_system_unit( + cost_component, + power_units, + base_power, + device_base_power, + ) + time_steps = get_time_steps(container) + pwl_cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) + sos_val = _get_sos_value(container, V, component) + for t in time_steps + # Run checks in every time step because each time step has a PWL cost function + data = cost_data[t] + compact_status = validate_compact_pwl_data(component, data, base_power) + if !uses_compact_power(component, V()) && compact_status == COMPACT_PWL_STATUS.VALID + error( + "The data provided is not compatible with formulation $V. Use a formulation compatible with Compact Cost Functions", + ) + # data = _convert_to_full_variable_cost(data, component) + elseif uses_compact_power(component, V()) && + compact_status != COMPACT_PWL_STATUS.VALID + @warn( + "The cost data provided is not in compact form. Will attempt to convert. Errors may occur." + ) + data = convert_to_compact_variable_cost(data) + else + @debug uses_compact_power(component, V()) compact_status name T V + end + cost_is_convex = PSY.is_convex(data) + break_points = PSY.get_x_coords(data) ./ base_power # TODO should this be get_x_lengths/get_breakpoint_upper_bounds? + _add_pwl_variables!(container, T, name, t, data) + _add_pwl_constraint!(container, component, U(), break_points, sos_val, t) + if !cost_is_convex + _add_pwl_sos_constraint!(container, component, U(), break_points, sos_val, t) + end + pwl_cost = + _get_pwl_cost_expression(container, component, t, data, multiplier * dt) + pwl_cost_expressions[t] = pwl_cost + end + return pwl_cost_expressions +end + +function _add_pwl_term!( + container::OptimizationContainer, + component::T, + cost_data::AbstractVector{PSY.PiecewiseStepData}, + ::U, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractServiceFormulation} + multiplier = objective_function_multiplier(U(), V()) + resolution = get_resolution(container) + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR + base_power = get_base_power(container) + # Re-scale breakpoints by Basepower + name = PSY.get_name(component) + time_steps = get_time_steps(container) + pwl_cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) + sos_val = _get_sos_value(container, V, component) + for t in time_steps + data = cost_data[t] + break_points = PSY.get_x_coords(data) ./ base_power + _add_pwl_variables!(container, T, name, t, data) + _add_pwl_constraint!(container, component, U(), break_points, sos_val, t) + _add_pwl_sos_constraint!(container, component, U(), break_points, sos_val, t) + pwl_cost = _get_pwl_cost_expression(container, component, t, data, multiplier * dt) + pwl_cost_expressions[t] = pwl_cost + end + return pwl_cost_expressions +end +=# +############################################################ +######## MarketBidCost: PiecewiseIncrementalCurve ########## +############################################################ + +""" +Creates piecewise linear market bid function using a sum of variables and expression for market participants. +Decremental offers are not accepted for most components, except Storage systems and loads. + +# Arguments + + - container::OptimizationContainer : the optimization_container model built in PowerSimulations + - var_key::VariableKey: The variable name + - component_name::String: The component_name of the variable container + - cost_function::MarketBidCost : container for market bid cost +""" +function _add_variable_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + cost_function::PSY.MarketBidCost, + ::U, +) where {T <: VariableType, U <: AbstractDeviceFormulation} + component_name = PSY.get_name(component) + @debug "Market Bid" _group = LOG_GROUP_COST_FUNCTIONS component_name + time_steps = get_time_steps(container) + initial_time = get_initial_time(container) + incremental_cost_curves = PSY.get_incremental_offer_curves(cost_function) + decremental_cost_curves = PSY.get_decremental_offer_curves(cost_function) + if isnothing(decremental_cost_curves) + error("Component $(component_name) is not allowed to participate as a demand.") + end + #= + variable_cost_forecast = PSY.get_variable_cost( + component, + op_cost; + start_time = initial_time, + len = length(time_steps), + ) + variable_cost_forecast_values = TimeSeries.values(variable_cost_forecast) + parameter_container = _get_cost_function_parameter_container( + container, + CostFunctionParameter(), + component, + T(), + U(), + eltype(variable_cost_forecast_values), + ) + =# + pwl_cost_expressions = + _add_pwl_term!( + container, + component, + cost_function, + incremental_cost_curves, + T(), + U(), + ) + jump_model = get_jump_model(container) + for t in time_steps + #= + set_multiplier!( + parameter_container, + # Using 1.0 here since we want to reuse the existing code that adds the mulitpler + # of base power times the time delta. + 1.0, + component_name, + t, + ) + set_parameter!( + parameter_container, + jump_model, + variable_cost_forecast_values[t], + component_name, + t, + ) + =# + add_to_expression!( + container, + ProductionCostExpression, + pwl_cost_expressions[t], + component, + t, + ) + add_to_objective_variant_expression!(container, pwl_cost_expressions[t]) + end + + # Service Cost Bid + #= + ancillary_services = PSY.get_ancillary_service_offers(op_cost) + for service in ancillary_services + _add_service_bid_cost!(container, component, service) + end + =# + return +end + +function _add_service_bid_cost!( + container::OptimizationContainer, + component::PSY.Component, + service::T, +) where {T <: PSY.Reserve{<:PSY.ReserveDirection}} + time_steps = get_time_steps(container) + initial_time = get_initial_time(container) + base_power = get_base_power(container) + forecast_data = PSY.get_services_bid( + component, + PSY.get_operation_cost(component), + service; + start_time = initial_time, + len = length(time_steps), + ) + forecast_data_values = PSY.get_cost.(TimeSeries.values(forecast_data)) + # Single Price Bid + if eltype(forecast_data_values) == Float64 + data_values = forecast_data_values + # Single Price/Quantity Bid + elseif eltype(forecast_data_values) == Vector{NTuple{2, Float64}} + data_values = [v[1][1] for v in forecast_data_values] + else + error("$(eltype(forecast_data_values)) not supported for MarketBidCost") + end + + reserve_variable = + get_variable(container, ActivePowerReserveVariable(), T, PSY.get_name(service)) + component_name = PSY.get_name(component) + for t in time_steps + add_to_objective_invariant_expression!( + container, + data_values[t] * base_power * reserve_variable[component_name, t], + ) + end + return +end + +function _add_service_bid_cost!(::OptimizationContainer, ::PSY.Component, ::PSY.Service) end diff --git a/src/devices_models/devices/common/objective_function/piecewise_linear.jl b/src/devices_models/devices/common/objective_function/piecewise_linear.jl new file mode 100644 index 0000000000..497584ade9 --- /dev/null +++ b/src/devices_models/devices/common/objective_function/piecewise_linear.jl @@ -0,0 +1,521 @@ +################################################## +################# SOS Methods #################### +################################################## + +function _get_sos_value( + container::OptimizationContainer, + ::Type{V}, + component::T, +) where {T <: PSY.Component, V <: AbstractDeviceFormulation} + if has_container_key(container, OnStatusParameter, T) + sos_val = SOSStatusVariable.PARAMETER + else + sos_val = sos_status(component, V()) + end + return sos_val +end + +function _get_sos_value( + container::OptimizationContainer, + ::Type{V}, + component::T, +) where {T <: PSY.Component, V <: AbstractServiceFormulation} + return SOSStatusVariable.NO_VARIABLE +end + +################################################## +################# PWL Variables ################## +################################################## + +# This cases bounds the data by 1 - 0 +function _add_pwl_variables!( + container::OptimizationContainer, + ::Type{T}, + component_name::String, + time_period::Int, + cost_data::PSY.PiecewiseLinearData, +) where {T <: PSY.Component} + var_container = lazy_container_addition!(container, PieceWiseLinearCostVariable(), T) + # length(PiecewiseStepData) gets number of segments, here we want number of points + pwlvars = Array{JuMP.VariableRef}(undef, length(cost_data) + 1) + for i in 1:(length(cost_data) + 1) + pwlvars[i] = + var_container[(component_name, i, time_period)] = JuMP.@variable( + get_jump_model(container), + base_name = "PieceWiseLinearCostVariable_$(component_name)_{pwl_$(i), $time_period}", + lower_bound = 0.0, + upper_bound = 1.0 + ) + end + return pwlvars +end + +################################################## +################# PWL Constraints ################ +################################################## + +""" +Implement the constraints for PWL variables. That is: + +```math +\\sum_{k\\in\\mathcal{K}} P_k^{max} \\delta_{k,t} = p_t \\\\ +\\sum_{k\\in\\mathcal{K}} \\delta_{k,t} = on_t +``` +""" +function _add_pwl_constraint!( + container::OptimizationContainer, + component::T, + ::U, + break_points::Vector{Float64}, + sos_status::SOSStatusVariable, + period::Int, +) where {T <: PSY.Component, U <: VariableType} + variables = get_variable(container, U(), T) + const_container = lazy_container_addition!( + container, + PieceWiseLinearCostConstraint(), + T, + axes(variables)..., + ) + len_cost_data = length(break_points) + jump_model = get_jump_model(container) + pwl_vars = get_variable(container, PieceWiseLinearCostVariable(), T) + name = PSY.get_name(component) + const_container[name, period] = JuMP.@constraint( + jump_model, + variables[name, period] == + sum(pwl_vars[name, ix, period] * break_points[ix] for ix in 1:len_cost_data) + ) + + if sos_status == SOSStatusVariable.NO_VARIABLE + bin = 1.0 + @debug "Using Piecewise Linear cost function but no variable/parameter ref for ON status is passed. Default status will be set to online (1.0)" _group = + LOG_GROUP_COST_FUNCTIONS + + elseif sos_status == SOSStatusVariable.PARAMETER + param = get_default_on_parameter(component) + bin = get_parameter(container, param, T).parameter_array[name, period] + @debug "Using Piecewise Linear cost function with parameter OnStatusParameter, $T" _group = + LOG_GROUP_COST_FUNCTIONS + elseif sos_status == SOSStatusVariable.VARIABLE + var = get_default_on_variable(component) + bin = get_variable(container, var, T)[name, period] + @debug "Using Piecewise Linear cost function with variable OnVariable $T" _group = + LOG_GROUP_COST_FUNCTIONS + else + @assert false + end + + const_normalization_container = lazy_container_addition!( + container, + PieceWiseLinearCostConstraint(), + T, + axes(variables)...; + meta = "normalization", + ) + + const_normalization_container[name, period] = JuMP.@constraint( + jump_model, + sum(pwl_vars[name, i, period] for i in 1:len_cost_data) == bin + ) + return +end + +""" +Implement the SOS for PWL variables. That is: + +```math +\\{\\delta_{i,t}, ..., \\delta_{k,t}\\} \\in \\text{SOS}_2 +``` +""" +function _add_pwl_sos_constraint!( + container::OptimizationContainer, + component::T, + ::U, + break_points::Vector{Float64}, + sos_status::SOSStatusVariable, + period::Int, +) where {T <: PSY.Component, U <: VariableType} + name = PSY.get_name(component) + @warn( + "The cost function provided for $(name) is not compatible with a linear PWL cost function. + An SOS-2 formulation will be added to the model. This will result in additional binary variables." + ) + + jump_model = get_jump_model(container) + pwl_vars = get_variable(container, PieceWiseLinearCostVariable(), T) + bp_count = length(break_points) + pwl_vars_subset = [pwl_vars[name, i, period] for i in 1:bp_count] + JuMP.@constraint(jump_model, pwl_vars_subset in MOI.SOS2(collect(1:bp_count))) + return +end + +################################################## +################ PWL Expressions ################# +################################################## + +function _get_pwl_cost_expression( + container::OptimizationContainer, + component::T, + time_period::Int, + cost_data::PSY.PiecewiseLinearData, + multiplier::Float64, +) where {T <: PSY.Component} + name = PSY.get_name(component) + pwl_var_container = get_variable(container, PieceWiseLinearCostVariable(), T) + gen_cost = JuMP.AffExpr(0.0) + cost_data = PSY.get_y_coords(cost_data) + for (i, cost) in enumerate(cost_data) + JuMP.add_to_expression!( + gen_cost, + cost * multiplier * pwl_var_container[(name, i, time_period)], + ) + end + return gen_cost +end + +function _get_pwl_cost_expression( + container::OptimizationContainer, + component::T, + time_period::Int, + cost_function::PSY.CostCurve{PSY.PiecewisePointCurve}, + ::U, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} + value_curve = PSY.get_value_curve(cost_function) + power_units = PSY.get_power_units(cost_function) + cost_component = PSY.get_function_data(value_curve) + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + cost_data_normalized = get_piecewise_pointcurve_per_system_unit( + cost_component, + power_units, + base_power, + device_base_power, + ) + multiplier = objective_function_multiplier(U(), V()) + resolution = get_resolution(container) + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR + return _get_pwl_cost_expression( + container, + component, + time_period, + cost_data_normalized, + multiplier * dt, + ) +end + +function _get_pwl_cost_expression( + container::OptimizationContainer, + component::T, + time_period::Int, + cost_function::PSY.FuelCurve{PSY.PiecewisePointCurve}, + ::U, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} + value_curve = PSY.get_value_curve(cost_function) + power_units = PSY.get_power_units(cost_function) + cost_component = PSY.get_function_data(value_curve) + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + cost_data_normalized = get_piecewise_pointcurve_per_system_unit( + cost_component, + power_units, + base_power, + device_base_power, + ) + fuel_cost = PSY.get_fuel_cost(cost_function) + fuel_cost_value = _get_fuel_cost_value( + container, + fuel_cost, + time_period, + ) + # Multiplier is not necessary here. There is no negative cost for fuel curves. + resolution = get_resolution(container) + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR + return _get_pwl_cost_expression( + container, + component, + time_period, + cost_data_normalized, + dt * fuel_cost_value, + ) +end + +################################################## +######## CostCurve: PiecewisePointCurve ########## +################################################## + +""" +Add PWL cost terms for data coming from a PiecewisePointCurve +""" +function _add_pwl_term!( + container::OptimizationContainer, + component::T, + cost_function::Union{ + PSY.CostCurve{PSY.PiecewisePointCurve}, + PSY.FuelCurve{PSY.PiecewisePointCurve}, + }, + ::U, + ::V, +) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} + # multiplier = objective_function_multiplier(U(), V()) + name = PSY.get_name(component) + value_curve = PSY.get_value_curve(cost_function) + cost_component = PSY.get_function_data(value_curve) + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + power_units = PSY.get_power_units(cost_function) + + # Normalize data + data = get_piecewise_pointcurve_per_system_unit( + cost_component, + power_units, + base_power, + device_base_power, + ) + + if all(iszero.((point -> point.y).(PSY.get_points(data)))) # TODO I think this should have been first. before? + @debug "All cost terms for component $(name) are 0.0" _group = + LOG_GROUP_COST_FUNCTIONS + return + end + + compact_status = validate_compact_pwl_data(component, data, base_power) + if !uses_compact_power(component, V()) && compact_status == COMPACT_PWL_STATUS.VALID + error( + "The data provided is not compatible with formulation $V. Use a formulation compatible with Compact Cost Functions", + ) + # data = _convert_to_full_variable_cost(data, component) + elseif uses_compact_power(component, V()) && compact_status != COMPACT_PWL_STATUS.VALID + @warn( + "The cost data provided is not in compact form. Will attempt to convert. Errors may occur." + ) + data = convert_to_compact_variable_cost(data) + else + @debug uses_compact_power(component, V()) compact_status name T V + end + + cost_is_convex = PSY.is_convex(data) + break_points = PSY.get_x_coords(data) + time_steps = get_time_steps(container) + pwl_cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) + sos_val = _get_sos_value(container, V, component) + for t in time_steps + _add_pwl_variables!(container, T, name, t, data) + _add_pwl_constraint!(container, component, U(), break_points, sos_val, t) + if !cost_is_convex + _add_pwl_sos_constraint!(container, component, U(), break_points, sos_val, t) + end + pwl_cost = + _get_pwl_cost_expression(container, component, t, cost_function, U(), V()) + pwl_cost_expressions[t] = pwl_cost + end + return pwl_cost_expressions +end + +""" +Add PWL cost terms for data coming from a PiecewisePointCurve for ThermalDispatchNoMin formulation +""" +function _add_pwl_term!( + container::OptimizationContainer, + component::T, + cost_function::Union{ + PSY.CostCurve{PSY.PiecewisePointCurve}, + PSY.FuelCurve{PSY.PiecewisePointCurve}, + }, + ::U, + ::V, +) where {T <: PSY.ThermalGen, U <: VariableType, V <: ThermalDispatchNoMin} + name = PSY.get_name(component) + value_curve = PSY.get_value_curve(cost_function) + cost_component = PSY.get_function_data(value_curve) + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + power_units = PSY.get_power_units(cost_function) + + # Normalize data + data = get_piecewise_pointcurve_per_system_unit( + cost_component, + power_units, + base_power, + device_base_power, + ) + @debug "PWL cost function detected for device $(name) using $V" + slopes = PSY.get_slopes(data) + if any(slopes .< 0) || !PSY.is_convex(data) + throw( + IS.InvalidValue( + "The PWL cost data provided for generator $(name) is not compatible with $U.", + ), + ) + end + + if validate_compact_pwl_data(component, data, base_power) == COMPACT_PWL_STATUS.VALID + error("The data provided is not compatible with formulation $V. \\ + Use a formulation compatible with Compact Cost Functions") + end + + if slopes[1] != 0.0 + @debug "PWL has no 0.0 intercept for generator $(component_name)" + # adds a first intercept a x = 0.0 and y below the intercept of the first tuple to make convex equivalent + intercept_point = (x = 0.0, y = first(data).y - COST_EPSILON) + data = PSY.PiecewiseLinearData(vcat(intercept_point, get_points(data))) + @assert PSY.is_convex(slopes) + end + + time_steps = get_time_steps(container) + pwl_cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) + break_points = PSY.get_x_coords(data) + sos_val = _get_sos_value(container, V, component) + for t in time_steps + _add_pwl_variables!(container, T, component_name, t, data) + _add_pwl_constraint!(container, component, U(), break_points, sos_val, t) + pwl_cost = + _get_pwl_cost_expression(container, component, t, cost_function, U(), V()) + pwl_cost_expressions[t] = pwl_cost + end + return pwl_cost_expressions +end + +""" +Creates piecewise linear cost function using a sum of variables and expression with sign and time step included. + +# Arguments + + - container::OptimizationContainer : the optimization_container model built in PowerSimulations + - var_key::VariableKey: The variable name + - component_name::String: The component_name of the variable container + - cost_function::PSY.CostCurve{PSY.PiecewisePointCurve}: container for piecewise linear cost +""" +function _add_variable_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + cost_function::Union{ + PSY.CostCurve{PSY.PiecewisePointCurve}, + PSY.FuelCurve{PSY.PiecewisePointCurve}, + }, + ::U, +) where {T <: VariableType, U <: AbstractDeviceFormulation} + component_name = PSY.get_name(component) + @debug "PWL Variable Cost" _group = LOG_GROUP_COST_FUNCTIONS component_name + # If array is full of tuples with zeros return 0.0 + value_curve = PSY.get_value_curve(cost_function) + cost_component = PSY.get_function_data(value_curve) + if all(iszero.((point -> point.y).(PSY.get_points(cost_component)))) # TODO I think this should have been first. before? + @debug "All cost terms for component $(component_name) are 0.0" _group = + LOG_GROUP_COST_FUNCTIONS + return + end + pwl_cost_expressions = + _add_pwl_term!(container, component, cost_function, T(), U()) + for t in get_time_steps(container) + add_to_expression!( + container, + ProductionCostExpression, + pwl_cost_expressions[t], + component, + t, + ) + add_to_objective_invariant_expression!(container, pwl_cost_expressions[t]) + end + return +end + +################################################## +###### CostCurve: PiecewiseIncrementalCurve ###### +######### and PiecewiseAverageCurve ############## +################################################## + +""" +Creates piecewise linear cost function using a sum of variables and expression with sign and time step included. + +# Arguments + + - container::OptimizationContainer : the optimization_container model built in PowerSimulations + - var_key::VariableKey: The variable name + - component_name::String: The component_name of the variable container + - cost_function::PSY.Union{PSY.CostCurve{PSY.PiecewiseIncrementalCurve}, PSY.CostCurve{PSY.PiecewiseAverageCurve}}: container for piecewise linear cost +""" +function _add_variable_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + cost_function::V, + ::U, +) where { + T <: VariableType, + V <: Union{ + PSY.CostCurve{PSY.PiecewiseIncrementalCurve}, + PSY.CostCurve{PSY.PiecewiseAverageCurve}, + }, + U <: AbstractDeviceFormulation, +} + # Create new PiecewisePointCurve + value_curve = PSY.get_value_curve(cost_function) + power_units = PSY.get_power_units(cost_function) + pointbased_value_curve = PSY.InputOutputCurve(value_curve) + pointbased_cost_function = + PSY.CostCurve(; value_curve = pointbased_value_curve, power_units = power_units) + # Call method for PiecewisePointCurve + _add_variable_cost_to_objective!( + container, + T(), + component, + pointbased_cost_function, + U(), + ) + return +end + +################################################## +###### FuelCurve: PiecewiseIncrementalCurve ###### +######### and PiecewiseAverageCurve ############## +################################################## + +""" +Creates piecewise linear fuel cost function using a sum of variables and expression with sign and time step included. + +# Arguments + + - container::OptimizationContainer : the optimization_container model built in PowerSimulations + - var_key::VariableKey: The variable name + - component_name::String: The component_name of the variable container + - cost_function::PSY.Union{PSY.FuelCurve{PSY.PiecewiseIncrementalCurve}, PSY.FuelCurve{PSY.PiecewiseAverageCurve}}: container for piecewise linear cost +""" +function _add_variable_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + cost_function::V, + ::U, +) where { + T <: VariableType, + V <: Union{ + PSY.FuelCurve{PSY.PiecewiseIncrementalCurve}, + PSY.FuelCurve{PSY.PiecewiseAverageCurve}, + }, + U <: AbstractDeviceFormulation, +} + # Create new PiecewisePointCurve + value_curve = PSY.get_value_curve(cost_function) + power_units = PSY.get_power_units(cost_function) + fuel_cost = PSY.get_fuel_cost(cost_function) + pointbased_value_curve = PSY.InputOutputCurve(value_curve) + pointbased_cost_function = + PSY.FuelCurve(; + value_curve = pointbased_value_curve, + power_units = power_units, + fuel_cost = fuel_cost, + ) + # Call method for PiecewisePointCurve + _add_variable_cost_to_objective!( + container, + T(), + component, + pointbased_cost_function, + U(), + ) + return +end diff --git a/src/devices_models/devices/common/objective_function/quadratic_curve.jl b/src/devices_models/devices/common/objective_function/quadratic_curve.jl new file mode 100644 index 0000000000..4f5a8aea09 --- /dev/null +++ b/src/devices_models/devices/common/objective_function/quadratic_curve.jl @@ -0,0 +1,252 @@ +# Add proportional terms to objective function and expression +function _add_quadraticcurve_variable_term_to_model!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + proportional_term_per_unit::Float64, + quadratic_term_per_unit::Float64, + time_period::Int, +) where {T <: VariableType} + resolution = get_resolution(container) + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR + if quadratic_term_per_unit >= eps() + cost_term = _add_quadratic_term!( + container, + T(), + component, + (quadratic_term_per_unit, proportional_term_per_unit), + dt, + time_period, + ) + else + cost_term = _add_proportional_term!( + container, + T(), + component, + proportional_term_per_unit * dt, + time_period, + ) + end + add_to_expression!( + container, + ProductionCostExpression, + cost_term, + component, + time_period, + ) + return +end + +# Dispatch for vector proportional/quadratic terms +function _add_quadraticcurve_variable_cost!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + proportional_term_per_unit::Vector{Float64}, + quadratic_term_per_unit::Vector{Float64}, +) where {T <: VariableType} + for t in get_time_steps(container) + _add_quadraticcurve_variable_term_to_model!( + container, + T(), + component, + proportional_term_per_unit[t], + quadratic_term_per_unit[t], + t, + ) + end + return +end + +# Dispatch for scalar proportional/quadratic terms +function _add_quadraticcurve_variable_cost!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + proportional_term_per_unit::Float64, + quadratic_term_per_unit::Float64, +) where {T <: VariableType} + for t in get_time_steps(container) + _add_quadraticcurve_variable_term_to_model!( + container, + T(), + component, + proportional_term_per_unit, + quadratic_term_per_unit, + t, + ) + end + return +end + +@doc raw""" +Adds to the cost function cost terms for sum of variables with common factor to be used for cost expression for optimization_container model. + +# Equation + +``` gen_cost = dt*sign*(sum(variable.^2)*cost_data[1] + sum(variable)*cost_data[2]) ``` + +# LaTeX + +`` cost = dt\times sign (sum_{i\in I} c_1 v_i^2 + sum_{i\in I} c_2 v_i ) `` + +for quadratic factor large enough. If the first term of the quadratic objective is 0.0, adds a +linear cost term `sum(variable)*cost_data[2]` + +# Arguments + +* container::OptimizationContainer : the optimization_container model built in PowerSimulations +* var_key::VariableKey: The variable name +* component_name::String: The component_name of the variable container +* cost_component::PSY.CostCurve{PSY.QuadraticCurve} : container for quadratic factors +""" +function _add_variable_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + cost_function::PSY.CostCurve{PSY.QuadraticCurve}, + ::U, +) where {T <: VariableType, U <: AbstractDeviceFormulation} + multiplier = objective_function_multiplier(T(), U()) + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + value_curve = PSY.get_value_curve(cost_function) + power_units = PSY.get_power_units(cost_function) + cost_component = PSY.get_function_data(value_curve) + quadratic_term = PSY.get_quadratic_term(cost_component) + proportional_term = PSY.get_proportional_term(cost_component) + proportional_term_per_unit = get_proportional_cost_per_system_unit( + proportional_term, + power_units, + base_power, + device_base_power, + ) + quadratic_term_per_unit = get_quadratic_cost_per_system_unit( + quadratic_term, + power_units, + base_power, + device_base_power, + ) + _add_quadraticcurve_variable_cost!( + container, + T(), + component, + multiplier * proportional_term_per_unit, + multiplier * quadratic_term_per_unit, + ) + return +end + +function _add_variable_cost_to_objective!( + ::OptimizationContainer, + ::T, + component::PSY.Component, + cost_function::PSY.CostCurve{PSY.QuadraticCurve}, + ::U, +) where { + T <: PowerAboveMinimumVariable, + U <: Union{AbstractCompactUnitCommitment, ThermalCompactDispatch}, +} + throw( + IS.ConflictingInputsError( + "Quadratic Cost Curves are not allowed for Compact formulations", + ), + ) + return +end + +function _add_fuel_quadratic_variable_cost!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + proportional_fuel_curve::Float64, + quadratic_fuel_curve::Float64, + fuel_cost::Float64, +) where {T <: VariableType} + _add_quadraticcurve_variable_cost!( + container, + T(), + component, + proportional_fuel_curve * fuel_cost, + quadratic_fuel_curve * fuel_cost, + ) +end + +function _add_fuel_quadratic_variable_cost!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + proportional_fuel_curve::Float64, + quadratic_fuel_curve::Float64, + fuel_cost::IS.TimeSeriesKey, +) where {T <: VariableType} + error("Not implemented yet") + _add_quadraticcurve_variable_cost!( + container, + T(), + component, + proportional_fuel_curve, + quadratic_fuel_curve, + ) +end + +@doc raw""" +Adds to the cost function cost terms for sum of variables with common factor to be used for cost expression for optimization_container model. + +# Equation + +``` gen_cost = dt*(sum(variable.^2)*cost_data[1]*fuel_cost + sum(variable)*cost_data[2]*fuel_cost) ``` + +# LaTeX + +`` cost = dt\times (sum_{i\in I} c_f c_1 v_i^2 + sum_{i\in I} c_f c_2 v_i ) `` + +for quadratic factor large enough. If the first term of the quadratic objective is 0.0, adds a +linear cost term `sum(variable)*cost_data[2]` + +# Arguments + +* container::OptimizationContainer : the optimization_container model built in PowerSimulations +* var_key::VariableKey: The variable name +* component_name::String: The component_name of the variable container +* cost_component::PSY.FuelCurve{PSY.QuadraticCurve} : container for quadratic factors +""" +function _add_variable_cost_to_objective!( + container::OptimizationContainer, + ::T, + component::PSY.Component, + cost_function::PSY.FuelCurve{PSY.QuadraticCurve}, + ::U, +) where {T <: VariableType, U <: AbstractDeviceFormulation} + multiplier = objective_function_multiplier(T(), U()) + base_power = get_base_power(container) + device_base_power = PSY.get_base_power(component) + value_curve = PSY.get_value_curve(cost_function) + power_units = PSY.get_power_units(cost_function) + cost_component = PSY.get_function_data(value_curve) + quadratic_term = PSY.get_quadratic_term(cost_component) + proportional_term = PSY.get_proportional_term(cost_component) + proportional_term_per_unit = get_proportional_cost_per_system_unit( + proportional_term, + power_units, + base_power, + device_base_power, + ) + quadratic_term_per_unit = get_quadratic_cost_per_system_unit( + quadratic_term, + power_units, + base_power, + device_base_power, + ) + fuel_cost = PSY.get_fuel_cost(cost_function) + # Multiplier is not necessary here. There is no negative cost for fuel curves. + _add_fuel_quadratic_variable_cost!( + container, + T(), + component, + multiplier * proportional_term_per_unit, + multiplier * quadratic_term_per_unit, + fuel_cost, + ) + return +end diff --git a/src/devices_models/devices/common/objective_functions.jl b/src/devices_models/devices/common/objective_functions.jl deleted file mode 100644 index b21004b6b9..0000000000 --- a/src/devices_models/devices/common/objective_functions.jl +++ /dev/null @@ -1,788 +0,0 @@ -function add_variable_cost!( - container::OptimizationContainer, - ::U, - devices::IS.FlattenIteratorWrapper{T}, - ::V, -) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} - for d in devices - op_cost_data = PSY.get_operation_cost(d) - _add_variable_cost_to_objective!(container, U(), d, op_cost_data, V()) - end - return -end - -function add_shut_down_cost!( - container::OptimizationContainer, - ::U, - devices::IS.FlattenIteratorWrapper{T}, - ::V, -) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} - multiplier = objective_function_multiplier(U(), V()) - for d in devices - op_cost_data = PSY.get_operation_cost(d) - cost_term = shut_down_cost(op_cost_data, d, V()) - iszero(cost_term) && continue - for t in get_time_steps(container) - _add_proportional_term!(container, U(), d, cost_term * multiplier, t) - end - end - return -end - -function add_proportional_cost!( - container::OptimizationContainer, - ::U, - devices::IS.FlattenIteratorWrapper{T}, - ::V, -) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} - multiplier = objective_function_multiplier(U(), V()) - for d in devices - op_cost_data = PSY.get_operation_cost(d) - cost_term = proportional_cost(op_cost_data, U(), d, V()) - iszero(cost_term) && continue - for t in get_time_steps(container) - _add_proportional_term!(container, U(), d, cost_term * multiplier, t) - end - end - return -end - -function add_proportional_cost!( - container::OptimizationContainer, - ::U, - devices::IS.FlattenIteratorWrapper{T}, - ::V, -) where {T <: PSY.ThermalGen, U <: OnVariable, V <: AbstractCompactUnitCommitment} - multiplier = objective_function_multiplier(U(), V()) - for d in devices - op_cost_data = PSY.get_operation_cost(d) - cost_term = proportional_cost(op_cost_data, U(), d, V()) - iszero(cost_term) && continue - for t in get_time_steps(container) - exp = _add_proportional_term!(container, U(), d, cost_term * multiplier, t) - add_to_expression!(container, ProductionCostExpression, exp, d, t) - end - end - return -end - -function _add_variable_cost_to_objective!( - container::OptimizationContainer, - ::T, - component::PSY.Component, - op_cost::PSY.OperationalCost, - ::U, -) where {T <: VariableType, U <: AbstractDeviceFormulation} - variable_cost_data = variable_cost(op_cost, T(), component, U()) - _add_variable_cost_to_objective!(container, T(), component, variable_cost_data, U()) - return -end - -function add_start_up_cost!( - container::OptimizationContainer, - ::U, - devices::IS.FlattenIteratorWrapper{T}, - ::V, -) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} - for d in devices - op_cost_data = PSY.get_operation_cost(d) - _add_start_up_cost_to_objective!(container, U(), d, op_cost_data, V()) - end - return -end - -function _add_start_up_cost_to_objective!( - container::OptimizationContainer, - ::T, - component::PSY.ThermalGen, - op_cost::Union{PSY.ThermalGenerationCost, PSY.MarketBidCost}, - ::U, -) where {T <: VariableType, U <: AbstractDeviceFormulation} - cost_term = start_up_cost(op_cost, component, U()) - iszero(cost_term) && return - multiplier = objective_function_multiplier(T(), U()) - for t in get_time_steps(container) - _add_proportional_term!(container, T(), component, cost_term * multiplier, t) - end - return -end - -const MULTI_START_COST_MAP = Dict{DataType, Int}( - HotStartVariable => 1, - WarmStartVariable => 2, - ColdStartVariable => 3, -) - -function _add_start_up_cost_to_objective!( - container::OptimizationContainer, - ::T, - component::PSY.ThermalMultiStart, - op_cost::PSY.ThermalGenerationCost, - ::U, -) where {T <: VariableType, U <: ThermalMultiStartUnitCommitment} - cost_terms = start_up_cost(op_cost, component, U()) - cost_term = cost_terms[MULTI_START_COST_MAP[T]] - iszero(cost_term) && return - multiplier = objective_function_multiplier(T(), U()) - for t in get_time_steps(container) - _add_proportional_term!(container, T(), component, cost_term * multiplier, t) - end - return -end - -function _get_cost_function_parameter_container( - container::OptimizationContainer, - ::S, - component::T, - ::U, - ::V, - cost_type::Type{W}, -) where { - S <: ObjectiveFunctionParameter, - T <: PSY.Component, - U <: VariableType, - V <: Union{AbstractDeviceFormulation, AbstractServiceFormulation}, - W, -} - if has_container_key(container, S, T) - return get_parameter(container, S(), T) - else - container_axes = axes(get_variable(container, U(), T)) - if has_container_key(container, OnStatusParameter, T) - sos_val = SOSStatusVariable.PARAMETER - else - sos_val = sos_status(component, V()) - end - return add_param_container!( - container, - S(), - T, - U, - sos_val, - uses_compact_power(component, V()), - W, - container_axes..., - ) - end -end - -""" -Adds to the cost function cost terms for sum of variables with common factor to be used for cost expression for optimization_container model. - -# Arguments - - - container::OptimizationContainer : the optimization_container model built in PowerSimulations - - var_key::VariableKey: The variable name - - component_name::String: The component_name of the variable container - - cost_component::PSY.CostCurve{PSY.LinearFunctionData} : container for cost to be associated with variable -""" -function _add_variable_cost_to_objective!( - container::OptimizationContainer, - ::T, - component::PSY.Component, - cost_function::PSY.CostCurve{PSY.LinearCurve}, - ::U, -) where {T <: VariableType, U <: AbstractDeviceFormulation} - multiplier = objective_function_multiplier(T(), U()) - base_power = get_base_power(container) - value_curve = PSY.get_value_curve(cost_function) - cost_component = PSY.get_function_data(value_curve) - proportional_term = PSY.get_proportional_term(cost_component) - resolution = get_resolution(container) - dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR - for time_period in get_time_steps(container) - linear_cost = _add_proportional_term!( - container, - T(), - component, - proportional_term * multiplier * base_power * dt, - time_period, - ) - add_to_expression!( - container, - ProductionCostExpression, - linear_cost, - component, - time_period, - ) - end - return -end - -@doc raw""" -Adds to the cost function cost terms for sum of variables with common factor to be used for cost expression for optimization_container model. - -# Equation - -``` gen_cost = dt*sign*(sum(variable.^2)*cost_data[1] + sum(variable)*cost_data[2]) ``` - -# LaTeX - -`` cost = dt\times sign (sum_{i\in I} c_1 v_i^2 + sum_{i\in I} c_2 v_i ) `` - -for quadratic factor large enough. If the first term of the quadratic objective is 0.0, adds a -linear cost term `sum(variable)*cost_data[2]` - -# Arguments - -* container::OptimizationContainer : the optimization_container model built in PowerSimulations -* var_key::VariableKey: The variable name -* component_name::String: The component_name of the variable container -* cost_component::PSY.CostCurve{PSY.QuadraticCurve} : container for quadratic factors -""" -function _add_variable_cost_to_objective!( - container::OptimizationContainer, - ::T, - component::PSY.Component, - cost_function::PSY.CostCurve{PSY.QuadraticCurve}, - ::U, -) where {T <: VariableType, U <: AbstractDeviceFormulation} - multiplier = objective_function_multiplier(T(), U()) - base_power = get_base_power(container) - value_curve = PSY.get_value_curve(cost_function) - cost_component = PSY.get_function_data(value_curve) - quadratic_term = PSY.get_quadratic_term(cost_component) - proportional_term = PSY.get_proportional_term(cost_component) - constant_term = PSY.get_constant_term(cost_component) - (constant_term == 0) || - throw(ArgumentError("Not yet implemented for nonzero constant term")) - resolution = get_resolution(container) - dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR - for time_period in get_time_steps(container) - if quadratic_term >= eps() - cost_term = _add_quadratic_term!( - container, - T(), - component, - (quadratic_term, proportional_term), - base_power, - multiplier * dt, - time_period, - ) - else - cost_term = _add_proportional_term!( - container, - T(), - component, - proportional_term * multiplier * base_power * dt, - time_period, - ) - end - add_to_expression!( - container, - ProductionCostExpression, - cost_term, - component, - time_period, - ) - end - return -end - -""" -Creates piecewise linear cost function using a sum of variables and expression with sign and time step included. - -# Arguments - - - container::OptimizationContainer : the optimization_container model built in PowerSimulations - - var_key::VariableKey: The variable name - - component_name::String: The component_name of the variable container - - cost_function::PSY.CostCurve{PSY.PiecewisePointCurve}: container for piecewise linear cost -""" -function _add_variable_cost_to_objective!( - container::OptimizationContainer, - ::T, - component::PSY.Component, - cost_function::PSY.CostCurve{PSY.PiecewisePointCurve}, - ::U, -) where {T <: VariableType, U <: AbstractDeviceFormulation} - component_name = PSY.get_name(component) - @debug "PWL Variable Cost" _group = LOG_GROUP_COST_FUNCTIONS component_name - # If array is full of tuples with zeros return 0.0 - value_curve = PSY.get_value_curve(cost_function) - cost_component = PSY.get_function_data(value_curve) - if all(iszero.((point -> point.y).(PSY.get_points(cost_component)))) # TODO I think this should have been first. before? - @debug "All cost terms for component $(component_name) are 0.0" _group = - LOG_GROUP_COST_FUNCTIONS - return - end - pwl_cost_expressions = _add_pwl_term!(container, component, cost_component, T(), U()) - for t in get_time_steps(container) - add_to_expression!( - container, - ProductionCostExpression, - pwl_cost_expressions[t], - component, - t, - ) - add_to_objective_invariant_expression!(container, pwl_cost_expressions[t]) - end - return -end - -function _get_sos_value( - container::OptimizationContainer, - ::Type{V}, - component::T, -) where {T <: PSY.Component, V <: AbstractDeviceFormulation} - if has_container_key(container, OnStatusParameter, T) - sos_val = SOSStatusVariable.PARAMETER - else - sos_val = sos_status(component, V()) - end - return sos_val -end - -function _get_sos_value( - container::OptimizationContainer, - ::Type{V}, - component::T, -) where {T <: PSY.Component, V <: AbstractServiceFormulation} - return SOSStatusVariable.NO_VARIABLE -end - -function _add_pwl_term!( - container::OptimizationContainer, - component::T, - cost_data::AbstractVector{PSY.LinearFunctionData}, - ::U, - ::V, -) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} - multiplier = objective_function_multiplier(U(), V()) - resolution = get_resolution(container) - dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR - base_power = get_base_power(container) - # Re-scale breakpoints by Basepower - time_steps = get_time_steps(container) - cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) - for t in time_steps - proportional_value = - PSY.get_proportional_term(cost_data[t]) * multiplier * base_power * dt - cost_expressions[t] = - _add_proportional_term!(container, U(), component, proportional_value, t) - end - return cost_expressions -end - -""" -Add PWL cost terms for data coming from the MarketBidCost -""" -function _add_pwl_term!( - container::OptimizationContainer, - component::T, - cost_data::AbstractVector{PSY.PiecewiseStepData}, - ::U, - ::V, -) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} - multiplier = objective_function_multiplier(U(), V()) - resolution = get_resolution(container) - dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR - base_power = get_base_power(container) - # Re-scale breakpoints by Basepower - name = PSY.get_name(component) - time_steps = get_time_steps(container) - pwl_cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) - sos_val = _get_sos_value(container, V, component) - for t in time_steps - # Run checks in every time step because each time step has a PWL cost function - data = cost_data[t] - compact_status = validate_compact_pwl_data(component, data, base_power) - if !uses_compact_power(component, V()) && compact_status == COMPACT_PWL_STATUS.VALID - error( - "The data provided is not compatible with formulation $V. Use a formulation compatible with Compact Cost Functions", - ) - # data = _convert_to_full_variable_cost(data, component) - elseif uses_compact_power(component, V()) && - compact_status != COMPACT_PWL_STATUS.VALID - @warn( - "The cost data provided is not in compact form. Will attempt to convert. Errors may occur." - ) - data = _convert_to_compact_variable_cost(data) - else - @debug uses_compact_power(component, V()) compact_status name T V - end - cost_is_convex = PSY.is_convex(data) - break_points = PSY.get_x_coords(data) ./ base_power # TODO should this be get_x_lengths/get_breakpoint_upper_bounds? - _add_pwl_variables!(container, T, name, t, data) - _add_pwl_constraint!(container, component, U(), break_points, sos_val, t) - if !cost_is_convex - _add_pwl_sos_constraint!(container, component, U(), break_points, sos_val, t) - end - pwl_cost = - _get_pwl_cost_expression(container, component, t, data, multiplier * dt) - pwl_cost_expressions[t] = pwl_cost - end - return pwl_cost_expressions -end - -function _add_pwl_term!( - container::OptimizationContainer, - component::T, - cost_data::AbstractVector{PSY.PiecewiseStepData}, - ::U, - ::V, -) where {T <: PSY.Component, U <: VariableType, V <: AbstractServiceFormulation} - multiplier = objective_function_multiplier(U(), V()) - resolution = get_resolution(container) - dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR - base_power = get_base_power(container) - # Re-scale breakpoints by Basepower - name = PSY.get_name(component) - time_steps = get_time_steps(container) - pwl_cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) - sos_val = _get_sos_value(container, V, component) - for t in time_steps - data = cost_data[t] - break_points = PSY.get_x_coords(data) ./ base_power - _add_pwl_variables!(container, T, name, t, data) - _add_pwl_constraint!(container, component, U(), break_points, sos_val, t) - _add_pwl_sos_constraint!(container, component, U(), break_points, sos_val, t) - pwl_cost = _get_pwl_cost_expression(container, component, t, data, multiplier * dt) - pwl_cost_expressions[t] = pwl_cost - end - return pwl_cost_expressions -end - -""" -Add PWL cost terms for data coming from a constant PWL cost function -""" -function _add_pwl_term!( - container::OptimizationContainer, - component::T, - data::PSY.PiecewiseLinearData, - ::U, - ::V, -) where {T <: PSY.Component, U <: VariableType, V <: AbstractDeviceFormulation} - multiplier = objective_function_multiplier(U(), V()) - resolution = get_resolution(container) - dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR - base_power = get_base_power(container) - # Re-scale breakpoints by Basepower - name = PSY.get_name(component) - - compact_status = validate_compact_pwl_data(component, data, base_power) - if !uses_compact_power(component, V()) && compact_status == COMPACT_PWL_STATUS.VALID - error( - "The data provided is not compatible with formulation $V. Use a formulation compatible with Compact Cost Functions", - ) - # data = _convert_to_full_variable_cost(data, component) - elseif uses_compact_power(component, V()) && compact_status != COMPACT_PWL_STATUS.VALID - @warn( - "The cost data provided is not in compact form. Will attempt to convert. Errors may occur." - ) - data = _convert_to_compact_variable_cost(data) - else - @debug uses_compact_power(component, V()) compact_status name T V - end - - cost_is_convex = PSY.is_convex(data) - break_points = PSY.get_x_coords(data) ./ base_power - time_steps = get_time_steps(container) - pwl_cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) - sos_val = _get_sos_value(container, V, component) - for t in time_steps - _add_pwl_variables!(container, T, name, t, data) - _add_pwl_constraint!(container, component, U(), break_points, sos_val, t) - if !cost_is_convex - _add_pwl_sos_constraint!(container, component, U(), break_points, sos_val, t) - end - pwl_cost = _get_pwl_cost_expression(container, component, t, data, multiplier * dt) - pwl_cost_expressions[t] = pwl_cost - end - return pwl_cost_expressions -end - -function _add_pwl_term!( - container::OptimizationContainer, - component::T, - data::PSY.PiecewiseLinearData, - ::U, - ::V, -) where {T <: PSY.ThermalGen, U <: VariableType, V <: ThermalDispatchNoMin} - multiplier = objective_function_multiplier(U(), V()) - resolution = get_resolution(container) - dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR - component_name = PSY.get_name(component) - @debug "PWL cost function detected for device $(component_name) using $V" - base_power = get_base_power(container) - slopes = PSY.get_slopes(data) - if any(slopes .< 0) || !PSY.is_convex(data) - throw( - IS.InvalidValue( - "The PWL cost data provided for generator $(component_name) is not compatible with $U.", - ), - ) - end - - if validate_compact_pwl_data(component, data, base_power) == COMPACT_PWL_STATUS.VALID - error("The data provided is not compatible with formulation $V. \\ - Use a formulation compatible with Compact Cost Functions") - end - - if slopes[1] != 0.0 - @debug "PWL has no 0.0 intercept for generator $(component_name)" - # adds a first intercept a x = 0.0 and y below the intercept of the first tuple to make convex equivalent - intercept_point = (x = 0.0, y = first(data).y - COST_EPSILON) - data = PSY.PiecewiseLinearData(vcat(intercept_point, get_points(data))) - @assert PSY.is_convex(slopes) - end - - time_steps = get_time_steps(container) - pwl_cost_expressions = Vector{JuMP.AffExpr}(undef, time_steps[end]) - break_points = PSY.get_x_coords(data) ./ base_power - sos_val = _get_sos_value(container, V, component) - for t in time_steps - _add_pwl_variables!(container, T, component_name, t, data) - _add_pwl_constraint!(container, component, U(), break_points, sos_val, t) - pwl_cost = _get_pwl_cost_expression(container, component, t, data, multiplier * dt) - pwl_cost_expressions[t] = pwl_cost - end - return pwl_cost_expressions -end - -# This cases bounds the data by 1 - 0 -function _add_pwl_variables!( - container::OptimizationContainer, - ::Type{T}, - component_name::String, - time_period::Int, - cost_data::PSY.PiecewiseLinearData, -) where {T <: PSY.Component} - var_container = lazy_container_addition!(container, PieceWiseLinearCostVariable(), T) - # length(PiecewiseStepData) gets number of segments, here we want number of points - pwlvars = Array{JuMP.VariableRef}(undef, length(cost_data) + 1) - for i in 1:(length(cost_data) + 1) - pwlvars[i] = - var_container[(component_name, i, time_period)] = JuMP.@variable( - get_jump_model(container), - base_name = "PieceWiseLinearCostVariable_$(component_name)_{pwl_$(i), $time_period}", - lower_bound = 0.0, - upper_bound = 1.0 - ) - end - return pwlvars -end - -function _add_pwl_variables!( - container::OptimizationContainer, - ::Type{T}, - component_name::String, - time_period::Int, - cost_data::PSY.PiecewiseStepData, -) where {T <: PSY.Component} - var_container = lazy_container_addition!(container, PieceWiseLinearCostVariable(), T) - # length(PiecewiseStepData) gets number of segments, here we want number of points - pwlvars = Array{JuMP.VariableRef}(undef, length(cost_data) + 1) - for i in 1:(length(cost_data) + 1) - pwlvars[i] = - var_container[(component_name, i, time_period)] = JuMP.@variable( - get_jump_model(container), - base_name = "PieceWiseLinearCostVariable_$(component_name)_{pwl_$(i), $time_period}", - ) - end - return pwlvars -end - -function _add_pwl_constraint!( - container::OptimizationContainer, - component::T, - ::U, - break_points::Vector{Float64}, - sos_status::SOSStatusVariable, - period::Int, -) where {T <: PSY.Component, U <: VariableType} - variables = get_variable(container, U(), T) - const_container = lazy_container_addition!( - container, - PieceWiseLinearCostConstraint(), - T, - axes(variables)..., - ) - len_cost_data = length(break_points) - jump_model = get_jump_model(container) - pwl_vars = get_variable(container, PieceWiseLinearCostVariable(), T) - name = PSY.get_name(component) - const_container[name, period] = JuMP.@constraint( - jump_model, - variables[name, period] == - sum(pwl_vars[name, ix, period] * break_points[ix] for ix in 1:len_cost_data) - ) - - if sos_status == SOSStatusVariable.NO_VARIABLE - bin = 1.0 - @debug "Using Piecewise Linear cost function but no variable/parameter ref for ON status is passed. Default status will be set to online (1.0)" _group = - LOG_GROUP_COST_FUNCTIONS - - elseif sos_status == SOSStatusVariable.PARAMETER - param = get_default_on_parameter(component) - bin = get_parameter(container, param, T).parameter_array[name, period] - @debug "Using Piecewise Linear cost function with parameter OnStatusParameter, $T" _group = - LOG_GROUP_COST_FUNCTIONS - elseif sos_status == SOSStatusVariable.VARIABLE - var = get_default_on_variable(component) - bin = get_variable(container, var, T)[name, period] - @debug "Using Piecewise Linear cost function with variable OnVariable $T" _group = - LOG_GROUP_COST_FUNCTIONS - else - @assert false - end - - JuMP.@constraint( - jump_model, - sum(pwl_vars[name, i, period] for i in 1:len_cost_data) == bin - ) - return -end - -function _add_pwl_sos_constraint!( - container::OptimizationContainer, - component::T, - ::U, - break_points::Vector{Float64}, - sos_status::SOSStatusVariable, - period::Int, -) where {T <: PSY.Component, U <: VariableType} - name = PSY.get_name(component) - @warn( - "The cost function provided for $(name) is not compatible with a linear PWL cost function. - An SOS-2 formulation will be added to the model. This will result in additional binary variables." - ) - - jump_model = get_jump_model(container) - pwl_vars = get_variable(container, PieceWiseLinearCostVariable(), T) - bp_count = length(break_points) - pwl_vars_subset = [pwl_vars[name, i, period] for i in 1:bp_count] - JuMP.@constraint(jump_model, pwl_vars_subset in MOI.SOS2(collect(1:bp_count))) - return -end - -function _get_pwl_cost_expression( - container::OptimizationContainer, - component::T, - time_period::Int, - cost_data::PSY.PiecewiseLinearData, - multiplier::Float64, -) where {T <: PSY.Component} - name = PSY.get_name(component) - pwl_var_container = get_variable(container, PieceWiseLinearCostVariable(), T) - gen_cost = JuMP.AffExpr(0.0) - cost_data = PSY.get_y_coords(cost_data) - for i in 1:length(cost_data) - JuMP.add_to_expression!( - gen_cost, - cost_data[i] * multiplier * pwl_var_container[(name, i, time_period)], - ) - end - return gen_cost -end - -function _get_pwl_cost_expression( - container::OptimizationContainer, - component::T, - time_period::Int, - cost_data::PSY.PiecewiseStepData, - multiplier::Float64, -) where {T <: PSY.Component} - # TODO: This functions needs to be reimplemented for the new model. The code is repeated - # because the internals will be different - name = PSY.get_name(component) - pwl_var_container = get_variable(container, PieceWiseLinearCostVariable(), T) - gen_cost = JuMP.AffExpr(0.0) - cost_data = PSY.get_y_coords(cost_data) - for i in 1:length(cost_data) - JuMP.add_to_expression!( - gen_cost, - cost_data[i] * multiplier * pwl_var_container[(name, i, time_period)], - ) - end - return gen_cost -end - -# These conversions are not properly done for the new models -function _convert_to_compact_variable_cost( - var_cost::PSY.PiecewiseLinearData, - p_min::Float64, - no_load_cost::Float64, -) - points = PSY.get_points(var_cost) - new_points = [(pp - p_min, c - no_load_cost) for (pp, c) in points] - return PSY.PiecewiseLinearData(new_points) -end - -# These conversions are not properly done for the new models -function _convert_to_compact_variable_cost( - var_cost::PSY.PiecewiseStepData, - p_min::Float64, - no_load_cost::Float64, -) - x = PSY.get_x_coords(var_cost) - y = vcat(PSY.get_y_coords(var_cost), PSY.get_y_coords(var_cost)[end]) - points = [(x[i], y[i]) for i in length(x)] - new_points = [(x = pp - p_min, y = c - no_load_cost) for (pp, c) in points] - return PSY.PiecewiseLinearData(new_points) -end - -# TODO: This method needs to be corrected to account for actual StepData. The TestData is point wise -function _convert_to_compact_variable_cost(var_cost::PSY.PiecewiseStepData) - p_min, no_load_cost = (PSY.get_x_coords(var_cost)[1], PSY.get_y_coords(var_cost)[1]) - return _convert_to_compact_variable_cost(var_cost, p_min, no_load_cost) -end - -function _convert_to_compact_variable_cost(var_cost::PSY.PiecewiseLinearData) - p_min, no_load_cost = first(PSY.get_points(var_cost)) - return _convert_to_compact_variable_cost(var_cost, p_min, no_load_cost) -end - -function _add_proportional_term!( - container::OptimizationContainer, - ::T, - component::U, - linear_term::Float64, - time_period::Int, -) where {T <: VariableType, U <: PSY.Component} - component_name = PSY.get_name(component) - @debug "Linear Variable Cost" _group = LOG_GROUP_COST_FUNCTIONS component_name - variable = get_variable(container, T(), U)[component_name, time_period] - lin_cost = variable * linear_term - add_to_objective_invariant_expression!(container, lin_cost) - return lin_cost -end - -function _add_quadratic_term!( - container::OptimizationContainer, - ::T, - component::U, - q_terms::NTuple{2, Float64}, - var_multiplier::Float64, - expression_multiplier::Float64, - time_period::Int, -) where {T <: VariableType, U <: PSY.Component} - component_name = PSY.get_name(component) - @debug "$component_name Quadratic Variable Cost" _group = LOG_GROUP_COST_FUNCTIONS component_name - var = get_variable(container, T(), U)[component_name, time_period] - q_cost_ = (var * var_multiplier) .^ 2 * q_terms[1] + var * var_multiplier * q_terms[2] - q_cost = q_cost_ * expression_multiplier - add_to_objective_invariant_expression!(container, q_cost) - return q_cost -end - -function _add_quadratic_term!( - container::OptimizationContainer, - ::T, - component::U, - q_terms::NTuple{2, Float64}, - var_multiplier::Float64, - expression_multiplier::Float64, - time_period::Int, -) where {T <: PowerAboveMinimumVariable, U <: PSY.ThermalGen} - component_name = PSY.get_name(component) - p_min = PSY.get_active_power_limits(component).min - @debug "$component_name Quadratic Variable Cost" _group = LOG_GROUP_COST_FUNCTIONS component_name - var = get_variable(container, T(), U)[component_name, time_period] - q_cost_ = - (var * var_multiplier) .^ 2 * q_terms[1] + - var * var_multiplier * (q_terms[2] + 2 * q_terms[1] * p_min) - q_cost = q_cost_ * expression_multiplier - add_to_objective_invariant_expression!(container, q_cost) - return q_cost -end diff --git a/src/devices_models/devices/common/rateofchange_constraints.jl b/src/devices_models/devices/common/rateofchange_constraints.jl index 7f9ec2a5d5..b5558b56b1 100644 --- a/src/devices_models/devices/common/rateofchange_constraints.jl +++ b/src/devices_models/devices/common/rateofchange_constraints.jl @@ -89,14 +89,14 @@ function add_linear_ramp_constraints!( name ∉ set_name && continue ramp_limits = PSY.get_ramp_limits(get_component(ic)) ic_power = get_value(ic) - @error "add rate_of_change_constraint" name ic_power + @debug "add rate_of_change_constraint" name ic_power con_up[name, 1] = JuMP.@constraint( get_jump_model(container), expr_up[name, 1] - ic_power <= ramp_limits.up * minutes_per_period ) con_down[name, 1] = JuMP.@constraint( get_jump_model(container), - ic_power - expr_dn[name, 1] <= ramp_limits.down * minutes_per_period + ic_power - expr_dn[name, 1] >= -1 * ramp_limits.down * minutes_per_period ) for t in time_steps[2:end] con_up[name, t] = JuMP.@constraint( @@ -106,8 +106,8 @@ function add_linear_ramp_constraints!( ) con_down[name, t] = JuMP.@constraint( get_jump_model(container), - variable[name, t - 1] - expr_dn[name, t] <= - ramp_limits.down * minutes_per_period + variable[name, t - 1] - expr_dn[name, t] >= + -1 * ramp_limits.down * minutes_per_period ) end end diff --git a/src/devices_models/devices/thermal_generation.jl b/src/devices_models/devices/thermal_generation.jl index 8b8fa437b3..28b12488da 100644 --- a/src/devices_models/devices/thermal_generation.jl +++ b/src/devices_models/devices/thermal_generation.jl @@ -74,7 +74,9 @@ initial_condition_default(::InitialTimeDurationOff, d::PSY.ThermalGen, ::Abstrac initial_condition_variable(::InitialTimeDurationOff, d::PSY.ThermalGen, ::AbstractThermalFormulation) = OnVariable() ########################Objective Function################################################## -proportional_cost(cost::PSY.ThermalGenerationCost, S::OnVariable, T::PSY.ThermalGen, U::AbstractThermalFormulation) = no_load_cost(cost, S, T, U) +# TODO: Decide what is the cost for OnVariable, if fixed or constant term in variable +#proportional_cost(cost::PSY.ThermalGenerationCost, S::OnVariable, T::PSY.ThermalGen, U::AbstractThermalFormulation) = no_load_cost(cost, S, T, U) +proportional_cost(cost::PSY.ThermalGenerationCost, S::OnVariable, T::PSY.ThermalGen, U::AbstractThermalFormulation) = PSY.get_fixed(cost) proportional_cost(cost::PSY.MarketBidCost, ::OnVariable, ::PSY.ThermalGen, ::AbstractThermalFormulation) = PSY.get_no_load_cost(cost) has_multistart_variables(::PSY.ThermalGen, ::AbstractThermalFormulation)=false @@ -104,32 +106,44 @@ uses_compact_power(::PSY.ThermalGen, ::ThermalCompactDispatch)=true variable_cost(cost::PSY.OperationalCost, ::ActivePowerVariable, ::PSY.ThermalGen, ::AbstractThermalFormulation)=PSY.get_variable(cost) variable_cost(cost::PSY.OperationalCost, ::PowerAboveMinimumVariable, ::PSY.ThermalGen, ::AbstractThermalFormulation)=PSY.get_variable(cost) -function no_load_cost(cost::PSY.ThermalGenerationCost, S::OnVariable, T::PSY.ThermalGen, U::AbstractThermalFormulation) - return no_load_cost(PSY.get_variable(cost), S, T, U) + PSY.get_fixed(cost) +""" +Theoretical Cost at power output zero. Mathematically is the intercept with the y-axis +""" +function no_load_cost(cost::PSY.ThermalGenerationCost, S::OnVariable, d::PSY.ThermalGen, U::AbstractThermalFormulation) + return _no_load_cost(PSY.get_variable(cost), d) end -# TODO given the old implementations, these functions seem to get the cost at *minimum* load, not *zero* load. Is that correct? -function no_load_cost(cost_function::PSY.CostCurve{PSY.PiecewisePointCurve}, ::OnVariable, ::PSY.ThermalGen, ::AbstractThermalFormulation) +function _no_load_cost(cost_function::PSY.CostCurve{PSY.PiecewisePointCurve}, d::PSY.ThermalGen) value_curve = PSY.get_value_curve(cost_function) cost = PSY.get_function_data(value_curve) return last(first(PSY.get_points(cost))) end -function no_load_cost(cost_function::PSY.CostCurve{PSY.LinearCurve}, ::OnVariable, d::PSY.ThermalGen, ::AbstractThermalFormulation) +function _no_load_cost(cost_function::Union{PSY.CostCurve{PSY.LinearCurve}, PSY.CostCurve{PSY.QuadraticCurve}}, d::PSY.ThermalGen) value_curve = PSY.get_value_curve(cost_function) - cost = PSY.get_function_data(value_curve) - return PSY.get_proportional_term(cost) * PSY.get_active_power_limits(d).min * PSY.get_system_base_power(d) + cost_component = PSY.get_function_data(value_curve) + # Always in \$/h + constant_term = PSY.get_constant_term(cost_component) + return constant_term +end + +function _no_load_cost(cost_function::PSY.FuelCurve{PSY.PiecewisePointCurve}, d::PSY.ThermalGen) + # value_curve = PSY.get_value_curve(cost_function) + # cost = PSY.get_function_data(value_curve) + return 0.0 end -function no_load_cost(cost_function::PSY.CostCurve{PSY.QuadraticCurve}, ::OnVariable, d::PSY.ThermalGen, ::AbstractThermalFormulation) - min_power = PSY.get_active_power_limits(d).min +function _no_load_cost(cost_function::Union{PSY.FuelCurve{PSY.LinearCurve}, PSY.FuelCurve{PSY.QuadraticCurve}}, d::PSY.ThermalGen) value_curve = PSY.get_value_curve(cost_function) - cost = PSY.get_function_data(value_curve) - evaluated = LinearAlgebra.dot( - [PSY.get_quadratic_term(cost), PSY.get_proportional_term(cost), PSY.get_constant_term(cost)], - [min_power^2, min_power, 1] - ) - return evaluated * PSY.get_system_base_power(d) + cost_component = PSY.get_function_data(value_curve) + # In Unit/h (unit typically in ) + constant_term = PSY.get_constant_term(cost_component) + fuel_cost = PSY.get_fuel_cost(cost_function) + if typeof(fuel_cost) <: Float64 + return constant_term * fuel_cost + else + error("Time series not implemented yet") + end end #! format: on diff --git a/src/feedforward/feedforward_constraints.jl b/src/feedforward/feedforward_constraints.jl index 52ac9c7d92..8c9ba13366 100644 --- a/src/feedforward/feedforward_constraints.jl +++ b/src/feedforward/feedforward_constraints.jl @@ -83,7 +83,7 @@ function _add_sc_feedforward_constraints!( devices::IS.FlattenIteratorWrapper{V}, model::DeviceModel{V, W}, ) where { - T <: FeedforwardSemiContinousConstraint, + T <: FeedforwardSemiContinuousConstraint, P <: OnStatusParameter, U <: Union{ActivePowerVariable, PowerAboveMinimumVariable}, V <: PSY.Component, @@ -128,7 +128,7 @@ function _add_sc_feedforward_constraints!( devices::IS.FlattenIteratorWrapper{V}, model::DeviceModel{V, W}, ) where { - T <: FeedforwardSemiContinousConstraint, + T <: FeedforwardSemiContinuousConstraint, P <: ParameterType, U <: VariableType, V <: PSY.Component, @@ -230,7 +230,7 @@ function add_feedforward_constraints!( end _add_sc_feedforward_constraints!( container, - FeedforwardSemiContinousConstraint, + FeedforwardSemiContinuousConstraint, parameter_type(), var, devices, diff --git a/src/feedforward/feedforwards.jl b/src/feedforward/feedforwards.jl index f0f04e1b54..a69b00a2ef 100644 --- a/src/feedforward/feedforwards.jl +++ b/src/feedforward/feedforwards.jl @@ -50,7 +50,22 @@ function get_feedforward_meta(ff::AbstractAffectFeedforward) end """ -Adds an upper bound constraint to a variable. + UpperBoundFeedforward( + component_type::Type{<:PSY.Component}, + source::Type{T}, + affected_values::Vector{DataType}, + add_slacks::Bool = false, + meta = CONTAINER_KEY_EMPTY_META + ) where {T} + +Constructs a parameterized upper bound constraint to implement feedforward from other models. + +# Arguments: +* `component_type::Type{<:PSY.Component}` : Specify the type of component on which the Feedforward will be applied +* `source::Type{T}` : Specify the VariableType, ParameterType or AuxVariableType as the source of values for the Feedforward +* `affected_values::Vector{DataType}` : Specify the variable on which the upper bound will be applied using the source values +* `add_slacks::Bool = false` : Add slacks variables to relax the upper bound constraint. + """ struct UpperBoundFeedforward <: AbstractAffectFeedforward optimization_container_key::OptimizationContainerKey @@ -87,7 +102,22 @@ get_optimization_container_key(ff::UpperBoundFeedforward) = ff.optimization_cont get_slacks(ff::UpperBoundFeedforward) = ff.add_slacks """ -Adds a lower bound constraint to a variable. + LowerBoundFeedforward( + component_type::Type{<:PSY.Component}, + source::Type{T}, + affected_values::Vector{DataType}, + add_slacks::Bool = false, + meta = CONTAINER_KEY_EMPTY_META + ) where {T} + +Constructs a parameterized lower bound constraint to implement feedforward from other models. + +# Arguments: +* `component_type::Type{<:PSY.Component}` : Specify the type of component on which the Feedforward will be applied +* `source::Type{T}` : Specify the VariableType, ParameterType or AuxVariableType as the source of values for the Feedforward +* `affected_values::Vector{DataType}` : Specify the variable on which the lower bound will be applied using the source values +* `add_slacks::Bool = false` : Add slacks variables to relax the lower bound constraint. + """ struct LowerBoundFeedforward <: AbstractAffectFeedforward optimization_container_key::OptimizationContainerKey @@ -149,7 +179,21 @@ function attach_feedforward!( end """ -Adds a constraint to make the bounds of a variable 0.0. Effectively allows to "turn off" a value. + SemiContinuousFeedforward( + component_type::Type{<:PSY.Component}, + source::Type{T}, + affected_values::Vector{DataType}, + meta = CONTAINER_KEY_EMPTY_META + ) where {T} + +It allows to enable/disable bounds to 0.0 for a specified variable. Commonly used to limit the +`ActivePowerVariable` in an Economic Dispatch problem by the commitment decision taken in +an another problem (typically a Unit Commitment problem). + +# Arguments: +* `component_type::Type{<:PSY.Component}` : Specify the type of component on which the Feedforward will be applied +* `source::Type{T}` : Specify the VariableType, ParameterType or AuxVariableType as the source of values for the Feedforward +* `affected_values::Vector{DataType}` : Specify the variable on which the semicontinuous limit will be applied using the source values """ struct SemiContinuousFeedforward <: AbstractAffectFeedforward optimization_container_key::OptimizationContainerKey @@ -203,8 +247,20 @@ function has_semicontinuous_feedforward( end """ -Fixes a Variable or Parameter Value in the model. Is the only Feed Forward that can be used + FixValueFeedforward( + component_type::Type{<:PSY.Component}, + source::Type{T}, + affected_values::Vector{DataType}, + meta = CONTAINER_KEY_EMPTY_META + ) where {T} + +Fixes a Variable or Parameter Value in the model from another problem. Is the only Feed Forward that can be used with a Parameter or a Variable as the affected value. + +# Arguments: +* `component_type::Type{<:PSY.Component}` : Specify the type of component on which the Feedforward will be applied +* `source::Type{T}` : Specify the VariableType, ParameterType or AuxVariableType as the source of values for the Feedforward +* `affected_values::Vector{DataType}` : Specify the variable on which the fix value will be applied using the source values """ struct FixValueFeedforward <: AbstractAffectFeedforward optimization_container_key::OptimizationContainerKey diff --git a/src/operation/decision_model.jl b/src/operation/decision_model.jl index c0959cbb11..091e9833c5 100644 --- a/src/operation/decision_model.jl +++ b/src/operation/decision_model.jl @@ -278,7 +278,7 @@ end function validate_time_series(model::DecisionModel{<:DefaultDecisionProblem}) sys = get_system(model) settings = get_settings(model) - available_resolutions = PSY.list_time_series_resolutions(sys) + available_resolutions = PSY.get_time_series_resolutions(sys) if get_resolution(settings) == UNSET_RESOLUTION && length(available_resolutions) != 1 throw( @@ -299,9 +299,7 @@ function validate_time_series(model::DecisionModel{<:DefaultDecisionProblem}) end if get_horizon(settings) == UNSET_HORIZON - # TODO: forecast horizon needs to return a TimePeriod value - resolution = get_resolution(settings) - set_horizon!(settings, PSY.get_forecast_horizon(sys) * resolution) + set_horizon!(settings, PSY.get_forecast_horizon(sys)) end counts = PSY.get_time_series_counts(sys) diff --git a/src/operation/emulation_model.jl b/src/operation/emulation_model.jl index 6b5d3fb16e..22d7d82258 100644 --- a/src/operation/emulation_model.jl +++ b/src/operation/emulation_model.jl @@ -279,7 +279,7 @@ function validate_time_series(model::EmulationModel{<:DefaultEmulationProblem}) end settings = get_settings(model) - available_resolutions = PSY.list_time_series_resolutions(sys) + available_resolutions = PSY.get_time_series_resolutions(sys) if get_resolution(settings) == UNSET_RESOLUTION && length(available_resolutions) != 1 throw( diff --git a/src/parameters/update_parameters.jl b/src/parameters/update_parameters.jl index dcb9f2a73c..be4ca5cdb0 100644 --- a/src/parameters/update_parameters.jl +++ b/src/parameters/update_parameters.jl @@ -531,7 +531,7 @@ function _update_pwl_cost_expression( ) where {T <: PSY.Component} pwl_var_container = get_variable(container, PieceWiseLinearCostVariable(), T) resolution = get_resolution(container) - dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR gen_cost = JuMP.AffExpr(0.0) slopes = PSY.get_slopes(cost_data) upb = get_breakpoint_upper_bounds(cost_data) @@ -553,7 +553,7 @@ function update_variable_cost!( time_period::Int, ) where {T <: PSY.Component} resolution = get_resolution(container) - dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR + dt = Dates.value(resolution) / MILLISECONDS_IN_HOUR base_power = get_base_power(container) component_name = PSY.get_name(component) cost_data = parameter_array[component_name, time_period] # TODO is this a new-style cost? diff --git a/src/simulation/simulation.jl b/src/simulation/simulation.jl index f668e47be0..cfc540d472 100644 --- a/src/simulation/simulation.jl +++ b/src/simulation/simulation.jl @@ -1040,7 +1040,7 @@ Solves the simulation model for sequential Simulations. - `sim::Simulation=sim`: simulation object created by Simulation() The optional keyword argument `exports` controls exporting of results to CSV files as -the simulation runs. Refer to [`export_results`](@ref) for a description of this argument. +the simulation runs. # Example ```julia diff --git a/src/utils/powersystems_utils.jl b/src/utils/powersystems_utils.jl index b6ff621ff1..7cf36f874c 100644 --- a/src/utils/powersystems_utils.jl +++ b/src/utils/powersystems_utils.jl @@ -126,6 +126,255 @@ function check_hvdc_line_limits_unidirectional(d::PSY.TwoTerminalHVDCLine) return end +################################################## +########### Cost Function Utilities ############## +################################################## + +""" +Obtain proportional (marginal or slope) cost data in system base per unit +depending on the specified power units +""" +function get_proportional_cost_per_system_unit( + cost_term::Float64, + unit_system::PSY.UnitSystem, + system_base_power::Float64, + device_base_power::Float64, +) + return _get_proportional_cost_per_system_unit( + cost_term, + Val{unit_system}(), + system_base_power, + device_base_power, + ) +end + +function _get_proportional_cost_per_system_unit( + cost_term::Float64, + ::Val{PSY.UnitSystem.SYSTEM_BASE}, + system_base_power::Float64, + device_base_power::Float64, +) + return cost_term +end + +function _get_proportional_cost_per_system_unit( + cost_term::Float64, + ::Val{PSY.UnitSystem.DEVICE_BASE}, + system_base_power::Float64, + device_base_power::Float64, +) + return cost_term * (system_base_power / device_base_power) +end + +function _get_proportional_cost_per_system_unit( + cost_term::Float64, + ::Val{PSY.UnitSystem.NATURAL_UNITS}, + system_base_power::Float64, + device_base_power::Float64, +) + return cost_term * system_base_power +end + +""" +Obtain quadratic cost data in system base per unit +depending on the specified power units +""" +function get_quadratic_cost_per_system_unit( + cost_term::Float64, + unit_system::PSY.UnitSystem, + system_base_power::Float64, + device_base_power::Float64, +) + return _get_quadratic_cost_per_system_unit( + cost_term, + Val{unit_system}(), + system_base_power, + device_base_power, + ) +end + +function _get_quadratic_cost_per_system_unit( + cost_term::Float64, + ::Val{PSY.UnitSystem.SYSTEM_BASE}, # SystemBase Unit + system_base_power::Float64, + device_base_power::Float64, +) + return cost_term +end + +function _get_quadratic_cost_per_system_unit( + cost_term::Float64, + ::Val{PSY.UnitSystem.DEVICE_BASE}, # DeviceBase Unit + system_base_power::Float64, + device_base_power::Float64, +) + return cost_term * (system_base_power / device_base_power)^2 +end + +function _get_quadratic_cost_per_system_unit( + cost_term::Float64, + ::Val{PSY.UnitSystem.NATURAL_UNITS}, # Natural Units + system_base_power::Float64, + device_base_power::Float64, +) + return cost_term * system_base_power^2 +end + +""" +Obtain the normalized PiecewiseLinear cost data in system base per unit +depending on the specified power units. + +Note that the costs (y-axis) are always in \$/h so +they do not require transformation +""" +function get_piecewise_pointcurve_per_system_unit( + cost_component::PSY.PiecewiseLinearData, + unit_system::PSY.UnitSystem, + system_base_power::Float64, + device_base_power::Float64, +) + return _get_piecewise_pointcurve_per_system_unit( + cost_component, + Val{unit_system}(), + system_base_power, + device_base_power, + ) +end + +function _get_piecewise_pointcurve_per_system_unit( + cost_component::PSY.PiecewiseLinearData, + ::Val{PSY.UnitSystem.SYSTEM_BASE}, + system_base_power::Float64, + device_base_power::Float64, +) + return cost_component +end + +function _get_piecewise_pointcurve_per_system_unit( + cost_component::PSY.PiecewiseLinearData, + ::Val{PSY.UnitSystem.DEVICE_BASE}, + system_base_power::Float64, + device_base_power::Float64, +) + points = cost_component.points + points_normalized = Vector{NamedTuple{(:x, :y)}}(undef, length(points)) + for (ix, point) in enumerate(points) + points_normalized[ix] = + (x = point.x * (device_base_power / system_base_power), y = point.y) + end + return PSY.PiecewiseLinearData(points_normalized) +end + +function _get_piecewise_pointcurve_per_system_unit( + cost_component::PSY.PiecewiseLinearData, + ::Val{PSY.UnitSystem.NATURAL_UNITS}, + system_base_power::Float64, + device_base_power::Float64, +) + points = cost_component.points + points_normalized = Vector{NamedTuple{(:x, :y)}}(undef, length(points)) + for (ix, point) in enumerate(points) + points_normalized[ix] = (x = point.x / system_base_power, y = point.y) + end + return PSY.PiecewiseLinearData(points_normalized) +end + +""" +Obtain the normalized PiecewiseStep cost data in system base per unit +depending on the specified power units. + +Note that the costs (y-axis) are in \$/MWh, \$/(sys pu h) or \$/(device pu h), +so they also require transformation. +""" +function get_piecewise_incrementalcurve_per_system_unit( + cost_component::PSY.PiecewiseStepData, + unit_system::PSY.UnitSystem, + system_base_power::Float64, + device_base_power::Float64, +) + return _get_piecewise_incrementalcurve_per_system_unit( + cost_component, + Val{unit_system}(), + system_base_power, + device_base_power, + ) +end + +function _get_piecewise_incrementalcurve_per_system_unit( + cost_component::PSY.PiecewiseStepData, + ::Val{PSY.UnitSystem.SYSTEM_BASE}, + system_base_power::Float64, + device_base_power::Float64, +) + return cost_component +end + +function _get_piecewise_incrementalcurve_per_system_unit( + cost_component::PSY.PiecewiseStepData, + ::Val{PSY.UnitSystem.DEVICE_BASE}, + system_base_power::Float64, + device_base_power::Float64, +) + x_coords = PSY.get_x_coords(cost_component) + y_coords = PSY.get_y_coords(cost_component) + ratio = device_base_power / system_base_power + x_coords_normalized = x_coords .* ratio + y_coords_normalized = y_coords ./ ratio + return PSY.PiecewiseStepData(x_coords_normalized, y_coords_normalized) +end + +function _get_piecewise_incrementalcurve_per_system_unit( + cost_component::PSY.PiecewiseStepData, + ::Val{PSY.UnitSystem.NATURAL_UNITS}, + system_base_power::Float64, + device_base_power::Float64, +) + x_coords = PSY.get_x_coords(cost_component) + y_coords = PSY.get_y_coords(cost_component) + x_coords_normalized = x_coords ./ system_base_power + y_coords_normalized = y_coords .* system_base_power + return PSY.PiecewiseStepData(x_coords_normalized, y_coords_normalized) +end + +################################################## +############### Auxiliary Methods ################ +################################################## + +# These conversions are not properly done for the new models +function convert_to_compact_variable_cost( + var_cost::PSY.PiecewiseLinearData, + p_min::Float64, + no_load_cost::Float64, +) + points = PSY.get_points(var_cost) + new_points = [(pp - p_min, c - no_load_cost) for (pp, c) in points] + return PSY.PiecewiseLinearData(new_points) +end + +# These conversions are not properly done for the new models +function convert_to_compact_variable_cost( + var_cost::PSY.PiecewiseStepData, + p_min::Float64, + no_load_cost::Float64, +) + x = PSY.get_x_coords(var_cost) + y = vcat(PSY.get_y_coords(var_cost), PSY.get_y_coords(var_cost)[end]) + points = [(x[i], y[i]) for i in length(x)] + new_points = [(x = pp - p_min, y = c - no_load_cost) for (pp, c) in points] + return PSY.PiecewiseLinearData(new_points) +end + +# TODO: This method needs to be corrected to account for actual StepData. The TestData is point wise +function convert_to_compact_variable_cost(var_cost::PSY.PiecewiseStepData) + p_min, no_load_cost = (PSY.get_x_coords(var_cost)[1], PSY.get_y_coords(var_cost)[1]) + return convert_to_compact_variable_cost(var_cost, p_min, no_load_cost) +end + +function convert_to_compact_variable_cost(var_cost::PSY.PiecewiseLinearData) + p_min, no_load_cost = first(PSY.get_points(var_cost)) + return convert_to_compact_variable_cost(var_cost, p_min, no_load_cost) +end + function _validate_compact_pwl_data( min::Float64, max::Float64, diff --git a/test/test_device_branch_constructors.jl b/test/test_device_branch_constructors.jl index fa4ed35940..8c37601938 100644 --- a/test/test_device_branch_constructors.jl +++ b/test/test_device_branch_constructors.jl @@ -243,7 +243,7 @@ end NetworkModel(PTDFPowerModel), ) - set_device_model!(template_uc, ThermalStandard, ThermalCompactUnitCommitment) + set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment) set_device_model!(template_uc, RenewableDispatch, FixedOutput) set_device_model!(template_uc, PowerLoad, StaticPowerLoad) set_device_model!(template_uc, DeviceModel(Line, StaticBranch)) diff --git a/test/test_device_hvdc.jl b/test/test_device_hvdc.jl index 8e0c6d45c0..2829182c3a 100644 --- a/test/test_device_hvdc.jl +++ b/test/test_device_hvdc.jl @@ -7,7 +7,7 @@ #duals=[CopperPlateBalanceConstraint], )) - set_device_model!(template_uc, ThermalStandard, ThermalCompactUnitCommitment) + set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment) set_device_model!(template_uc, RenewableDispatch, RenewableFullDispatch) set_device_model!(template_uc, PowerLoad, StaticPowerLoad) set_device_model!(template_uc, DeviceModel(Line, StaticBranch)) @@ -25,7 +25,7 @@ #duals=[CopperPlateBalanceConstraint], )) - set_device_model!(template_uc, ThermalStandard, ThermalCompactUnitCommitment) + set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment) set_device_model!(template_uc, RenewableDispatch, RenewableFullDispatch) set_device_model!(template_uc, PowerLoad, StaticPowerLoad) set_device_model!(template_uc, DeviceModel(Line, StaticBranch)) diff --git a/test/test_device_thermal_generation_constructors.jl b/test/test_device_thermal_generation_constructors.jl index 6cd7a00cf5..52b49fb744 100644 --- a/test/test_device_thermal_generation_constructors.jl +++ b/test/test_device_thermal_generation_constructors.jl @@ -1,4 +1,101 @@ test_path = mktempdir() + +@testset "Test Thermal Generation Cost Functions " begin + test_cases = [ + ("linear_cost_test", 4664.88, ThermalBasicUnitCommitment), + ("linear_fuel_test", 4664.88, ThermalBasicUnitCommitment), + ("quadratic_cost_test", 3301.81, ThermalDispatchNoMin), + ("quadratic_fuel_test", 3331.12, ThermalDispatchNoMin), + ("pwl_io_cost_test", 3421.64, ThermalBasicUnitCommitment), + ("pwl_io_fuel_test", 3421.64, ThermalBasicUnitCommitment), + ("pwl_incremental_cost_test", 3424.43, ThermalBasicUnitCommitment), + ("pwl_incremental_fuel_test", 3424.43, ThermalBasicUnitCommitment), + ("non_convex_io_pwl_cost_test", 3047.14, ThermalBasicUnitCommitment), + ("fixed_market_bid_cost", 3047.14, ThermalBasicUnitCommitment) + ] + for (i, cost_reference, thermal_formulation) in test_cases + @testset "$i" begin + sys = build_system(PSITestSystems, "c_$(i)") + template = ProblemTemplate(NetworkModel(CopperPlatePowerModel)) + set_device_model!(template, ThermalStandard, thermal_formulation) + set_device_model!(template, PowerLoad, StaticPowerLoad) + model = DecisionModel( + template, + sys; + name = "UC_$(i)", + optimizer = HiGHS_optimizer, + system_to_file = false, + optimizer_solve_log_print = true, + ) + @test build!(model; output_dir = test_path) == PSI.ModelBuildStatus.BUILT + @test solve!(model) == PSI.RunStatus.SUCCESSFULLY_FINALIZED + results = OptimizationProblemResults(model) + expr = read_expression(results, "ProductionCostExpression__ThermalStandard") + var_unit_cost = sum(expr[!, "Test Unit"]) + @test isapprox(var_unit_cost, cost_reference; atol = 1) + @test expr[!, "Test Unit"][end] == 0.0 + end + end +end + +@testset "Test Thermal Generation Cost Functions Fuel Cost time series" begin + test_cases = [ + "linear_fuel_test_ts", + "quadratic_fuel_test_ts", + "pwl_io_fuel_test_ts", + "pwl_incremental_fuel_test_ts", + "market_bid_cost" + ] + for i in test_cases + @testset "$i" begin + sys = build_system(PSITestSystems, "c_$(i)") + template = ProblemTemplate(NetworkModel(CopperPlatePowerModel)) + set_device_model!(template, ThermalStandard, ThermalBasicUnitCommitment) + #= + model = DecisionModel( + template, + sys; + name = "UC_$(i)", + optimizer = HiGHS_optimizer, + system_to_file = false, + ) + @test build!(model; output_dir = test_path) == PSI.ModelBuildStatus.BUILT + @test solve!(model) == PSI.RunStatus.SUCCESSFULLY_FINALIZED + =# + end + end +end + +@testset "Test Thermal Generation MarketBidCost models" begin + test_cases = [ + ("fixed_market_bid_cost", 20532.76), + #"market_bid_cost", + ] + for (i, cost_reference) in test_cases + @testset "$i" begin + sys = build_system(PSITestSystems, "c_$(i)") + template = ProblemTemplate(NetworkModel(CopperPlatePowerModel)) + set_device_model!(template, ThermalStandard, ThermalBasicUnitCommitment) + set_device_model!(template, PowerLoad, StaticPowerLoad) + model = DecisionModel( + template, + sys; + name = "UC_$(i)", + optimizer = HiGHS_optimizer, + system_to_file = false, + optimizer_solve_log_print = true, + ) + @test build!(model; output_dir = test_path) == PSI.ModelBuildStatus.BUILT + @test solve!(model) == PSI.RunStatus.SUCCESSFULLY_FINALIZED + results = OptimizationProblemResults(model) + expr = read_expression(results, "ProductionCostExpression__ThermalStandard") + var_unit_cost = sum(expr[!, "Test Unit1"]) + @test isapprox(var_unit_cost, cost_reference; atol = 1) + @test expr[!, "Test Unit1"][end] == 0.0 + end + end +end + ################################### Unit Commitment tests ################################## @testset "Thermal UC With DC - PF" begin bin_variable_keys = [ @@ -606,38 +703,7 @@ end build!(UC; output_dir = mktempdir(; cleanup = true)) @test build!(UC; output_dir = mktempdir(; cleanup = true)) == PSI.ModelBuildStatus.BUILT moi_tests(UC, 56, 0, 56, 14, 21, true) - psi_checksolve_test(UC, [MOI.OPTIMAL], 13143.5) -end - -## PWL linear Cost implementation test -@testset "Solving UC with CopperPlate testing Convex PWL" begin - template = get_thermal_standard_uc_template() - UC = DecisionModel( - UnitCommitmentProblem, - template, - PSB.build_system(PSITestSystems, "c_linear_pwl_test"); - optimizer = HiGHS_optimizer, - initialize_model = false, - ) - @test build!(UC; output_dir = mktempdir(; cleanup = true)) == PSI.ModelBuildStatus.BUILT - moi_tests(UC, 32, 0, 8, 4, 14, true) - psi_checksolve_test(UC, [MOI.OPTIMAL], 13046.32, 0.01) -end - -@testset "Solving UC with CopperPlate testing PWL-SOS2 implementation" begin - template = get_thermal_standard_uc_template() - UC = DecisionModel( - UnitCommitmentProblem, - template, - PSB.build_system(PSITestSystems, "c_sos_pwl_test"); - optimizer = cbc_optimizer, - initialize_model = false, - ) - @test build!(UC; output_dir = mktempdir(; cleanup = true)) == PSI.ModelBuildStatus.BUILT - moi_tests(UC, 32, 0, 8, 4, 14, true) - # Cbc can have reliability issues with SoS. The objective function target in the this - # test was calculated with CPLEX do not change if Cbc gets a bad result - psi_checksolve_test(UC, [MOI.OPTIMAL], 13746.13, 10.0) + psi_checksolve_test(UC, [MOI.OPTIMAL], 8223.50) end #= Test disabled due to inconsistency between the models and the data @@ -815,7 +881,7 @@ end sys_5 = build_system(PSITestSystems, "c_sys5_uc") template_uc = ProblemTemplate(NetworkModel(PTDFPowerModel; PTDF_matrix = PTDF(sys_5))) - set_device_model!(template_uc, ThermalStandard, ThermalCompactUnitCommitment) + set_device_model!(template_uc, ThermalStandard, ThermalBasicUnitCommitment) set_device_model!(template_uc, RenewableDispatch, FixedOutput) set_device_model!(template_uc, PowerLoad, StaticPowerLoad) set_device_model!(template_uc, DeviceModel(Line, StaticBranch)) @@ -837,32 +903,3 @@ end on_sundance = on[!, "Sundance"] @test all(isapprox.(on_sundance, 1.0)) end - -#= -# NOTE not a comprehensive test, should expand as part of the cost refactor -@testset "Test no_load_cost" begin - sys = build_system(PSITestSystems, "c_sys5_uc") - comp = get_component(ThermalStandard, sys, "Sundance") - sys_base_power = get_base_power(sys) - set_base_power!(comp, 123.4) - min_limit = PSY.get_active_power_limits(comp).min - @test isapprox( - PSI.no_load_cost( - PSY.LinearFunctionData(5.0), - OnVariable(), - comp, - ThermalBasicUnitCommitment(), - ), - 5.0 * min_limit * sys_base_power, - ) - @test isapprox( - PSI.no_load_cost( - QuadraticFunctionData(3.0, 5.0, 0.0), - OnVariable(), - comp, - ThermalBasicUnitCommitment(), - ), - (3.0 * min_limit^2 + 5.0 * min_limit) * sys_base_power, - ) -end -=# diff --git a/test/test_simulation_build.jl b/test/test_simulation_build.jl index 2cc10818fa..0b3f2e4d7c 100644 --- a/test/test_simulation_build.jl +++ b/test/test_simulation_build.jl @@ -185,28 +185,28 @@ end ac_power_model = PSI.get_simulation_model(PSI.get_models(sim), :ED) c = PSI.get_constraint( PSI.get_optimization_container(ac_power_model), - FeedforwardSemiContinousConstraint(), + FeedforwardSemiContinuousConstraint(), ThermalStandard, "ActivePowerVariable_ub", ) @test !isempty(c) c = PSI.get_constraint( PSI.get_optimization_container(ac_power_model), - FeedforwardSemiContinousConstraint(), + FeedforwardSemiContinuousConstraint(), ThermalStandard, "ActivePowerVariable_lb", ) @test !isempty(c) c = PSI.get_constraint( PSI.get_optimization_container(ac_power_model), - FeedforwardSemiContinousConstraint(), + FeedforwardSemiContinuousConstraint(), ThermalStandard, "ReactivePowerVariable_ub", ) @test !isempty(c) c = PSI.get_constraint( PSI.get_optimization_container(ac_power_model), - FeedforwardSemiContinousConstraint(), + FeedforwardSemiContinuousConstraint(), ThermalStandard, "ReactivePowerVariable_lb", ) diff --git a/test/test_simulation_execute.jl b/test/test_simulation_execute.jl index ef5eb645f6..6bc4581769 100644 --- a/test/test_simulation_execute.jl +++ b/test/test_simulation_execute.jl @@ -59,9 +59,7 @@ function test_2_stage_decision_models_with_feedforwards(in_memory) template_uc, c_sys5_hy_uc; name = "UC", - optimizer = HiGHS_optimizer, - - ), + optimizer = HiGHS_optimizer), DecisionModel( template_ed, c_sys5_hy_ed; diff --git a/test/test_utils/mock_operation_models.jl b/test/test_utils/mock_operation_models.jl index ed50617e3f..a7c53f10ea 100644 --- a/test/test_utils/mock_operation_models.jl +++ b/test/test_utils/mock_operation_models.jl @@ -12,7 +12,7 @@ function PSI.DecisionModel( kwargs..., ) where {T <: PM.AbstractPowerModel} settings = PSI.Settings(sys; kwargs...) - available_resolutions = PSY.list_time_series_resolutions(sys) + available_resolutions = PSY.get_time_series_resolutions(sys) if length(available_resolutions) == 1 PSI.set_resolution!(settings, first(available_resolutions)) else