From 41040b0619396d4882386e99612f5a310c031a4b Mon Sep 17 00:00:00 2001 From: tillh Date: Fri, 3 Jul 2026 17:56:22 +0200 Subject: [PATCH 1/6] Add LevelDependentRateTES storage node Introduce the advanced thermal energy storage node LevelDependentRateTES, whose charge and discharge rate limits depend on the storage level via one to three c-rate anchor points per direction. Includes the struct and constructors, input checks, the state-of-charge dependent capacity constraints (piecewise-linear with a tight per-constraint big-M), and the region-selection binary variables. Also declares the visualize_c_rates stub that the Plots extension fills in. --- src/EnergyModelsHeat.jl | 6 +- src/checks.jl | 101 ++++++++++ src/constraint_functions.jl | 358 +++++++++++++++++++++++++++++++++++- src/model.jl | 25 +++ src/structures/node.jl | 166 ++++++++++++++++- src/visualization.jl | 28 +++ 6 files changed, 677 insertions(+), 7 deletions(-) create mode 100755 src/visualization.jl 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..f226a89 100644 --- a/src/constraint_functions.jl +++ b/src/constraint_functions.jl @@ -1,4 +1,4 @@ -""" + """ EMB.constraints_capacity(m, n::HeatPump, ๐’ฏ::TimeStructure, modeltype::EnergyModel) Method for creating the constraints on the maximum capacity of a [`HeatPump`](@ref). @@ -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 From 2fce2cc529946c13cfe9125647bebdf0404d1498 Mon Sep 17 00:00:00 2001 From: tillh Date: Fri, 3 Jul 2026 17:56:22 +0200 Subject: [PATCH 2/6] Add Plots extension for visualize_c_rates Provide the visualize_c_rates implementation as a Plots package extension so the core package keeps no hard plotting dependency. Declare Plots under [weakdeps]/[extensions] with a [compat] entry. --- Project.toml | 7 ++ ext/PlotsExt.jl | 191 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 198 insertions(+) create mode 100644 ext/PlotsExt.jl 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/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 From e7faa1f9e4814d747a6b96b3db725fdc1dfd49cd Mon Sep 17 00:00:00 2001 From: tillh Date: Fri, 3 Jul 2026 17:56:22 +0200 Subject: [PATCH 3/6] Add tests for LevelDependentRateTES Cover the input checks, the constructors, the level and capacity constraints, and a coupling check asserting the used charge/discharge rates never exceed the theoretical c-rate curve evaluated at the previous storage level. --- test/test_checks.jl | 78 ++++++++++++++++ test/test_heat_storage.jl | 182 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 260 insertions(+) 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 From 6b1a164d5307026d455d9332166c55dcb1393513 Mon Sep 17 00:00:00 2001 From: tillh Date: Fri, 3 Jul 2026 17:56:22 +0200 Subject: [PATCH 4/6] Document LevelDependentRateTES and c-rate visualization Add the node documentation page (with the c-rate curve plot and a worked example of the anchor points), register the node and visualize_c_rates in the public/internal API references, and add the NEWS entry. --- NEWS.md | 16 ++ docs/make.jl | 1 + docs/src/library/internals/methods-EMB.md | 1 + docs/src/library/internals/methods-EMH.md | 2 + docs/src/library/public.md | 7 + docs/src/nodes/crate_curves_example.png | Bin 0 -> 41685 bytes docs/src/nodes/leveldependentratetes.md | 169 ++++++++++++++++++++++ 7 files changed, 196 insertions(+) create mode 100644 docs/src/nodes/crate_curves_example.png create mode 100644 docs/src/nodes/leveldependentratetes.md 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/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 0000000000000000000000000000000000000000..f6dd2f030af64a8494ae6c89886eb3ccbc08fb10 GIT binary patch literal 41685 zcmYJb2Rzm7|37{Z*^!w|vS%`~J1DaE$liPJosJ|Sgoy0DiI7!7R+JSPm60ubh2QIR ze?R~Kx$nn)SMSsNeZ8;ib-iBC@w!B7sNEyNr@==c5JXCf@>&Q4783%2v4@KZe{!E? zmLL9yW2t;k9&v&GC%d&M34vfjD9PW^@y}SF@sk~0y14v@y$Qo_UE<)bb#(EQpdTFf zJ)CC5R=Ip8bSud8>3znx$HZ4>^YkLcHSXAOX$PDRrFd!I_O z_r@08^V2^Iy=2c=l`g?k5s1|*YJUHI5uSpB&x2ScyF=nO+k%IO7yaZ32M{IqoOSn)f=wKjSxU9Iic=m0J z?O?iLoh@foR(@gO$nl*vXp;HtQ zZQ%(QJv}`;d(rXk@1mljuseyhwY7Qo;%8ff4wZ@X)Nb6kQKnxdA|#}(>SSd#zqaur^d%Su5i17DNmW=PZ z!N$RzM>xzVEp6?!ka%BTUtH!UW3>=v68WAucshN7mP6|n0is!ewj}~9n*x`ZnBlgk zwTDMdbaZdVZ7*kMXL5Xn`}dz_WzE&v^#A(x>-_KTP`vBMuCA_6pKzI3$?@lZ{UW9p z_uiat%FfP)Uj^&EjEKNwWMouVA7Yk!9Q<>;MIDL!@SIuh=FOY_lCyg&BOB}M0))!y z>KDJW!h%CWvZO;!{`~9=oThL6{95YtS67UX-DiO`6_m!%YpDXQ5>E#$EiF^i*A*2N zk&%(LU&hB%9%t`cjFsvJ1O(LC4OM=bEzxFAQ|bQcP9u(ZY_T?8K`Y_A9!b=HQXBy+}?j;Y++<5Osx-LDJt{YBSE4 zkB=`fFwnB*s8#v*K-wdc6S#pRzhj(~0VO4+3R|7vp+Vl|7*?e;694RZX!!NJSKbj5)tH$_ELX+(51HMe&bl}o=24{Ov{6*M-AF){tv z`2LoTxaeLyLuF+pEGjBLL50C(wMGv3BS#w?fDWn<@@))j*O_P zsHmu_>gp@2s`B#l`}z9n>FeLa?Ooh?`I4G1b!}lmM^A5QYb)#3tM7Go>J0Jvy1I)C z3m+FfovamARLVwSQCT_Y=nu1x*drXGz%8}9y1FdMfW)k< zq518T{k4e(m+w<+iUnWU)nBu+v9X}&=;%%kx6G~^XlW&CaCm!oSO`$ZGf3T6R2&{3 zKR-X+Dk<=6`nxhxw7tE3Gi~7W=g*y;k4{LLP+Z*HNl8g11YUh}l*CaZl_nA>jr;fS zD=0+9#1On1)SEGJa9Ek|iLa}z9UUEo<%6deoscE;T3c9HI5;>QA0Ol3;Iszrd%C+H zEeCa3-(**LSsiw94pScB?;q?Pf=I_mCy+#>H^JV(dGyc%hWM4rE#LPFrBWg%Cnqj5 zVW&=%JFj6as>`zE9IvdbY-_4sx#`1Hy_3^L+Qr_8`lwa}6Ks~`tg08VtzM+0jFYOf zC9NzkGodgsF{J}}AbjPw);?Sl_vS!yFf%_qNud#+DNs&zny7>!a(eNZG?O!X|M~Hr zh95q^J)0yM)F;by3=MAw>{yR|O?q$6W74csZ^kMn*4X_-5T$W(ezqo!Qi!2sMxCFm zRb8q0!zpC|o2zsCqCKYg!CP0@5T#BK7(;17_v%+a2snSvxGmhx9En0{ghbsRA0LPP zt8G#QNf0*&*65Af+-swrzP_hMcs`HMx|%pVd;^)6;e1cOiu=1 ziB1%U2;aK(LS`ItpjJCq>1AP$Umt0t-qibyq1$LeTF3>AM{(4lK^aIlKXz?ZyxC2_ehcuj7|*IM(f z&E6|JKg{}mZ9y`+i(BtBu?fkw%B-CjAD>3Tmj`J-oPFogonQ`Ae}6F~1Y}t^`>&u7 zy@xjT_Ma9#ORnSN;|K06C>a~mA&X(Eg-t{M9(kaKVP}=NEqWe9(32ovzkYr3cURrg zGB+#B46-yS627ThY31tr+4|;Nx4HAxu`<{i!^6W^mx%{5CEjkzW#hCiE-i)HN5)ec zK-LMCz8KDtBPAv6GcP^8fonNhtcjb`7IHe<8&AXZ=~Wm-TD*Mx__5<1bqa%&loWVKi~CQcyZfJHZhghXXHTAJ=<7p3yb^K2 zzj7r$I(o9+i7K^!Y;;sbO^r>dHa~yi{iDke;gBnizhAw273tF8`>;Ov*c0_UGO`QO z5*Znpm8B)0S(~Kg$0x z5bHA61L<2h+|tqvodUG9uW0=K1_TDOpctfsk8(^jgp(j8718zEQ5fLn_^geO86n=WHlrfz|w%&$JJ1P}_f2P@|Q(*mbrbIoN+A9&&$(2Ym7I2Vr(at?EPheHrAH%y8p*|;#P^BDrgx#Uf$FVEv-Zv~gRlsG%Ts7NO0pxor0FZ^3z1SL3z#Hh|TcJ&5QRdvu| zq%gYcS5G|Sb(FV{4=EB5lZ}l{7DGvPb~qDmNXQw#BrED)gqI$M8W|Xf-M&phL4g;E z9Q%t7y;k{8wsm%c=>-}b#W6Pck-xUuFKh*wjIn^P*QR{C_49~%d*t1)^+sz z>(tc4rb&c-g!yv`gqebz92PSacZ`#!!^lrZLO_? z&jdsHgoNs7hkd&#;<`PEl&`@dNXPJFrrw1F5 zF^enj##$HcK%nSCu)calOG`WP?L|{l0KA2fQTl3JMrNk2eqOUkf!0c?ZkZ#>qHYo( zlCG}qsQLj}*T;`cDCt1ukC4=0W2IzdWW0J+CjM5z!(;oNUrtd`6v;`%DNcwpB6L5VQAUrqTCd7f%Hx;(>KY&vm9(2Sd&^%SIqDIQ=y8xkTDx=S zPD4Y3RrfPN0fFqtWn16h^2^AW@ug0fm}}DC4h;>hce5KA87Z6Z7GI@?U&sVL%UWxt zR|u94TOQ5W2EOK7|$Gw)<{OL9K`C z;}ThRd8E&9?@XR92psW0vLSVlrD>T59!(x3#s=in=ar`$PJw zhkZFZDvvt@nH4YtAcwqn*>S~puW%y<$;k~6!O4JInNXHiRs{I?I%Ne@06zc)=S^w# z^!x=V1s_qDsW({(bSjgq^p9 zg(Hy`5E3(t-E_ouBZUy=E&zJ`VXYtuWJv{y`2YC!l3BKj#NFS42W8&&T5_6aL7lfQKJ^gB!A3c3$&c%MtkbfU!s@T|A zx9P^5;JdUwXHL_N9;IWtYHDh+Q4CtUVxyY^*?gOlVdsDAh)I|MN&P|FmK7-WW`lJx- zw1|iZyP?-1_NtIbl$HCZr;STH3m2 z`;5}T#p0Jv7kZO)iu&0tZ>9m1qKxT+9MtT)K@;CYD`<&Y#TT+2ki|VdIl0f8NK;i& zSy^~8B+sI`k^Asuz3PlhH{|r-225$^sBH0Dr~VY)dWi$pt9ZiRcs<(yl*lgCc>Puz zqwrn-#+=%JFP5f2Jik@M&hEGCOmhK%Cn4j}>}H56?FgpT<@8+tprCoN(1lM4S4?;c zLuCLgL(#>%vpRe8hPb#9OG4wjcVWe5WZ|}tE5B1=*ej`JiYx5X?r+VZp$()807Afh z{i;k`Nsviws->rkTbaO^dRJE+SK9*|Z9jecW%*Iv*3^GLw-d&AjuI3H>X^#kygi7-8@JLGjg>c{fJv1>kRbLb8 z@J>}#6%c6{0n{JJ2Yr?`{*tCMWEdKG^{~!=awyt^9{Ne_9-suaP*rkai- z0ahXAxu6Jm4WMkL$-4&Y-m4)YA&a+&u3fv<+}uq4nAfyLOj1&kg~f9FR~PkT2PY>; z<(#5wsv{af2Poids;a7VJz+;bqX7AKb$8=jxud-6K7Xz}t!Lx{I6vJQ zB#JCn&ypT}dwRH~pV>Rp(edMFXOzzE0wZ&UyLaIgxq>iX2sI_K1|0qQIVIOG;j^k^ zHCu8RxW76E!Q(Xj?UQ9qcPwq{%a_H*yfMEeG;it*sgVh zGA5PJ#K=l7kS*!9`;QKQdx@-_ot;jMC7GF-o-Bh`5)@e>lR`~nOb^~)Bj8D7wYfQo zW?_UvJsvGEMo>Qb#!TCSR~8oXo8$8GtmNfUtKBj>0NbvSk;$pR)n_nWU-^P@#)dynftg4LtJrI3L^L9!MO-s_mq>ap*8XK2aS4+=nu!hSm z7JB%Gun^KS-`~2y*r2)?@TAn&`#{7<+nf~t{TrO1Xl-w==tEkP*lRu7B@M8k)VyzS z5O6Cm7gz0@H|ullGEtqqy}b$d3;HNtvq_nbb=amITS*BLQ~ysGfNu21$c^_oBq_f;hJN$T&yC z+SazTq+~4gMXZ*GhoCS?Y25p}1qB7~-@l)noHXW}fC@V{rn|4Glo-6fs_W?J2ns-D zr7P5>>FH@l@+VvMp=T$+@_;pAW(+pj`C#X0 zP(zUJw?bB$0|kYN<}Pgn6_ z08}vzlZPpnKnMCKR9@SI8^UK)h~s~M0%!yM*)PY)Qdr*FGBYv~2W-6hAd$7o^bC9@V>i z7Xg|G>6`NbVB^7of#Rq4)zzzj)2hFgba-Ndh5qVr<@sgg4K_B@)<9uMO#%WMlP^RK6m0xF4Ap-S4$iV}#88$EZz$M5E7{17vpzjia5B7zK) z$wn(uo{;@gcrh7fRE%I-HNIty0i#sbsMkXUl8+DpKrA1e;+0>z!_ltGL7sK0KoKF6~pENe4hT{aQ{(XX%&C#xKxhFd|$OqtAJP{}dSJt{_>L1HAO zWs)T0gw;9xx8L-fQ#eqH$^eRfTwL5Nppb)&@7S1kG&R3u-+A;oOU4b@{FF`dt5*`i zM}Oj7W8a$(qzOg!FNPXYQ^x4dUn|I_$PFa|WbKZ$A5*M7`f93RA>wN8`e(#61?R<7 zg^7(#K_J#ekM0n53MiX9eddaapHe`NbepRzEBjugp0%^H^W@1B$oIfM74F@;0*oQ_ z-!io5AJ1t2I?~7}D-d=Ot~`&tc!TxFaJ8ueGOMer>*xgJDD)L z&ryqiO-YHv*(MFFPq8XFg~2s8@<{L(a3iG8wFQq1igf3?tQ?+ z(4hH=8V(;R**y*(Tg7KySVe!P+=I|+{EH{0`0uk`DV3nUyF4#ERcxLQAIE;JfQvw6 zW3y0r!llKC;D-97*&_w-9*hd&f8lSGUAq=_qr`*Ig(0ll=oy~q^N4O7Y)S+oJmv1u zx|(h^W4Gn5otydcpOrW)qwD|n&Kr)zH%B|K?vWsyt*EMb>6hUO=YLn|o+cvS_YUSG zD}~+5@A|!TCmlxn0tf`}CC()L^Jx{_{!x1IHNsa)essg7cOGKH;1rm8zY@GFnB-K! z5}Mr`!;QKkCR>D)ak*^e)8JsM*Kgt(PJ|RB+oax)?Rz>-pbLTqJnqJdkd-SA_BSM_ zcv9x{$Ih*zDINVHm1#-eu%+s{`@za>AVpfx9e0*Gx~43m5F-beGq8I!ua2{vfB3e z7&5b8G*f@CISDq-YLX zQO`I1c1Ne0hv_eer$hKdCK%x0x*rpNns~FsHp_yh^)dR=p*_XOtv7@aT-eyy0BY0I zx7dccd3ZpTZK$dFx$SROx7gR$2iqX!#fxv>zIm^W$97rj>gbdk)vp3T-Y5O54bUKw zHJaqtwm)c%)P#E2h?Xm`cxtXXFDT0DJ;tw{U5BCj%PhZl6oK740Wv!9I?e}_e8%4@ z9*BalP&K&@ijIwoORH5+TzGi6S?F=smoI|c+-1$pC&1y%=#O{3w`N1$M_>d0U5A1U zsS1X*YS0AgE-ZB9S@Hq|x$y6d(aCdN!@%fQZJfn(hsIwWFg0gasiEj%7%0rD&m& zjh&tBQciQi^73+dIgd?@3L6N9thg-TI>^W{GP5Qq0uPOei2=JrgF~m1^J+Bw9;AF} z@9=wZbby*HKrhhXkn~-r9#bTd_w@8MHqO{F04}i7{?Sr!jDsBC%F1euuOERJ%oO-z ztGg`qtfCTfVP%Ojtg$pMe)-bJ_wVoM>x=U65HfoM%oMchE-ow6j={vj0=K0RpcptV zzH8&EOe9Pw(-!|P4K5E>mH6TD9)5nQz1F3)@N@$JMm@2o7bVB=eH(}4O-?eerGd2D zUWKWt6E6j=UaUc$3IYW3=<*eCqF8l%XPUfvW@cJJLn6eM(y6#hUmqVo^pt`d1O!?Vl5i67(7*nH zfhWN7y5ku^W^q9c1L*)KEbMrZii#=~z;5u&cS3}HE^F@>o40Y-3FlRf`fk=0NDw1@ z9L6%sB&?z{1@V8)bBl0!5v{jTjlA4FF96;KmWDTsT|HCW2QVs#@fsT@&E86CYK06A zl75@h@ItQNs&mZRL#3sq0T<5#-2wUss^aa83ujRgk@3O7jaCUjbZ4iBXy?wyr#?1T zRY&JVS=ky2j|^t%5eW;f<&TbsAkfgErD+a6zD*$b@I)kT)4O*dwEiv13Wajcz`!ud z0Qv)v1hcka2WMwvP0gRbe%*vhh296>9FWKe9j@5uXgpk8g9hhoV@E*fVNy}UKwF^( zQPR`5f^MN*q_#2LwEN1f?G7%;N@w7pczAe#%LVlXkPfwo^UK)SL8{g>L^#`h(+qa8 zdJ0nfjrU1nm79>@eq1;5d6=ws2cOkDs#u|%iMt?Q@b{Tq3+P(w%gbW0>`)P)x^aCU z8yg$SxJ^b*4t652i}H$!k^|7JaImrW!J`VhIL%5*(vGd@UGx-rbv8@?Q%sFQ9sDei z7*!nv^3ZZU0so-l?$E$_V&wTW%{H@P?tb&fOWWoJJTLE z_wnNbD3OICuGZG%G&DxFjcgM;IFuQgOdR;6+3i_NVq zH&<7+4RYp>zjv>*ut09hzQKHBEP{*{PxD4kK~clREp&qQsA%puF=Om*s3dIG9?kET zoxXBm;O*@V)|8CT>L_@{B@J%h8#a2A*uhK@5)o+&JqxtA=NOB#ppNN+^Vz52uqy-kK>>5a>F`6)?gtxnVVOGq`r8>4zasR0Lesy18BldEgp=Tt>l9H0XgG>oOK-V*EAwh5*^!x7ai?ntBsW%R=&7MU^ zuY9YyJ_dfCS>6rhre(zc^#TN!)TF}h7n6{{f9 zZv`DVui|qevLj25emnc@X#eY84&nKu_KibPyj)IK}#`4kzj21Oia85iU{`goBDd7+~5d0nww|W z-Yb0!Fu1~~-thwp50nnrKU*etu!@j3K<^5GRORd2(AEZ~*2gELAhsXw?4;jvNK|11 z@@LfKSx{SA<^rhwo!{2M;UUEOX}YizD9DN=A7TF|bE8WXtTZG5x?A?cFe1R=x50_o zYCEljLIzX>R4C_>!a?@-JVbb5nfQ`J*I5j1Pryi+#p$8R{_O7;NOWF41E&uy zY8UMVw4~~A6;#qiK7I_UBqYsor{?a$iVAS=%UinQbE~T6n!M~7;wRp^nfdta9~~V5 zBj7b_o1UF*1LhzoCq0336dW$&8exW&6#%-Q&Z{h z%?*rwD zvVCUE2lL-pSXijSMk^gG_Sgcnqpx4S1e~4JE3(pnI}atB`q_nXmawq!#%!xZYQF+* z&Ez^fcwTH4W)xNiP_AIw01SHa`}cQEO%~KkKYv<)5s=ys$Q-QjLbc33+d8tVS6Q$2 zJ|VphPRm+>s4@pRgH)h5D061*oU#&&co?^?*S>Y=>IqIjR!VphzwPe`{FbfKZGAK6 zE$uCU#Cc1fUDI#fy7lzw(_-~$ehouIB?>Z#rNPgiXTd~oEyUVNef`>0OY7J2GBqP3 zNLEJ|XIo*=_2C_Lf&^t2QXOl@b!}5;CBeTpPHQH4EhWy9qo3Y>EZM!NK2sO*MZo^ue^f| z_`}^@0Quy-`Ae|vCMSibDlrg5rCfg>{+XnI*0nS=-7!2?1f|yu6X!&j7Kq~h-X6#j z)wMRgi8sK17#*={Q<0a4U2`22{dgz|u{6Sc0|V-V4x{@&b_w*K^-<_c!%`)}VGtE;QO{K7!AJ-c=0GnqU`ifo=6w0j_d z-W>dwqyWu{&INmc?EuxpikgiC7nrYv_j0zek)Gbc$-Wa2Q7j}V&9|ly)|S+^=H}rf zEJ`<}q-Fr=`Fv1Ekz@_EdoJSd`^d_6e5Nttc;VaGpfgbCW*>Ob7-AC^;ap zlxngi0dod}5|md>IXSteB}jE3^OroJ;)Ie1pc$Uhcf8o&*SFw#JE}AK*)wln-=%GT zNG#FO(W7O!5kGd1N7Er%;Ol^{`%-B5Oj|%B1FY%d?ELra?A|l#ha^{{g+x~?{D*dox zE;|_l;&x44U0!7jx>$pcGv6Jn{rTJZ@|Qep0!lG)aY{f{FJ4Fn?tK~_7Lbw(fwdZy zG=R!;3I!@lCd{}nay zfYo$#MjJimL1Tp+JJ%AhyE?{=($Li%gMOW|%1Q<@vfUZKxs923RLp4a<zks;3J}3>Yu071?F{e?IeyoW%*B9h9q4td zOG_elpEE#kgFY`91XPr$^lKCd+3Q+t6dB+0i+ObG{Kn{SzuFW7KS|HTT=%b+8ucF4 z45#CV+eSDszJlJfa~5(^>rHTfpb*j?zfiRu&@G@J94xv2ybYQ&=pn;_;N z$7nJDq3=Lcg8r{C@2Pm`$O1LctbP#$M`*_BBS^S(hY1Adgvrty)6Q73W^9s+STZfAh&@qy1BW@lJE)O z)@sT-&&B&79l_9E4Soia^tGVBSToPjPZ;~sr5|GFz^5QrfW%D6Z0+j$`o#-*A~x7U z!vHv;H{}iZG%(pFZNX_YRS@?z0|NtKq6$(!gB~oX$behMmcL^mZqqH ze%z=ua@pN+fFnp9r6(xwqT*;B4#lQDIc0PQ4v37^)zC?{1k$QAnDx-E6$<4RkRs?# zJNPrOPIT1NV&)y;P_j9GV^&}#^_B|?r`Oea1qc6q<1lhJj_$VoFcV10EiEmOKfw{x z;O1bYa&dA}rcej3SdxLCoxP{Gx1_Q02vQDA4@-g~2RR!Xo7S4YzyBY9NpvR>zEUDg z6D~7fYFWnkPsp#JTpTLlad+o3a|V_@(N3N@exQX~)M8g#({PylWWC9v;A zeGEuJlX0adZo3MbvYMJ+>nuLP-Z~jWU!I@Q%8~Be+00n(Ul3J55Us?el}TCCu%wUc zlr@S|h4`mf0JsOqj#ImE;FYL)eF#W`&`W~0I(LsDETQc&Z%Qfgg(53%j(UdX2E-SL znWel?v0>loR~k>YhqVKb0(JmR0ERqz8}}$Z#po6%pdlnjg9(L0L~{$QPPFqixegrz zoLN^YD$X`9-6uh;#~-zPagYW8(wf zDfCrBu+9w(!0vX2JW_84TJZk6 z4d(X&FeEFoe);kR7#OfF3=E9XWNCod0N^hEbYAh@p&=42ec=0om)uCf;R3fO5| zv8B-cxV|n4?TD|Xg24AUZf&(U=JP}Mnnn@R0{{XcAti+}*&Q?zh!iwdi;9ZU;Mm*U zJ%yGk9f@8?gLl3az~DeIh2evN55*4z2zajUoG1(zAU42wtFEUMFs9dH9`LQ8fPkj< z_ApTGXWS$Z*~}zkAASMN@}MZPGTXRA1`3~S>q9^y)Dphs#(Y&@8Xd<WKkOpZ&~Osc0tpf{VRQ?ZYXGRsqoIJ+f&m(U z3kJMK;}L)dn45~yQmsl!$o`Oppf2;8G>?K)3_ZqZZySpG1t4uma4^^FQ}t!tOx~c2 zL+caQngYmA{i*!WQW*=)MECBoi-~>v^2Oulhe+tB%hQ8W)D7SB^zv#l=F_DD?jI*G-|+dOq5l3O;ooy)|RV6LL1uuc@@`HF1hUn za;YTOn{L2$(L0M_j8(zjL=(dz#?4v1q7Epj51p?uuo(^rjldUa1cmI^FH11dflJ&< zSj9&SXG8MQm2Q44`10@8#4R-+Y+0)aOoCxIuobaz2(u-y5!u_HJZE2+&c>5`y>^zD zZgLO3S1*1TqCjY+NzSc(bJWPsnQZS(n?Gqs4O*Q&$$v7VWIoc=b5&o#ko0?D=H@>c zI6Fu7p9s7C!bE~KfUIX&_{WfOl%OFf4gb58fSWZ_iND?R&GE0{+BKs<^d*e?lslV{ z>p-G{0omI^tLI7_Lbjh-DTcp2fZXEc?8h6<1QIZ=;c6p=hl`7gx)60nM)&T*N!%H=m#|9fXy zio?(PJ<>g6O@FVgx%&{IhotzQfT~eCkL23Syhq*rs^O`PGf2`1nk)ZQaT-)PQ`z2M zw<^Nrhj%H!Nf8=UG$hLJlMreWsQC-?0Qy9Ut65B8S4 zmsMc`)A9d}J(yWbi&v~Bngs7QpH#h6D>}81Vg7p*KQb{wNW%I!4QkBni!9Q|U+8he z*)IK4yUlVxY&{c`CZ8By8k!9+@qkyB9r^FC4e8OJw*2~e10{4cS7p1SioTSD_3w^D z*MnZ^d2WizCpPDWYSgXCqsPJW-#Bn`OOfuLLYx#{?})y7Gosfxo#VetTw2Uj8Phye z3pFlNE3sz({ubYVmr84CFOpiKt(I4J?aDh}@}Psp``;w~NW#sUD+vsxvZ?`9TJX*O zzj6HMn8YU)_N2~pr(eME;ijhJ5EsX6+Fg#X0`U%rlU zu;q1T9(}*nTp?#R^nPXdx0Zh_SPcU`s@6AX3m;x{o8m%s7U3YWE}iSZHV@Y5dZe6N zEQ?;SdWPTj=rFYX&obE3KjV7y;ezmk!1)_p9xHD2voiiY>xUGW8U?>YG-x6K3sM!OXb?P|N}}%}e(>eQRuVmo(Eo;Ec1hFIH>zcsw|^kzcPpI*^~}Tn z?;0bdH!61dAL-2}b)x3ZKIn(>{5Mk<(JYu;S<@ZVD(3_Bgi>;Ov`Qr}zmLh4T~MRv z-*PPM;TWhh9UZbw#-L}T`6RlyZa^&t5LHrIY6R37YN7Q_!`3Y@G-XY2BmYSnSuPeu zqM@WS4hZ|+)%Ffg(9vsak+_Br9vrSs@L#`<#Y}$fT3SjgoK2^#+H(jg7XwWm2h)jydam5kN{E1 z5gS^f78Vx3JrDKP2E2YX8g*g>h^zPSXh#XfLV_4uQpHj3b^E6RsAd}TP96^t!VD}# za9inc7;4DXibo`D9Op4Hx|+OkC(<&9JD7lJ&dWOpZNcj5SD1B-j7oF1FWzE6_b5>5 zU~dU%c9@(rJl+7=1EiIs=xd=dmS+hc()(U^lBAHP_dXNyiBNDH%8sq1JYRxT>IVwJVHSXdUU0h_1D~o9CTu-e zCGfptNYcS_0N@4!MZ+?w2_L7Y(`2Y)pd0FAHgWefUQke&&dPBu zq@UegYrX=CPbQ4!1OO#AT{!LRU)J-rURh{in)u84AN9IaeGO9%f9L)Aj&)FG7iu?9e6K2y(0uDV^NhOUvIULOU&nQGtTZ%T@(f zz5fXRJV>41FZ_Gxk7UY`U_t!gM#JWexk$*V1M!Z&^k`I4hJ&r=l;3}Mcj_l;9`L1? zG2iD;RymAh`G_@Nnmkwa&u>h{vigZ(A4*P9(A5*7R3e`@K{7boi7xz3`gPFkR{At1 zhAgJ*O+H3!Io6f$6OXoiZ?3J167mzMsAY3Ie0+QHse-?EqFpCn7nzAIyDQ7J5FsPV z-S>2KN^C7uG#)Y(A{i6i4hi|YwXo0@7~oAIAyOo0Hyd9I2=0*Vbz#LmjLoYNY0P35 z@?!B)rn>vUh3DPVX*z{BtvrnW`k{Qu#VI57qzH|oTBoS|m8xsv!O%-Vjm0K*DLixw zzZ_I9Vn*rd4~#gf{Nt~^r+T-Bxh}Ytwz|&{w127grnksRM z+?c9=0OcP#t{e!9{-rnI)a~nwklsM#+^5AG)#H@Y|N11USKj^oySW12=5H#OgB>4C zKbLd6g+lAz(341r^JtSP+?IC*r9n@{BG!aw@$x()c2&(}C2l06`$WgTu2^tkmEzg_rAf>d;U9&M*4RX;_-cPreUvYj@+g`z~foM)R) z5Vy3^kh!QuEY-%fOw>XUP!qzZWUV`CXUDU~puDJP!rnnv*>d2Hj9&*?_suLKyU(5- z!iLL+(zmZGLar7$EkWKkCqVp(?sGk3R(oP@zqT?o`=n+Jk`AJnV5Fz#GBY^?1BfZ< znVA*FOq4Nsxw&wx1G?RpV=fy-%-%I{yWZW#@T_d*W3TluL+JpPXt4x&u~kco>`on9P6=43!kb z(jsm;LQl7Z#Qt9B#!4e+0S5-DLPu27VK-?d$K!5`rLu@+WrR^+m_?F*4Nu{67~A5_ zR=gPDz`!r`Mb{pYVP@Z&H9f`_-aaf7qrqs?`%I}gVqdhzzuiD8h9uT($bqDXEo zWtaFK5e7{wy_Md&joRhRzBNY{k(-w5hd9h^v7;~3p5-=omC7$L*tOroSWU>>}$3)nYt7jpq z@st3x3esMa+~?n{z<<~C;~EXGRJhaXFn!VFfA1*?N9Djjc=hrn^dL{2>jN-<#D)m} zIzuAyH7Y?O&8JW%FAI+?F>izubxPmVse9Z&&7wc*{mux(o+&xLQa5s@b5smy6%r(% z-JQw~anNz0=LXdCX>pBKrs99S00ElP>#zUb#bLg&9EU5LN5ZaI z8|qWSbwGsxPlKUN25p+^REU#OF~h&MsNsT65E(>8#*))e^tNF_r+NMjojVfoqWhiv%aoL*&CNIx z%eqN8U!^k-+4Na#EO_zXv*Dya7M7eT<+m|6YB%<}gMIVHKh-RPQ>Rbu!%Sx9-HyL} z^uuaE*ejv23=rt6p2mR=VG@z(;qz}-cEB+s@|doFM-`yCD{XydM8$wdPe}>SWd(7R z+9qx^%@9xD`HZ?NV-2-dt7*V|JTDZvf<5_${V98Z>XEAfLYhD+C1ZxOf5E{J7frrcz;Lx z93ckb^kYwQuTTkNbs~ina{{E|2th5UI(fXO5@xsk!@?+?PuNAFx-TeeD>A)%c(Z|k zWHbB#m=0DyH2;A00%|3bQf+RoCCnxK7DVJUQMcD+W!y8&7)Q>T`!$bRsR8oQWyq!BHVi|odB&N#}+T0$> zK2J!f0sjK(g|&?hyacpW_Km*yteN5-O%g1%`czg;IFu_57>Bh^A}%wX977>GNxcBF zF@bLs#TJAlk#5)a-Dm2yq~-Y=i`I}bD2e$1m+K(KQxj< zzYpY6=V5xR;g{vNYKSsxd@7If&@{wD=uzWEO3Q`bdu?qEtDh3v77fg@SPp@~$FbKz zd8ye&m$@77;e^~2oDV{gJhlKI=ep9@i3zGw?qQan)H?#5e1zr@!u0E1%jNOOio*~6 zKkOH_?5ZBWpO7y>Z9d5s)_g%pne|PDYXPKCNAD;+K=YdJ=ubSl*9h# z-stdO-=lubfn~t~h0kJUL<)GHBH!b~`YFXzT$%k~;6ME1pKeLO%%knMT>hU7hS*=1 z@+19qblHvfQg-CZr_O^z>z2%cX36J29paSmyjO?>znRT0)lQ3Q-2rZqeV$;w5;^mA z-5M>Ioj6{yGeH>9J?mejp{3QH;~;+uXAQu~@$vClqyav-Hd)8c`xB0~9tCCz7Ye_7 zDw(YE`aav&b9w@qj!UmB;^S+&Zx7%Q*^}=kKKL6Mj^zS&9zyI0MSuN{ik@K*^t82XZm-^19B&Co0ge**{gr9VmYblbP zx6_FIsY0Oqwxs|KY?^aGS+taJqnD!K$5vG%yUqNGzzhG7hwY=Tc&Dv>fJ(N6S{as_b7uoUn z*m?Kne47=YLp*J9-nXkeMH-RE3X?~O$S{!-`{s++;e}5~@M7X!m_F)e#D%}~%r9Xh6y)m_aA*)C3n{}UGxG<2Zr1v3*e`&q5)$lR?IT{AGdt1ee={4|d3 zs?>SP+jhRW)2AF)fhXM}s?8Yv%6a+6 zQFbr z#}DW&sd!pz*$`Xz`6VMdU)xJ=Exo@Q(*@ex_TizZfB**<*L@WgZ+CY%;hU41s(8UAkZtj={Z^Rqi*4S#Y&1WsWMVh0 zyl#G~ZQOAj#f1`#_K3uV<=O{5q>oorO}iVa(S1wVHLcv$;6Fk$gCjrmZ^5bQSX#XB z6s4p{+o}1X$3kR1f?=ECv|bv>YfUUF79N@+s}FzqQf`OA&qFdNw>i1|=P+|JeO6A< z5R;XpHb?HN!+lNp#Nr|*e@a=NVvfLZ?fBhI$)QtD21*i)8cdc9JP{HNhgZYQw3y7Y zB9WKAu8pT19`^7l#HCyN`HIdBncaN1DTyh4yxB6;VntTKM=Bm&(NKFZBKm~mC_Zad z$o|M3fwc)XW>X6n!fl-W`MHyrnTYCZsfR;f-BHh-AH8hS zeDc61#3I$5yzUP}Tv}|s@ eNJ-*joY}y3^!;wpVXXa(Vl;(=e;|9AK6s$aNpZ;n zeuk0XY=9l~T7Zz7&r90D-7A(@^0O_oBULr&LLylmKR4NBp+e4bPPidq&80u)BW2m> z4WxvC-rF5|C!1)Vn<5JjQ*omz+!4(jJtipBb7<1qP1ZLhRY(vl)ZpNuRCab=hsK?( zM3%`0|Hc>+@z0N3X})G!^v6FfQ-3n3hs51ZO5)s}9*~}XqA2;MlB}_i`PSJ{z{q3% z+l=i`>?r#_G)tU%ySGXfGly$8Wiw}7F+e?s*?=D{ktybZGv|a3$Fvbe4(pNr(})AzUNg4RAW9f^_7j{vxP0x*XpQ? z!c~cnCBGvMwi_QK)PAjoFS?RcXU|K%&0=t}&LWS_qNsbScNaIo!`hNCGwqXV-28%p zUvubZSKRT^ud9Apl&UUw-pZ`|&1Zg=AoWK57P+)+T>1Sn0b*Je$rNHm?bdkj;EcCc z=A#!K^IUSRgeaaMcSTqq8L?>NmVO?eqO((_$qCh$(Ip{gx~~+7=Uy02U}f)V!4H)| zcP!g|raSHQLUeTb-)HOk>ejuzII>9a)!#)x%8J$PcTDfxx36B#jEtnO1PJ(VYocGO zA^pCBGxE=iqWnr6#B_vJv;aH%@A?T-iGZDjOuyNZFZMJvG(g_9jE&LmIOxa1PUFzQ zde-p}3tC<(d5D!!!ydD5S*pv@L$p{cCHX{*fN9^myZQ9DEt_Xv9>?RKRUH9DIscEU zuYjs@-J;zfNQs0J5+X>0pn$Z1q|(w2Qqq#rozjv50wN;a9g+f)(xHGf(gG5aZ*4qx zy!Xy=N`|f}LU#y&SuDKLEt*YXcD`I!kG{0g5@h{vQ&}=BC7CEW8L(Q8OGBz}1 z$VCR62;fj|N|71?*yj`qxZE#iPC>&8c)owa*c)B`JQFOqDPq{xw))%lJXOWZ`gP4| z&XFbE$;R8G$p@v*h}1V0|Lfsn zd{QGwb^@Y-g*}(Vhp@@X7WVYS)VGh?=zLwq^oZ}4loz9svo*0{V_AO>(VzXdfzHGk zrB$!mh2QqA!8v)Nj<6+=Y9_mH&>T$V1bYC`4-XBE`6P*` zA&ZQ8{lIZwqtUSLRCx8&2VcvYW_ztYuCi)wvZ_5W()Zpk)?SNTU4DJD0#*9!*xm2w zocIv4uP8V*yib1olnf_*^Y$$-=trDHp8)L<2q`nJE2+S}KR$NzZ~5Jd!~R5@1)-yw z6VofVyiD)lZaJ|_vy*RV=lt-8?mfX+hSLlTo^m|ryQXsjVrO~i31voQB{>5a`t%fI zDTrCa6&|ujfT<^c*??rK`?ugzP)q?n4Oz*5$NQ{x{=Ppt<#@Xkn=GoPI*peStqoth zxFja^qq9~2*HiwlxvIkYTAdBW%w|rV7LJ!xawA>bdQ#((!s&ya0*(T@ADCT@$P z@yos#bu6Qd^79X4pPWvww8RhKH}^Y!JK`8lx)_&OCu0=XG-%x z*5MP~>xi0%okRsB`W|+GoQH%Xx{3^88(~1_`hlu?`}*X$fUyF-1GRKnb*G9^kKUQ0S>Z7D=W*h)mQXaIf`m^#dx2$IM=E+*k(WqT%u(0QSU-$cOEKZ|oF{wY= zopB7$>Arg@@~c-pSM>5kZGr!QY9DjF1u$6yq5(&U^Yb|u8C?zzLgz_WZ0@64_%lRY zHBSl`EX~Sv?+~kan_;1D_}1$w%16AaB)4&Fe$;2V_VC0d1_~X$&}4jthUSwm&?M(i z-l_}srp7#nQe4t27in@W6+-EK@kRaOjosL7Olf#{dzIPED(k5=dZN-bkFJ%!_1|)9 zeYv)x-3`$@O#}M|8yama}>Z2u-vYFtW2Z{Vw%~vc+{QtT!LHbz1!7eb@_9 zX__MMu@g2i5?jg!;D1Qea_vR8Ie6kg77mHlBid>35j6T58pZaoM zm44>gR zcM%}_2uQ+IEfjQt^c%VZ>NojqbwAZxrktlb{rp^p1qUIw4>^ysRPu%Lo9e4V` zFvZwYer|Skxb<%1rgOm@oy`eS2Bt=RLro#~BJvu=P#GTGL0Z*Z3KUCVe?e)TF5nt4 z^0Sm4?T;w)r`iBB1dy#&+d|j~Th*?XKbtCz9?2flE?)K&M5LPTN)6Wja^>9zE009& zVt=f#cC0nJ{cv0Rh>V_E60w3KtPt;0DK;0PyiYSlGfF0{fS`79=zVe6#gIN_+lK3X zwyuvvM1qxlHiV7R-1e2vUJB3Vtgmdyyv+Oh!i_8D5-jcuK5vgIH#@ei(XR|Wm^&GEnMWsuDAh7uYs z=S|zq@>{cm0}9g#a7dN`W8Hf@(Mlyvg_P5!_K%eQCuZFc%c{kjFG=2U(ne*w>xL-pyOhbRVA z?zCTi-jCp!iLOig*+K+P&g1$19RY>%C5icQhA*Zp)B#65U2$nsz2aLzU*&7wJ9PyO zRJkipCtgUoYEHcKW_HhVJY}6-zI1Z}wLm;S{`;42UlxH-evkd25hjll^1MnX8&ljK6T=@HL67j9%Niirmst zxi)A96`w|emTcG-i{q_-;l*~Wa}XGP^>lpldja_KEO%k^%8zaD)(^6Hm82=vha@6#$rOX$I8-6lV>atQt_di*@{p9_$C)=Rj#{Uz>aOyofq?(&g zer$hfG3znILPG(F{OZ{RGaK9aPk;4?Kgsu|$zm0I_xKSCBQ7IY1O!^%>Y7zp05?8PySI#UY$ z2fF7rWa)&W`Etg-dLu=ZS*AY{ZemB|UlJx0CY@(DO3%ZfPqi z0EIamY&zZT0TCpjlMPPMuRw%k^H`LJUa^EkyQ*-DY5C{U+Wg{w$`NO?fIf5*9> z(g)gm{P|a38gZv!Q6IH_cREtuHA8jf^4%z;OD7)wEY%%E)gRV1O0ZC>CLuB-S3G5$ zAfZ@3V&caIH@e6l>oHrReoe+G9Um%StUrTX(S2a(;OMh3G4X*m5&Wy~c3V(}kxk*= z-1{aB0wnYSipTZlZXFYq*J^KwLeAHZ&YIFU7AcK*-$Ni_`SkhEHrGA8Z9)r{%}yom#t_KRca}9yd_Y0GYa=Kz?SBVJgER8<6n{R@GB*Cm7z&{#^yaYe z@X~pA>JE9?7#ZJknGXRw91!fZCvV_30ouXgCH%cE;l-04<1(HZm&|ES4dO6~+%K8h z_56?V@7oVb+%F$=G%kLYZpq&uV`W`NwRnm;+Lg;}oy=76?TTE>vT`IMjKu6v;9KzPWI5O<8A3&1sBE`YNGPkyVn9>J;V-UgtB+>Z( zYakZ_mBD?lk7xVy;fwhVJS<&e;Qv2Y@Hrf2~8)TMp z=&l_eskKLU%K6g{>2RXdvOnfmtajgi8nEyc9XDfCk~=PJU1rib%P^vAhx$3Qy*6uN zZbtPN+eVJp=Q7QIk!E=G#9;GhlQ+cVSN$JnDFI2t3{rbv-*HbA-TYB*E-s{ZKhTe) z5%=-+W$CqJt2{4*$}e#ZEpnJ_72eX!>A$s;vugS7PUL8P)-%*b~x zC+Rv?n&l6N7ca%|Nhe*Cmq@>CD7`8e9K=!;yfvO`YEcu;%a)h=8DokBu8r%WWF~N^ zJhiod)arwEcL0Ltf%i;>d8}vVAaE9pGN0Mv{!uQb!CqeYE9sf`G~H4_Ku3Z|S~Z=T z**cG(Zxj35=URI)-l66)T-I!r-bZcX?pm3~W~6;X;P>)yJ9aY?MUZhXdI4nOAfE!# z*OO65?ON~wuR}xKcDopU*b(HH|8wnbD1#%WY*R)fE6Lc-udVML_&0c$F!QNq$tmyO zl3h}7BL783-g@hTv9gh1>BkD{kL_Z6;0b2D9^ZU+v348@&KHWNZ7! zAm4&73OJM0v@|4LoQn+ftkM$`og7B<3s zu;qTPfA3#&BNt7HIMtzECY^9LnZK2}U=XrY+g2Vx>0F;-Y(m)_eY6N_k zcMpt$Kw1fi2pG6zKv@d$K}X29qopj>%;bzG{i&_u@4OT$i$>d#naHIA09tY9GQ954S(2!-L`61;Lz_4 zh~fP`Ie}L7)VWU@lAv8$Y6qcQj2b3FuT@LWYGj3MpWpJYT_<&<`oUK5*r(iz)v}U& z;prCKToc;6kv$=IX%?JS%c1IOWDra~q2AeE9C-tsVakkZbZ9^$+ycdsuIUHH`iCDZ z-90^)78X~pQ!qZWM%F%2Y{+7WZ&*Pd1Q)sVjzW&Gk?@mv3)P<)7S_qkuDmM?A1n+SrDKfQqB%kp~R_fp& z!ytl^vhs33*n%LahKh<~m-CceCna=#(Pfh2Zq zLAD=~mPj=sA5XnXo)qtFd7QzeEXKJ|^Z2=qH+vQa@&UaD@I$3C(qFNa7EU_4R6i@` zB-(mz+IhqH=>2^s`(#*=r29U*9Of(y6wWr7PNutJNpRiCyW5w`1@JnKhpyp*>4c3y zm#pN+LYfyggUZTE5dA2JrgmZBJz830P#M)K9Lr^BB>uPGsjX_-=0k8uYnK;}_u5nk zKZErOeLJ(&LTu>Nu|QcS;FH_|<9;)-SCW>`OV<6n?mcF5bVVLL&-F z2*B%C%T*z1w1p*nmMIA+-z7l9ys|)7_xbm@h%aOGQ0S8MEL#0K154RLU~#fw>6TTq zX*>FG*Dq+mz%Jb4Oog61m2x!nAPWqbOxyHLiTEQ+uZPGmJk+P%P6wh{=$~9P9)O!; z*C=^kR1_)w2OaBjA^&t+J_eWa%?diJrK<87?3~^7Ofyz1ua<_<3@c-Wa3fL(!s+(-qzWuZDxV{8sO+ia>p?p`Xn#SH#xJXMn zIS3|Xt0TtnyfKise*oVY;`l`OzH=ocR2>Be!B0@G)hlrl0Y>F<(SXWJ;je#>pqrFi zxeC|hU0$lbsR;9d6vg%H3oCf^S@*0wnR>WGnpDg5Stg6)TRCgR$wB9cj5zYyJw##* z)6GOr5r3sz+JK6xs=AWfh+oErma(OGzI_}|irg?Fjy-$NN~;tfQ#AI8$vq4I`p@Oh zkKa*wHX#3un^6nLhfQ8RFiOfPesSz}(1G%xJJ|S(ohF%$m2+Mz6|c`;SAm)(2w~Id zTUuGs2Gs1ye?-{wRjpwgHTJ_3WUwjQRO6TpZ>Tb@t}Uv#Ix6JipuNCSiN35^zE)lI z&WHO#rZo|BcRuFshvJRX9JRM=Yw1@Jn|Duuuqi4k`jPDlg8zz9LX0RZ^vKG0`NYYn zXGuCqyiYEF^qz0V5HV&E?aJ39Gk&3D{kQVhNB9riCe+);;zQP@F%<#{-f)Svf%yv% zDWF1i2@bX*{g}Y{arBdiT_Fn`^S=Vjvabx$s`)x=o(}0%o*!4?Nqm+x|M|AXJTO{NkOtDk*PH^j{n{h98QH+ zpU(OPuv7j7@$o_(&}hT>I1(ToNaS%;2nF%>f7colEs6waOfu ze()3(S*sJ@=rZZ(@zPCBtVn-PRP%%&h_3%&Do~E@{t7~*ew9F3pJIoo7us75(5VG{ zC2((m!VkjNP@w$`9RvCQ{HI`PI=0;YCqJD`2x}_RH6iXTE7>}b5Ha!Kuxmp1JPZyR z5W1&UHn)+Zwh8vlmI%2&Tz{6$`a3fvD+NobeK>3Q4N{G`$bc0BD6-G+JRq?nm9JnX zi+tfi%w+>JCtWN58howH>LA#*WFhcw(kX1)#0*ZpX}Gu<{%Dd_s-^o`hE=|lFYQlf z6??uOZNce!qEMqg>Q9Es8*hxVPVwv}BudygIM66u9pPs~nzjDvlU&ARg`02x@WXY; z&#j=!v=8L4E8(Y;i+|UBf+3^*eq2QCmpJ-9eB1P>tZ@K->$4cLEz(hGZ2O-}V$Lry zXPJwxYk*M$L=gySE`tBLb^A7euVz5Ev{+v_Tg?_7c@DNF6_O^+_{+v~zCEOzp{wGa z#IyJEU0dXe=T;qkfvL)vRXA|+o|;!I{!IZ-2w&Fxr_jAjEi7LymPptfavG*>Q_Y@7 z!=Bgp;3;Pa9Bo_Zd@jP2DMNLh{ZYH|E*tKG(BPu$wQNk6OOqD z%-3);HOaukSGL$xd0-AeJv{%@r{x&;2O%Gzy>WiXi_jtDB!?U!+%0wWw++NZUeSkm zYINCPAL3s7=ju>f2O`*FZ;HqI{ld9Y!zo$TsM5=1BCH-#pEp@#G55{=Niib%mK5j1 zsT>;{o;>scyGS2M2dV)H5vB*@o_Il)4W!e0@=<@^7Zr=hz4PsD(=en4lLc^xzVe>X zycZo#R6Ira7O{xRneNY-?o(THpRZs)=xlCexe>)vitbh=!XCuj^-%*}gQ3IL%Rbi< z{TFSdkJC@j!#>crM+$Yq90JeRFuW=jPw=-0bRJ>S0jqf+Vyj2XZYt*+=c*F~rJc4! zolNc)st0pfj$+;Oxc2mD>RL174W${az z<&RURvH~Hx8gKDp=29WwuZ2vf_6Yw5#XbEQIU{{El`s1(w(;2_V!!W5l23Ysnvj3l z&zd)fJiEkGg(Xrxd_wS$3WN$h ztsnTC6m2=U&=hDn%JP%pn<%abOwAsdXU%?%ICypsH{x^uvwLW)$WDM81e@ca^WVgt zuaXBBg5a(sjCp}V;dpM=DYB(pW4WCzTapxpm$d5y@!)`>^}o!F8t!7=Pupwiu3YdS z^{GYD(eHN0P!YIulA}gJN$Mu|6bA)oLx-d!0i-JoS^#?f0fZT0J{K@XemT}btXsd` zSPuZ6aeZ$}&_05sI!-!Khzb`h=~mvE3Rl%Tsv?fc?|8eu8F`B;P|@LdekXpe!|tU| z(CkAC*~Jr%&K$KcBjuYPb{xVjB3^&k?TNCWHt;>VP-C^8PqNhAcHH<)H*%}({1#@>{C!rc&43Ka32vHMmTGExpb$=;2|9oClHABOBu+cC~jtJk@EBp{NjFx!& z&vV`Z5u;FzMJ@udNmKPO<@&Y_U(K{uMS2otwcFm~=Sa7`eNI}0`dbJ<4TASlvzOL- zVb8$gz$L-3nDo^;IQDwbA(fW#y!y}pQay=a=TX!e)J{Sqn-NObn+pq^Cf#w6?*Ysj z@*+4iD{1aqh1D~8`T1e5U;kapFSY?$ctEkxiJD<(TW3LC4oYTG>#rFn@VYbICrMC< z9(;bFGx{)t+ELTg_f*N8YSfUZFcSR>f-i}*$He?s;RBr#nU_CNt7nHiOF)t6Pd7_3 zJO?Ux^?x zHNTFZ3(X*4gOTT4h3=v#j-tTfbR*Wz*pVmsRwsmc4Z@9{PDDjG8qg*NxBzHT1Vn>A ziN3x*mcu9?riQIauf+o@0WLGbL5M{q+>91cFV-97xEtB_rGxXlq?uhw`ev%3&_YQQ~4X(z+*sEU1nIl4A>!zgrb;PyKiUK;v6Bt7yOhs zf{yC5f~6E~{<+5dbBF{oo!Hn{9P1_~DQge>(DZJKunl$q91dEhmhB@7^%8Cvx&qz? zndGap^Bgk;vS^5qmXk()cZjdd8}SO?kuBMrt9;-JO z5IYQdfe}@SX5kwBfoCt|mDdQK>-jtE`bmUO`L*kr{%uEp{aC9 zm6bYp5N2egt!c&+&n*XK$Uh@`ygE)1w&*w;1Nc6U@2R=8zA7w1;CZG&&kdd`beH;; zoepg6kJP*YBf(IaYI9pLOby+D!n)SRHow}%-?ekh60&|!PAQa&b&lHp8hqEl7kNth z%`ny8SQXlbTinv35Ek(?Sxg`W&;BxOsg=!(8N22#8Y>#>cKh3`j{^`xt~pf=n?N!_ zuTFnW(Q+6fAEEmd%KiS}ZD4`<>!~SnoHq7LiaT~PE;%|B)g+mLd0W zkp+41&6_WX%OPLpl7`e%US9QJh_GK11#T^JWqF&Zxiq^JObq3jPBLT*oGbVVbLnAa zKLtf@_vaS7N8Wtb>V#qmNdgF#xzGAyfK(dti5unLmY0@10jQz?h=%m@KiMH6nAya3 za9h**U)A+CVwy2A(H<}=-l4w~jgn!AJkzDTtK`pKO$aX4`MdcF9t?a83Jm0{4C#5?y!+&X z85ns6k^6=q9C1T7iIkP+-qD5qh?!#H{+609p@t7{G zo3y<0fh;`*(l;S&KGfrReq@I~^cuWM;njL%f8lW|`J&m-}Y79i*uE&0bow>hY9r9})scL1#_?`QgVqtW>Y) zyST~jv9eBqgr%e7A3%UG#&6&eq4ZGxKkdQwQ`>m`_x%5LheU~N&O%(Bx3XuY;pPf{ zq*keV1F0<3YCS@e&^lWJM3RZ5IJp2*9B$!+x&ie?x?35gV=B=vM%H^fucuKJVgKKb zUK)qa62aL`9Rw~)NcO`3HS$SCHt7ExnixhbtQOG9O4x)gaXkE!K-y1|skOPZ<}(VD zDb&kn_B??DcbZaK8_9oenT5zSrW)PUX#L$PBhzJWBqYw zS@(Y4L)Leu<~6;%$rRyBZI<(OJ`X{F8AJ}6j^@LlISMW5R%bvm@zQ?D5)e!C_gPZu zwc;Y$LUtP4>f5VKPagNPKPi+`H6K>>ysBI~-<>{7(C|K86ykTpye{rs#h_4LJ9jE! z1SA7`C|sP#wENlg{X8c+@c&48GH6Xh_`CiOxA_8}%px59zl@+V6=|V_mr0A}c})}wk$tw?_aJvh4BQ5wLw!)au9$eodVdLpKeE?GKKZ|1 zo6%;bmR&`Cl)9()D+d6(K{sd%0PN+U>zrHFP)DoNvc~xAKFJFCK0maepyvT8vr1kg zq|D&^nODf~So^Zg#_MZjP$J!0R*>QdG7HOA8A&zuI}4%{4GI!?g-pXnV-r{sFul+J zz7P#z3}yfPM+sE10Zqy7R>5gnA{fu8(sa@ASKRkZ-9rQ#ph%ZVe|6l4J zQ~f2T2$B7BdxfJM{`K z7Y6AgDbNH8tr61W2IpI6ARH*OHV@_AI(xDE-@krRhCrf_m!3oNmGM`lyLIz;=mP4! z29w1Gx-DQ-0C&ghzbT8v;bUdB9G*WNbq)a@0mA3Y2uBP-r^%N$_;jx8(MPjt&o3xB zRMo?i*=Wdxy&6Uyio7r?@ES!mtvvcPjYQ%hNZdRzn;PSYxB(4?y@zr4IUmj;q|~6K zz3ilVUjN}eRPH``ea~Q!B)r%33x-;Fg)@TO8zn_zWepn z<7Pb9m01uN*=pG9much%;b)jm2Fk0oHS}h}9b}DqUwb9f9~+G@@jTImX)86kJi$A5 zUXIQWwcj@#5x5R8)uh1*jh4fcYQDKgot{Sdmz<|n6fDcjA{=#z*Jl^R=AZAek zyLxhiH#4GlCeo?3zcW?Pb=ST)sc#w7ihuvMRVPN<9Q)mOVyGPk>1cWjkF+j`47vDR zeB2vs`~&gYGiqXNfm?~Ou>i|Rl4(3xO4KKd=m4eNI(R#bh4aD=F#+Cqh2L?tm%lNM zTJbS)#4IK?`v#%3w(YHaPP_2ZA2#)yAYe-Jxda-{&`^iEEKr=Km2T+>lvd_?Z3YKK zH}?h~y)VgcqypI>C_B}&q*<25NI7FFjnAtZ=c^Ugoov;g{8{C+7|8*Bo=W6YP2^|? z;G@C>1$LNkyM{ftMs{AWVPy^#-QRLes%u*CRc#H23tsTBS8sAz$yuua?`>pl%`gW1#F0(4*-E@6goOL!QKAUV7 zVYh2+BNAk&I=9GLmpt5LJ%000i$A^olP*d{MmJ)wc3u;RDuBd-;j7=Rsz78PT5ce+ zZCEm`9P+8^xR``c6qWuF5*&=!+!H~q28^^apSxjpi75rD83V63og#fK2nz0R{vcWH z_UBq=B_hisvhzHZ`pM>W54N&Jh^*_Rr;CT7-HEOyv*dJRbH0%M?xrWkk69qTZ4T4zt;&QE>h0lhqQBp zOw~YAn%rDapt2aw$=+?&n^jb(2`lF3uJM&+;>5Qp@9#uio}-zdY~xi%z5X|N*TV$l zjWSA7@k<25E6YCqT~Q4fnl<0c(fRihfycoEiC1lD?v__!)9Z9tB=+NJz}#$?<)jyI zsXpD}G}UoAg_;}S=AdLw!D%WBT$JqeY8it9uYK>g27BZWG$)l=oF2l0RGtflj9(qj zjMykXi7xq`xpzfPLLv}k>_ILGwD#@ic3?&$45Dj-NS%{&4~o1bXs5oq)yw*0TGH5m zXG5C9lm7jB*_JmZnXYLu7d<40BL*42>zurr>o%)?=)tLjS~L#c^2!PoC=da{7tqAt z1oy%xt$;`Gke&X$e$pa8+3VO3pHR?*P4$&8g9Hk~=VMUtaCXHzIa+v2wT&&WR2?pV zRSgZ8@o464#hF#%eQ^$Y2}OLBiOAQg?={?|6nY%brLtxLx?nO#Kkyr!5sy?lXZQmG zoIG_Z0+orgLvSvPfCmc#nZ?I1b5-(SV0PvAQVJM0UVTbX@INd7d#tX03&|(` zZ*(K+=HHL#sA-RiI!6hfDJ3pEcS%!cj{lV{(vU?h0^Yo$&HCKf4%~92h^qYdzwgLOn}GO{wJ2=XBU0Jy;drg!y{KD zpCS#OBj1Jml98wR)~%`%1Ww?19&A0pnE*Qs#S%ICVdVeF%fo~L-PcwJGdsP{_ozB1 zKn2x?Bgv+ADM4PiQr+Um=Ps4#p4>CPEiWac@p;dJ9Ua?#*H!c@i+jL4J-a5IV-*}u z_+@GLpSOYzFbJvB%f9`c*mTadJ<|>Y>3;FF%dm{K)p53%zMycB1-LCDRA0Lar`Zdq_lOGw-{cQ{{*{QG z!qKPPngluI*X-e`UN20kY&^IMAG~p`)nIIZ+`M`iT(c$DD8zRDiK(N2-`C^zfr_V$ywK_{`$mW89Nis*+Bcu5gsau!yJ?Kh zI_2KfA;(y;%rR>6hpO)wuM+I6V6)!TyN1%EtHBakeOec(NM`>{kXaqc+Cem{{=xcK ze0}~FHQq}Stpa#mQ+Qp}8&*#GpZm~KF5U{<+&TNmh5;uaS@35%`Aun=x^nu|s1z&pByyKE*4a!!CLRusR5x>ZWVUC+LG< z9;q;*ODOP5XvcPcc;92dIcgsGpXD34EuF?n-nSRm;Ho&K?jI;B3SVZDUubbU$0_>c zcN_y>WjZq6LsTehp?R{tE`JTH?1vlWE7=@FXtbaI(}XzoTXc#_m${NU#90H>0O3CT zoeli4AAg@f<$lP#?}q^!fAim~Xwy3($WRX97?=KzmAqY9^t?=g(drk>Oo5N#C$Et? zBoK=FZM*sH;r!b!`TJD@!1_x<=djGb^J+6B(1$J}=;3yl^*tz%@4p|sG~EN{Gfx*X zJ#)9NX{$_Rkzf;)&<<2>s~erm%WskhniTze_`5?kEdxVN#fH=dwaF9T zM@no6;5t}xQ;{-XE^X?@D_ou|hH)#pD68$dAIOhswSC8HPmnZe75O=ml6838K`<$Rz!Eca%>NSRVx~PdsaH2m>IIY!o#eLywZqw zTW91o2tL3-dO4K!n81}Gc|byG)^~=yf`@_;?T-EgOHg;bE)(8>h^$Hfd46G>)Zrg0_v*`4 ztGP~R{MT$oT~dgh*x{P48S+I}S%dE0`(cmSiSd&&YL_Fo2ZoW^=DK+Edt*W*VHG@K z@TQJT{c*EBn`!x2v5tuNp!bGu=mgs32cw%gwYGwuxaCEBZm7>^y-8!hr@r@%rBL;% zrnqeLy~E3E@599_O&m2#M4_p}Xj8aGh};CeH~S-9%X0cRq{TkmDC!pp9iLf{9WvQ7Lxld%dl(&>8!{o0YP}OCt{<>c+bLk@6h!^lo~amoa>@& zOiArausGVx!TZa1eWK8Ti~#tgX8K8-)!qLhKc+m9$Aa|iiXgPYlm{^*a#OWQZ%U_+ z*`q@If<;_BQ-OlggGHLJU+bfJud(f=+CC-C2i~avX^+CULM(|f377j}OUggG7 zzlZ2Ii~uy#z67MDK^%TacmuIoWQ_L)3HI9OSFphOu~<$o$i~?j0??-S(V%~Ec`R*S zIPx4n(eI~;Nzz7dk-1BZs>YCiaV*?x_2(9s{e0#aR zHI;VLtDabwerq{!JpfZ9_0wF+;g>2?p#2Dmb^O~c42|dxsl2Y@_>D}_m)eZt&1)1L zXW9|Luk4xAXPFfKedq<|vzSAhSuA9h(XI7)jRW}t_}VFCWA|uaIfTAo1pZsy&KiNV z`f`kMW4M2Zc`GY-44nLeF4%~Z*W*7cdknvfzg#!Gy%tO%ji+9$M*271zrj|iNhM)B zS8q;=E``qKD%_CHLR#JCHb>>yYGv@}!~Am<&(Eok2exe>~o$4f)U$ zWH9f1F5;^ByW$0A@Xx34se|F8LGNozSkzJ5_-lPi>ch`N&^Nz;FM7C5&~eHzb@#tu z9&vkn#|ygK72QwH<5qT%HgI#~Pg{1?a27PMFLspwSsHTNI9AH`e3!C@CsT@BlgPVJ z8Jty+^LPA$FcPQsP?I_K-@+w6x%RObcNw0$1>9hGE()f`cm!6}aXb!EPoP$eU|j0b z*fexIFip+g*Ba_QZ|L+#9+D#My-%ZK$_yqM)2KF|6)vq!Lqz!9^zh{+3fE$hDZQm! zBM^hy6K#!j)gT6_nY)R^o0rK~kK~JsY%8wjTB9Cv%0ZRYAYY5OOkIBGH zb1N{ubUN$-N0SIx?)4_UCzx2hZ<%l+F8waz4v>|rEo&rkz@FC&3uBD(k!NBm-50C(ODs+KWc*;sFSK_ znNXM9Bi39;TK*RKL{gl$M^c)`1(&9#dAoP32dK0!!ogRbuz_RHE!e4Ff_|H;BYPA1 zcHHH)uOX!~?s_t#Am7MJT#jo~rx@b1FOU!e{%HJ0BwiBt|NN1U%&$3x3+@1hfZj=J zZjFyc%S!kq{YGjs)7{m;-wscwiCX=4@ohC7Vr&0kA{|mm4<+(4^jm)aowezh)S&(D zqHp-}fFOzjGNJaQ*A*KT6{0O)J6#0}cYQ4M#=WKjy11(fo}t}~mqjLo`kR~L`nSGv zike#AKn8?+?wA@NcPE}_3>t8iX%y zboCbcjpDv?-fl$+D-9|p;V2(`@jGHIjvOgHfk5yt>bLz(JO3kcPy2nM&~vA0pV8!M zpyso$kmG~Hxx6=*04EtqxQ62g*nzl&>^ENx`p|KWJ~U;g17G-;6p%@m#(Dq?uug{z zwx+&kXSNeq2DJ z#E`v<)p+vn7W(wx{6Gtj42WLY4T$ZF(d}5XN0G6zw$Q15i%;ls(&nM~_VRd&F!Y2h zfa-;oHQd53V(faIW($q#CuT*^Km0BCS>n=UHZwsTtxeT>8 zwg+oDcz@tM*#!4gPf1$^wfs@T(kJk> zy)*b9i%>s8_HD8g2RRE1Z&{AwlvhyQd+HiSBtrjeiQO9s9s~k@gGo;Uj5@5&%j=$G z*d%&yb$)TM^(ZIzt`6d8IZ$MTmP?QZGv~oP;&B?PpN}X}w4G1hP1n`?hFPkm@=nM} zoRr#?{2Bt4RYO`Ieb`5|P?mC)@b3kp&Kt~zMiXN0Vn!R>;X91_>_9&*Rfx?;*7Q@q z-o1bFZ2o>Hd&(pCQ2D!L5r@yoB22gj+aMZ5%tkC(N&eUp9VayFe8lNzE;!h>%G6Q^ zdH``u%9onjplD}UTJ!$Bfno57MRlzPVRC7g7eTBkzuX4M-9w3R!E>wIc7E_&s%u04 zXh%FMD_i2HQtk-M2`tb8>7cA361|7~{B@>%Wb}I%eibDp_^;2``I+|)nm*UYs;~}w zB$;m1Z|hN;8j#WwK80y@zklB!^+hnl1kpgz9m*-}8CNLJ!9XZxBBGb_C@+KVC%5aG zerZ0ftdCs63~5wf7iT2~F#YV$Z-dHyS%YT3_QL?E%rsX@c=Dh9^I zLWpDEFTBDw&KKHI2cw^@Dw?BzyEMR%Nazj3u8APxJ32a`B7VjDs!!(vQ@r$%j7w+x zd&8$PmYr57-WA11PngobdC&2=Ntd9Z`s%PTNL8ok@$bW~B7V(1u4E)Fb zYYx&Htd!0I9Q9C3i?0p%)Z9?N1dD&9DWzFqjbwMN}0jv_G^kH5UZ@VYU@~ z?gRvpVM-K0V}tOauI2=@;@~adkYQk8fNI7Qr;i4Pt}ZSEW)%RDgr5TJY0#v^qaK-@ zTwPo1UfIL)gK3CZFb%Y#;v>J(eP5X9tpqdTVX`I6od*~kuz?xq>1hH0@(p7u885K` zQj(H-HFP{KqiV)`p@t9LTaa1F!xm1UmtQFlOdZD<;@sph19WtlV=I1JOCIS*&jjSTIKO8Xn$)KDVr_ zEG%Ki(z1aMmfDLK41Ajk*~4%*&;PDaS%{$?tLy8hgH}#*a%euxScT&Q|L~8pJ<@UY zj|eLbA!=O2Q~QWVdjK75dW0w4?Y_)b6Yci}=Av9PV9#LgCxCn~sRKoyX=(AZ@nE!hsLI0$ad2=jFff330LUYN z_HrUIe?5hfzF;U}9GdMY%sd?cWxdz5G&D%T)oWK5VAcf)6$`<{;7U6`EGOV?f-Fy6 zZ7mEJ)@^XH0}UNe9T3BWL1r*k61eDgwNvnU7?az4QE;_5Bco+{T9-L~_pis4lD0NE z3CYU(dKf;?F<}Us5R72YdfQlIfFT#CzO&Gt2kQtn z6Eg$k_FX|=04P$h@i0@jzpKmj3~`OKWyqouCa$%RN6eqSY!;X=ZaSyoxha8ZSF3}E zhZiB~03(56a(toMI_T9cElE5N5W!PnRt1#+(8r#%0%K5ZL;UdmeOnN0fRT03Q*Nf$ z*!GKzybf)bdUBXT`wIr9!8A#!%t4S3o}QXwyoM+H4u15^j9%H=)>i+*4!-FB>FLa) znmp4e9&NyoK^sIVSR{g2aH+Vo7BEeQ7#5){q8Nsz1mhCK41$CrMg~YiGaL!%REDaQ zXakdqON{-h+`03&XKi3s926Dy9Q9_et2g13(g72`Ih&2?w9wu_a~FdFv8W@ z+S-coRWX-_1mHepP~+{}%{DK!jq4oO5N(7+b$PPn{clMD|LZKv2p|@T|CEzzYlkncq9QAzu`ahSwIL;R5e>?jj(t%0`hYQ~P^pNU zBna~$rZ69Nc1}NO_#jj*@k81212C2u%m6zRq59Af17}G&`5)jN6I`*I=I85MqtmsS z&EadD5EOtLlc;TKVsr3Dd>ldPz;N0&eCCr+jphFMe&=g#MaAzha3J;pQE8qZOO70Y zvomiS4k9nX7FpA@E< zyHJ%>tWY$e2LlBQS;~Ajlh*BsE)PCA7Aw4CXq55R@Qe4J2VdHG5CdCKqk7@D85x@L z@{=VcmO-%KuTFbU;q$PjRh{0bD(JR`M-^ii$NB|3{$#T{Mu?g z@$aAJlI-W>3Tgieb^yhF)0vy#pCX}{%jH6(j6Oq7wx4U^3oF4jq*Q1$mddzGU~LZ^ zs>MzuIiTxS*2nIiPRGz~GfkGJqI8M(>eX8)=moulxbM*j5qOhQ1@mQuhG5xqQ{BcI zKV%&*p?`VgWRSYIwpx?V@yO18r|Cop-rEWlg1+TZ;o!J&_Q_2>b}4X~}Hw8(L$cr{#URp2&)tqAFXsFa|< zKzI2fY94RxJ(H>Jn{qryR6a`T?C5|5ieF_^|~k-7*>_eo%e$iV6@}+&*x7DHICuC$+rVh4av89oV?-!S0Uag|+OTL$bF+1?=b|!}qM3GpGka@hG;Gw=l$$%}amyl-EO26#gs9z>mDG1M zOcuc?$p{NW)F?#t*Ku*fcIdu`PdNZDfBm{qiEEEphq<}nQFr2!8IAwCb<2%@JJJE* zrkISAknFayT|(hUG+KRc3e(!}xIppC?xqDt#UhaiG8cl)wC8@}H85t0WSsW)Q*{g$ z3t-y<^FB)!9cjR^Ylu^%vgQvm%jn=iE!0SeCtO|i54d28 z$HGIvTZ<>%fPjFSZ5}832L@P2#uH)1f(V9j#yCXH`$@0wz+%KPWt>E67Yuw@5iAx9 zSLW)~4T!5b7Zp-AjI)b@fz6w5z_~Tg))u4a6HrP0+w$0{*3rHB909VHa6(*N?^;I6 ztf>$?itoNFODz5TfK!Ac|8XQHWI`(rWn^UQv-3H){qQLCo8;w}SOfT|8SO*dGp((( zB}=C6t~fx{Izj?$`7`ZT2n&2~kTQW_6Rdgu(Gdem_4)bsq~T;KT2ulBy;pS)?Cboh zs&vd^;FSMZTAJKbxAk%wo^P|>zlb7jZfUt(SqV%pRuTg4Mn*>T({MlMF#oSd-hSj- Wb;hl@vaEnmBU=L_IG5Q8(*FX3p(_3W literal 0 HcmV?d00001 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. From 910a4dc01e179f6553edeee02e67b32da6c040f4 Mon Sep 17 00:00:00 2001 From: tillh Date: Fri, 3 Jul 2026 17:56:22 +0200 Subject: [PATCH 5/6] Add LevelDependentRateTES example Add a runnable example combining an electricity source, an electric boiler, a LevelDependentRateTES and a heat demand, plotting the theoretical c-rate curves and the realised rates. Add the example's Plots and DataFrames dependencies. --- examples/Project.toml | 2 + examples/advanced_example.jl | 204 +++++++++++++++++++++++++++++++++++ 2 files changed, 206 insertions(+) create mode 100755 examples/advanced_example.jl 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) From 34013917257c534af5ba2fa95405b10a95c4f528 Mon Sep 17 00:00:00 2001 From: tillh Date: Fri, 3 Jul 2026 18:44:49 +0200 Subject: [PATCH 6/6] formatting --- src/constraint_functions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constraint_functions.jl b/src/constraint_functions.jl index f226a89..8048c34 100644 --- a/src/constraint_functions.jl +++ b/src/constraint_functions.jl @@ -1,4 +1,4 @@ - """ +""" EMB.constraints_capacity(m, n::HeatPump, ๐’ฏ::TimeStructure, modeltype::EnergyModel) Method for creating the constraints on the maximum capacity of a [`HeatPump`](@ref).