diff --git a/NEWS.md b/NEWS.md index 723b23d..75617fc 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,21 @@ # Release Notes +## Unreleased + +### New thermal energy storage node + +* Introduced the `AbstractTES` node `LevelDependentRateTES`, the most advanced thermal energy + storage node, with state-of-charge dependent charge and discharge rate limits described by + one to three c-rate anchor points per direction (piecewise-linear with binary region + selection), including documentation, checks, and tests. + The node is restricted to `OperationalModel`. + +### Other changes + +* Added the function `visualize_c_rates` for plotting the theoretical c-rate curves, provided + through a `Plots` package extension (`ext/PlotsExt.jl`) so the core package keeps no hard + plotting dependency. + ## Version 0.2.0 (2026-04-14) ### Breaking changes diff --git a/Project.toml b/Project.toml index 85b17dd..db0e822 100644 --- a/Project.toml +++ b/Project.toml @@ -8,8 +8,15 @@ EnergyModelsBase = "5d7e687e-f956-46f3-9045-6f5a5fd49f50" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" TimeStruct = "f9ed5ce0-9f41-4eaa-96da-f38ab8df101c" +[weakdeps] +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" + +[extensions] +PlotsExt = "Plots" + [compat] EnergyModelsBase = "0.10.0" JuMP = "1.23.5" +Plots = "1" TimeStruct = "0.9.1" julia = "1.10" diff --git a/docs/make.jl b/docs/make.jl index 9be779e..f44e035 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -49,6 +49,7 @@ Documenter.makedocs( "Nodes" => [ "HeatPump" => "nodes/heatpump.md", "ThermalEnergyStorage" => "nodes/thermalenergystorage.md", + "LevelDependentRateTES" => "nodes/leveldependentratetes.md", "HeatExchanger" => "nodes/heatexchanger.md", ], "How to" => [ diff --git a/docs/src/library/internals/methods-EMB.md b/docs/src/library/internals/methods-EMB.md index 49cd3fa..e8b7297 100644 --- a/docs/src/library/internals/methods-EMB.md +++ b/docs/src/library/internals/methods-EMB.md @@ -11,6 +11,7 @@ Pages = ["methods-EMB.md"] ```@docs EMB.create_link EMB.variables_link +EMB.variables_node ``` ## [Constraint methods](@id lib-int-met_emb-con) diff --git a/docs/src/library/internals/methods-EMH.md b/docs/src/library/internals/methods-EMH.md index 4453c8e..b440098 100644 --- a/docs/src/library/internals/methods-EMH.md +++ b/docs/src/library/internals/methods-EMH.md @@ -13,6 +13,8 @@ Pages = ["methods-EMH.md"] ## [Utility methods](@id lib-int-met-util) ```@docs +EMH.c_rate_points_charge +EMH.c_rate_points_discharge EMH.dh_fraction EMH.dh_upgrade EMH.fraction_equal_mass diff --git a/docs/src/library/public.md b/docs/src/library/public.md index d86fac1..8181fb1 100644 --- a/docs/src/library/public.md +++ b/docs/src/library/public.md @@ -25,10 +25,17 @@ DHPipe HeatPump ThermalEnergyStorage BoundRateTES +LevelDependentRateTES HeatExchanger EnergyModelsHeat.DirectHeatUpgrade ``` +## [Functions](@id lib-pub-functions) + +```@docs +EnergyModelsHeat.visualize_c_rates +``` + !!! todo "Changes" Once exported, remove `EnergyModelsHeat.` in front of `DirectHeatUpgrade` diff --git a/docs/src/nodes/crate_curves_example.png b/docs/src/nodes/crate_curves_example.png new file mode 100644 index 0000000..f6dd2f0 Binary files /dev/null and b/docs/src/nodes/crate_curves_example.png differ diff --git a/docs/src/nodes/leveldependentratetes.md b/docs/src/nodes/leveldependentratetes.md new file mode 100644 index 0000000..d5d9c00 --- /dev/null +++ b/docs/src/nodes/leveldependentratetes.md @@ -0,0 +1,169 @@ +# [LevelDependentRateTES nodes](@id nodes-LDRTES) + +The [`LevelDependentRateTES`](@ref) is the most advanced thermal energy storage node in +`EnergyModelsHeat`. +Like [`ThermalEnergyStorage`](@ref) and [`BoundRateTES`](@ref) it behaves mostly like a +[`RefStorage`](@extref EnergyModelsBase.RefStorage) with thermal energy losses that are +proportional to the storage level (see *[ThermalEnergyStorage nodes](@ref nodes-TES)* for the +description of the heat loss). +Its distinguishing feature is that the maximum charge and discharge rates depend on the +current *state of charge* of the storage. + +In a real thermal energy storage the rate at which heat can be added or withdrawn is not +constant: charging typically slows down as the storage fills up, and discharging slows down as +it empties. +The [`LevelDependentRateTES`](@ref) captures this through user-provided **c-rate curves**: a +small number of `[level, rate]` anchor points that describe the maximum rate as a function of +the storage level at the end of the previous operational period. + +!!! danger "Only operational models are supported" + A [`LevelDependentRateTES`](@ref) can only be used together with an + [`OperationalModel`](@extref EnergyModelsBase.OperationalModel). + In an [`InvestmentModel`](@extref EnergyModelsBase.InvestmentModel) the installed storage, + charge and discharge capacities become decision variables. + The slopes and intercepts of the c-rate curves then depend on those capacity variables, so + the rate limit would multiply a capacity variable by the storage-level variable. + This product of two decision variables is a *bilinear* (nonconvex quadratic) term, which + the linear and mixed-integer-linear solvers used in this package (*e.g.*, `HiGHS`) cannot + handle. + +!!! danger "StorageBehavior for LevelDependentRateTES" + Like the other thermal energy storage nodes, a [`LevelDependentRateTES`](@ref) can only + use [`Cyclic`](@extref EnergyModelsBase.Cyclic) storage behaviors, and when using + `RepresentativePeriods` this is reduced to + [`CyclicRepresentative`](@extref EnergyModelsBase.CyclicRepresentative). + +## [Introduced type and its fields](@id nodes-LDRTES-fields) + +[`LevelDependentRateTES`](@ref) uses the same standard fields as the other thermal energy +storage nodes (see *[ThermalEnergyStorage nodes](@ref nodes-TES-fields-stand)*): `id`, +`charge`, `level`, `discharge`, `stor_res`, `input`, `output` and `data`. +In addition it introduces the heat loss factor and the c-rate anchor points. + +### [Additional fields](@id nodes-LDRTES-fields-new) + +- **`heat_loss_factor::Float64`** :\ + The heat lost relative to the storage level of the previous operational period, given as a + fraction in the range ``[0, 1]``. + It corresponds to the loss occurring between two operational periods for a given operational + duration of 1. + +- **`c_rate_points_charge::Vector{<:Vector{<:Real}}`** and + **`c_rate_points_discharge::Vector{<:Vector{<:Real}}`** :\ + One to three `[level, rate]` anchor points describing the maximum (dis-)charge rate as a + function of the storage level. + For **charging**, the curve starts at the point ``(0, capacity(charge(n)))`` (the full + installed charge rate at an empty storage) and runs through the provided points in *strictly + ascending* level order. + For **discharging**, the curve starts at the point ``(capacity(level(n)), capacity(discharge(n)))`` + (the full installed discharge rate at a full storage) and runs through the provided points in + *strictly descending* level order. + All rates must be non-negative, the charge levels must lie in ``(0, capacity(level(n))]`` and + the discharge levels in ``[0, capacity(level(n)))``. + The number of anchor points sets the number of linear segments: + - one point gives a single straight line, + - two or three points give a piecewise-linear curve, which introduces binary variables to + select the active segment. + +As an example, the storage in `examples/advanced_example.jl` has an installed charge rate of +30, discharge rate of 40 and level capacity of 100, combined with the three-point curves + +```julia +c_rate_points_charge = [[10.0, 30.0], [50.0, 20.0], [100.0, 10.0]] +c_rate_points_discharge = [[75.0, 40.0], [50.0, 20.0], [25.0, 15.0]] +``` + +The charge curve therefore stays at the full rate of 30 up to a level of 10 and then falls to +20 at level 50 and 10 at the full level of 100. +The discharge curve reaches its full rate of 40 for levels at or above 75 and falls to 20 at +level 50 and 15 at level 25 as the storage empties. +The resulting theoretical curves are shown below. +Near a full storage the charge rate is additionally limited by the remaining free capacity, and +near an empty storage the discharge rate is limited by the available energy, which is why both +curves bend towards zero at the respective ends. + +![Theoretical charge and discharge c-rate curves for the example anchor points.](crate_curves_example.png) + +!!! tip "Visualising the curves" + The plot above was produced with the function + [`visualize_c_rates`](@ref EnergyModelsHeat.visualize_c_rates), which plots the theoretical + (dis-)charge curves a given set of anchor points produces, so you can check your input + before solving a model. + It is provided through a package extension and becomes available once you run + `using Plots`. + +## [Mathematical description](@id nodes-LDRTES-math) + +In the following description we use the variable and function names from the model, as in the +*[ThermalEnergyStorage](@ref nodes-TES-math)* page. + +### [Variables](@id nodes-LDRTES-math-var) + +[`LevelDependentRateTES`](@ref) uses all standard [`RefStorage`](@extref EnergyModelsBase.RefStorage) +variables (see *[ThermalEnergyStorage variables](@ref nodes-TES-math-var)*). +In addition, when a c-rate curve uses more than one anchor point, the following auxiliary +binary variables are declared (in [`EMB.variables_node`](@ref EnergyModelsBase.variables_node)): + +- ``\texttt{bin\_region\_charge}[n, t, 1:2]`` selects the active segment of the charge curve. +- ``\texttt{bin\_region\_discharge}[n, t, 1:2]`` selects the active segment of the discharge curve. + +Two binaries per direction are sufficient to encode up to three segments. + +### [Constraints](@id nodes-LDRTES-math-con) + +[`LevelDependentRateTES`](@ref) uses the standard `Storage` constraints, with the heat-loss +adjustment to `constraints_level_iterate` shared by all thermal energy storage nodes (see +*[Level constraints](@ref nodes-TES-math-con-level)*). +It overrides `constraints_capacity` to add the state-of-charge dependent rate limits. + +The standard capacity bounds are kept: + +```math +\begin{aligned} +\texttt{stor\_level}[n, t] & ≤ \texttt{stor\_level\_inst}[n, t] \\ +\texttt{stor\_charge\_use}[n, t] & ≤ \texttt{stor\_charge\_inst}[n, t] \\ +\texttt{stor\_discharge\_use}[n, t] & ≤ \texttt{stor\_discharge\_inst}[n, t] +\end{aligned} +``` + +On top of these, the rate is limited by the c-rate curve evaluated at the storage level of the +*previous* operational period ``t_{prev}``. + +**Single anchor point.** +With a single charge anchor point ``(x_1, y_1)``, the charge limit is the straight line from +``(0, capacity(charge(n)))`` to ``(x_1, y_1)``: + +```math +\texttt{stor\_charge\_use}[n, t] ≤ +\frac{y_1 - capacity(charge(n), t)}{x_1} \times \texttt{stor\_level}[n, t_{prev}] + +capacity(charge(n), t) +``` + +Analogously, the single-segment discharge limit is the straight line from +``(capacity(level(n)), capacity(discharge(n)))`` to the discharge anchor point ``(x_1, y_1)``: + +```math +\texttt{stor\_discharge\_use}[n, t] ≤ +\frac{capacity(discharge(n), t) - y_1}{capacity(level(n), t_{prev}) - x_1} +\times \left(\texttt{stor\_level}[n, t_{prev}] - capacity(level(n), t_{prev})\right) + +capacity(discharge(n), t) +``` + +**Two or three anchor points.** +With more anchor points the curve becomes piecewise linear. +For each segment a linear limit of the same form as above is built, and the binary variables +``\texttt{bin\_region\_charge}`` / ``\texttt{bin\_region\_discharge}`` together with a +big-``M`` formulation determine which segment applies in period ``t``, based on where +``\texttt{stor\_level}[n, t_{prev}]`` falls relative to the anchor-point levels. +Only the limit of the active segment is binding; the inactive segments and the region bounds +are relaxed by big-``M`` terms. +Because a [`LevelDependentRateTES`](@ref) is restricted to an +[`OperationalModel`](@extref EnergyModelsBase.OperationalModel), every slope, intercept and +capacity is a fixed parameter, so each big-``M`` is computed from the geometry to be just large +enough to make its constraint non-binding: the region bounds are relaxed up to the installed +level capacity, and each inactive rate segment is lifted just above the installed (dis-)charge +capacity over the storage-level range on which it is inactive. +This avoids the artificial level cap that a single, too-small constant would impose, while +keeping the model's relaxation tight. +The result is that the realised (dis-)charge rate never exceeds the piecewise-linear c-rate +curve at the current state of charge. diff --git a/examples/Project.toml b/examples/Project.toml index c080c44..6ad738a 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -1,7 +1,9 @@ [deps] +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" EnergyModelsBase = "5d7e687e-f956-46f3-9045-6f5a5fd49f50" EnergyModelsHeat = "ad1b8b27-e232-4da9-b498-bea9c19a30d7" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" TimeStruct = "f9ed5ce0-9f41-4eaa-96da-f38ab8df101c" diff --git a/examples/advanced_example.jl b/examples/advanced_example.jl new file mode 100755 index 0000000..b89ff18 --- /dev/null +++ b/examples/advanced_example.jl @@ -0,0 +1,204 @@ +using Pkg +# Activate the local environment including EnergyModelsHeat and its dependencies +Pkg.activate(@__DIR__) +# Use dev version if run as part of tests +haskey(ENV, "EMX_TEST") && Pkg.develop(path = joinpath(@__DIR__, "..")) +# Install the dependencies. +Pkg.instantiate() + +# Import the required packages +using EnergyModelsBase +using EnergyModelsHeat +using JuMP +using TimeStruct +using HiGHS +using DataFrames +# Loading `Plots` activates the `visualize_c_rates` plotting extension of EnergyModelsHeat +using Plots + +const EMB = EnergyModelsBase +const EMH = EnergyModelsHeat +const TS = TimeStruct + +""" + generate_level_dependent_tes_example() + +Generate the data for an example consisting of an electricity source, an electric boiler, a +[`LevelDependentRateTES`](@ref) and a heat demand. The boiler converts electricity to heat, +which can either be delivered directly to the demand or stored in the thermal energy storage. + +The storage uses state-of-charge dependent charge and discharge rate limits: charging gets +slower as the storage fills up, and discharging gets slower as it empties. This example +demonstrates the flexibility this provides when the electricity price and the heat demand +vary over the day. +""" +function generate_level_dependent_tes_example() + @info "Generate case data - LevelDependentRateTES example" + + # Define the different resources and their emission intensity in t CO₂/MWh + Power = ResourceCarrier("Power", 0.0) + Heat = ResourceCarrier("Heat", 0.0) + CO2 = ResourceEmit("CO2", 1.0) + products = [Power, Heat, CO2] + + # Creation of the time structure with 24 hourly operational periods in a single year + op_number = 24 + operational_periods = SimpleTimes(op_number, 1) + T = TwoLevel(1, 1, operational_periods; op_per_strat = op_number) + + # Creation of the model type with global data + model = OperationalModel( + Dict(CO2 => FixedProfile(1e6)), # Emission cap for CO₂ in t/a + Dict(CO2 => FixedProfile(0)), # Emission price for CO₂ in €/t + CO2, # CO₂ instance + ) + + # A deterministic daily heat demand (MW) with a morning and an evening peak + heat_demand = OperationalProfile([ + 10, 8, 6, 6, 8, 15, 30, 40, 35, 25, 20, 18, + 18, 20, 22, 25, 30, 40, 45, 40, 30, 22, 16, 12, + ]) + # A deterministic electricity price (€/MWh) that is cheap at night and expensive at peaks + el_price = OperationalProfile([ + 5, 4, 4, 4, 5, 8, 15, 25, 20, 12, 10, 10, + 10, 10, 12, 15, 20, 30, 35, 28, 18, 12, 8, 6, + ]) + + ### Thermal energy storage parameters ### + charge_capacity = 30 + level_capacity = 100 + discharge_capacity = 40 + # Anchor points [storage level, maximum rate]: charging slows down as the storage fills, + # in ascending level order + c_rate_points_charge = [[10.0, 30.0], [50.0, 20.0], [100.0, 10.0]] + # Discharging slows down as the storage empties, in descending level order + c_rate_points_discharge = [[75.0, 40.0], [50.0, 20.0], [25.0, 15.0]] + + # Create the individual nodes: an electricity source (1), an electric boiler (2), + # a heat demand (3) and the thermal energy storage (4). + nodes = [ + RefSource( + "electricity source", # Node id + FixedProfile(50), # Installed capacity in MW + el_price, # Variable OPEX in €/MWh + FixedProfile(0), # Fixed OPEX in €/MW/a + Dict(Power => 1), # Output from the node, in this case Power + ), + RefNetworkNode( + "electric boiler", + FixedProfile(50), # Installed capacity in MW + FixedProfile(0), # Variable OPEX in €/MWh + FixedProfile(0), # Fixed OPEX in €/MW/a + Dict(Power => 1), # Input resource and corresponding input ratio + Dict(Heat => 1), # Output resource and corresponding output ratio + ), + RefSink( + "heat demand", # Node id + heat_demand, # Required demand in MW + Dict(:surplus => FixedProfile(0), :deficit => FixedProfile(1000)), + Dict(Heat => 1), # Energy carrier and corresponding ratio to demand + ), + LevelDependentRateTES{CyclicRepresentative}( + "thermal energy storage", + StorCap(FixedProfile(charge_capacity)), # Installed charge rate in MW + StorCap(FixedProfile(level_capacity)), # Installed storage capacity in MWh + StorCap(FixedProfile(discharge_capacity)), # Installed discharge rate in MW + Heat, # Stored resource + 0.01, # Heat loss factor (fraction per period) + c_rate_points_charge, + c_rate_points_discharge, + Dict(Heat => 1), # Input resource and input ratio + Dict(Heat => 1), # Output resource and output ratio + ), + ] + + # Connect all nodes for the overall energy/mass balance + links = [ + Direct("source-boiler", nodes[1], nodes[2], Linear()), + Direct("boiler-demand", nodes[2], nodes[3], Linear()), + Direct("boiler-storage", nodes[2], nodes[4], Linear()), + Direct("storage-demand", nodes[4], nodes[3], Linear()), + ] + + ### c-rate curves visualization ### + # `visualize_c_rates` becomes available through the Plots extension of EnergyModelsHeat. + # It shows the theoretical (dis-)charge curves the chosen anchor points produce. + visualize_c_rates( + charge_capacity, + discharge_capacity, + level_capacity, + c_rate_points_charge, + c_rate_points_discharge, + ) + + # Input data structure + case = Case(T, products, [nodes, links], [[get_nodes, get_links]]) + return case, model +end + +""" + process_results(m, case) + +Collect the storage results into a `DataFrame` and plot the realised charge and discharge +rates against the state of charge, illustrating the state-of-charge dependent rate limits. +""" +function process_results(m, case) + # Extract the nodes and resources from the data + _, _, _, tes = get_nodes(case) + Heat = get_products(case)[2] + 𝒯 = get_time_struct(case) + periods = collect(𝒯) + + # The rate limit in period t depends on the level at the end of the previous period + level = [value(m[:stor_level][tes, t]) for t ∈ periods] + prev_level = vcat(0.0, level[1:(end-1)]) + charge = [value(m[:flow_in][tes, t, Heat]) for t ∈ periods] + discharge = [value(m[:flow_out][tes, t, Heat]) for t ∈ periods] + + results = DataFrame( + period = 1:length(periods), + prev_level = prev_level, + charge = charge, + discharge = discharge, + ) + + p = plot(layout = (1, 2), size = (1200, 500), legend = :topright) + scatter!( + p[1], + results.prev_level, + results.charge, + label = "Charge", + markercolor = :red, + xlabel = "State of charge", + ylabel = "Rate", + title = "Charge rate vs state of charge", + ) + scatter!( + p[2], + results.prev_level, + results.discharge, + label = "Discharge", + markercolor = :blue, + xlabel = "State of charge", + ylabel = "Rate", + title = "Discharge rate vs state of charge", + ) + return results, p +end + +# Generate the case and model data and run the model +case, model = generate_level_dependent_tes_example() +m = EMB.create_model(case, model) + +@info "Solving the model..." +optimizer = optimizer_with_attributes( + HiGHS.Optimizer, + MOI.Silent() => true, + "mip_rel_gap" => 0.01, +) +set_optimizer(m, optimizer) +optimize!(m) + +results, combined_plot = process_results(m, case) +# Skip the interactive display when the example is run as part of the automated tests +haskey(ENV, "EMX_TEST") || display(combined_plot) diff --git a/ext/PlotsExt.jl b/ext/PlotsExt.jl new file mode 100644 index 0000000..6d8cf8e --- /dev/null +++ b/ext/PlotsExt.jl @@ -0,0 +1,191 @@ +""" +Package extension of `EnergyModelsHeat` providing the implementation of +[`visualize_c_rates`](@ref EnergyModelsHeat.visualize_c_rates). + +The extension is only loaded once the `Plots` package is available (declared as a weak +dependency of `EnergyModelsHeat`). This keeps plotting out of the core package, consistent +with the convention in the `EnergyModelsX` ecosystem of not adding plotting dependencies to +the technology packages. +""" +module PlotsExt + +using EnergyModelsHeat +using Plots + +# ========================================================= +# Helper functions for computing charge/discharge curves +# ========================================================= + +""" + compute_charge_curve(levels, cap_charge, c_rate_points) + +Return the maximum charge rate at each storage level in `levels`. + +- `levels`: array of storage levels to evaluate. +- `cap_charge`: installed charge rate (curve value at an empty storage). +- `c_rate_points`: one to three `[level, rate]` pairs in ascending storage-level order. + +The formulas correspond to the piecewise-linear constraints built for a +[`LevelDependentRateTES`](@ref EnergyModelsHeat.LevelDependentRateTES). +""" +function compute_charge_curve(levels, cap_charge, c_rate_points) + if length(c_rate_points) == 1 + x1, y1 = c_rate_points[1] + + m = (y1 - cap_charge) / x1 + b = cap_charge + + return m .* levels .+ b + + elseif length(c_rate_points) == 2 + (x1, y1), (x2, y2) = c_rate_points + + # First segment + m1 = (y1 - cap_charge) / x1 + b1 = cap_charge + + # Second segment + m2 = (y2 - y1) / (x2 - x1) + b2 = y2 - m2 * x2 + + return [level <= x1 ? m1 * level + b1 : m2 * level + b2 for level ∈ levels] + + elseif length(c_rate_points) == 3 + (x1, y1), (x2, y2), (x3, y3) = c_rate_points + + m1 = (y1 - cap_charge) / x1 + b1 = cap_charge + + m2 = (y2 - y1) / (x2 - x1) + b2 = y2 - m2 * x2 + + m3 = (y3 - y2) / (x3 - x2) + b3 = y3 - m3 * x3 + + return [ + level <= x1 ? m1 * level + b1 : + (level <= x2 ? m2 * level + b2 : m3 * level + b3) for level ∈ levels + ] + else + error("c_rate_points_charge must contain between 1 and 3 coordinate pairs.") + end +end + +""" + compute_discharge_curve(levels, cap_discharge, cap_level, c_rate_points) + +Return the maximum discharge rate at each storage level in `levels`. + +- `levels`: array of storage levels to evaluate. +- `cap_discharge`: installed discharge rate (curve value at a full storage). +- `cap_level`: installed storage level capacity. +- `c_rate_points`: one to three `[level, rate]` pairs in descending storage-level order. + +The formulas correspond to the piecewise-linear constraints built for a +[`LevelDependentRateTES`](@ref EnergyModelsHeat.LevelDependentRateTES). +""" +function compute_discharge_curve(levels, cap_discharge, cap_level, c_rate_points) + if length(c_rate_points) == 1 + x1, y1 = c_rate_points[1] + + m = (cap_discharge - y1) / (cap_level - x1) + b = cap_discharge - m * cap_level + + return m .* levels .+ b + + elseif length(c_rate_points) == 2 + (x1, y1), (x2, y2) = c_rate_points + + # First segment (region 1) + m1 = (cap_discharge - y1) / (cap_level - x1) + b1 = cap_discharge - m1 * cap_level + + # Second segment (region 2) + m2 = (y1 - y2) / (x1 - x2) + b2 = y2 - m2 * x2 + + return [level >= x1 ? m1 * level + b1 : m2 * level + b2 for level ∈ levels] + + elseif length(c_rate_points) == 3 + (x1, y1), (x2, y2), (x3, y3) = c_rate_points + + # Region 1: highest SOC (cap_level down to x1) + m1 = (cap_discharge - y1) / (cap_level - x1) + b1 = cap_discharge - m1 * cap_level + + # Region 2: between x1 and x2 + m2 = (y1 - y2) / (x1 - x2) + b2 = y2 - m2 * x2 + + # Region 3: between x2 and x3 (lowest SOC) + m3 = (y2 - y3) / (x2 - x3) + b3 = y3 - m3 * x3 + + return [ + level >= x1 ? m1 * level + b1 : + (level >= x2 ? m2 * level + b2 : m3 * level + b3) for level ∈ levels + ] + else + error("c_rate_points_discharge must contain between 1 and 3 coordinate pairs.") + end +end + +function EnergyModelsHeat.visualize_c_rates( + charge_capacity, + discharge_capacity, + level_capacity, + c_rate_points_charge, + c_rate_points_discharge, +) + # Storage level range to evaluate + levels = range(0, level_capacity, length = 201) + + # ========================================================= + # Compute curves + # ========================================================= + + charge_curve = compute_charge_curve(levels, charge_capacity, c_rate_points_charge) + discharge_curve = compute_discharge_curve( + levels, + discharge_capacity, + level_capacity, + c_rate_points_discharge, + ) + maximum_charge_curve = [-level + level_capacity for level ∈ levels] + maximum_discharge_curve = [level for level ∈ levels] + + combined_charge = Vector{Float64}(undef, length(charge_curve)) + for i ∈ eachindex(charge_curve, maximum_charge_curve) + combined_charge[i] = min(charge_curve[i], maximum_charge_curve[i]) + end + + combined_discharge = Vector{Float64}(undef, length(discharge_curve)) + for i ∈ eachindex(discharge_curve, maximum_discharge_curve) + combined_discharge[i] = min(discharge_curve[i], maximum_discharge_curve[i]) + end + + # ========================================================= + # Plot + # ========================================================= + + p = plot( + levels, + combined_charge, + label = "Max Charge Rate", + xlabel = "Storage Level", + ylabel = "Charge/Discharge Rate", + linewidth = 2, + ylims = (0, max(charge_capacity, discharge_capacity)), + size = (800, 500), + ) + + plot!(levels, combined_discharge, label = "Max Discharge Rate", linewidth = 2) + + title!("C-Rate Curves (Theoretical)") + + display(p) + @info "Theoretical c-rate curves plotted" + return p +end + +end diff --git a/src/EnergyModelsHeat.jl b/src/EnergyModelsHeat.jl index 09d5236..96a60e0 100644 --- a/src/EnergyModelsHeat.jl +++ b/src/EnergyModelsHeat.jl @@ -2,7 +2,7 @@ Main module for `EnergyModelsHeat`: a framework for energy system models with thermal components. -It exports the resource `ResourceHeat` and structures for DH pipe, heat pump and heat exchanger. +It exports the resource `ResourceHeat` and structures for DH pipe, heat pump, thermal energy storage and heat exchanger. """ module EnergyModelsHeat @@ -22,6 +22,7 @@ include("model.jl") include("constraint_functions.jl") include("utils.jl") include("legacy_constructor.jl") +include("visualization.jl") # Custom input validation include("checks.jl") @@ -34,6 +35,7 @@ export DHPipe export HeatPump export PinchData export HeatExchanger -export AbstractTES, ThermalEnergyStorage, BoundRateTES +export AbstractTES, ThermalEnergyStorage, BoundRateTES, LevelDependentRateTES +export visualize_c_rates end diff --git a/src/checks.jl b/src/checks.jl index 3c74cc4..6a37c7f 100644 --- a/src/checks.jl +++ b/src/checks.jl @@ -196,6 +196,107 @@ function EMB.check_node( end end +""" + EMB.check_node(n::LevelDependentRateTES{T}, 𝒯, modeltype::EnergyModel, check_timeprofiles::Bool) where {T<:EMB.StorageBehavior} + +This method checks that the *[`LevelDependentRateTES`](@ref)* node is valid. + +It reuses the standard checks of a `Storage` node through calling the function +[`EMB.check_node_default`](@extref EnergyModelsBase.check_node_default), and adds checks on the +heat loss factor as well as on the c-rate anchor points that shape the state-of-charge +dependent charge and discharge limits. + +## Checks +- The `TimeProfile` of the field `capacity` in the type in the field `charge` is required + to be non-negative if the chosen composite type has the field `capacity`. +- The `TimeProfile` of the field `capacity` in the type in the field `level` is required + to be non-negative. +- The `TimeProfile` of the field `capacity` in the type in the field `discharge` is required + to be non-negative if the chosen composite type has the field `capacity`. +- The `TimeProfile` of the field `fixed_opex` is required to be non-negative and + accessible through a `StrategicPeriod` as outlined in the function + `check_fixed_opex(n, 𝒯ᴵⁿᵛ, check_timeprofiles)` for the chosen composite type. +- The values of the dictionary `input` are required to be non-negative. +- The values of the dictionary `output` are required to be non-negative. +- The value of the field `heat_loss_factor` is required to be in the range ``[0, 1]``. +- The fields `c_rate_points_charge` and `c_rate_points_discharge` are each required to hold + between 1 and 3 points, where each point is a 2-element `[level, rate]` pair with a + non-negative rate. +- The `level` values of `c_rate_points_charge` must be strictly ascending and lie in the + range ``(0, capacity(level(n), t)]`` for all ``t ∈ \\mathcal{T}``. +- The `level` values of `c_rate_points_discharge` must be strictly descending and lie in the + range ``[0, capacity(level(n), t))`` for all ``t ∈ \\mathcal{T}``. + +## Warnings +- The `StorageBehavior` should not be `CyclicStrategic` when using `RepresentativePeriods`. +""" +function EMB.check_node( + n::LevelDependentRateTES{T}, + 𝒯, + modeltype::EnergyModel, + check_timeprofiles::Bool, +) where {T<:EMB.StorageBehavior} + EMB.check_node_default(n, 𝒯, modeltype, check_timeprofiles) + + @assert_or_log( + heat_loss_factor(n) ≥ 0, + "The heat_loss_factor field must be non-negative." + ) + @assert_or_log( + heat_loss_factor(n) ≤ 1, + "The heat_loss_factor field must be less or equal to 1." + ) + + # The smallest installed level capacity bounds the admissible anchor-point levels + level_cap = minimum(capacity(level(n), t) for t ∈ 𝒯) + + # The expected format is repeated in every message so a user that hits one error + # immediately sees the structure all anchor points must follow. + format_charge = + "The field c_rate_points_charge must contain 1 to 3 points, each a 2-element " * + "[level, rate] pair with the levels in strictly ascending order and within " * + "(0, installed level capacity], e.g. [[10.0, 30.0], [50.0, 20.0], [100.0, 10.0]]." + format_discharge = + "The field c_rate_points_discharge must contain 1 to 3 points, each a 2-element " * + "[level, rate] pair with the levels in strictly descending order and within " * + "[0, installed level capacity), e.g. [[75.0, 40.0], [50.0, 20.0], [25.0, 15.0]]." + + points_charge = c_rate_points_charge(n) + points_discharge = c_rate_points_discharge(n) + + # --- Charge anchor points --- + @assert_or_log(1 ≤ length(points_charge) ≤ 3, format_charge) + @assert_or_log(all(length(p) == 2 for p ∈ points_charge), format_charge) + @assert_or_log(all(p[2] ≥ 0 for p ∈ points_charge), format_charge) + @assert_or_log(all(0 < p[1] ≤ level_cap for p ∈ points_charge), format_charge) + @assert_or_log( + issorted(points_charge; by = p -> p[1], lt = <) && + allunique(p[1] for p ∈ points_charge), + format_charge + ) + + # --- Discharge anchor points --- + @assert_or_log(1 ≤ length(points_discharge) ≤ 3, format_discharge) + @assert_or_log(all(length(p) == 2 for p ∈ points_discharge), format_discharge) + @assert_or_log(all(p[2] ≥ 0 for p ∈ points_discharge), format_discharge) + @assert_or_log(all(0 ≤ p[1] < level_cap for p ∈ points_discharge), format_discharge) + @assert_or_log( + issorted(points_discharge; by = p -> p[1], rev = true, lt = <) && + allunique(p[1] for p ∈ points_discharge), + format_discharge + ) + + if (T <: CyclicStrategic) && + isa(𝒯, TwoLevel{S,T,U} where {S,T,U<:RepresentativePeriods}) + @warn( + "Using `CyclicStrategic` with a `LevelDependentRateTES` and `RepresentativePeriods` " * + "results in errors for the calculation of the heat loss. It is not advised " * + "to utilize this `StorageBehavior`. Use instead `CyclicRepresentative`.", + maxlog = 1 + ) + end +end + """ EMB.check_link(l::DHPipe, 𝒯, modeltype::EnergyModel, check_timeprofiles::Bool) diff --git a/src/constraint_functions.jl b/src/constraint_functions.jl index 8244477..8048c34 100644 --- a/src/constraint_functions.jl +++ b/src/constraint_functions.jl @@ -218,3 +218,359 @@ function EMB.constraints_capacity( constraints_capacity_installed(m, n, 𝒯, modeltype) end + +""" + EMB.constraints_capacity(m, n::LevelDependentRateTES, 𝒯::TimeStructure, modeltype::OperationalModel) + +Method for creating the constraints on the maximum capacity of a [`LevelDependentRateTES`](@ref). + +It adjusts the capacity constraints to account for state-of-charge dependent charge and +discharge limits. The limits are drawn from the c-rate anchor points on the node (one to three +per direction); the rate allowed in operational period `t` depends on the storage level at the +end of the previous operational period `t_prev`. + +With a single anchor point the limit is a single line. With two or three anchor points the +curve is piecewise linear, and the binary variables `bin_region_charge`/`bin_region_discharge` +(declared in [`EMB.variables_node`](@ref)) together with a big-``M`` formulation select the +active region. Because the node is restricted to an `OperationalModel`, all slopes and +intercepts are fixed parameters, so every big-``M`` is computed tightly from the geometry (the +installed level capacity for the region-selection bounds, and the smallest lift that keeps each +inactive rate segment above the installed rate) instead of a single loose constant. This keeps +the inactive constraints non-binding without weakening the model's relaxation. + +!!! warning "Only operational models are supported" + This method only dispatches on an [`OperationalModel`](@extref EnergyModelsBase.OperationalModel). + With an [`InvestmentModel`](@extref EnergyModelsBase.InvestmentModel) the installed storage, + charge and discharge capacities become decision variables instead of fixed parameters. The + slopes and intercepts of the c-rate curves then depend on those capacity variables, so the + rate limit would multiply a capacity variable by the storage-level variable + ``\\texttt{stor\\_level}[n, t_{prev}]``. This product of two decision variables is a + bilinear (nonconvex quadratic) term, which the linear/mixed-integer-linear solvers used in + this package (*e.g.*, `HiGHS`) cannot handle. A [`LevelDependentRateTES`](@ref) must + therefore be used with an [`OperationalModel`](@extref EnergyModelsBase.OperationalModel). +""" +function EMB.constraints_capacity( + m, + n::LevelDependentRateTES, + 𝒯::TimeStructure, + modeltype::OperationalModel, +) + @constraint(m, [t ∈ 𝒯], m[:stor_level][n, t] <= m[:stor_level_inst][n, t]) + @constraint(m, [t ∈ 𝒯], m[:stor_charge_use][n, t] <= m[:stor_charge_inst][n, t]) + @constraint(m, [t ∈ 𝒯], m[:stor_discharge_use][n, t] <= m[:stor_discharge_inst][n, t]) + + constraints_capacity_installed(m, n, 𝒯, modeltype) + + ### CHARGING ### + + x_1_c = c_rate_points_charge(n)[1][1] + y_1_c = c_rate_points_charge(n)[1][2] + + if length(c_rate_points_charge(n)) == 1 + @info "Using linear C_rate gradient for charging of $n." + + for (t_prev, t) ∈ withprev(𝒯) + if t_prev === nothing + continue + end + + m_charge = (y_1_c - capacity(charge(n), t)) / x_1_c + b_charge = capacity(charge(n), t) + + @constraint( + m, + m[:stor_charge_use][n, t] <= + m_charge * m[:stor_level][n, t_prev] + b_charge + ) + end + elseif length(c_rate_points_charge(n)) == 2 + x_2_c = c_rate_points_charge(n)[2][1] + y_2_c = c_rate_points_charge(n)[2][2] + + @warn "Using piecewise linear C-rate gradients for charging of $n. This introduces binary variables." + + for (t_prev, t) ∈ withprev(𝒯) + if t_prev === nothing + continue + end + + cap_charge = capacity(charge(n), t) + level_cap = capacity(level(n), t_prev) + region_bin = m[:bin_region_charge][n, t, 1] + + @constraint( + m, + m[:stor_level][n, t_prev] <= x_1_c + (level_cap - x_1_c) * region_bin + ) + @constraint(m, m[:stor_level][n, t_prev] >= x_1_c - x_1_c * (1 - region_bin)) + + m_1_charge = (y_1_c - cap_charge) / x_1_c + b_1_charge = cap_charge + m_2_charge = (y_2_c - y_1_c) / (x_2_c - x_1_c) + b_2_charge = y_2_c - m_2_charge * x_2_c + + M_1 = max( + 0.0, + cap_charge - min( + m_1_charge * x_1_c + b_1_charge, + m_1_charge * level_cap + b_1_charge, + ), + ) + M_2 = max(0.0, cap_charge - min(b_2_charge, m_2_charge * x_1_c + b_2_charge)) + + @constraint( + m, + m[:stor_charge_use][n, t] <= + m_1_charge * m[:stor_level][n, t_prev] + b_1_charge + region_bin * M_1 + ) + @constraint( + m, + m[:stor_charge_use][n, t] <= + m_2_charge * m[:stor_level][n, t_prev] + b_2_charge + + (1 - region_bin) * M_2 + ) + end + elseif length(c_rate_points_charge(n)) == 3 + x_2_c = c_rate_points_charge(n)[2][1] + y_2_c = c_rate_points_charge(n)[2][2] + x_3_c = c_rate_points_charge(n)[3][1] + y_3_c = c_rate_points_charge(n)[3][2] + + @warn "Using piecewise linear C-rate gradients for charging of $n. This introduces binary variables." + + for (t_prev, t) ∈ withprev(𝒯) + if t_prev === nothing + continue + end + + cap_charge = capacity(charge(n), t) + level_cap = capacity(level(n), t_prev) + lower_split_bin = m[:bin_region_charge][n, t, 1] + upper_split_bin = m[:bin_region_charge][n, t, 2] + + # Enforce a valid region encoding (the upper split implies the lower split) + @constraint(m, upper_split_bin <= lower_split_bin) + + @constraint( + m, + m[:stor_level][n, t_prev] <= x_1_c + (level_cap - x_1_c) * lower_split_bin + ) + @constraint( + m, + m[:stor_level][n, t_prev] >= x_1_c - x_1_c * (1 - lower_split_bin) + ) + @constraint( + m, + m[:stor_level][n, t_prev] <= x_2_c + (level_cap - x_2_c) * upper_split_bin + ) + @constraint( + m, + m[:stor_level][n, t_prev] >= x_2_c - x_2_c * (1 - upper_split_bin) + ) + + m_1_charge = (y_1_c - cap_charge) / x_1_c + b_1_charge = cap_charge + m_2_charge = (y_2_c - y_1_c) / (x_2_c - x_1_c) + b_2_charge = y_2_c - m_2_charge * x_2_c + m_3_charge = (y_3_c - y_2_c) / (x_3_c - x_2_c) + b_3_charge = y_3_c - m_3_charge * x_3_c + + M_1 = max( + 0.0, + cap_charge - min( + m_1_charge * x_1_c + b_1_charge, + m_1_charge * level_cap + b_1_charge, + ), + ) + M_2 = + max(0.0, cap_charge - min(b_2_charge, m_2_charge * level_cap + b_2_charge)) + M_3 = max(0.0, cap_charge - min(b_3_charge, m_3_charge * x_2_c + b_3_charge)) + + # Region selection: + # - lower_split_bin = 0, upper_split_bin = 0 => region 1 (lowest SOC) + # - lower_split_bin = 1, upper_split_bin = 0 => region 2 (middle SOC) + # - lower_split_bin = 1, upper_split_bin = 1 => region 3 (highest SOC) + @constraint( + m, + m[:stor_charge_use][n, t] <= + m_1_charge * m[:stor_level][n, t_prev] + b_1_charge + lower_split_bin * M_1 + ) + @constraint( + m, + m[:stor_charge_use][n, t] <= + m_2_charge * m[:stor_level][n, t_prev] + b_2_charge + + (1 - lower_split_bin) * M_2 + + upper_split_bin * M_2 + ) + @constraint( + m, + m[:stor_charge_use][n, t] <= + m_3_charge * m[:stor_level][n, t_prev] + b_3_charge + + (1 - upper_split_bin) * M_3 + ) + end + end + + ### DISCHARGING ### + + x_1_d = c_rate_points_discharge(n)[1][1] + y_1_d = c_rate_points_discharge(n)[1][2] + + if length(c_rate_points_discharge(n)) == 1 + @info "Using linear C_rate gradient for discharging of $n." + + for (t_prev, t) ∈ withprev(𝒯) + if t_prev === nothing + continue + end + + m_discharge = + (capacity(discharge(n), t) - y_1_d) / (capacity(level(n), t_prev) - x_1_d) + b_discharge = + capacity(discharge(n), t) - m_discharge * capacity(level(n), t_prev) + + @constraint( + m, + m[:stor_discharge_use][n, t] <= + m_discharge * m[:stor_level][n, t_prev] + b_discharge + ) + end + elseif length(c_rate_points_discharge(n)) == 2 + x_2_d = c_rate_points_discharge(n)[2][1] + y_2_d = c_rate_points_discharge(n)[2][2] + + @warn "Using piecewise linear C-rate gradients for discharging of $n. This introduces binary variables." + + for (t_prev, t) ∈ withprev(𝒯) + if t_prev === nothing + continue + end + + cap_discharge = capacity(discharge(n), t) + level_cap = capacity(level(n), t_prev) + region_bin = m[:bin_region_discharge][n, t, 1] + + @constraint( + m, + m[:stor_level][n, t_prev] <= x_1_d + (level_cap - x_1_d) * region_bin + ) + @constraint(m, m[:stor_level][n, t_prev] >= x_1_d - x_1_d * (1 - region_bin)) + + m_1_discharge = (cap_discharge - y_1_d) / (level_cap - x_1_d) + b_1_discharge = cap_discharge - m_1_discharge * level_cap + m_2_discharge = (y_1_d - y_2_d) / (x_1_d - x_2_d) + b_2_discharge = y_2_d - m_2_discharge * x_2_d + + M_1 = max( + 0.0, + cap_discharge - min(b_1_discharge, m_1_discharge * x_1_d + b_1_discharge), + ) + M_2 = max( + 0.0, + cap_discharge - min( + m_2_discharge * x_1_d + b_2_discharge, + m_2_discharge * level_cap + b_2_discharge, + ), + ) + + @constraint( + m, + m[:stor_discharge_use][n, t] <= + m_1_discharge * m[:stor_level][n, t_prev] + b_1_discharge + + (1 - region_bin) * M_1 + ) + @constraint( + m, + m[:stor_discharge_use][n, t] <= + m_2_discharge * m[:stor_level][n, t_prev] + b_2_discharge + + region_bin * M_2 + ) + end + elseif length(c_rate_points_discharge(n)) == 3 + x_2_d = c_rate_points_discharge(n)[2][1] + y_2_d = c_rate_points_discharge(n)[2][2] + x_3_d = c_rate_points_discharge(n)[3][1] + y_3_d = c_rate_points_discharge(n)[3][2] + + @warn "Using piecewise linear C-rate gradients for discharging of $n. This introduces binary variables." + + for (t_prev, t) ∈ withprev(𝒯) + if t_prev === nothing + continue + end + + cap_discharge = capacity(discharge(n), t) + level_cap = capacity(level(n), t_prev) + upper_split_bin = m[:bin_region_discharge][n, t, 1] + lower_split_bin = m[:bin_region_discharge][n, t, 2] + + # Being in the upper region implies being in the middle/upper combined region + @constraint(m, upper_split_bin <= lower_split_bin) + + @constraint( + m, + m[:stor_level][n, t_prev] <= x_1_d + (level_cap - x_1_d) * upper_split_bin + ) + @constraint( + m, + m[:stor_level][n, t_prev] >= x_1_d - x_1_d * (1 - upper_split_bin) + ) + @constraint( + m, + m[:stor_level][n, t_prev] <= x_2_d + (level_cap - x_2_d) * lower_split_bin + ) + @constraint( + m, + m[:stor_level][n, t_prev] >= x_2_d - x_2_d * (1 - lower_split_bin) + ) + + m_1_discharge = (cap_discharge - y_1_d) / (level_cap - x_1_d) + b_1_discharge = cap_discharge - m_1_discharge * level_cap + m_2_discharge = (y_1_d - y_2_d) / (x_1_d - x_2_d) + b_2_discharge = y_2_d - m_2_discharge * x_2_d + m_3_discharge = (y_2_d - y_3_d) / (x_2_d - x_3_d) + b_3_discharge = y_3_d - m_3_discharge * x_3_d + + M_1 = max( + 0.0, + cap_discharge - min(b_1_discharge, m_1_discharge * x_1_d + b_1_discharge), + ) + M_2 = max( + 0.0, + cap_discharge - + min(b_2_discharge, m_2_discharge * level_cap + b_2_discharge), + ) + M_3 = max( + 0.0, + cap_discharge - min( + m_3_discharge * x_2_d + b_3_discharge, + m_3_discharge * level_cap + b_3_discharge, + ), + ) + + # Region selection: + # - upper_split_bin = 1, lower_split_bin = 1 => region 1 (highest SOC) + # - upper_split_bin = 0, lower_split_bin = 1 => region 2 (middle SOC) + # - upper_split_bin = 0, lower_split_bin = 0 => region 3 (lowest SOC) + @constraint( + m, + m[:stor_discharge_use][n, t] <= + m_1_discharge * m[:stor_level][n, t_prev] + b_1_discharge + + (1 - upper_split_bin) * M_1 + ) + @constraint( + m, + m[:stor_discharge_use][n, t] <= + m_2_discharge * m[:stor_level][n, t_prev] + b_2_discharge + + upper_split_bin * M_2 + + (1 - lower_split_bin) * M_2 + ) + @constraint( + m, + m[:stor_discharge_use][n, t] <= + m_3_discharge * m[:stor_level][n, t_prev] + b_3_discharge + + lower_split_bin * M_3 + ) + end + end +end diff --git a/src/model.jl b/src/model.jl index 4e3f02e..0ba8a02 100644 --- a/src/model.jl +++ b/src/model.jl @@ -40,3 +40,28 @@ function EMB.create_link( @constraint(m, [t ∈ 𝒯, p ∈ inputs(l)], m[:link_in][l, t, p] ≤ m[:link_cap_inst][l, t]) constraints_capacity_installed(m, l, 𝒯, modeltype) end + +""" + EMB.variables_node(m, 𝒩::Vector{<:LevelDependentRateTES}, 𝒯, modeltype::OperationalModel) + +Declare the auxiliary binary variables required by [`LevelDependentRateTES`](@ref) nodes. + +For each node and operational period, two binaries per direction encode which of up to three +piecewise-linear regions of the state-of-charge dependent c-rate curve is active: +- `bin_region_charge[n, t, 1:2]` selects the active region of the charge curve. +- `bin_region_discharge[n, t, 1:2]` selects the active region of the discharge curve. + +The variables are only declared for an [`OperationalModel`](@extref EnergyModelsBase.OperationalModel), +see the note on investment models in [`EMB.constraints_capacity`](@ref). +""" +function EMB.variables_node( + m, + 𝒩::Vector{<:LevelDependentRateTES}, + 𝒯, + modeltype::OperationalModel, +) + # Two binaries encode the active region of the piecewise charge curve (3 regions max) + @variable(m, bin_region_charge[𝒩, 𝒯, 1:2], Bin) + # Two binaries encode the active region of the piecewise discharge curve (3 regions max) + @variable(m, bin_region_discharge[𝒩, 𝒯, 1:2], Bin) +end diff --git a/src/structures/node.jl b/src/structures/node.jl index 02ddc7e..21088bf 100644 --- a/src/structures/node.jl +++ b/src/structures/node.jl @@ -246,7 +246,7 @@ heat losses do not occur while charging or discharging, *i.e.*, they are proport storage level. !!! warning "StorageBehavior" - `BoundRateTES` in its current implementation only supports + `ThermalEnergyStorage` in its current implementation only supports [`CyclicRepresentative`](@extref EnergyModelsBase.CyclicRepresentative) as storage behavior when using [`RepresentativePeriods`](@extref TimeStruct.RepresentativePeriods). This input is not a required input due to the inclusion of a constructor. @@ -262,7 +262,8 @@ storage level. `ThermalEnergyStorage` node. Depending on the chosen type, the discharge parameters can include variable OPEX, fixed OPEX, and/or a capacity. - **`stor_res::Resource`** is the stored [`Resource`](@extref EnergyModelsBase.Resource). -- **`heat_loss_factor::Float64`** are the relative heat losses in percent. +- **`heat_loss_factor::Float64`** are the relative heat losses as a fraction (in the range + ``[0, 1]``) of the storage level of the previous time period. - **`input::Dict{<:Resource,<:Real}`** are the input [`Resource`](@extref EnergyModelsBase.Resource)s with conversion value `Real`. - **`output::Dict{<:Resource,<:Real}`** are the generated [`Resource`](@extref EnergyModelsBase.Resource)s @@ -370,7 +371,8 @@ ratio between the (dis-)charge rate and the installed storage capacity. - **`level::AbstractStorageParameters`** are the level parameters of the `BoundRateTES`. Depending on the chosen type, the level parameters can include variable OPEX and/or fixed OPEX. - **`stor_res::Resource`** is the stored [`Resource`](@extref EnergyModelsBase.Resource). -- **`heat_loss_factor::Float64`** are the relative heat losses in percent. +- **`heat_loss_factor::Float64`** are the relative heat losses as a fraction (in the range + ``[0, 1]``) of the storage level of the previous time period. - **`level_discharge::Float64`** is the ratio of maximum discharge rate and installed storage level. - **`level_charge::Float64`** is the ratio of maximum charge rate and installed storage level. - **`input::Dict{<:Resource,<:Real}`** are the input [`Resource`](@extref EnergyModelsBase.Resource)s @@ -462,7 +464,163 @@ function BoundRateTES( end """ - heat_loss_factor(n::ThermalEnergyStorage) + LevelDependentRateTES{T} <: AbstractTES{T} + +A `LevelDependentRateTES` that functions mostly like a [`RefStorage`](@extref EnergyModelsBase.RefStorage) +with the additional option to include thermal losses and state-of-charge dependent charge and +discharge rates. Heat losses are quantified through a heat loss factor that describes the +amount of thermal energy that is lost in relation to the storage level from the previous +time period. + +Charge and discharge limits can follow c-rate curves defined by one to three anchor points +(`c_rate_points_charge` and `c_rate_points_discharge`). For charging, the curve starts at +`(0, capacity(charge))` and interpolates through the provided points in ascending order of +storage level. For discharging, the curve starts at `(capacity(level), capacity(discharge))` +and interpolates through the provided points in descending order of storage level. With more +than one anchor point, piecewise-linear constraints and binaries are introduced to select the +active region. + +!!! warning "StorageBehavior" + `LevelDependentRateTES` in its current implementation only supports + [`CyclicRepresentative`](@extref EnergyModelsBase.CyclicRepresentative) as storage behavior. + This input is not a required input due to the utilization of an inner constructor. + +# Fields +- **`id`** is the name/identifier of the node. +- **`charge::AbstractStorageParameters`** are the charging parameters of the + `LevelDependentRateTES` node. Depending on the chosen type, the charge parameters can + include variable OPEX, fixed OPEX, and/or a capacity. +- **`level::AbstractStorageParameters`** are the level parameters of the `LevelDependentRateTES`. + Depending on the chosen type, the charge parameters can include variable OPEX and/or fixed OPEX. +- **`discharge::AbstractStorageParameters`** are the discharging parameters of the + `LevelDependentRateTES` node. Depending on the chosen type, the discharge parameters can + include variable OPEX, fixed OPEX, and/or a capacity. +- **`stor_res::Resource`** is the stored [`Resource`](@extref EnergyModelsBase.Resource). +- **`heat_loss_factor::Float64`** is the relative heat loss per operational period duration. +- **`c_rate_points_charge::Vector{<:Vector{<:Real}}`** are one to three `(level, rate)` pairs + in ascending storage-level order describing the maximum charge rate at the previous level. +- **`c_rate_points_discharge::Vector{<:Vector{<:Real}}`** are one to three `(level, rate)` pairs + in descending storage-level order describing the maximum discharge rate at the previous level. +- **`input::Dict{<:Resource,<:Real}`** are the input [`Resource`](@extref EnergyModelsBase.Resource)s + with conversion value `Real`. +- **`output::Dict{<:Resource,<:Real}`** are the generated [`Resource`](@extref EnergyModelsBase.Resource)s + with conversion value `Real`. Only relevant for linking and the stored + [`Resource`](@extref EnergyModelsBase.Resource) as the output value is not utilized in + the calculations. +- **`data::Vector{<:ExtensionData}`** is the additional data (*e.g.*, for investments). The + field `data` is conditional through usage of a constructor. +""" +struct LevelDependentRateTES{T} <: AbstractTES{T} + id::Any + charge::EMB.AbstractStorageParameters + level::EMB.AbstractStorageParameters + discharge::EMB.AbstractStorageParameters + stor_res::Resource + heat_loss_factor::Float64 + c_rate_points_charge::Vector{<:Vector{<:Real}} + c_rate_points_discharge::Vector{<:Vector{<:Real}} + input::Dict{<:Resource,<:Real} + output::Dict{<:Resource,<:Real} + data::Vector{<:ExtensionData} +end + +function LevelDependentRateTES{T}( + id, + charge::EMB.AbstractStorageParameters, + level::EMB.AbstractStorageParameters, + discharge::EMB.AbstractStorageParameters, + stor_res::Resource, + heat_loss_factor::Float64, + c_rate_points_charge::Vector{<:Vector{<:Real}}, + c_rate_points_discharge::Vector{<:Vector{<:Real}}, + input::Dict{<:Resource,<:Real}, + output::Dict{<:Resource,<:Real}, +) where {T<:EMB.StorageBehavior} + return LevelDependentRateTES{T}( + id, + charge, + level, + discharge, + stor_res, + heat_loss_factor, + c_rate_points_charge, + c_rate_points_discharge, + input, + output, + ExtensionData[], + ) +end + +function LevelDependentRateTES( + id::Any, + charge::EMB.AbstractStorageParameters, + level::EMB.AbstractStorageParameters, + discharge::EMB.AbstractStorageParameters, + stor_res::Resource, + heat_loss_factor::Float64, + c_rate_points_charge::Vector{<:Vector{<:Real}}, + c_rate_points_discharge::Vector{<:Vector{<:Real}}, + input::Dict{<:Resource,<:Real}, + output::Dict{<:Resource,<:Real}, + data::Vector{<:ExtensionData}, +) + return LevelDependentRateTES{CyclicRepresentative}( + id, + charge, + level, + discharge, + stor_res, + heat_loss_factor, + c_rate_points_charge, + c_rate_points_discharge, + input, + output, + data, + ) +end + +function LevelDependentRateTES( + id::Any, + charge::EMB.AbstractStorageParameters, + level::EMB.AbstractStorageParameters, + discharge::EMB.AbstractStorageParameters, + stor_res::Resource, + heat_loss_factor::Float64, + c_rate_points_charge::Vector{<:Vector{<:Real}}, + c_rate_points_discharge::Vector{<:Vector{<:Real}}, + input::Dict{<:Resource,<:Real}, + output::Dict{<:Resource,<:Real}, +) + return LevelDependentRateTES{CyclicRepresentative}( + id, + charge, + level, + discharge, + stor_res, + heat_loss_factor, + c_rate_points_charge, + c_rate_points_discharge, + input, + output, + ExtensionData[], + ) +end + +""" + c_rate_points_charge(n::LevelDependentRateTES) + +Return the `(level, rate)` breakpoints that shape the charge c-rate curve. +""" +c_rate_points_charge(n::LevelDependentRateTES) = n.c_rate_points_charge +""" + c_rate_points_discharge(n::LevelDependentRateTES) + +Return the `(level, rate)` breakpoints that shape the discharge c-rate curve. +""" +c_rate_points_discharge(n::LevelDependentRateTES) = n.c_rate_points_discharge + +""" + heat_loss_factor(n::AbstractTES) Returns the heat loss factor for storage `n`. """ diff --git a/src/visualization.jl b/src/visualization.jl new file mode 100755 index 0000000..b1684bb --- /dev/null +++ b/src/visualization.jl @@ -0,0 +1,28 @@ +""" + visualize_c_rates( + charge_capacity, + discharge_capacity, + level_capacity, + c_rate_points_charge, + c_rate_points_discharge, + ) + +Plot the *theoretical* maximum charge and discharge rate curves of a +[`LevelDependentRateTES`](@ref) as a function of the storage level, so the user can inspect +the (dis-)charge curves their input produces before solving a model. + +The arguments mirror the corresponding fields of a [`LevelDependentRateTES`](@ref): +- **`charge_capacity`** is the installed charge rate (the curve's value at an empty storage). +- **`discharge_capacity`** is the installed discharge rate (the curve's value at a full storage). +- **`level_capacity`** is the installed storage level capacity. +- **`c_rate_points_charge`** are the one to three `[level, rate]` anchor points of the charge + curve, in ascending storage-level order. +- **`c_rate_points_discharge`** are the one to three `[level, rate]` anchor points of the + discharge curve, in descending storage-level order. + +!!! note "Plotting extension" + `visualize_c_rates` is provided through a package extension that is only loaded once the + `Plots` package is available. Run `using Plots` before calling this function. Without + `Plots` loaded, calling it raises an informative error. +""" +function visualize_c_rates end diff --git a/test/test_checks.jl b/test/test_checks.jl index 4b816a2..c9d062b 100644 --- a/test/test_checks.jl +++ b/test/test_checks.jl @@ -166,6 +166,84 @@ end EnergyModelsBase.TEST_ENV = false end +# Test that the fields of a LevelDependentRateTES are correctly checked +# - EMB.check_node(n::LevelDependentRateTES{T}, 𝒯, modeltype::EnergyModel, check_timeprofiles::Bool) +@testitem "Checks - LevelDependentRateTES" setup = [TESTestData] begin + using EnergyModelsBase + using TimeStruct + + # Set the global to true to suppress the error message + EnergyModelsBase.TEST_ENV = true + + # Function for setting up the system for testing a `LevelDependentRateTES` node + check_graph = TESTestData.tes_test_case + type = LevelDependentRateTES + + # Resources used in the analysis + heat = ResourceCarrier("Heat", 0.0) + + # Test that a wrong capacity is caught by the checks + @test_throws AssertionError check_graph(; type, charge_cap = FixedProfile(-25)) + @test_throws AssertionError check_graph(; type, level_cap = FixedProfile(-25)) + @test_throws AssertionError check_graph(; type, discharge_cap = FixedProfile(-25)) + + # Test that a wrong heat loss factor is caught by the checks + @test_throws AssertionError check_graph(; type, heat_loss_factor = 1.2) + @test_throws AssertionError check_graph(; type, heat_loss_factor = -0.4) + + # Test that a wrong input or output is caught by the checks + @test_throws AssertionError check_graph(; type, input = Dict(heat => -0.5)) + @test_throws AssertionError check_graph(; type, output = Dict(heat => -0.5)) + + # Test that invalid charge anchor points are caught by the checks + # Too many points + @test_throws AssertionError check_graph(; + type, + c_rate_points_charge = [[2.0, 10.0], [6.0, 8.0], [10.0, 6.0], [18.0, 4.0]], + ) + # Levels not strictly ascending + @test_throws AssertionError check_graph(; + type, + c_rate_points_charge = [[10.0, 7.0], [2.0, 10.0]], + ) + # Level above the installed level capacity (20) + @test_throws AssertionError check_graph(; + type, + c_rate_points_charge = [[2.0, 10.0], [25.0, 4.0]], + ) + # Negative rate + @test_throws AssertionError check_graph(; + type, + c_rate_points_charge = [[2.0, -10.0]], + ) + + # Test that invalid discharge anchor points are caught by the checks + # Levels not strictly descending + @test_throws AssertionError check_graph(; + type, + c_rate_points_discharge = [[2.0, 1.0], [18.0, 5.0]], + ) + # Level not below the installed level capacity (20) + @test_throws AssertionError check_graph(; + type, + c_rate_points_discharge = [[20.0, 5.0], [10.0, 3.0]], + ) + + # Test that the warning regarding the time structure is thrown + msg = + "Using `CyclicStrategic` with a `LevelDependentRateTES` and `RepresentativePeriods` " * + "results in errors for the calculation of the heat loss. It is not advised " * + "to utilize this `StorageBehavior`. Use instead `CyclicRepresentative`." + type_cs = LevelDependentRateTES{CyclicStrategic} + oper = RepresentativePeriods(8760, [0.5, 0.5], SimpleTimes(4, 1)) + # `match_mode = :any` because building the node also emits the piecewise-region warning + # from the constraint builder; we only assert the storage-behavior warning appears. + @test_logs (:warn, msg) match_mode = :any check_graph(; type = type_cs, oper) + + # Set the global again to false + EnergyModelsBase.TEST_ENV = false +end + # Test that the fields of a DHPipe are correctly checked # - EMB.check_link(l::DHPipe, 𝒯, modeltype::EnergyModel, check_timeprofiles::Bool) @testitem "Checks - DHPipe" setup = [DHPipeTestData] begin diff --git a/test/test_heat_storage.jl b/test/test_heat_storage.jl index 6c7c60a..9d09899 100644 --- a/test/test_heat_storage.jl +++ b/test/test_heat_storage.jl @@ -25,6 +25,8 @@ output = Dict(heat_use => 1), level_charge = 0.125, level_discharge = 0.25, + c_rate_points_charge = [[2.0, 10.0], [10.0, 7.0], [18.0, 4.0]], + c_rate_points_discharge = [[18.0, 5.0], [10.0, 3.0], [2.0, 1.0]], supply_cap = FixedProfile(0.7), oper = SimpleTimes(4, 2), ) @@ -61,6 +63,19 @@ input, output, ) + elseif type <: LevelDependentRateTES + tes = type( + "TES", + StorCap(charge_cap), + StorCap(level_cap), + StorCap(discharge_cap), + stor_res, + heat_loss_factor, + c_rate_points_charge, + c_rate_points_discharge, + input, + output, + ) end # Create the test nodes @@ -357,3 +372,170 @@ end @test sum(value.(m[:stor_discharge_use][tes, t]) ≈ 0.5 for t ∈ 𝒯) == 2 end end + +@testitem "LevelDependentRateTES" setup = [TESTestData] begin + using JuMP + using EnergyModelsBase + using TimeStruct + const EMH = EnergyModelsHeat + + # Create the case and modeltype + type = LevelDependentRateTES + m, case, modeltype = TESTestData.tes_test_case(; type) + optimize!(m) + + # Extract the individual elements and resources + tes = get_nodes(case)[2] + heat_use = get_products(case)[2] + 𝒯 = get_time_struct(case) + 𝒯ᴵⁿᵛ = strategic_periods(𝒯) + lvl = value.(m[:stor_level][tes, :]) + + @testset "LevelDependentRateTES - Utility functions" begin + # Test the EMB extraction functions + @test charge(tes) == StorCap(FixedProfile(10)) + @test level(tes) == StorCap(FixedProfile(20)) + @test discharge(tes) == StorCap(FixedProfile(5)) + @test storage_resource(tes) == heat_use + @test inputs(tes) == [heat_use] + @test outputs(tes) == [heat_use] + @test node_data(tes) == ExtensionData[] + + # Test the EMHEat extraction functions + @test EMH.heat_loss_factor(tes) == 0.05 + @test EMH.c_rate_points_charge(tes) == [[2.0, 10.0], [10.0, 7.0], [18.0, 4.0]] + @test EMH.c_rate_points_discharge(tes) == [[18.0, 5.0], [10.0, 3.0], [2.0, 1.0]] + end + + @testset "LevelDependentRateTES - Constructor" begin + c_charge = [[2.0, 10.0], [10.0, 7.0], [18.0, 4.0]] + c_discharge = [[18.0, 5.0], [10.0, 3.0], [2.0, 1.0]] + tes_1 = LevelDependentRateTES{CyclicRepresentative}( + "TES", + StorCap(FixedProfile(10)), + StorCap(FixedProfile(20)), + StorCap(FixedProfile(5)), + heat_use, + 0.05, + c_charge, + c_discharge, + Dict(heat_use => 1), + Dict(heat_use => 1), + ExtensionData[], + ) + tes_2 = LevelDependentRateTES{CyclicRepresentative}( + "TES", + StorCap(FixedProfile(10)), + StorCap(FixedProfile(20)), + StorCap(FixedProfile(5)), + heat_use, + 0.05, + c_charge, + c_discharge, + Dict(heat_use => 1), + Dict(heat_use => 1), + ) + tes_3 = LevelDependentRateTES( + "TES", + StorCap(FixedProfile(10)), + StorCap(FixedProfile(20)), + StorCap(FixedProfile(5)), + heat_use, + 0.05, + c_charge, + c_discharge, + Dict(heat_use => 1), + Dict(heat_use => 1), + ExtensionData[], + ) + for field ∈ fieldnames(LevelDependentRateTES) + @test getproperty(tes_1, field) == getproperty(tes_2, field) + @test getproperty(tes_1, field) == getproperty(tes_3, field) + end + end + + @testset "LevelDependentRateTES - Constraints-level" begin + # Test that the loss is correctly included in the first periods + @test all( + isapprox( + lvl[t], + lvl[t_prev] + + value.(m[:stor_level_Δ_op][tes, t]) * duration(t) - + lvl[t_prev] * EMH.heat_loss_factor(tes) * duration(t); + atol = 1e-6) for + t_inv ∈ 𝒯ᴵⁿᵛ for (t_prev, t) ∈ withprev(t_inv) if !isnothing(t_prev) + ) + + # Test that the total loss is correct + heat_stored = sum(lvl[t] * duration(t) for t ∈ 𝒯) + heat_in = sum(value.(m[:flow_in][tes, t, heat_use]) * duration(t) for t ∈ 𝒯) + heat_out = sum(value.(m[:flow_out][tes, t, heat_use]) * duration(t) for t ∈ 𝒯) + @test heat_stored * EMH.heat_loss_factor(tes) ≈ heat_in - heat_out + end + + @testset "LevelDependentRateTES - Constraints-capacity" begin + # The auxiliary region-selection binaries are declared + @test haskey(object_dictionary(m), :bin_region_charge) + @test haskey(object_dictionary(m), :bin_region_discharge) + + # The installed capacities are respected + @test all(value.(m[:stor_level][tes, t]) ≤ capacity(level(tes), t) + 1e-6 for t ∈ 𝒯) + @test all( + value.(m[:stor_charge_use][tes, t]) ≤ capacity(charge(tes), t) + 1e-6 for + t ∈ 𝒯 + ) + @test all( + value.(m[:stor_discharge_use][tes, t]) ≤ capacity(discharge(tes), t) + 1e-6 for + t ∈ 𝒯 + ) + + # The used charge and discharge rates never exceed the theoretical c-rate curve + # evaluated at the storage level of the previous operational period. The curves are + # recomputed here independently from the anchor points, so this checks that the + # constraints reproduce the intended state-of-charge dependent limits. + pts_c = EMH.c_rate_points_charge(tes) + pts_d = EMH.c_rate_points_discharge(tes) + + # Maximum charge rate at level `x`: piecewise line from (0, cap) through the ascending + # anchor points. + function charge_curve(x, cap) + x1, y1 = pts_c[1] + seg1 = (y1 - cap) / x1 * x + cap + length(pts_c) == 1 && return seg1 + x2, y2 = pts_c[2] + seg2 = (y2 - y1) / (x2 - x1) * (x - x2) + y2 + length(pts_c) == 2 && return x ≤ x1 ? seg1 : seg2 + x3, y3 = pts_c[3] + seg3 = (y3 - y2) / (x3 - x2) * (x - x3) + y3 + return x ≤ x1 ? seg1 : (x ≤ x2 ? seg2 : seg3) + end + + # Maximum discharge rate at level `x`: piecewise line from (cap_level, cap) through the + # descending anchor points. + function discharge_curve(x, cap, cap_level) + x1, y1 = pts_d[1] + seg1 = (cap - y1) / (cap_level - x1) * (x - cap_level) + cap + length(pts_d) == 1 && return seg1 + x2, y2 = pts_d[2] + seg2 = (y1 - y2) / (x1 - x2) * (x - x2) + y2 + length(pts_d) == 2 && return x ≥ x1 ? seg1 : seg2 + x3, y3 = pts_d[3] + seg3 = (y2 - y3) / (x2 - x3) * (x - x3) + y3 + return x ≥ x1 ? seg1 : (x ≥ x2 ? seg2 : seg3) + end + + @test all( + value(m[:stor_charge_use][tes, t]) ≤ + charge_curve(lvl[t_prev], capacity(charge(tes), t)) + 1e-6 for + (t_prev, t) ∈ withprev(𝒯) if !isnothing(t_prev) + ) + @test all( + value(m[:stor_discharge_use][tes, t]) ≤ + discharge_curve( + lvl[t_prev], + capacity(discharge(tes), t), + capacity(level(tes), t_prev), + ) + 1e-6 for (t_prev, t) ∈ withprev(𝒯) if !isnothing(t_prev) + ) + end +end