Rolling horizon (#1327)

* Change rep_period profiles to handle rolling horizon

Implement ProfileWithRollingHorizon and make profiles.rep_period use
that instead of Vector{Float64}.
Change _profile_aggregate to dispatch on Vector{Float64} and
ProfileWithRollingHorizon.
These should allow the rolling horizon variables to be aggregated
the same way as the normal profiles.

Part of #1365

* Implement high-level run_rolling_horizon function with placeholders

Implement run_rolling_horizon with the main idea. Use placeholders
to specify places where a more complicated function needs to be called.
Update the EnergyProblem structure to hold an inner EnergyProblem
for the rolling horizon model.

Part of #1365

* Create rolling horizon parameters

Add ParametricOptInterface.
Implement adding the rolling horizon parameters in src/rolling-horizon/create.jl.
Update create_model to handle the rolling horizon parameters and wrap the
solver in POI.Optimizer if it's rolling horizon.

Part of #1365

* Implement update of rolling horizon parameters

Implement the update functions for rolling horizon and add them
to the run_rolling_horizon function.

Part of #1365

* Test adding rolling horizon parameters

Part of #1365

* Add functions for validation and data preparation for rolling horizon

Add validation of input in the rolling horizon function.
Add function to prepare the input to get in the rolling horizon and to get
out.
Add function to save the solution in the relevant tables.

Part of #1365

* Add unit tests for rolling horizon

Add tests of the various parts of rolling horizon.

Part of #1365

* Update the storage constraint to use the rolling initial storage level

Update src/constraints/storage.jl to use the initial_storage_level stored
in the rolling horizon.

Add test for the rolling horizon objective values to control future changes.

Part of #1365

Author: abelsiqueira

Project.toml

@@ -13,6 +13,7 @@ JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
@@ -28,6 +29,7 @@ JSON = "0.21.4, 1"
 JuMP = "1"
 MathOptInterface = "1"
 OrderedCollections = "1"
+ParametricOptInterface = "0.12.1"
 SparseArrays = "1"
 Statistics = "1"
 TOML = "1"

src/TulipaEnergyModel.jl

@@ -15,6 +15,7 @@ using TulipaIO: TulipaIO
 using HiGHS: HiGHS
 using JuMP: JuMP, @constraint, @expression, @objective, @variable
 using MathOptInterface: MathOptInterface
+using ParametricOptInterface: ParametricOptInterface as POI
 
 ## Others
 using LinearAlgebra: LinearAlgebra
@@ -25,10 +26,10 @@ using TimerOutputs: TimerOutput, @timeit
 const to = TimerOutput()
 
 # Definitions and auxiliary files
-include("utils.jl")
 include("run-scenario.jl")
 include("model-parameters.jl")
 include("structures.jl")
+include("utils.jl")
 
 # Data
 include("input-schemas.jl")
@@ -39,6 +40,9 @@ include("data-preparation.jl")
 # Data massage and model preparation
 include("model-preparation.jl")
 
+# Rolling horizon
+include("rolling-horizon/rolling-horizon.jl")
+
 # Model creation
 for folder_name in ["variables", "constraints", "expressions"]
     folder_path = joinpath(@__DIR__, folder_name)

src/constraints/storage.jl

@@ -5,10 +5,27 @@ export add_storage_constraints!
 
 Adds the storage asset constraints to the model.
 """
-function add_storage_constraints!(connection, model, variables, expressions, constraints, profiles)
+function add_storage_constraints!(
+    connection,
+    model,
+    variables,
+    expressions,
+    constraints,
+    profiles;
+    rolling_horizon = false,
+)
     var_storage_level_rep_period = variables[:storage_level_rep_period]
     var_storage_level_over_clustered_year = variables[:storage_level_over_clustered_year]
 
+    rolling_horizon_lookup = if rolling_horizon
+        Dict{Int,Int}(
+            row.var_storage_id::Int => row.id::Int for
+            row in DuckDB.query(connection, "FROM param_initial_storage_level")
+        )
+    else
+        Dict{Int,Int}()
+    end
+
     ## REP-PERIOD CONSTRAINTS (within a representative period)
     # - Balance constraint (using the lowest temporal resolution)
     let table_name = :balance_storage_rep_period, cons = constraints[table_name]
@@ -27,7 +44,12 @@ function add_storage_constraints!(connection, model, variables, expressions, con
                         sum,
                         0.0,
                     )
-                    initial_storage_level = row.initial_storage_level::Union{Float64,Missing}
+                    initial_storage_level = if rolling_horizon && row.time_block_start == 1
+                        rolling_horizon_id = rolling_horizon_lookup[row.id]
+                        variables[:param_initial_storage_level].container[rolling_horizon_id]
+                    else
+                        row.initial_storage_level::Union{Float64,Missing}
+                    end
                     storage_charging_efficiency = row.storage_charging_efficiency::Float64
                     storage_discharging_efficiency = row.storage_discharging_efficiency::Float64
 

src/create-model.jl

@@ -69,9 +69,19 @@ function create_model(
     model_file_name = "",
     enable_names = true,
     direct_model = false,
+    rolling_horizon = false,
+    rolling_horizon_window_length = 0,
 )
+    if rolling_horizon
+        @assert rolling_horizon_window_length > 0
+    end
+
     ## Optimizer
-    optimizer_with_attributes = JuMP.optimizer_with_attributes(optimizer, optimizer_parameters...)
+    optimizer_with_attributes = if rolling_horizon
+        JuMP.optimizer_with_attributes(() -> POI.Optimizer(optimizer()), optimizer_parameters...)
+    else
+        JuMP.optimizer_with_attributes(optimizer, optimizer_parameters...)
+    end
 
     ## Model
     if direct_model
@@ -113,6 +123,17 @@ function create_model(
         constraints,
     )
 
+    ## Rolling Horizon Parameters
+    if rolling_horizon
+        add_rolling_horizon_parameters!(
+            connection,
+            model,
+            variables,
+            profiles,
+            rolling_horizon_window_length,
+        )
+    end
+
     ## Expressions
     expressions = Dict{Symbol,TulipaExpression}()
 
@@ -165,7 +186,8 @@ function create_model(
         variables,
         expressions,
         constraints,
-        profiles,
+        profiles;
+        rolling_horizon,
     )
 
     @timeit to "add_hub_constraints!" add_hub_constraints!(model, constraints)

src/model-preparation.jl

@@ -572,19 +572,22 @@ function add_expressions_to_constraints!(connection, variables, constraints)
 end
 
 function prepare_profiles_structure(connection)
+    # Independent of being rolling horizon or not, these are complete
     rep_period = Dict(
-        (row.profile_name, row.year, row.rep_period) => [
-            row.value for row in DuckDB.query(
-                connection,
-                "SELECT profile.value
-                FROM profiles_rep_periods AS profile
-                WHERE
-                    profile.profile_name = '$(row.profile_name)'
-                    AND profile.year = $(row.year)
-                    AND profile.rep_period = $(row.rep_period)
-                ",
-            )
-        ] for row in DuckDB.query(
+        (row.profile_name, row.year, row.rep_period) => ProfileWithRollingHorizon(
+            Float64[
+                row.value for row in DuckDB.query(
+                    connection,
+                    "SELECT profile.value
+                    FROM profiles_rep_periods AS profile
+                    WHERE
+                        profile.profile_name = '$(row.profile_name)'
+                        AND profile.year = $(row.year)
+                        AND profile.rep_period = $(row.rep_period)
+                    ",
+                )
+            ],
+        ) for row in DuckDB.query(
             connection,
             "SELECT DISTINCT
                 profiles.profile_name,
@@ -596,7 +599,7 @@ function prepare_profiles_structure(connection)
     )
 
     over_clustered_year = Dict(
-        (row.profile_name, row.year) => [
+        (row.profile_name, row.year) => Float64[
             row.value for row in DuckDB.query(
                 connection,
                 "SELECT profile.value

src/rolling-horizon/create.jl

@@ -0,0 +1,77 @@
+"""
+    add_scalar_rolling_horizon_parameter!(model, variables, connection, param_table_name)
+
+Create the rolling horizon Parameters for the table `param_table_name`.
+This assumes that the table named `param_table_name` has a column
+`original_value` and uses that as initial value for the variable.
+
+The container is saved in the `model[key]` and the `TulipaVariable` for the
+parameter is stored in `variables[key]`, where `key = Symbol(param_table_name)`.
+"""
+function add_scalar_rolling_horizon_parameter!(model, variables, connection, param_table_name)
+    initial_value =
+        [row.original_value::Float64 for row in DuckDB.query(connection, "FROM $param_table_name")]
+    num_rows = length(initial_value)
+    param = TulipaVariable(connection, param_table_name)
+    param.container = @variable(model, [1:num_rows] in JuMP.Parameter.(initial_value))
+
+    key = Symbol(param_table_name)
+    model[key] = param.container
+    variables[key] = param
+    return param
+end
+
+"""
+    add_rolling_horizon_parameters!(connection, model, variables, profiles, window_length)
+
+Create Parameters to handle rolling horizon.
+
+The profile parameters are attached to `profiles.rep_period`.
+
+The other parameters are the ones that have initial value (currently only initial_storage_level).
+These must be filtered from the corresponding indices table when time_block_start = 1.
+The corresponding parameters is saved in the variables and in the model.
+"""
+function add_rolling_horizon_parameters!(connection, model, variables, profiles, window_length)
+    # Profiles
+    for (_, profile_object) in profiles.rep_period
+        profile_object.rolling_horizon_variables =
+            @variable(model, [1:window_length] in JuMP.Parameter(0.0))
+    end
+
+    # Scalar rolling horizon parameters
+    # These need a table filtering where time_block_start = 1, and the current
+    # strategy is to create a table `param_NAME` and call
+    # add_scalar_rolling_horizon_parameter!.
+
+    ## initial_storage_level
+    DuckDB.query(
+        connection,
+        """
+        DROP SEQUENCE IF EXISTS id;
+        CREATE SEQUENCE id START 1;
+        CREATE OR REPLACE TABLE param_initial_storage_level AS
+        SELECT
+            nextval('id') as id,
+            var.asset,
+            var.year,
+            var.rep_period,
+            var.id as var_storage_id,
+            asset_milestone.initial_storage_level as original_value
+        FROM var_storage_level_rep_period as var
+        LEFT JOIN asset_milestone
+            ON var.asset = asset_milestone.asset
+            AND var.year = asset_milestone.milestone_year
+        WHERE time_block_start = 1;
+        DROP SEQUENCE id;
+        """,
+    )
+    add_scalar_rolling_horizon_parameter!(
+        model,
+        variables,
+        connection,
+        "param_initial_storage_level",
+    )
+
+    return
+end