Identifying the strategic period of a period.

[JulStraus]

Problem statement

As outlined in #73, we currently allow for identifying whether a TimePeriod, AbstractOperationalScenario, or AbstractRepresentativePeriod is within an AbstractStrategicPeriod using the function _sp_period. We require it for the discounting approach to identify which strategic period has to be utilized when the input to discount or objective_weight is not a strategic period. The current system works but has a major issue as outlined in #73: It only checks that the field sp is the same (or sp and branch in the case of TwoLevelTree).
This implies that when working with multiple time structures, e.g., when using Plasmo, we can experience problems.

In addition, we only utilize the current system for strategic periods as we do not need it for a different purpose. It can be easily implemented as well for other types.

Solution

A new implementation for strategic periods is given by

_check_period(sp::AbstractStrategicPeriod, t::TimePeriod) = t in sp
_check_period(sp::AbstractStrategicPeriod, t::AbstractOperationalScenario) = t in opscenarios(sp)
_check_period(sp::AbstractStrategicPeriod, t::AbstractRepresentativePeriod) = t in repr_periods(sp)

_identify_sp(sp::AbstractStrategicPeriod, ts::TimeStructure) = sp
function _identify_sp(t::Union{TimePeriod, AbstractOperationalScenario, AbstractRepresentativePeriod}, ts::TimeStructure)
    sps = collect(strat_periods(ts))
    filter_fun(sp) = _check_period(sp, t)
    per = findfirst(filter_fun, sps)
    isnothing(per) && throw(ErrorException("Period not part of any strategic period"))
    return sps[per]
end

This implementation can be extended for both representative periods

_check_period(rp::AbstractRepresentativePeriod, t::TimePeriod) = t in rp
_check_period(rp::AbstractRepresentativePeriod, t::AbstractOperationalScenario) = t in opscenarios(rp)

_identify_rp(rp::AbstractRepresentativePeriod, ts::TimeStructure) = rp
function _identify_rp(t::Union{TimePeriod, AbstractOperationalScenario}, ts::TimeStructure)
    rps = collect(repr_periods(ts))
    filter_fun(sp) = _check_period(sp, t)
    per = findfirst(filter_fun, rps)
    isnothing(per) && throw(ErrorException("Period not part of any representative period"))
    return rps[per]
end

and operational scenarios

_check_period(osc::AbstractOperationalScenario, t::TimePeriod) = t in osc

_identify_osc(osc::AbstractOperationalScenario, ts::TimeStructure) = osc
function _identify_osc(t::TimePeriod, ts::TimeStructure)
    oscs = collect(opscenarios(ts))
    filter_fun(sp) = _check_period(sp, t)
    per = findfirst(filter_fun, oscs)
    isnothing(per) && throw(ErrorException("Period not part of any operational scenario"))
    return oscs[per]
end

Discussion

The new approach is in practice slower than the old approach as we do not only compare the field sp, but the complete OperationalPeriod. Due to the application of findfirst, its penalty is dependent on which period we are interested in (in which strategic period it is located) and how the time structure looks like. The penalty can be high if we have a lot of strategic periods and operational periods.

As an example, consider the following time structure:

# 15 strategic periods with 4 weeks each
ts = TwoLevel(15, 5, SimpleTimes(168*4, 1); op_per_strat=8760)

The required time through specifying the functions

function test_old(t, ts)
    for k ∈ range(1,1e5)
        TimeStruct._sp_period(t, ts)
    end
end
function test_new(t, ts)
    for k ∈ range(1,1e5)
        _identify_sp(t, ts)
    end
end

is dependent on the chosen time period (which strategic period it is in). If it is in the first strategic period, it requires on my PC 0.11 s compared to 0.08 s for the old method. If it is in the 15th, that changes to 0.52 s compared to 0.07 s.

This leaves the question whether we want to include it or not. Do we see any problems by just comparing a subset of parameters (_strat_per and comparable for the new methods) of a period? What are your thoughts @hellemo and @trulsf?

[trulsf]

I think the safe approach as suggested has some efficiency challenges with the current implementation for iterating time structures.
The following example shows an illustration (where I have used _identify_sp() instead of _sp_period):

scenarios = OperationalScenarios([SimpleTimes(24,1) for i in 1:5])
repr = RepresentativePeriods(4, 1000, scenarios)
periods = TwoLevel(10, 5, repr; op_per_strat = 8760)

model = Model()
@variable(model, x[periods] >= 0)
@time @objective(model, Min, sum(objective_weight(t, periods, 0.05; type="start") * x[t] for t in periods))

This will take several seconds even for this moderately sized time structure. Increasing to a week for each scenario, it takes a lot of time.

The problem seems to reside with the calculation of _multiple_adj which involves the _total_duration function that currently has an inefficient linear implementation. Thus the overall complexity will be of order n^3.

If we are to use this approach, there needs to be a more efficient implementation of the _total_duration function, e.g., through the use of some caching mechanism.

The other challenge is that it is not guaranteed to catch all errors with regard to having multiple time structures. The check t in sp will only check that t == tt for some OperationalPeriod tt in the strategic period. The following example shows that this may not be enough to discern that the time structure is the same

periods1 = TwoLevel(3, [1, 1, 1], [day, day, day], 8760.0)
periods2 = TwoLevel(3, [1, 2, 1], [day, day, day], 8760.0)
periods3 = TwoLevel(3, [2, 1, 1], [day, day, day], 8760.0)

ts1 = collect(periods1)

val1 = discount(ts1[34], periods1, 0.05; type="start") # val1 = 0.9523809523809523
val2 = discount(ts1[34], periods2, 0.05; type="start") # Errors with period not part of any strategic period
val3 = discount(ts1[34], periods3, 0.05; type="start") # val3 = 0.9070294784580498

Thus, subtle errors may creep in any anyway. Maybe a first step is to improve documentation to note that users should be aware. It is also possible to hide the discount and multiple functions, and instead just expose functions like present_value or scale that works on iterators/collections:

present_value(f, periods, discounter) = sum(discount(discounter, t) * f(t) for t in periods)
scale(f, periods) = map(t -> multiple(t) * f(t), periods)

# Ex
pv = present_value(t -> capex[t] * x[t], periods, discounter)

[JulStraus]

Your example is interesting. I am always a bit confused in which case equality with custom types exists within Julia. While your example seems to be understandable, it is quite often that one has to point to the exact instance when comparing equality.

I am always in favor of improving documentation, but I would not necessarily change to objective_value or discount. We could make them internal and create a new function, but that should be accompanied with other changes as it is only a minor change, but breaking.