Add LevelDependentRateTES node with state-of-charge dependent charge/discharge rates

[tillho]

IMPORTANT This is not a pressing issue. Don't rush, enjoy your summer break and have a look at this whenever it suits you or when you're back :) God sommer!

Summary

Adds a new thermal energy storage node: LevelDependentRateTES. Like ThermalEnergyStorage and BoundRateTES it behaves like a RefStorage with heat losses proportional to the storage level, but its distinguishing feature is that the maximum charge and discharge rates depend on the current state of charge, described by user-provided c-rate curves.

Design notes

• c-rate curves: for charging, the curve starts at (0, capacity(charge)) and runs through the anchor points in ascending level order; for discharging, it starts at (capacity(level), capacity(discharge)) and runs in descending order. One point gives a single line; two or three give a piecewise-linear curve whose active segment is selected by big-M binaries.
• Tight big-M: because the node is restricted to OperationalModel, all slopes and intercepts are constants, so each big-M is computed from the geometry.
• OperationalModel only: in an InvestmentModel the installed capacities become decision variables, which would turn the rate limit into a bilinear (capacity × storage-level) term that the LP/MILP solvers used here cannot handle. This is enforced by dispatch and documented.
• SOS2: was considered but not used because I think the explicit formulation is more transparent for users, which is something that can be discussed.

[JulStraus]

While it is ok that the node is only valid without investments, the restriction to create the methods only for OperationalModel has a few problems:

1. It should be possible to use it both with InvestmentModel and RecHorOperationalModel, but without StorageInvData. This is in the current situation not the case.
2. The node falls back to a standard AbstractTES in the case of using an InvestmentModel, even if it does not include StorageInvData. That should not be the default, at least in my opinion.

The approach should be instead implemented similar to the approach chosen in EnergyModelsHydrogen.

There are also a few open questions for me regarding the break points. While I understand why they are formed the way they are, I think we can improve on the concept. specifically, we do not need binary variables if the bounds are convex.

I will post a first draft for how we can implement that more elegant.

[JulStraus]

Please do not add a newline to fit it within ca 90 characters. It is a Markdown file where our standard is to use a new line for each new sentence. That makes maintenance significantly easier.

That must be adjusted.

In addition, you should follow the standard approach for describing a node. I do prefer that all information for a given node is on a single page.

[JulStraus]

The @ref does not do anything in the example here, so I would suggest to change it.

Suggested change

	[`LevelDependentRateTES`](@ref) and a heat demand. The boiler converts electricity to heat,
	`LevelDependentRateTES` and a heat demand. The boiler converts electricity to heat,

[JulStraus]

I would move these parts into the normal code. That makes it easier for the reader to understand.

[JulStraus]

I would suggest to have the functions in subfiles, e.g., a file called calculate.jl for compute_charge_curve and compute_discharge_curve (also, I think a more descriptive name would be calculate_charge_curve and calculate_discharge_curve and visualize.jl for visualize_c_rates

[JulStraus]

This function should be based on the node and not the individual parameters, i.e., EnergyModelsHeat.visualize_c_rates(n::LevelDependentRateTES) where the individual parameters are extracted from the node through

charge_capacity = capacity(charge(n))
discharge_capacity = capacity(discharge(n))
level_capacity =  = capacity(level(n))
c_rate_points_charge = EnergyModelsHeat.c_rate_points_charge(n)
c_rate_points_discharge = EnergyModelsHeat.c_rate_points_discharge(n)

This is still a bit problematic as the capacities are TimeProfiles, but that should be relevant anyhow as we can have changing capacities per strategic period, even with this setup.

[JulStraus]

Use the same format as for ThermalEnergyStorage`.

You can skip the level constraints if you provide a link to the thermal energy storage page.

[JulStraus]

Would it be possible to add to the figures the bounds? That would improve the understanding. In addition, the figures xlabel and the charge y label gets cropped in my version.

[JulStraus]

General comments for the implementation:

1. No @warn or @info in the code.
2. The function is rather long with if statements. This can be circumvented through creating the limits on the rate actual subtypes and not using a Vector{Vector}, see the comment in node.jl.
3. You do not have any constraints on the first operational period of any substructure. This is not a good idea in my opinion as it can lead to entirely different behavior in the first period. That is especially a problem if you have a TwoLevel{RepresentativePeriods{OperationalScenarios{SimpleTimes}}} time structure. This should be fixed similar to the approach chosen for EnergyModelsHydrogen. Through introduction of subtypes, this may be fixed for your type on the lowest level outlined above.
4. I am not really certain that the implementation is easier to read or more transparent compared to a SOS2 implementation ;)

[JulStraus]

I would strongly advise you to create specific types for this behavior. There are quite a few reasons for that,

1. it reduces the issues for variable creation as binaries are not always required,
2. it simplifies the actual implementation within constraints_capacity for calling a subfunction, depending on the type,
3. it allows you to have investments possible if (and only if) the bounds are convex as you may in this situation allow for inequality constraints that are not piecewise linear, and
4. it simplifies the checks and avoid potential problems,

A potential structure is posted at the end.

[JulStraus]

I would strongly advise you to create specific types for this behavior. There are quite a few reasons for that,

1. it reduces the issues for variable creation as binaries are not always required (convex bounds do not require binaries),
2. it simplifies the actual implementation within constraints_capacity for calling a subfunction, depending on the type,
3. it allows you to have investments possible if (and only if) the bounds are convex AND the ratio between the charge and level capacity is constant in the investments (could be enforced in the investments, if we want to) as you may in this situation allow for inequality constraints that are not piecewise linear, and
4. it simplifies the checks and avoid potential problems,

I propose a potential structure in a different comment.

[JulStraus]

The following approach could be utilized for incorporating the bounds:

# Type for specifying whether it is convex or concave
abstract type StorBoundsType end
struct ConvexBounds <: StorBoundsType end
struct ConcaveBounds <: StorBoundsType end

# Type with all slopes and intercepts
struct StorBounds{T}
    slopes::Vector
    intercepts::Vector
    levels::Vector
end

# Constructor for calculating whether the system is convex or concave as well as for calculating
# the slopes and intercepts
function StorBounds(x::Vector{Float64}, y::Vector{Float64})
    bool = true
    intercepts = []
    slopes = []
    levels = []
    for k ∈ 2:length(x)
        push!(levels, x[k])
        push!(slopes, (y[k] - y[k-1]) / (x[k] - x[k-1]))
        push!(intercepts, y[k] - x[k] * slopes[k-1])
        if k ≥ 3 && slopes[k-1] > slopes[k-2]
            bool = false
        end
    end
    idx =  findall(isinf, slopes)
    deleteat!(slopes, idx)
    deleteat!(intercepts, idx)
    deleteat!(levels, idx)
    if bool
        StorBounds{ConvexBounds}(slopes, intercepts, levels)
    else
        StorBounds{ConcaveBounds}(slopes, intercepts, levels)
    end
end


# Functions for creating specific type
function ChargeStorBounds(
    level::Vector{<:Number},
    charge::Vector{<:Number},
)
    insert!(level, 1, 0.0)
    push!(level, 1.0)
    insert!(charge, 1, 1.0)
    push!(charge, 0.0)
    return StorBounds(level, charge)
end
function DischargeStorBounds(
    level::Vector{<:Number},
    discharge::Vector{<:Number},
)
    insert!(level, 1, 0.0)
    push!(level, 1.0)
    insert!(discharge, 1, 0.0)
    push!(discharge, 1.0)
    return StorBounds(level, discharge)
end

# Testing
y = [1.0, 0.9,]
x = [0.2, 0.5,]
s = ChargeStorBounds(x, y)

y = [0.2, 0.8]
x = [0.0, 0.5]
s = DischargeStorBounds(x, y)

The advantage is that we

1. know directly whether the system is convex or concave and can use the information for creating binaries, or not,
2. allow for more potential break points than the hard-coded number,
3. reduce the number of code lines through reutilization of information.

The number of binaries can be identified through the following function:

number_bin(bounds::StorBounds{ConvexBounds}) = 0
number_bin(bounds::StorBounds{ConcaveBounds}) = length(bounds.intercepts)

and the constraints could be implemented in the convex case as

function constraints_cap_bounds(
    m,
    n::LevelDependentRateTES,
    t::TS.OperationalPeriod,
    t_prev::TS.OperationalPeriod,
    bounds::StorBounds{ConvexBounds},
)
    ratio = capacity(charge(n, t))/capacity(level(n, t))
    for (a, b) ∈ zip(bounds.slopes, bounds.intercepts)
        @constraint(m,
            m[:stor_charge_use][n, t] ≤
                a * ratio * m[:stor_level][n, t_prev] + b * capacity(charge(n))
        )
    end
end