From 4896fb653fe4e2d13bbd5a4547b793d8ddf7b243 Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 30 Aug 2024 11:22:37 +1200 Subject: [PATCH] [docs] add Rolling horizon tutorial Co-authored-by: datejada --- docs/Project.toml | 1 + .../tutorials/algorithms/rolling_horizon.csv | 169 +++++++++ .../tutorials/algorithms/rolling_horizon.jl | 331 ++++++++++++++++++ 3 files changed, 501 insertions(+) create mode 100644 docs/src/tutorials/algorithms/rolling_horizon.csv create mode 100644 docs/src/tutorials/algorithms/rolling_horizon.jl diff --git a/docs/Project.toml b/docs/Project.toml index 255889e1268..c4ad9ab2dfe 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -26,6 +26,7 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" MultiObjectiveAlgorithms = "0327d340-17cd-11ea-3e99-2fd5d98cecda" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" PATHSolver = "f5f7c340-0bb3-5c69-969a-41884d311d1b" +ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/docs/src/tutorials/algorithms/rolling_horizon.csv b/docs/src/tutorials/algorithms/rolling_horizon.csv new file mode 100644 index 00000000000..dfae1396ab3 --- /dev/null +++ b/docs/src/tutorials/algorithms/rolling_horizon.csv @@ -0,0 +1,169 @@ +t,demand_MW,solar_pu +2000-01-01T00:00:00,51.6,0 +2000-01-01T01:00:00,49.2,0 +2000-01-01T02:00:00,46.5,0 +2000-01-01T03:00:00,44.3,0 +2000-01-01T04:00:00,43.3,0 +2000-01-01T05:00:00,42.1,0 +2000-01-01T06:00:00,39.8,0 +2000-01-01T07:00:00,40.2,0 +2000-01-01T08:00:00,41.3,0.212560386 +2000-01-01T09:00:00,45,0.608695652 +2000-01-01T10:00:00,49.3,0.845410628 +2000-01-01T11:00:00,54.3,0.995169082 +2000-01-01T12:00:00,56,1 +2000-01-01T13:00:00,54.9,0.763285024 +2000-01-01T14:00:00,53.3,0.309178744 +2000-01-01T15:00:00,53.5,0.009661836 +2000-01-01T16:00:00,57.5,0 +2000-01-01T17:00:00,65,0 +2000-01-01T18:00:00,66.2,0 +2000-01-01T19:00:00,64.5,0 +2000-01-01T20:00:00,61,0 +2000-01-01T21:00:00,59,0 +2000-01-01T22:00:00,58.7,0 +2000-01-01T23:00:00,54.1,0 +2000-01-02T00:00:00,49.7,0 +2000-01-02T01:00:00,46.5,0 +2000-01-02T02:00:00,44.8,0 +2000-01-02T03:00:00,44.5,0 +2000-01-02T04:00:00,46,0 +2000-01-02T05:00:00,48.6,0 +2000-01-02T06:00:00,52.6,0 +2000-01-02T07:00:00,59,0 +2000-01-02T08:00:00,65.1,0.096618357 +2000-01-02T09:00:00,70.1,0.256038647 +2000-01-02T10:00:00,73.5,0.391304348 +2000-01-02T11:00:00,76.2,0.47826087 +2000-01-02T12:00:00,76.8,0.531400966 +2000-01-02T13:00:00,75.1,0.434782609 +2000-01-02T14:00:00,73.2,0.202898551 +2000-01-02T15:00:00,72.5,0.014492754 +2000-01-02T16:00:00,75.2,0 +2000-01-02T17:00:00,80.7,0 +2000-01-02T18:00:00,80.7,0 +2000-01-02T19:00:00,77.5,0 +2000-01-02T20:00:00,71.3,0 +2000-01-02T21:00:00,67.6,0 +2000-01-02T22:00:00,65.8,0 +2000-01-02T23:00:00,60.4,0 +2000-01-03T00:00:00,54.7,0 +2000-01-03T01:00:00,50.9,0 +2000-01-03T02:00:00,48.5,0 +2000-01-03T03:00:00,47.7,0 +2000-01-03T04:00:00,48.2,0 +2000-01-03T05:00:00,48.5,0 +2000-01-03T06:00:00,49.1,0 +2000-01-03T07:00:00,53.3,0 +2000-01-03T08:00:00,58.9,0.09178744 +2000-01-03T09:00:00,64.6,0.265700483 +2000-01-03T10:00:00,68.8,0.367149758 +2000-01-03T11:00:00,72,0.400966184 +2000-01-03T12:00:00,72.4,0.347826087 +2000-01-03T13:00:00,70.9,0.251207729 +2000-01-03T14:00:00,69.5,0.111111111 +2000-01-03T15:00:00,69.5,0.009661836 +2000-01-03T16:00:00,72.5,0 +2000-01-03T17:00:00,77.3,0 +2000-01-03T18:00:00,77.4,0 +2000-01-03T19:00:00,73.9,0 +2000-01-03T20:00:00,68,0 +2000-01-03T21:00:00,64.1,0 +2000-01-03T22:00:00,62.8,0 +2000-01-03T23:00:00,58.1,0 +2000-01-04T00:00:00,52.8,0 +2000-01-04T01:00:00,49.1,0 +2000-01-04T02:00:00,47,0 +2000-01-04T03:00:00,45.9,0 +2000-01-04T04:00:00,46.1,0 +2000-01-04T05:00:00,45.5,0 +2000-01-04T06:00:00,44.1,0 +2000-01-04T07:00:00,46.5,0.004830918 +2000-01-04T08:00:00,50.3,0.256038647 +2000-01-04T09:00:00,55.6,0.700483092 +2000-01-04T10:00:00,60.3,0.888888889 +2000-01-04T11:00:00,65.6,0.93236715 +2000-01-04T12:00:00,65.9,0.787439614 +2000-01-04T13:00:00,63.2,0.550724638 +2000-01-04T14:00:00,60.7,0.275362319 +2000-01-04T15:00:00,60.1,0.019323671 +2000-01-04T16:00:00,63.4,0 +2000-01-04T17:00:00,71.3,0 +2000-01-04T18:00:00,73.1,0 +2000-01-04T19:00:00,70.9,0 +2000-01-04T20:00:00,66.8,0 +2000-01-04T21:00:00,64.2,0 +2000-01-04T22:00:00,63.9,0 +2000-01-04T23:00:00,58.9,0 +2000-01-05T00:00:00,54,0 +2000-01-05T01:00:00,50.7,0 +2000-01-05T02:00:00,49.4,0 +2000-01-05T03:00:00,49.6,0 +2000-01-05T04:00:00,51.7,0 +2000-01-05T05:00:00,56.9,0 +2000-01-05T06:00:00,66.2,0 +2000-01-05T07:00:00,76.3,0.009661836 +2000-01-05T08:00:00,82,0.29468599 +2000-01-05T09:00:00,83.8,0.628019324 +2000-01-05T10:00:00,85.9,0.777777778 +2000-01-05T11:00:00,87.7,0.893719807 +2000-01-05T12:00:00,87.7,0.874396135 +2000-01-05T13:00:00,86.2,0.743961353 +2000-01-05T14:00:00,84.7,0.444444444 +2000-01-05T15:00:00,83.9,0.057971014 +2000-01-05T16:00:00,85.9,0 +2000-01-05T17:00:00,92,0 +2000-01-05T18:00:00,92,0 +2000-01-05T19:00:00,89,0 +2000-01-05T20:00:00,82,0 +2000-01-05T21:00:00,77.2,0 +2000-01-05T22:00:00,74.1,0 +2000-01-05T23:00:00,67,0 +2000-01-06T00:00:00,61.8,0 +2000-01-06T01:00:00,58,0 +2000-01-06T02:00:00,56.3,0 +2000-01-06T03:00:00,56.4,0 +2000-01-06T04:00:00,57.7,0 +2000-01-06T05:00:00,60.6,0 +2000-01-06T06:00:00,67.4,0 +2000-01-06T07:00:00,75.7,0.009661836 +2000-01-06T08:00:00,79.7,0.256038647 +2000-01-06T09:00:00,81.7,0.584541063 +2000-01-06T10:00:00,84.2,0.821256039 +2000-01-06T11:00:00,86.3,0.942028986 +2000-01-06T12:00:00,86,0.884057971 +2000-01-06T13:00:00,83.8,0.661835749 +2000-01-06T14:00:00,81.5,0.328502415 +2000-01-06T15:00:00,80.9,0.028985507 +2000-01-06T16:00:00,83.8,0 +2000-01-06T17:00:00,90.7,0 +2000-01-06T18:00:00,90.7,0 +2000-01-06T19:00:00,88.2,0 +2000-01-06T20:00:00,82.1,0 +2000-01-06T21:00:00,77.2,0 +2000-01-06T22:00:00,73.9,0 +2000-01-06T23:00:00,67.5,0 +2000-01-07T00:00:00,61.8,0 +2000-01-07T01:00:00,57.9,0 +2000-01-07T02:00:00,56.9,0 +2000-01-07T03:00:00,57.5,0 +2000-01-07T04:00:00,59.2,0 +2000-01-07T05:00:00,64.8,0 +2000-01-07T06:00:00,77.9,0 +2000-01-07T07:00:00,89.3,0.004830918 +2000-01-07T08:00:00,94.1,0.154589372 +2000-01-07T09:00:00,94.4,0.434782609 +2000-01-07T10:00:00,95.9,0.589371981 +2000-01-07T11:00:00,97.3,0.70531401 +2000-01-07T12:00:00,96.7,0.647342995 +2000-01-07T13:00:00,95.6,0.531400966 +2000-01-07T14:00:00,93.7,0.265700483 +2000-01-07T15:00:00,92.7,0.028985507 +2000-01-07T16:00:00,94,0 +2000-01-07T17:00:00,100,0 +2000-01-07T18:00:00,99.2,0 +2000-01-07T19:00:00,95.8,0 +2000-01-07T20:00:00,88.9,0 +2000-01-07T21:00:00,83.3,0 +2000-01-07T22:00:00,79.2,0 +2000-01-07T23:00:00,71.3,0 diff --git a/docs/src/tutorials/algorithms/rolling_horizon.jl b/docs/src/tutorials/algorithms/rolling_horizon.jl new file mode 100644 index 00000000000..535f4f07bc2 --- /dev/null +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -0,0 +1,331 @@ +# Copyright (c) 2024 Diego Tejada and contributors #src +# #src +# Permission is hereby granted, free of charge, to any person obtaining a copy #src +# of this software and associated documentation files (the "Software"), to deal #src +# in the Software without restriction, including without limitation the rights #src +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell #src +# copies of the Software, and to permit persons to whom the Software is #src +# furnished to do so, subject to the following conditions: #src +# #src +# The above copyright notice and this permission notice shall be included in all #src +# copies or substantial portions of the Software. #src +# #src +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR #src +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, #src +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE #src +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER #src +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, #src +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE #src +# SOFTWARE. #src + +# # Rolling horizon problems +# +# **This tutorial was originally contributed by Diego Tejada.** +# +# The term "rolling horizon" refers to solving a time-dependent model +# repeatedly, where the planning interval is shifted forward in time during each +# solution step. +# +# With the features provided in [ParametricOptInterface.jl](@ref), setting up +# such a model is quite straightforward. This tutorial explains the necessary +# steps to implement a basic model with a rolling horizon in a simple generation +# expansion problem (GEP). + +# ## Required packages +# +# This tutorial uses the following packages + +using JuMP +import CSV +import DataFrames +import Dates +import HiGHS +import ParametricOptInterface as POI +import Plots +import StatsPlots + +# ## The optimization model +# +# The model is a simplified GEP problem in which we decide the new capacity in +# renewables for a power system with a given thermal and storage capacity. +# +# ### Variables +# +# - Investment: $i \geq 0$ +# - Renewable production: $r_t \geq 0$ +# - Thermal production: $0 \leq p_t \leq \overline{P}$ +# - Storage level: $0 \leq s_t \leq \overline{S}$ +# - Storage charging: $0 \leq c_t \leq \overline{C}$ +# - Storage discharging: $0 \leq d_t \leq \overline{D}$ +# +# ### Parameters that will change each window +# +# - Demand at time $t$: $D_t$ +# - Renewable availability at time $t$: $A_t$ +# - Initial storage: $S_0$ +# +# ### Constraints +# +# 1. **Balance Constraint:** +# +# $p_t + r_t + d_t = D_t + c_t, \quad \forall t$ +# +# 2. **Storage Dynamics for $t \geq 2$:** +# +# $s_t = s_{t-1} + \eta^c \cdot c_t - \frac{d_t}{\eta^d}, \quad \forall t \in \{2, \ldots, \mathcal{T}\}$ +# +# 3. **Initial Storage:** +# +# $s_1 = S_0 +\eta^c \cdot c_1 - \frac{d_1}{\eta^d}$ +# +# 4. **Maximum Renewable Availability:** +# +# $r_t \leq A_t \cdot i, \quad \forall t$ +# +# ### Objective Function +# +# The objective function to minimize is the total cost: +# +# $\min \left(I \cdot i + \sum_{t} O \cdot p_t\right)$ + +# In large-scale optimization problems, this model is solved in two steps: +# +# 1. *Investment decisions step*: This involves a simplified version of the +# model, for example, without integer variables and representative periods. +# 2. *Operation decisions step*: After determining the values of the investments +# from the previous step, this step involves solving an operational problem +# to decide on production, storage levels, charging, and discharging. +# +# The second step is also computationally intensive. A common practice is to use +# the rolling horizon idea to solve multiple identical problems of a smaller +# size. These problems differ only in parameters such as demand, renewable +# profiles, or initial conditions. +# +# This example focuses on the second step, aiming to determine the operational +# variables for a given investment using a rolling horizon strategy. +# +# ## Parameter definition and input data +# +# There are two main parameters for a rolling horizon basic implementation: the +# optimization window and the move forward. +# +# **Optimization Window** (optimization_window): It defines how many periods +# (e.g., hours) we will optimize each time. For this example, we set the default +# value in 48h, meaning we will optimize two days each time. + +optimization_window = 48 + +# **Move Forward** (move_forward): It defines how many periods (e.g., hours) we +# will move forward to optimize the next optimization window. For this example, +# we set the default value in 24h, meaning we will move 1 day ahead each time. + +move_forward = 24 + +# Note that the move forward parameter must be lower or equal to the +# optimization window parameter to work correctly. + +@assert optimization_window >= move_forward "optimization_window must be greater or equal to move_forward" + +# Let's explore the input data in file [rolling_horizon.csv](rolling_horizon.csv). +# We have a total time horizon of a week (i.e., 168h), an electricity demand, +# and a solar production profile. + +filename = joinpath(@__DIR__, "rolling_horizon.csv") +time_series = CSV.read(filename, DataFrames.DataFrame); + +# We define the solar investment (e.g., 150 MW) to determine the solar +# production during the operation optimization step. +# +# In addition, we can determine some basic information about the rolling +# horizon, such as the number of windows that we are going to optimize given the +# problem's time horizon. +# +# Finally, we can see a plot representing the first two optimization windows and +# the move forward parameter to have a better idea of how the rolling horizon +# works. + +## Scale the solar profile +solar_investment = 150 +time_series.solar_MW = solar_investment * time_series.solar_pu + +## input data calculation for the Rolling Horizon +total_time_length = size(time_series, 1) +number_of_windows = ceil(Int, total_time_length / move_forward) +println("number of windows:", number_of_windows) + +#- +p1 = Plots.plot( + time_series.t, + [time_series.demand_MW, time_series.solar_MW]; + ylabel = "MW", + label = ["demand" "solar"], + color = [:ivory4 :darkorange1], + linewidth = 3, + ylims = (0, solar_investment), + xticks = 0:12:total_time_length, +) +Plots.vspan!(p1, [1, optimization_window]; alpha = 0.25, label = "") +Plots.annotate!( + p1, + 18, + time_series[12, :solar_MW], + Plots.text("optimization\n window 1", :top, :left, 8), +) +p2 = Plots.plot( + time_series.t, + [time_series.demand_MW, time_series.solar_MW]; + xlabel = "Hours", + ylabel = "MW", + label = ["" ""], + color = [:ivory4 :darkorange1], + linewidth = 3, + ylims = (0, solar_investment), + xticks = 0:12:total_time_length, +) +Plots.vspan!( + p2, + [move_forward, move_forward + optimization_window]; + alpha = 0.25, + label = "", +) +Plots.annotate!( + p2, + 42, + time_series[12, :solar_MW], + Plots.text("optimization\n window 2", :top, :left, 8), +) +Plots.plot!( + p2, + [1; move_forward], + [100; 100]; + arrow = 2, + label = "", + c = :black, +) +Plots.annotate!(p2, 5, 130, Plots.text(" move\nforward", :top, :left, 8)) +Plots.plot(p1, p2; layout = (2, 1)) + +# ## Rolling horizon first window +# +# We first sample the initial input data and get the parameter values for the +# first optimization window. +# +# We also create a helper index `t_minus_1` to get easy access to the previous +# hour using the function `mod1`. + +## Create data of the first window +time_series_filter = 1:optimization_window +availability = time_series.solar_pu[time_series_filter] +demand = time_series.demand_MW[time_series_filter] +## Create index t and t-1 +t = 1:optimization_window +t_minus_1 = + mod1.(optimization_window:(2*optimization_window-1), optimization_window) + +# We have all the information we need to create and optimize the first window in +# the model. + +model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) +set_silent(model) +@variable(model, 0 <= i) +@variable(model, 0 <= r[t]) +@variable(model, 0 <= p[t] <= 150) +@variable(model, 0 <= s[t] <= 40) +@variable(model, 0 <= c[t] <= 10) +@variable(model, 0 <= d[t] <= 10) +@variable(model, D[t] in Parameter.(demand[t])) +@variable(model, A[t] in Parameter.(availability[t])) +@variable(model, So in Parameter(0.0)) +@constraint(model, balance, p[t] .+ r[t] .+ d[t] .== D[t] .+ c[t]) +@constraint( + model, + storage[t in 2:optimization_window], + s[t] == s[t-1] + 0.9 * c[t] - d[t] / 0.9 +) +@constraint(model, init_storage, s[1] == So + 0.9 * c[1] - d[1] / 0.9) +@constraint(model, max_ava, r[t] .<= A[t] * i) +@objective(model, Min, 100 * i + sum(50 * p[t])) +fix(i, solar_investment; force = true) +optimize!(model) + +# After the optimization, we can store the results in vectors. It's important to +# note that despite optimizing for 48 hours (the default value), we only store +# the values for the "move forward" parameter (e.g., 24 hours or one day using +# the default value). This approach ensures that there is a buffer of additional +# periods or hours beyond the "move forward" parameter to prevent the storage +# from depleting entirely at the end of the specified hours. + +objective_function_per_window = zeros(number_of_windows) +renewable_production = zeros(total_time_length) +storage_level = zeros(total_time_length) + +# Store results from the first window +renewable_production[1:move_forward] = value.(model[:r])[1:move_forward] +storage_level[1:move_forward] = value.(model[:s])[1:move_forward] +objective_function_per_window[1] = objective_value(model) +println("Objective function first window: ", objective_function_per_window[1]) + +# ### Rolling horizon for the following windows +# +# For the following windows on the horizon, we: +# +# 1. Update the parameter values from the input data for that window +# 2. Update the parameters in the models using the ParametricOptInterface.jl +# 3. Solve the model for that window +# 4. Store the results +# +# Although this is a small problem, the benefits of using +# ParametricOptInterface.jl can be seen in the simplex iterations in the first +# window compared to the ones in the subsequent ones. + +for window in 2:number_of_windows + # Update window data + window_start = Int(1 + (window - 1) * move_forward) + window_end = Int(window_start + optimization_window - 1) + time_series_filter = mod1.(window_start:window_end, total_time_length) + availability = time_series.solar_pu[time_series_filter] + demand = time_series.demand_MW[time_series_filter] + initial_storage = storage_level[window_start-1] + # Update parameters in the model + MOI.set.(model, POI.ParameterValue(), model[:D], demand) + MOI.set.(model, POI.ParameterValue(), model[:A], availability) + MOI.set(model, POI.ParameterValue(), model[:So], initial_storage) + # Optimize again + optimize!(model) + # Store results for each window + window_end_output = + minimum([Int(window_start + move_forward - 1) total_time_length]) + output_filter = window_start:window_end_output + last_output_value = + minimum([move_forward (window_end_output - window_start + 1)]) + renewable_production[output_filter] = value.(model[:r])[1:last_output_value] + storage_level[output_filter] = value.(model[:s])[1:last_output_value] + objective_function_per_window[window] = objective_value(model) +end + +# We can explore the outputs in the following graphs: + +Plots.plot( + objective_function_per_window; + label = false, + linewidth = 3, + xlabel = "Window", + ylabel = "\$", + title = "Objective Function per Window", +) + +#- + +Plots.plot( + [time_series.demand_MW, renewable_production, storage_level]; + label = ["demand" "solar" "battery"], + linewidth = 3, + xlabel = "Hours", + ylabel = "MW", +) + +# **Final remark**: [ParametricOptInterface.jl](@ref) offers an easy way to +# update the parameter of an optimization problem that will be solved several +# times, as in the rolling horizon implementation. It has the benefit of +# avoiding rebuilding the model each time we want to solve it with new +# information in a new window.