Add functions to compute time resolution periods and matrices

src/TulipaEnergyModel.jl

@@ -7,8 +7,9 @@ using Graphs
 using HiGHS
 using JuMP
 
+include("input_tables.jl")
 include("io.jl")
 include("model.jl")
-include("input_tables.jl")
+include("time-resolution.jl")
 
 end

src/time-resolution.jl

@@ -0,0 +1,132 @@
+export resolution_matrix, compute_rp_periods
+
+"""
+    M = resolution_matrix(rp_periods, time_steps)
+
+Computes the resolution balance matrix using the array of `rp_periods` and the array of `time_steps`.
+The elements in these arrays must be ranges.
+
+## Examples
+
+The following two examples are for two flows/assets with resolutions of 3h and 4h, so that the representative period has 4h periods.
+
+```jldoctest
+rp_periods = [1:4, 5:8, 9:12]
+time_steps = [1:4, 5:8, 9:12]
+resolution_matrix(rp_periods, time_steps)
+
+# output
+
+3×3 Matrix{Float64}:
+ 1.0  0.0  0.0
+ 0.0  1.0  0.0
+ 0.0  0.0  1.0
+```
+
+```jldoctest
+rp_periods = [1:4, 5:8, 9:12]
+time_steps = [1:3, 4:6, 7:9, 10:12]
+resolution_matrix(rp_periods, time_steps)
+
+# output
+
+3×4 Matrix{Float64}:
+ 1.0  0.333333  0.0       0.0
+ 0.0  0.666667  0.666667  0.0
+ 0.0  0.0       0.333333  1.0
+```
+"""
+function resolution_matrix(
+    rp_periods::AbstractVector{<:UnitRange{<:Integer}},
+    time_steps::AbstractVector{<:UnitRange{<:Integer}},
+)
+    matrix = [
+        length(period ∩ time_step) / length(time_step) for period in rp_periods,
+        time_step in time_steps
+    ]
+
+    return matrix
+end
+
+"""
+    rp_periods = compute_rp_periods(array_time_steps)
+
+Given the time steps of various flows/assets in the `array_time_steps` input, compute the representative period splits.
+Each element of `array_time_steps` is an array of ranges with the following assumptions:
+
+  - An element is of the form `V = [r₁, r₂, …, rₘ]`, where each `rᵢ` is a range `a:b`.
+  - `r₁` starts at 1.
+  - `rᵢ₊₁` starts at the end of `rᵢ` plus 1.
+  - `rₘ` ends at some value `N`, that is the same for all elements of `array_time_steps`.
+
+Notice that this implies that they form a disjunct partition of `1:N`.
+
+The output will also be an array of ranges with the conditions above.
+The output is constructed greedily, i.e., it selects the next largest breakpoint following the algorithm below:
+
+ 0. Input: `Vᴵ₁, …, Vᴵₚ`, a list of time step ranges. Each element of `Vᴵⱼ` is a range `r = r.start:r.end`. Output: `V`.
+ 1. Compute the end of the representative period `N` (all `Vᴵⱼ` should have the same end)
+ 2. Start with an empty `V = []`
+ 3. Define the beginning of the range `s = 1`
+ 4. Define an array with all the next breakpoints `B` such that `Bⱼ` is the first `r.end` such that `r.end ≥ s` for each `r ∈ Vᴵⱼ`.
+ 5. The end of the range will be the `e = max Bⱼ`.
+ 6. Define `r = s:e` and add `r` to the end of `V`.
+ 7. If `e = N`, then END
+ 8. Otherwise, define `s = e + 1` and go to step 4.
+
+## Examples
+
+```jldoctest
+time_steps1 = [1:4, 5:8, 9:12]
+time_steps2 = [1:3, 4:6, 7:9, 10:12]
+compute_rp_periods([time_steps1, time_steps2])
+
+# output
+
+3-element Vector{UnitRange{Int64}}:
+ 1:4
+ 5:8
+ 9:12
+```
+
+```jldoctest
+time_steps1 = [1:1, 2:3, 4:6, 7:10, 11:12]
+time_steps2 = [1:2, 3:4, 5:5, 6:7, 8:9, 10:12]
+compute_rp_periods([time_steps1, time_steps2])
+
+# output
+
+5-element Vector{UnitRange{Int64}}:
+ 1:2
+ 3:4
+ 5:6
+ 7:10
+ 11:12
+```
+"""
+function compute_rp_periods(
+    array_time_steps::AbstractVector{<:AbstractVector{<:UnitRange{<:Integer}}},
+)
+    # Get Vᴵ₁, the last range of it, the last element of the range
+    representative_period_end = array_time_steps[1][end][end]
+    for time_steps in array_time_steps
+        # Assumption: All start at 1 and end at N
+        @assert time_steps[1][1] == 1
+        @assert representative_period_end == time_steps[end][end]
+    end
+    rp_periods = UnitRange{Int}[] # List of ranges
+    period_start = 1
+
+    while period_start < representative_period_end
+        # The first range end larger than period_start for each range in each time_steps.
+        breakpoints = (
+            first(r[end] for r in time_steps if r[end] ≥ period_start) for
+            time_steps in array_time_steps
+        )
+        period_end = maximum(breakpoints)
+        @assert period_end ≥ period_start
+        push!(rp_periods, period_start:period_end)
+        period_start = period_end + 1
+    end
+    return rp_periods
+end

test/runtests.jl

@@ -70,3 +70,50 @@ end
         )
     end
 end
+
+@testset "Time resolution" begin
+    @testset "resolution_matrix" begin
+        rp_periods = [1:4, 5:8, 9:12]
+        time_steps = [1:4, 5:8, 9:12]
+        expected = [
+            1.0 0.0 0.0
+            0.0 1.0 0.0
+            0.0 0.0 1.0
+        ]
+        @test resolution_matrix(rp_periods, time_steps) == expected
+
+        time_steps = [1:3, 4:6, 7:9, 10:12]
+        expected = [
+            1.0 1/3 0.0 0.0
+            0.0 2/3 2/3 0.0
+            0.0 0.0 1/3 1.0
+        ]
+        @test resolution_matrix(rp_periods, time_steps) == expected
+
+        time_steps = [1:6, 7:9, 10:10, 11:11, 12:12]
+        expected = [
+            2/3 0.0 0.0 0.0 0.0
+            1/3 2/3 0.0 0.0 0.0
+            0.0 1/3 1.0 1.0 1.0
+        ]
+        @test resolution_matrix(rp_periods, time_steps) == expected
+    end
+
+    @testset "compute_rp_periods" begin
+        # regular
+        time_steps1 = [1:4, 5:8, 9:12] # every 4 hours
+        time_steps2 = [1:3, 4:6, 7:9, 10:12] # every 3 hours
+        time_steps3 = [i:i for i = 1:12] # hourly
+
+        @test compute_rp_periods([time_steps1, time_steps2]) == time_steps1
+        @test compute_rp_periods([time_steps1, time_steps2, time_steps3]) == time_steps1
+        @test compute_rp_periods([time_steps2, time_steps3]) == time_steps2
+
+        # Irregular
+        time_steps4 = [1:6, 7:9, 10:11, 12:12]
+        time_steps5 = [1:2, 3:4, 5:12]
+        @test compute_rp_periods([time_steps1, time_steps4]) == [1:6, 7:9, 10:12]
+        @test compute_rp_periods([time_steps1, time_steps5]) == [1:4, 5:12]
+        @test compute_rp_periods([time_steps4, time_steps5]) == [1:6, 7:12]
+    end
+end