diff --git a/NEWS.md b/NEWS.md index 6a5f875..ecc4ab4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,14 @@ # Release notes +## Unversioned + +### Breaking changes + +#### Rework of `PeriodDemandSink` + +* Rewrote `PeriodDemandSink` with `PeriodPartition` (introduced in `TimeStruct` 0.9.12) to increase flexibility of node with respect to the time structure. +* Rewriting changed input arguments as well as behavior of the node. + ## Version 0.3.0 (2026-04-16) ### Breaking changes @@ -9,7 +18,7 @@ * Early retirement changes the model behavior. * Model worked without adjustments except for compatibility updates. * Breaking change still included to maintain the possibility to do bug fixes in version 0.1.x for existing models with `EnergyModelsBase` v0.9.x. - * UpdMoved from `Data` to `ExtensionData` in fields and tests. + * Moved from `Data` to `ExtensionData` in fields and tests. * Improved `MinUpDownTimeNode`: * Rewrote the function `constraints_capacity` based on core functionality. * Updated the test set to identify all potential problems with the node. diff --git a/Project.toml b/Project.toml index b447233..a38acc4 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,7 @@ JuMP = "4076af6c-e467-56ae-b986-b466b2749572" TimeStruct = "f9ed5ce0-9f41-4eaa-96da-f38ab8df101c" [compat] -EnergyModelsBase = "0.10" +EnergyModelsBase = "0.10.5" EnergyModelsRenewableProducers = "0.7" JuMP = "1.23" TimeStruct = "0.9" diff --git a/README.md b/README.md index c2ca56b..c176a85 100644 --- a/README.md +++ b/README.md @@ -13,9 +13,9 @@ > The different node types are partly experimental. > They have furthermore some limitations with respect to the chosen `TimeStructure` or whether they are able to handle investments. > As a consequence, it is advised to read the documentation for each node to identify their usefulness. -> Is is planned to removesome nodes and rewrite the behaviour of other nodes to improve their flexibility. +> Is is planned to remove some nodes and rewrite the behaviour of other nodes to improve their flexibility. > -> Version 0.2.8 is based on [`EnergyModelsBase` v0.8](https://energymodelsx.github.io/EnergyModelsBase.jl/v0.8/) while version 0.2.9 is adding support for [`EnergyModelsBase` v0.9](https://energymodelsx.github.io/EnergyModelsBase.jl/v0.9/). +> Among others, using `PeriodDemandSink` in combination with `EnergyModelsGUI` results in errors when trying to access field values with `PartitionProfile`. ## Usage @@ -26,7 +26,7 @@ Please refer to the *[documentation](https://energymodelsx.github.io/EnergyModel ## Cite -If you find `EnergyModelsFlex` useful in your work, we kindly request that you cite the following [publication](https://doi.org/10.69997/sct.119164): +If you find `EnergyModelsFlex` useful in your work, we kindly request that you cite the following *[publication](https://doi.org/10.69997/sct.119164)*: ```bibtex @article{venas2025energymodelsflex, diff --git a/docs/make.jl b/docs/make.jl index 73cd771..a86e54d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -74,8 +74,10 @@ makedocs( "Links" => Any[ "CapacityCostLink"=>"links/capacitycostlink.md", ], - "How-to" => - Any["Contribute"=>"how-to/contribute.md"], + "How-to" => Any[ + "Contribute"=>"how-to/contribute.md", + "Update models"=>"how-to/update-models.md", + ], "Library" => Any[ "Public"=>"library/public.md", "Internals"=>String[ diff --git a/docs/src/how-to/update-models.md b/docs/src/how-to/update-models.md new file mode 100644 index 0000000..36228a3 --- /dev/null +++ b/docs/src/how-to/update-models.md @@ -0,0 +1,33 @@ +# [Update your model to the latest versions](@id how_to-update) + +`EnergyModelsFlex` is still in a pre-release version. +Hence, there are frequently breaking changes occuring, although we plan to keep backwards compatibility. +This document is designed to provide users with information regarding how they have to adjust their nodes to keep compatibility to the latest changes. + +## [Adjustments from 0.3.0](@id how_to-update-03) + +The introduction of `PartitionProfile` in *[`TimeStruct` v0.9.12](https://github.com/sintefore/TimeStruct.jl/releases/tag/v0.9.12)* allowed a rewrite of `PeriodDemandSink`: + +```julia +# The previous nodal description for a `PeriodDemandSink` was given by: +PeriodDemandSink( + id::Any, + period_length::Int, + period_demand::Array{<:Real}, + cap::TimeProfile, + penalty::Dict{Symbol,<:TimeProfile}, + input::Dict{<:Resource,<:Real}, + data::Vector{<:ExtensionData} = ExtensionData[], +) + +# This translates to the following new version +PeriodDemandSink( + id, + cap, + FixedProfile(period_length), + PartitionProfile(period_demand), + penalty, + input, + data +) +``` diff --git a/docs/src/index.md b/docs/src/index.md index 27530f0..1ac2444 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -64,8 +64,8 @@ Depth = 1 ```@contents Pages = [ + "how-to/update-models.md", "how-to/contribute.md", - "how-to/utilize.md", ] Depth = 1 ``` diff --git a/docs/src/library/internals/methods-EMF.md b/docs/src/library/internals/methods-EMF.md index a181447..5be9842 100644 --- a/docs/src/library/internals/methods-EMF.md +++ b/docs/src/library/internals/methods-EMF.md @@ -9,7 +9,6 @@ Pages = ["methods-EMF.md"] ## [Check methods](@id lib-int-met-check) ```@docs -EMF.check_period_ts EMF.check_limits_default EMF.check_input EMF.check_cap_price_periods diff --git a/docs/src/library/internals/methods-fields.md b/docs/src/library/internals/methods-fields.md index 63ae180..5b6bf68 100644 --- a/docs/src/library/internals/methods-fields.md +++ b/docs/src/library/internals/methods-fields.md @@ -11,8 +11,8 @@ Pages = ["methods-fields.md"] ```@docs EnergyModelsFlex.period_demand -EnergyModelsFlex.period_length -EnergyModelsFlex.period_index +EnergyModelsFlex.period_duration +EnergyModelsFlex.periods EnergyModelsFlex.number_of_periods ``` @@ -35,4 +35,4 @@ EnergyModelsFlex.cap_price ```@docs EnergyModelsFlex.limits EnergyModelsFlex.heat_resource -``` \ No newline at end of file +``` diff --git a/docs/src/nodes/sink/perioddemand.md b/docs/src/nodes/sink/perioddemand.md index ef8fc6f..55db89e 100644 --- a/docs/src/nodes/sink/perioddemand.md +++ b/docs/src/nodes/sink/perioddemand.md @@ -10,7 +10,12 @@ This node can, *e.g.*, be combined with [`MinUpDownTimeNode`](@ref), to allow pr !!! warning "TimeStructure for node" This node is designed for **uniform or repetitive duration of operational periods**. - Irregular durations may cause misalignment of shifted loads, especially if the field `period_length` does not align with the chosen [`SimpleTimes`](@extref TimeStruct.SimpleTimes) structure representing the operational periods. + Irregular durations may cause misalignment of shifted loads, especially if the field `period_duration` does not align with the chosen [`SimpleTimes`](@extref TimeStruct.SimpleTimes) structure representing the operational periods. + +!!! warning "`PeriodDemandSink` and `EnergyModelsGUI`" + Some of the fields of this node cannot be represented in `EnergyModelsGUI`. + The reason for that limitation is that `EnergyModelsGUI` does not yet support partitions of `TimePeriod`s. + `EnergyModelsGUI` can still be utilized for all other fields. ## [Introduced type and its fields](@id nodes-perioddemandsink-fields) @@ -33,12 +38,11 @@ The standard fields are given as: - **`penalty::Dict{Symbol,<:TimeProfile}`**:\ The penalty dictionary is used for providing penalties for soft constraints to allow for both over and under delivering the demand.\ It must include the fields `:surplus` and `:deficit`. - In addition, it is crucial that the sum of both values is larger than 0 to avoid an unconstrained model. + In addition, it is crucial that the sum of both values in each demand period is larger than 0 to avoid an unconstrained model. !!! warning "Chosen values" - The current implementation does not represent the proper cost due to the summation. - Instead, you must consider the duration of an operational period and the field `period_length` when providing a value. - In this case, the value should be multiplied by ``1/period\_length(n) \times duration(t)``. + The implementation is relative to the chosen `period_duration` (see below). + If the period duration is ``24``, then the cost is for the unsatisfied demand within the ``24`` demand period, multiplied with the probability and the repetitons within a strategic period. - **`input::Dict{<:Resource,<:Real}`**:\ The field `input` includes [`Resource`](@extref EnergyModelsBase.Resource)s with their corresponding conversion factors as dictionaries.\ @@ -60,18 +64,27 @@ The standard fields are given as: [`AbstractPeriodDemandSink`](@ref EnergyModelsFlex.AbstractPeriodDemandSink)s require additional fields to specify both the periods and their respective demands: -- **`period_length::Int`**:\ - Defines how many operational periods are included in a single demand period.\ - For instance, if the duration of the operational periods is 1 hour and `period_length = 24`, then each demand period spans one day. - The demand of this node (for a given day, see below) must then be filled on a daily basis, without any restrictions on *when* during the day the demand must be filled. +- **`period_duration::TimeProfile`**:\ + Defines the total duration of a single demand period.\ + For instance, if the duration of 1 of the operational time structure is 1 hour and `period_duration = FixedProfile(24)`, then each demand period spans one day. + The demand of this node (for a given day, see below) must then be filled on a daily basis, without any restrictions on *when* during the day the demand must be filled given the available capacity.\ + Due to a constructor, it can either be specified as number (the same duration in all demand periods), as a vector (varying duration of each demand period), or as a time profile (*e.g.*, varying period durations due to varying operational time structures). + It cannot be specified as `OperationalProfile`. + +- **`period_demand::TimeProfile`**:\ + The total demand to be met during each demand period. + The length of this time profile should match the number of periods (*e.g.*, days) in the time structure. + If the time structure represents one year with hourly resolution and the demand periods correspond to a day, this time profile must then have 365 elements. -- **`period_demand::Vector{<:Real}`**:\ - The total demand to be met during each demand period. The length of this vector should match the number of periods (*e.g.*, days) in the time structure. If the time structure represents on year with hourly resolution, this vector must then have 365 elements. + It is best to utilize the [`PartitionProfile`](@extref TimeStruct.PartitionProfile) type if the demand is varying. + If it is constant, you can also utilize [`StrategicProfile`][@extref TimeStruct.StrategicProfile], [`RepresentativeProfile`][@extref TimeStruct.RepresentativeProfile], or [`ScenarioProfile`][@extref TimeStruct.ScenarioProfile], depending on your chosen time structure. + It cannot be specified as `OperationalProfile`. !!! warning "Time consistency" - Ensure that the `period_demand` vector length aligns with the total time horizon divided by `period_length`. Mismatches can lead to indexing errors or inconsistent demand enforcement. + Ensure that the `period_demand` time profile length aligns with the operational time horizon duration divided by `period_duration` + Mismatches can lead to indexing errors or inconsistent demand enforcement. -These fields are at the 2ⁿᵈ and 3ʳᵈ position below the field `id` as shown in [`PeriodDemandSink`](@ref). +These fields are at the 3ʳᵈ and 4ᵗʰ position below the field `cap` as shown in [`PeriodDemandSink`](@ref). ## [Mathematical description](@id nodes-perioddemandsink-math) @@ -105,12 +118,12 @@ The variables include: #### [Additional variables](@id nodes-perioddemandsink-math-add) [`AbstractPeriodDemandSink`](@ref EnergyModelsFlex.AbstractPeriodDemandSink) nodes declare in addition several variables through dispatching on the method [`EnergyModelsBase.variables_element()`](@ref) for including constraints for deficits and surplus for individual resources as well as what the fraction satisfied by each resource. -These variables are for a [`AbstractPeriodDemandSink`](@ref EnergyModelsFlex.AbstractPeriodDemandSink) node ``n`` in demand periods ``i``: +These variables are for a [`AbstractPeriodDemandSink`](@ref EnergyModelsFlex.AbstractPeriodDemandSink) node ``n`` in demand periods ``t_pd``: -- ``\texttt{demand\_sink\_surplus}[n, i]``:\ - Surplus of energy delivered beyond the required `period_demand` in demand period `i`. -- ``\texttt{demand\_sink\_deficit}[n, i]``:\ - Deficit of energy delivered relative to the `period_demand` in period `i`. +- ``\texttt{demand\_sink\_surplus}[n, t_pd]``:\ + Surplus of energy delivered beyond the required `period_demand` of demand period `t_pd` . +- ``\texttt{demand\_sink\_deficit}[n, t_pd]``:\ + Deficit of energy delivered relative to the `period_demand` of demand period `t_pd` . ### [Constraints](@id nodes-perioddemandsink-math-con) @@ -155,33 +168,37 @@ These standard constraints are: - `constraints_data`:\ This function is only called for specified additional data, see above. -The function `constraints_capacity` is extended with a new method to account for the calculation of the period demand deficit and surplus through: +The function `constraints_capacity` is extended with a new method to account for the calculation of the period demand deficit and surplus. + +The overall balance is modified as ```math -\texttt{cap\_use}[n, t] + \texttt{sink\_deficit}[n, t] = \texttt{cap\_inst}[n, t] + \texttt{sink\_surplus}[n, t] +\texttt{cap\_use}[n, t] + \texttt{sink\_deficit}[n, t] = \texttt{cap\_inst}[n, t] ``` +while operational period surplus ``\texttt{sink\_surplus}[n, t]`` is fixed to 0. + +The surplus and deficit of the demand period can then be calculated as + ```math \begin{aligned} -\texttt{demand\_sink\_deficit}[n, i] + & \sum_{t \in P_i} \texttt{​cap\_use}[n,t] = \\ -& \texttt{demand\_sink\_surplus}[n, i] + period\_demand(n, i) +\texttt{demand\_sink\_deficit}[n, t_{pd}] + & \sum_{t \in t_{pd}} \texttt{​cap\_use}[n, t] \times duration(t) = \\ +& \texttt{demand\_sink\_surplus}[n, t_{pd}] + period\_demand(n, t_{pd}) \end{aligned} ``` -where ``P_i`` is the set of operational periods in demand period ``i``. +where ``t_{pd}`` is the demand period consisting of a set of operational periods. As a consequence, `constraints_opex_var` requires as well a new method as we only consider the deficit within a complete period: ```math \begin{aligned} -\texttt{opex\_var}[n, t_{inv}] = \sum_{t ∈ t_{inv}}(& \texttt{demand\_sink\_surplus}[n, i_t] \times \texttt{surplus\_penalty}(n, t) + \\ -& \texttt{demand\_sink\_deficit}[n, i_t] \times \texttt{deficit\_penalty}(n, t)) \times \\ -& scale\_op\_sp(t_{inv}, t) +\texttt{opex\_var}[n, t_{inv}] = \sum_{t_{pd} ∈ periods(t_{inv})}(& \texttt{demand\_sink\_surplus}[n, t_{pd}] \times \texttt{surplus\_penalty}(n, t_{pd}) + \\ +& \texttt{demand\_sink\_deficit}[n, t_{pd}] \times \texttt{deficit\_penalty}(n, t_{pd})) \times \\ +& scale\_op\_sp(t_{inv}, first(t_{pd})) / duration(first(t_{pd})) \end{aligned} ``` -where ``i_t`` is the period index such that ``t \in P_{i_t}``. - -!!! tip "The function `scale_op_sp`" - The function [``scale\_op\_sp(t_{inv}, t)``](@extref EnergyModelsBase.scale_op_sp) calculates the scaling factor between operational and investment periods. - It also takes into account potential operational scenarios and their probability as well as representative periods. +!!! note "`scale_op_sp` and `duration`" + The function [`scale_op_sp(t_inv, t)`](@extref EnergyModelsBase.scale_op_sp) calculates the scaling factor between operational and investment periods including the `duration` of each operational period. + It hence must be divided by the `duration` of the first operational period to avoid including the duration of the operational period. diff --git a/examples/flexible_demand.jl b/examples/flexible_demand.jl index 546a5f2..54dc1ec 100644 --- a/examples/flexible_demand.jl +++ b/examples/flexible_demand.jl @@ -28,7 +28,7 @@ CO2 = ResourceEmit("CO2", 0) # This would be a huge incentive to produce during the weekend, if we allowed the # `PeriodDemandSink` capacity during the weekend. day = [1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 9, 8, 7, 6, 5, 4, 3, 2] -el_cost = [repeat(day, 5)..., fill(0, 2 * 24)...] +el_cost = vcat(repeat(day, 5), zeros(48)) grid = RefSource( "grid", @@ -41,16 +41,16 @@ grid = RefSource( # The production can only run between 6-20 on weekdays, with a capacity of 300 kW. # First, define the maximum capacity for a regular weekday (24 hours). # The capacity is 0 between 0 am and 6 am, 300 kW between 6 am and 8 pm, and 0 again between 8 pm and midnight. -weekday_prod = [fill(0, 6)..., fill(300, 14)..., fill(0, 4)...] -@assert length(weekday_prod) == 24 +weekday_prod = vcat(zeros(6), fill(300, 14), zeros(4)) + # Repeat a weekday 5 times, for a workweek, then no production on the weekends. -week_prod = [repeat(weekday_prod, 5)..., fill(0, 2 * 24)...] +week_prod = vcat(repeat(weekday_prod, 5), fill(0, 2 * 24)) demand = PeriodDemandSink( "demand_product", - 24, # 24 hours per day. - [fill(1500, 5)..., 0, 0], # Demand of 1500 units per day, and nothing (0) in the weekend. OperationalProfile(week_prod), # kW - installed capacity + 24, # 24 hours per day. + PartitionProfile(vcat(ones(5)*1500, [0, 0])), # Demand of 1500 units per day, and nothing (0) in the weekend. Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e8)), # € / Demand - Price for not delivering products Dict(Product => 1), ) diff --git a/src/EnergyModelsFlex.jl b/src/EnergyModelsFlex.jl index b5ed038..55c455d 100644 --- a/src/EnergyModelsFlex.jl +++ b/src/EnergyModelsFlex.jl @@ -25,6 +25,8 @@ for node_type ∈ ["source", "sink", "network", "storage", "link"] include("$node_type/checks.jl") end +include("legacy_constructors.jl") + export MinUpDownTimeNode, ActivationCostNode, ElectricBattery, LoadShiftingNode export PeriodDemandSink, MultipleInputSink export PayAsProducedPPA, StorageEfficiency, LimitedFlexibleInput, Combustion diff --git a/src/legacy_constructors.jl b/src/legacy_constructors.jl new file mode 100644 index 0000000..8c51571 --- /dev/null +++ b/src/legacy_constructors.jl @@ -0,0 +1,35 @@ +## Sinks + +function PeriodDemandSink( + id::Any, + period_length::Int, + period_demand::Array{<:Real}, + cap::TimeProfile, + penalty::Dict{Symbol,<:TimeProfile}, + input::Dict{<:Resource,<:Real}, + data::Vector{<:ExtensionData} = ExtensionData[], +) + @warn( + "The used implementation of a `PeriodDemandSink` will be discontinued in the near " * + "future. The new implementation uses the approach of `PeriodPartition` with " * + "a change in the field positions.\n" * + "In practice, three changes have to be incorporated: \n 1. the capacity is moved to " * + "the 2ⁿᵈ position, \n 2. `period_length` is renamed `period_duration`, moved to the " * + "3ʳᵈ position, and can accept as well a `Vector` or `TimeProfile`s as input, and\n" * + " 3. `period_demand` is moved to the 4ᵗʰ position and requires as input a " * + "`PartitionProfile` of the previously provided `Vector`.\n" * + "See the documentation (https://energymodelsx.github.io/EnergyModelsFlex.jl/stable/how-to/update-models/#Adjustments-from-0.3.0) " * + "on how to update your model to the latest version.", + maxlog = 1 + ) + + return PeriodDemandSink( + id, + cap, + FixedProfile(period_length), + PartitionProfile(period_demand), + penalty, + input, + data + ) +end diff --git a/src/sink/checks.jl b/src/sink/checks.jl index e493d52..3493275 100644 --- a/src/sink/checks.jl +++ b/src/sink/checks.jl @@ -11,11 +11,13 @@ additional checks on the data. - The field `cap` is required to be non-negative. - The values of the dictionary `input` are required to be non-negative. - The dictionary `penalty` is required to have the keys `:deficit` and `:surplus`. +- The values `:deficit` and `:surplus` of the dictionary `penalty` are required to be + indexable by a `PeriodPartition`. - The sum of the values `:deficit` and `:surplus` in the dictionary `penalty` has to be non-negative to avoid an infeasible model. -- The remainder of the divison of the lowest time structure by the period length must be 0. -- The length of the period demand must equal the length of the lowest period times the - parameter period length. +- The individual periods must all satisfy the specified duration(s). +- The field `period_demand` is required to be non-negative and indexable by a + `PeriodPartition`. """ function EMB.check_node( n::PeriodDemandSink, @@ -23,52 +25,59 @@ function EMB.check_node( modeltype::EnergyModel, check_timeprofiles::Bool, ) - EMB.check_node_default(n, 𝒯, modeltype, check_timeprofiles) + 𝒯ᵖᵈ = periods(n, 𝒯) + per_dur = period_duration(n) + bool = true - # Check that the period length and demand is working with the time structure - # The check will only be activated in 0.3 - # for (idx_sp, ts_oper) ∈ enumerate(𝒯.operational) - # sub_msg = "the operational time structure in strategic period $(idx_sp)" - # check_period_ts(ts_oper, n, sub_msg) - # end -end - -""" - check_period_ts(ts::RepresentativePeriods, n::PeriodDemandSink, msg::String) - check_period_ts(ts::OperationalScenarios, n::PeriodDemandSink, msg::String) - check_period_ts(ts::SimpleTimes, n::PeriodDemandSink, msg::String) - -Function for checking that the timestructure is valid in combination with the chosen period -structure in a [`PeriodDemandSink`(@ref). -""" -function check_period_ts(ts::RepresentativePeriods, n::PeriodDemandSink, msg::String) - for (idx, ts_oper) ∈ enumerate(ts.rep_periods) - sub_msg = msg * " in representative period $(idx)" - check_period_ts(ts_oper, n, sub_msg) - end -end -function check_period_ts(ts::OperationalScenarios, n::PeriodDemandSink, msg::String) - for (idx, ts_oper) ∈ enumerate(ts.scenarios) - sub_msg = msg * " in operational scenario $(idx)" - check_period_ts(ts_oper, n, sub_msg) - end -end -function check_period_ts(ts::SimpleTimes, n::PeriodDemandSink, msg::String) - len = period_length(n) - n_per = number_of_periods(n) @assert_or_log( - length(ts)%len == 0, - "The specified period length does not work with $(msg)." + all(capacity(n, t) ≥ 0 for t ∈ 𝒯), + "The capacity must be non-negative." ) @assert_or_log( - length(ts)/len ≤ n_per, - "The vector `period_demand` is shorter than the $(msg)." + all(inputs(n, p) ≥ 0 for p ∈ inputs(n)), + "The values for the Dictionary `input` must be non-negative." ) - if length(ts)%len == 0 & length(ts)/len < n_per - @warn( - "The vector `period_demand` is longer than required in $(msg). " * - "The last $(Int(n_per-length(ts)/len)) values will be omitted.", - maxlog=1, + @assert_or_log( + :surplus ∈ keys(n.penalty) && :deficit ∈ keys(n.penalty), + "The entries `:surplus` and `:deficit` are required in the dictionary `penalty`." + ) + if :surplus ∈ keys(n.penalty) + message = "are not allowed for the key `:surplus` in the dictionary `penalty`." + bool *= EMB.check_partition_profile(surplus_penalty(n), message) + else + bool = false + end + if :deficit ∈ keys(n.penalty) + message = "are not allowed for the key `:deficit` in the dictionary `penalty`." + bool *= EMB.check_partition_profile(deficit_penalty(n), message) + else + bool = false + end + + if bool + @assert_or_log( + all(surplus_penalty(n, t_pd) + deficit_penalty(n, t_pd) ≥ 0 for t_pd ∈ 𝒯ᵖᵈ), + "An inconsistent combination of `:surplus` and `:deficit` leads to an infeasible model." + ) + end + + message = "are not allowed for the field `:period_duration`." + bool = EMB.check_partition_profile(period_duration(n), message) + if bool + @assert_or_log( + all(sum(duration(t) for t ∈ t_pd) ≥ per_dur[t_pd] for t_pd ∈ 𝒯ᵖᵈ), + "The duration of the last period on the `SimpleTimes` level is shorter than " * + "specified. This is caused by inconsistently specified `period_duration` and" * + "time structure." + ) + end + + message = "are not allowed for the field `:period_demand`." + bool = EMB.check_partition_profile(period_demand(n), message) + if bool + @assert_or_log( + all(period_demand(n, t_pd) ≥ 0 for t_pd ∈ 𝒯ᵖᵈ), + "The period demand must be non-negative." ) end end diff --git a/src/sink/constraint_functions.jl b/src/sink/constraint_functions.jl index b2fc180..f0444b9 100644 --- a/src/sink/constraint_functions.jl +++ b/src/sink/constraint_functions.jl @@ -13,42 +13,28 @@ function EMB.constraints_capacity( 𝒯::TimeStructure, modeltype::EnergyModel, ) - @constraint( - m, - [t ∈ 𝒯], - m[:cap_use][n, t] + m[:sink_deficit][n, t] == - m[:cap_inst][n, t] + m[:sink_surplus][n, t] - ) + # Declaration of the required subsets. + 𝒯ᵖᵈ = periods(n, 𝒯) - # Need to constraint the used capacity to the installed capacity per - # operational period. Instead, the node may get input in operational periods - # when the cap field is 0. This is ok for regular sink nodes, but this node - # only penalizes surplus or deficit over a period. @constraint( m, [t ∈ 𝒯], - m[:cap_use][n, t] <= m[:cap_inst][n, t] + m[:cap_use][n, t] + m[:sink_deficit][n, t] == m[:cap_inst][n, t] ) - # Create a list mapping the demand period i to the operational periods it contains. - num_periods = number_of_periods(n, 𝒯) - period2op = [[] for i ∈ 1:num_periods] + # Fix the surplus to 0 for t ∈ 𝒯 - period_id = period_index(n, t) - push!(period2op[period_id], t) + fix(m[:sink_surplus][n, t], 0; force = true) end - for i ∈ 1:num_periods - # Sum all values inside period i. - period_total = sum(m[:cap_use][n, t] for t ∈ period2op[i]) - # Define the demand_sink_deficit as the difference between the period demand and - # the total capacity used. - @constraint( - m, - period_total + m[:demand_sink_deficit][n, i] == - period_demand(n, i) + m[:demand_sink_surplus][n, i] - ) - end + # Set the energy balance for the partition duration + @constraint( + m, + [t_pd ∈ 𝒯ᵖᵈ], + m[:demand_sink_deficit][n, t_pd] + + sum(m[:cap_use][n, t] * duration(t) for t ∈ t_pd) == + m[:demand_sink_surplus][n, t_pd] + period_demand(n, t_pd) + ) EMB.constraints_capacity_installed(m, n, 𝒯, modeltype) end @@ -58,8 +44,9 @@ end Function for creating the constraint on the variable OPEX of an [`AbstractPeriodDemandSink`](@ref). -The method is adjusted from the default method through utilizing the period demand surplus -and deficit instead of the operational period deficit or surplus. +The method is changed from the standard approach through calculating both the demand period +surplus and deficit in addition to the operational period deficit. The operational period +surplus is fixed to 0 to avoid problems in the calculations. """ function EMB.constraints_opex_var(m, n::AbstractPeriodDemandSink, 𝒯ᴵⁿᵛ, ::EnergyModel) # Only penalise the total surplus and deficit in each period, not in the @@ -68,10 +55,11 @@ function EMB.constraints_opex_var(m, n::AbstractPeriodDemandSink, 𝒯ᴵⁿᵛ, m, [t_inv ∈ 𝒯ᴵⁿᵛ], m[:opex_var][n, t_inv] == sum( - ( - m[:demand_sink_surplus][n, period_index(n, t)] * surplus_penalty(n, t) + - m[:demand_sink_deficit][n, period_index(n, t)] * deficit_penalty(n, t) - ) * scale_op_sp(t_inv, t) for t ∈ t_inv + ( + m[:demand_sink_surplus][n, t_pd] * surplus_penalty(n, t_pd) + + m[:demand_sink_deficit][n, t_pd] * deficit_penalty(n, t_pd) + ) * scale_op_sp(t_inv, first(t_pd)) / duration(first(t_pd)) + for t_pd ∈ periods(n, t_inv) ) ) end diff --git a/src/sink/datastructures.jl b/src/sink/datastructures.jl index cbc9442..6757cd6 100644 --- a/src/sink/datastructures.jl +++ b/src/sink/datastructures.jl @@ -26,77 +26,100 @@ abstract type AbstractMultipleInputSinkStrat <: AbstractMultipleInputSink end A `PeriodDemandSink` is a [`Sink`](@extref EnergyModelsBase.Sink) that has a demand that can be satisfied any time during a period of defined length. If the chosen time structure has -operational periods of a duration of 1 hour and the demand should be fulfilled daily, -`period_length` should be 24. The demand for each day is then set as an array as the +operational periods of a duration of 1 hour and the demand should be fulfilled daily, +`period_duration` should be 24. The demand for each day is then set as a time profile in the `period_demand` field. The `cap` field is the maximum capacity that can be fulfilled in each operational period. # Fields - **`id::Any`** is the name/identifier of the node. -- **`period_length::Int`** is the number of periods in which the period demand can be - satisfied. -- **`period_demand::Array{<:Real}`** is the demand within each of the periods. - **`cap::TimeProfile`** is the installed capacity. +- **`period_duration::TimeProfile`** is the sum of the durations of the individual + operational periods within a given demand period. Due to a constructor, it can either be + specified as number (the same duration in all demand periods), as a vector (varying + duration of each demand period), or as a time profile (*e.g.*, varying period durations + due to varying operational time structures). It cannot be specified as `OperationalProfile`. +- **`period_demand::TimeProfile`** is the demand within each of the periods as time profile. + It cannot be specified as `OperationalProfile`. - **`penalty::Dict{Symbol,<:TimeProfile}`** are penalties for surplus or deficits. The - dictionary requires the fields `:surplus` and `:deficit`. + dictionary requires the fields `:surplus` and `:deficit`. - **`input::Dict{<:Resource,<:Real}`** are the input [`Resource`](@extref EnergyModelsBase.Resource)s with conversion value `Real`. -- **`data::Vector{<:ExtensionData}`** is the additional data (*e.g.*, for investments). The field `data` - is conditional through usage of a constructor. +- **`data::Vector{<:ExtensionData}`** is the additional data (*e.g.*, for investments). The + field `data` is conditional through usage of a constructor. """ struct PeriodDemandSink <: AbstractPeriodDemandSink id::Any - period_length::Int - period_demand::Array{<:Real} cap::TimeProfile + period_duration::TimeProfile + period_demand::TimeProfile penalty::Dict{Symbol,<:TimeProfile} input::Dict{<:Resource,<:Real} data::Vector{<:ExtensionData} end function PeriodDemandSink( id, - period_length::Int, - period_demand::Vector{<:Real}, cap::TimeProfile, + period_duration::Union{Number, Vector{<:Number}}, + period_demand::TimeProfile, + penalty::Dict{Symbol,<:TimeProfile}, + input::Dict{<:Resource,<:Real}, + data::Vector{<:ExtensionData}, +) + if isa(period_duration, Number) + per_dur = FixedProfile(period_duration) + elseif isa(period_duration, Vector{<:Number}) + per_dur = PartitionProfile(period_duration) + end + return PeriodDemandSink(id, cap, per_dur, period_demand, penalty, input, data) +end +function PeriodDemandSink( + id, + cap::TimeProfile, + period_duration::Union{Number, Vector{<:Number}, TimeProfile}, + period_demand::TimeProfile, penalty::Dict{Symbol,<:TimeProfile}, input::Dict{<:Resource,<:Real}, ) - PeriodDemandSink(id, period_length, period_demand, cap, penalty, input, ExtensionData[]) + return PeriodDemandSink(id, cap, period_duration, period_demand, penalty, input, ExtensionData[]) end """ period_demand(n::AbstractPeriodDemandSink) - period_demand(n::AbstractPeriodDemandSink, i::Int) + period_demand(n::AbstractPeriodDemandSink, t_pd::TS.PeriodPartition) -Returns the period demand of `AbstractPeriodDemandSink` `n` as Array or in demand period `i`. +Returns the period demands of `AbstractPeriodDemandSink` `n` as a `TimeProfile` or in +demand period `t_pd`. """ period_demand(n::AbstractPeriodDemandSink) = n.period_demand -period_demand(n::AbstractPeriodDemandSink, i) = n.period_demand[i] +period_demand(n::AbstractPeriodDemandSink, t_pd::TS.PeriodPartition) = + n.period_demand[t_pd] """ - period_length(n::AbstractPeriodDemandSink) + period_duration(n::AbstractPeriodDemandSink) -Returns the length of the demand period of `AbstractPeriodDemandSink` `n`. +Returns the demand periods of `AbstractPeriodDemandSink` `n` as `TimeProfile` or in demand +period `t_pd`. """ -period_length(n::AbstractPeriodDemandSink) = n.period_length +period_duration(n::AbstractPeriodDemandSink) = n.period_duration +period_duration(n::AbstractPeriodDemandSink, t_pd::TS.PeriodPartition) = + n.period_duration[t_pd] """ - number_of_periods(n::AbstractPeriodDemandSink) - number_of_periods(n::AbstractPeriodDemandSink, 𝒯::TimeStructure) + periods(n::AbstractPeriodDemandSink, ts::TS.TimeStructure) -Returns the number of periods for a `PeriodDemandSink` `n`. If a `TimeStructure` is provided -it calculates it based on the chosen time structure. +Returns the demand periods for a `PeriodDemandSink` `n` for the given time structure. """ -number_of_periods(n::AbstractPeriodDemandSink) = length(period_demand(n)) -number_of_periods(n::AbstractPeriodDemandSink, 𝒯::TimeStructure) = - Int(length(𝒯) / period_length(n)) +periods(n::AbstractPeriodDemandSink, ts::TS.TimeStructure) = + partition_duration(ts, period_duration(n)) """ - period_index(n::AbstractPeriodDemandSink, t) + number_of_periods(n::AbstractPeriodDemandSink, ts::TS.TimeStructure) -Returns the index of the period (*e.g.*, day) that a operational period `t` belongs to. +Returns the number of demand periods for a `PeriodDemandSink` `n` for the given time structure. """ -period_index(n::AbstractPeriodDemandSink, t) = Int(ceil(t.period.op / period_length(n))) +number_of_periods(n::AbstractPeriodDemandSink, ts::TS.TimeStructure) = + length(periods(n, ts)) """ struct MultipleInputSink <: AbstractMultipleInputSink diff --git a/src/sink/model.jl b/src/sink/model.jl index 9272e44..7cdb35d 100644 --- a/src/sink/model.jl +++ b/src/sink/model.jl @@ -1,11 +1,12 @@ """ EMB.variables_element(m, 𝒩ˢⁱⁿᵏ::Vector{<:AbstractPeriodDemandSink}, 𝒯, ::EnergyModel) -Creates the following additional variables for **ALL** [`PeriodDemandSink`](@ref) nodes: -- `demand_sink_surplus[n, i]` is a non-negative variable indicating a surplus in demand in - each period `i`. -- `demand_sink_deficit[n, i]` is a non-negative variable indicating a deficit in demand for - each period `i`. +Creates the following additional variables for **ALL** [`AbstractPeriodDemandSink`](@ref) +nodes: +- `demand_sink_surplus[n, t_pd]` is a non-negative variable indicating a surplus in demand in + each demand period `t_pd`. +- `demand_sink_deficit[n, t_pd]` is a non-negative variable indicating a deficit in demand for + each demand period `t_pd`. !!! note "Definition of period" The period in the description above does not correspond to an operational period as known @@ -18,10 +19,11 @@ function EMB.variables_element( 𝒯, ::EnergyModel, ) - @variable(m, demand_sink_surplus[n ∈ 𝒩ˢⁱⁿᵏ, i=1:number_of_periods(n, 𝒯)] ≥ 0) - @variable(m, demand_sink_deficit[n ∈ 𝒩ˢⁱⁿᵏ, i=1:number_of_periods(n, 𝒯)] ≥ 0) + @variable(m, demand_sink_surplus[n ∈ 𝒩ˢⁱⁿᵏ, periods(n, 𝒯)] ≥ 0) + @variable(m, demand_sink_deficit[n ∈ 𝒩ˢⁱⁿᵏ, periods(n, 𝒯)] ≥ 0) end + """ EMB.variables_element(m, 𝒩::Vector{<:AbstractMultipleInputSinkStrat}, 𝒯, ::EnergyModel) diff --git a/test/network/test_MinUpDownTimeNode.jl b/test/network/test_MinUpDownTimeNode.jl index be18f4e..dcfce3e 100644 --- a/test/network/test_MinUpDownTimeNode.jl +++ b/test/network/test_MinUpDownTimeNode.jl @@ -18,7 +18,7 @@ function test_case_minupdown( 𝒫 = [power, product, CO2] day = [1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 9, 8, 7, 6, 5, 4, 3, 2] - el_cost = [repeat(day, 5)..., fill(0, 2 * 24)...] + el_cost = vcat(repeat(day, 5), zeros(48)) src = RefSource( "src", @@ -43,15 +43,15 @@ function test_case_minupdown( # The demand can only be satisfied between 6-20 on weekdays, with a capacity of 200. # No production on weekends. - weekday_prod = [fill(0, 6)..., fill(200, 14)..., fill(0, 4)...] - week_prod = [repeat(weekday_prod, 5)..., fill(0, 2 * 24)...] + weekday_prod = vcat(zeros(6), fill(300, 14), zeros(4)) + week_prod = vcat(repeat(weekday_prod, 5), fill(0, 2 * 24)) snk = PeriodDemandSink( "demand_product", - 24, # 24 hours per day - [fill(1500, 5)..., 0, 0], # Demand 1500 units per day, and nothing (0) in the weekend. - OperationalProfile(week_prod), # kW - installed capacity - Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e8)), # € / Demand - Price for not delivering products + OperationalProfile(week_prod), + 24, + PartitionProfile([fill(1500, 5)..., 0, 0]), + Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e4)), # € / Demand - Price for not delivering products Dict(product => 1), ) diff --git a/test/sink/test_PeriodDemandSink.jl b/test/sink/test_PeriodDemandSink.jl index ed7ba4a..67b8d4f 100644 --- a/test/sink/test_PeriodDemandSink.jl +++ b/test/sink/test_PeriodDemandSink.jl @@ -1,30 +1,74 @@ -using EnergyModelsBase -using EnergyModelsFlex - -using HiGHS -using JuMP -using Test -using TimeStruct - # Resources used in the tests Power = ResourceCarrier("Power", 0) -CO2 = ResourceEmit("CO2", 0) +CO2 = ResourceEmit("CO₂", 0) + +function per_dem_snk_case(; + snk = nothing, + repr = false, + 𝒯 = TwoLevel( + 1, 1, + SimpleTimes(repeat(vcat([2, 2, 2], ones(14), [4]), 7)), + op_per_strat=8760. + ), +) + + # The production can only run between 6-20 on weekdays, with capacity of 200. + # No production on weekends. + weekday_prod = vcat(zeros(3), ones(14)*200, [0]) + price_day = [1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 9, 8, 6.5, 6, 3.5] -function create_system(demand; 𝒯 = TwoLevel(1, 1, SimpleTimes(7 * 24, 1))) - day = [1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 9, 8, 7, 6, 5, 4, 3, 2] - el_cost = [repeat(day, 5)..., fill(0, 2 * 24)...] + # Demand for 1500 units per day, and nothing (0) in the weekend with a maximum production + # of 200 per hour in between 6:00 and 20:00 + if repr + el_cost = RepresentativeProfile([ + OperationalProfile(price_day), + OperationalProfile(price_day), + OperationalProfile(price_day), + OperationalProfile(price_day), + FixedProfile(1e9), + FixedProfile(0), + FixedProfile(0), + ]) + week_prod = RepresentativeProfile([ + OperationalProfile(weekday_prod), + OperationalProfile(weekday_prod), + OperationalProfile(weekday_prod), + OperationalProfile(weekday_prod), + OperationalProfile(weekday_prod), + FixedProfile(0), + FixedProfile(0), + ]) + per_dem = RepresentativeProfile([fill(1500, 5)..., 0, 0]) + snk_sur = RepresentativeProfile(vcat([-8], zeros(6))) + else + el_cost = OperationalProfile(vcat(repeat(price_day, 4), fill(1e9, 18), zeros(36))) + week_prod = OperationalProfile(vcat(repeat(weekday_prod, 5), zeros(36))) + per_dem = PartitionProfile([fill(1500, 5)..., 0, 0]) + snk_sur = PartitionProfile(vcat([-8], zeros(6))) + end - grid = RefSource( + src = RefSource( "grid", FixedProfile(500), # kW - OperationalProfile(el_cost), + el_cost, FixedProfile(0), Dict(Power => 1), ) + if isnothing(snk) + snk = PeriodDemandSink( + "demand_product", + week_prod, + 24, + per_dem, + Dict(:surplus => snk_sur, :deficit => FixedProfile(1e4)), + Dict(Power => 1), + ) + end + 𝒫 = [Power, CO2] - 𝒩 = [grid, demand] - ℒ = [Direct("grid-demand", grid, demand)] + 𝒩 = [src, snk] + ℒ = [Direct("grid-snk", src, snk)] case = Case(𝒯, 𝒫, [𝒩, ℒ]) modeltype = OperationalModel( @@ -33,115 +77,309 @@ function create_system(demand; 𝒯 = TwoLevel(1, 1, SimpleTimes(7 * 24, 1))) CO2, ) m = create_model(case, modeltype) - return case, modeltype, m + return m, case, modeltype end -function create_demand_node() - # The production can only run between 6-20 on weekdays, with capacity of 200. - # No production on weekends. - weekday_prod = [fill(0, 6)..., fill(200, 14)..., fill(0, 4)...] - @assert length(weekday_prod) == 24 - week_prod = [repeat(weekday_prod, 5)..., fill(0, 2 * 24)...] +# Test that the fields of a `PeriodDemandSink` are correctly checked +# - EMB.check_node(n::PeriodDemandSink, 𝒯, modeltype::EnergyModel, check_timeprofiles::Bool) +@testset "Check functions" begin + # Set the global to true to suppress the error message + EMB.TEST_ENV = true - demand = PeriodDemandSink( - "demand_product", - # 24 hours per day. - 24, - # Produce 1500 units per day, and nothing (0) in the weekend. - [fill(1500, 5)..., 0, 0], - OperationalProfile(week_prod), # kW - installed capacity - Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e8)), # € / Demand - Price for not delivering products - Dict(Power => 1), + function check_per_dem_sink(; + cap = FixedProfile(10), + per_dur = 24, + per_dem = PartitionProfile([fill(1500, 5)..., 0, 0]), + penalty = Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e4)), + input = Dict(Power => 1), + 𝒯 = TwoLevel(2, 1, SimpleTimes(repeat(vcat([2, 2, 2], ones(14), [4]), 7))), ) - return demand -end + snk = PeriodDemandSink( + "demand_product", + cap, + per_dur, + per_dem, + penalty, + input, + ) + + return per_dem_snk_case(; snk, 𝒯) + end + # Test that a wrong capacity is caught by the checks + @test_throws AssertionError check_per_dem_sink(; cap=FixedProfile(-25)) -function test_max_grid_production(force_max_production::Bool) - demand = create_demand_node() - case, model, m = create_system(demand) + # Test that a wrong input is caught by the checks + @test_throws AssertionError check_per_dem_sink(; input = Dict(Power => -1)) - # Force the grid node to produce at max capacity. - source = get_nodes(case)[1] - @show source + # Test that a wrong penalty dictionary is caught + penalties = [ + Dict(:surplus => FixedProfile(0)), + Dict(:deficit => FixedProfile(0)), + Dict(:surplus => OperationalProfile([0]), :deficit => FixedProfile(1e4)), + Dict(:surplus => FixedProfile(0), :deficit => OperationalProfile([1e4])), + Dict(:surplus => FixedProfile(-1e5), :deficit => FixedProfile(1e4)), + ] + for penalty ∈ penalties + @test_throws AssertionError check_per_dem_sink(; penalty) + end + + # Test that a wrong period duration is caught by the checks, including in other time + # structures + @test_throws AssertionError check_per_dem_sink(; per_dur=25) + @test_throws AssertionError check_per_dem_sink(; per_dur=StrategicProfile([25, 24])) + week = SimpleTimes(repeat(vcat([2, 2, 2], ones(14), [4]), 7)) + opscen = OperationalScenarios(2, [week, week], [0.5, 0.5]) + 𝒯 = TwoLevel(2, 1, opscen; op_per_strat=8760.) + @test_throws AssertionError check_per_dem_sink(; per_dur=25, 𝒯) + @test_throws AssertionError check_per_dem_sink(; per_dur=StrategicProfile([25, 24]), 𝒯) + rep = RepresentativePeriods(2, 8760., [.5, .5], [week, week]) + 𝒯 = TwoLevel(2, 1, rep; op_per_strat=8760.) + @test_throws AssertionError check_per_dem_sink(; per_dur=25, 𝒯) + @test_throws AssertionError check_per_dem_sink(; per_dur=StrategicProfile([25, 24]), 𝒯) - set_optimizer(m, OPTIMIZER) - set_optimizer_attribute(m, MOI.Silent(), true) - optimize!(m) + # Test that a wrong period demand is caught by the checks + @test_throws AssertionError check_per_dem_sink(; per_dem=OperationalProfile([25])) + @test_throws AssertionError check_per_dem_sink(; per_dem=FixedProfile(-10)) - # Test optimal solution - @test termination_status(m) == MOI.OPTIMAL + # Set the global again to false + EMB.TEST_ENV = false end -# Test that the fields of a `PeriodDemandSink` are correctly checked -# - EMB.check_node(n::PeriodDemandSink, 𝒯, modeltype::EnergyModel, check_timeprofiles::Bool) -# @testset "Checks" begin -# function create_check_case(; -# cap = FixedProfile(10), -# period_length = 24, -# period_demand = [fill(1500, 5)..., 0, 0], -# T = TwoLevel(1, 1, SimpleTimes(7 * 24, 1)), -# ) -# demand = PeriodDemandSink( -# "demand_product", -# period_length, -# period_demand, -# cap, -# Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e8)), -# Dict(Power => 1), -# ) - -# return create_system(demand; T) -# end -# # Test that a wrong capacity is caught by the checks -# # this implies that the default checks are working -# @test_throws AssertionError create_check_case(cap=FixedProfile(-25)) - -# # Test that a wrong period length is caught by the checks., including in other time -# # structures -# @test_throws AssertionError create_check_case(period_length=25) -# week = SimpleTimes(168, 1); -# opscen = OperationalScenarios(2, [week, week], [0.5, 0.5]); -# T = TwoLevel(1, 1, opscen; op_per_strat=8760.); -# @test_throws AssertionError create_check_case(;period_length=25, T) -# rep = RepresentativePeriods(2, 8760., [.5, .5], [week, week]); -# T = TwoLevel(1, 1, rep; op_per_strat=8760.); -# @test_throws AssertionError create_check_case(;period_length=25, T) - -# # Test that a wrong period demand is caught by the checks -# @test_throws AssertionError create_check_case(period_demand=[25]) - -# # Test that larger period demands are caught by the checks and print a warning -# msg = -# "The vector `period_demand` is longer than required in " * -# "the operational time structure in strategic period 1. " * -# "The last 23 values will be omitted." -# @test_logs (:warn, msg) create_check_case(period_demand=ones(30)); -# end - -@testset "run-production" begin - test_max_grid_production(false) +@testset "Utility functions" begin + # Create the node and time structure + cap = FixedProfile(10) + per_dur = 24 + per_dem = PartitionProfile([fill(1500, 5)..., 0, 0]) + snk = PeriodDemandSink( + "demand_product", + cap, + per_dur, + per_dem, + Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e4)), + Dict(Power => 0.5), + ) + 𝒯 = TwoLevel(1, 1, SimpleTimes(7 * 24, 1)) + 𝒯ᵖᵈ = EMF.periods(snk, 𝒯) + + @testset "Utility - constructor" begin + # Test that all constructor methods are working + snk_2 = PeriodDemandSink( + "demand_product", + cap, + per_dur, + per_dem, + Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e4)), + Dict(Power => 0.5), + ExtensionData[] + ) + snk_3 = PeriodDemandSink( + "demand_product", + cap, + FixedProfile(per_dur), + per_dem, + Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e4)), + Dict(Power => 0.5), + ) + snk_4 = PeriodDemandSink( + "demand_product", + cap, + PartitionProfile(ones(7)*24), + per_dem, + Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e4)), + Dict(Power => 0.5), + ) + snk_5 = PeriodDemandSink( + "demand_product", + per_dur, + [fill(1500, 5)..., 0, 0], + cap, + Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e4)), + Dict(Power => 0.5), + ) + snk_6 = PeriodDemandSink( + + "demand_product", + per_dur, + [fill(1500, 5)..., 0, 0], + cap, + Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1e4)), + Dict(Power => 0.5), + ExtensionData[] + ) + + for field ∈ fieldnames(PeriodDemandSink) + if field ≠ :period_duration && field ≠ :period_demand + @test getproperty(snk, field) == getproperty(snk_2, field) + @test getproperty(snk, field) == getproperty(snk_3, field) + @test getproperty(snk, field) == getproperty(snk_4, field) + @test getproperty(snk, field) == getproperty(snk_5, field) + @test getproperty(snk, field) == getproperty(snk_6, field) + else + @test all( + getproperty(snk, field)[t_pd] == getproperty(snk_2, field)[t_pd] + for t_pd ∈ 𝒯ᵖᵈ) + @test all( + getproperty(snk, field)[t_pd] == getproperty(snk_3, field)[t_pd] + for t_pd ∈ 𝒯ᵖᵈ) + @test all( + getproperty(snk, field)[t_pd] == getproperty(snk_4, field)[t_pd] + for t_pd ∈ 𝒯ᵖᵈ) + @test all( + getproperty(snk, field)[t_pd] == getproperty(snk_5, field)[t_pd] + for t_pd ∈ 𝒯ᵖᵈ) + @test all( + getproperty(snk, field)[t_pd] == getproperty(snk_6, field)[t_pd] + for t_pd ∈ 𝒯ᵖᵈ) + end + end + end + + @testset "Utility - Identification functions" begin + # Test that all identification functions are working + @test EMB.has_input(snk) + @test !EMB.has_emissions(snk) + @test !EMB.has_output(snk) + end + + @testset "Utility - Extraction functions" begin + # Test that all EMB extraction functions are working + @test capacity(snk) == FixedProfile(10) + @test all(capacity(snk, t) == 10 for t ∈ 𝒯) + @test inputs(snk) == [Power] + @test inputs(snk, Power) == 0.5 + @test surplus_penalty(snk) == FixedProfile(0) + @test all(surplus_penalty(snk, t) == 0 for t ∈ 𝒯) + @test deficit_penalty(snk) == FixedProfile(1e4) + @test all(deficit_penalty(snk, t) == 1e4 for t ∈ 𝒯) + @test node_data(snk) == ExtensionData[] + + # Test that all EMF extraction functions are working + @test EMF.period_demand(snk) == per_dem + @test EMF.periods(snk, 𝒯) == partition_duration(𝒯, per_dur) + @test all( + EMF.period_demand(snk, t_dp) == per_dem[t_dp] for t_dp ∈ EMF.periods(snk, 𝒯) + ) + end + + @testset "Utility - Other functions" begin + # Test that all other functions required for a PeriodDemandSink are working + @test EMF.number_of_periods(snk, 𝒯) == 7 + end end -@testset "test-daily-demand-fulfilled" begin - demand = create_demand_node() - case, model, m = create_system(demand) +@testset "Constraint implementation" begin + + # Create a test set for testing the invariants + function per_sink_tests(m, case; repr=false, oscs=1) + set_optimizer(m, OPTIMIZER) + optimize!(m) - set_optimizer(m, OPTIMIZER) - optimize!(m) + # Test optimal solution + general_tests(m) - # Test optimal solution - general_tests(m) + # Extract the required values from the case and node + 𝒯 = get_time_struct(case) + 𝒯ᴵⁿᵛ = strategic_periods(𝒯) + t_inv = first(𝒯ᴵⁿᵛ) + snk = get_nodes(case)[2] + pers = EMF.periods(snk, 𝒯) - demand = get_nodes(case)[2] - vals = get_values(m, :cap_use, demand, get_time_struct(case)) + # Adjust the variables based on the functions + if repr + main_day = OperationalProfile(vcat(zeros(3), ones(6)*200, zeros(6), [100, 200], [0])) + prod = RepresentativeProfile(vcat( + OperationalProfile(vcat(zeros(3), ones(8)*200, zeros(3), ones(3)*200, [0])), + main_day, + main_day, + main_day, + FixedProfile(0), + FixedProfile(0), + FixedProfile(0), + )) + demand_pd = RepresentativeProfile([2200, fill(1500, 3)..., 0, 0, 0]) + deficit_pd = RepresentativeProfile([0, 0, 0, 0, 1500, 0, 0]) + surplus_pd = RepresentativeProfile([700, 0, 0, 0, 0, 0, 0]) + else + prod = OperationalProfile(vcat( + vcat(zeros(3), ones(8)*200, zeros(3), ones(3)*200, [0]), + repeat(vcat(zeros(3), ones(6)*200, zeros(6), [100, 200], [0]), 3), + zeros(54), + )) + demand_pd = PartitionProfile([2200, fill(1500, 3)..., 0, 0, 0]) + deficit_pd = PartitionProfile([0, 0, 0, 0, 1500, 0, 0]) + surplus_pd = PartitionProfile([700, 0, 0, 0, 0, 0, 0]) + end - period_length = demand.period_length - num_periods = EnergyModelsFlex.number_of_periods(demand, get_time_struct(case)) + # Test the variable generation + @test length(m[:demand_sink_surplus][snk, :]) == 14 * oscs + @test length(m[:demand_sink_deficit][snk, :]) == 14 * oscs - # Test that the demand is fulfilled for each period. - for i ∈ 1:num_periods - period_values = vals[((i-1)*period_length+1):(i*period_length)] - period_total = sum(val for val ∈ period_values) - @test period_total == demand.period_demand[i] + # Tests for the capacity function + # EMB.constraints_capacity(m, n::AbstractPeriodDemandSink, 𝒯::TimeStructure, modeltype::EnergyModel) + + # Test that the individual deficits and surpluses are correctly calculated + @test all( + value.(m[:sink_deficit][snk, t]) + value.(m[:cap_use][snk, t]) ≈ + value.(m[:cap_inst][snk, t]) for t ∈ 𝒯, + atol = TEST_ATOL + ) + # Test that the surplus is fixed to 0 + @test all(is_fixed.(m[:sink_surplus][snk, t]) for t ∈ 𝒯) + @test all(value.(m[:sink_surplus][snk, t]) ≈ 0 for t ∈ 𝒯) + + # Test that the production is as planned based on the cost with no production in period + # 5 due to the prohibitive costs + cap = capacity(snk) + deficit = cap - prod + @test all(value.(m[:cap_use][snk, t]) ≈ prod[t] for t ∈ 𝒯) + @test all(value.(m[:sink_deficit][snk, t]) ≈ deficit[t] for t ∈ 𝒯) + + # Test that the demand is fulfilled for the first 4 periods and the last 2 + @test all(sum(value.(m[:cap_use][snk, t]) for t ∈ t_pd) ≈ demand_pd[t_pd] for t_pd ∈ pers) + @test all(value.(m[:demand_sink_deficit][snk, t_pd]) ≈ deficit_pd[t_pd] for t_pd ∈ pers) + @test all(value.(m[:demand_sink_surplus][snk, t_pd]) ≈ surplus_pd[t_pd] for t_pd ∈ pers) + + # Test the upper bound on the installed capacity and the value for the capacity + @test all(value.(m[:cap_use][snk, t]) ≲ value.(m[:cap_inst][snk, t]) for t ∈ 𝒯) + @test all(is_fixed.(m[:cap_inst][snk, t]) for t ∈ 𝒯) + @test all(value.(m[:cap_inst][snk, t]) ≈ capacity(snk, t) for t ∈ 𝒯) + + # Test that the fixed OPEX is set to 0 + # - EMB.constraints_opex_fixed(m, n::Sink, 𝒯ᴵⁿᵛ, modeltype::EnergyModel) + @test all(is_fixed.(m[:opex_fixed][snk, t_inv]) for t_inv ∈ 𝒯ᴵⁿᵛ) + @test all(value.(m[:opex_fixed][snk, t_inv]) ≈ 0 for t_inv ∈ 𝒯ᴵⁿᵛ) + + # Test that the variable OPEX is correctly calculated + # - EMB.constraints_opex_fixed(m, n::AbstractPeriodDemandSink, 𝒯ᴵⁿᵛ, modeltype::EnergyModel) + @test all( + value.(m[:opex_var][snk, t_inv]) ≈ + sum( + (value.(m[:demand_sink_deficit][snk, t_pd]) * deficit_penalty(snk, t_pd) + + value.(m[:demand_sink_surplus][snk, t_pd]) * surplus_penalty(snk, t_pd)) * + multiple_strat(t_inv, first(t_pd)) * probability(first(t_pd)) + for t_pd ∈ EMF.periods(snk, t_inv)) + for t_inv ∈ 𝒯ᴵⁿᵛ) + @test all( + value.(m[:opex_var][snk, t_inv]) ≈ 8760/24/7 * (1e4*1500 - 8 * 700) + for t_inv ∈ 𝒯ᴵⁿᵛ) + + return objective_value(m) end + + # Create and optimize the model + ops = vcat([2, 2, 2], ones(14), [4]) + 𝒯 = TwoLevel(2, 1, SimpleTimes(repeat(ops, 7)), op_per_strat=8760.) + m, case, modeltype = per_dem_snk_case(; 𝒯) + obj_1 = per_sink_tests(m, case) + + 𝒯 = TwoLevel(2, 1, OperationalScenarios(2, SimpleTimes(repeat(ops, 7))), op_per_strat=8760.) + m, case, modeltype = per_dem_snk_case(; 𝒯) + obj_2 = per_sink_tests(m, case; oscs=2) + + 𝒯 = TwoLevel(2, 1, RepresentativePeriods(7, 8760, SimpleTimes(ops)), op_per_strat=8760.) + m, case, modeltype = per_dem_snk_case(; 𝒯, repr=true) + obj_3 = per_sink_tests(m, case; repr=true) + + @test obj_1 ≈ obj_2 + @test obj_1 ≈ obj_3 end diff --git a/test/utils.jl b/test/utils.jl index e0540fc..b72129f 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,18 +1,16 @@ -const ATOL = 1e-7 - """ a ≲ b -Approximate ≤ comparison with absolute tolerance `ATOL`. +Approximate ≤ comparison with absolute tolerance `TEST_ATOL`. """ -≲(a::Real, b::Real) = a <= b + ATOL +≲(a::Real, b::Real) = a ≤ b + TEST_ATOL """ a ≳ b -Approximate ≥ comparison with absolute tolerance `ATOL`. +Approximate ≥ comparison with absolute tolerance `TEST_ATOL`. """ -≳(a::Real, b::Real) = a + ATOL >= b +≳(a::Real, b::Real) = a + TEST_ATOL ≥ b function get_values(m, variable, node, iterable) return [JuMP.value(m[variable][node, t]) for t ∈ iterable]