From f33c9774f08c576319dec926fb95c76419d3c422 Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 30 Aug 2024 11:22:37 +1200 Subject: [PATCH 01/11] [docs] add Rolling horizon tutorial Co-authored-by: datejada --- docs/Project.toml | 5 + docs/make.jl | 1 + .../tutorials/algorithms/rolling_horizon.csv | 169 +++++++++++ .../tutorials/algorithms/rolling_horizon.jl | 275 ++++++++++++++++++ 4 files changed, 450 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 c4a317f2078..5c594ba70dd 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" @@ -48,6 +49,7 @@ Documenter = "=1.6.0" DocumenterCitations = "1" Dualization = "0.5" Enzyme = "0.12.14" +ForwardDiff = "0.10" GLPK = "=1.2.1" HTTP = "1.5.4" HiGHS = "=1.9.2" @@ -56,11 +58,14 @@ Ipopt = "=1.6.6" JSON = "0.21" JSONSchema = "1" Literate = "2.8" +MarkdownAST = "0.1" MathOptInterface = "=1.31.1" MultiObjectiveAlgorithms = "=1.3.3" PATHSolver = "=1.7.7" +ParametricOptInterface = "0.8.1" Plots = "1" SCS = "=2.0.1" SQLite = "1" +SpecialFunctions = "2" StatsPlots = "0.15" Tables = "1" diff --git a/docs/make.jl b/docs/make.jl index c52062d7bd3..54ea801b60b 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -377,6 +377,7 @@ const _PAGES = [ "tutorials/algorithms/benders_decomposition.md", "tutorials/algorithms/cutting_stock_column_generation.md", "tutorials/algorithms/tsp_lazy_constraints.md", + "tutorials/algorithms/rolling_horizon.md", "tutorials/algorithms/parallelism.md", ], "Applications" => [ diff --git a/docs/src/tutorials/algorithms/rolling_horizon.csv b/docs/src/tutorials/algorithms/rolling_horizon.csv new file mode 100644 index 00000000000..d74486700ee --- /dev/null +++ b/docs/src/tutorials/algorithms/rolling_horizon.csv @@ -0,0 +1,169 @@ +day,hour,demand_MW,solar_pu +01,00,51.6,0 +01,01,49.2,0 +01,02,46.5,0 +01,03,44.3,0 +01,04,43.3,0 +01,05,42.1,0 +01,06,39.8,0 +01,07,40.2,0 +01,08,41.3,0.212560386 +01,09,45,0.608695652 +01,10,49.3,0.845410628 +01,11,54.3,0.995169082 +01,12,56,1 +01,13,54.9,0.763285024 +01,14,53.3,0.309178744 +01,15,53.5,0.009661836 +01,16,57.5,0 +01,17,65,0 +01,18,66.2,0 +01,19,64.5,0 +01,20,61,0 +01,21,59,0 +01,22,58.7,0 +01,23,54.1,0 +02,00,49.7,0 +02,01,46.5,0 +02,02,44.8,0 +02,03,44.5,0 +02,04,46,0 +02,05,48.6,0 +02,06,52.6,0 +02,07,59,0 +02,08,65.1,0.096618357 +02,09,70.1,0.256038647 +02,10,73.5,0.391304348 +02,11,76.2,0.47826087 +02,12,76.8,0.531400966 +02,13,75.1,0.434782609 +02,14,73.2,0.202898551 +02,15,72.5,0.014492754 +02,16,75.2,0 +02,17,80.7,0 +02,18,80.7,0 +02,19,77.5,0 +02,20,71.3,0 +02,21,67.6,0 +02,22,65.8,0 +02,23,60.4,0 +03,00,54.7,0 +03,01,50.9,0 +03,02,48.5,0 +03,03,47.7,0 +03,04,48.2,0 +03,05,48.5,0 +03,06,49.1,0 +03,07,53.3,0 +03,08,58.9,0.09178744 +03,09,64.6,0.265700483 +03,10,68.8,0.367149758 +03,11,72,0.400966184 +03,12,72.4,0.347826087 +03,13,70.9,0.251207729 +03,14,69.5,0.111111111 +03,15,69.5,0.009661836 +03,16,72.5,0 +03,17,77.3,0 +03,18,77.4,0 +03,19,73.9,0 +03,20,68,0 +03,21,64.1,0 +03,22,62.8,0 +03,23,58.1,0 +04,00,52.8,0 +04,01,49.1,0 +04,02,47,0 +04,03,45.9,0 +04,04,46.1,0 +04,05,45.5,0 +04,06,44.1,0 +04,07,46.5,0.004830918 +04,08,50.3,0.256038647 +04,09,55.6,0.700483092 +04,10,60.3,0.888888889 +04,11,65.6,0.93236715 +04,12,65.9,0.787439614 +04,13,63.2,0.550724638 +04,14,60.7,0.275362319 +04,15,60.1,0.019323671 +04,16,63.4,0 +04,17,71.3,0 +04,18,73.1,0 +04,19,70.9,0 +04,20,66.8,0 +04,21,64.2,0 +04,22,63.9,0 +04,23,58.9,0 +05,00,54,0 +05,01,50.7,0 +05,02,49.4,0 +05,03,49.6,0 +05,04,51.7,0 +05,05,56.9,0 +05,06,66.2,0 +05,07,76.3,0.009661836 +05,08,82,0.29468599 +05,09,83.8,0.628019324 +05,10,85.9,0.777777778 +05,11,87.7,0.893719807 +05,12,87.7,0.874396135 +05,13,86.2,0.743961353 +05,14,84.7,0.444444444 +05,15,83.9,0.057971014 +05,16,85.9,0 +05,17,92,0 +05,18,92,0 +05,19,89,0 +05,20,82,0 +05,21,77.2,0 +05,22,74.1,0 +05,23,67,0 +06,00,61.8,0 +06,01,58,0 +06,02,56.3,0 +06,03,56.4,0 +06,04,57.7,0 +06,05,60.6,0 +06,06,67.4,0 +06,07,75.7,0.009661836 +06,08,79.7,0.256038647 +06,09,81.7,0.584541063 +06,10,84.2,0.821256039 +06,11,86.3,0.942028986 +06,12,86,0.884057971 +06,13,83.8,0.661835749 +06,14,81.5,0.328502415 +06,15,80.9,0.028985507 +06,16,83.8,0 +06,17,90.7,0 +06,18,90.7,0 +06,19,88.2,0 +06,20,82.1,0 +06,21,77.2,0 +06,22,73.9,0 +06,23,67.5,0 +07,00,61.8,0 +07,01,57.9,0 +07,02,56.9,0 +07,03,57.5,0 +07,04,59.2,0 +07,05,64.8,0 +07,06,77.9,0 +07,07,89.3,0.004830918 +07,08,94.1,0.154589372 +07,09,94.4,0.434782609 +07,10,95.9,0.589371981 +07,11,97.3,0.70531401 +07,12,96.7,0.647342995 +07,13,95.6,0.531400966 +07,14,93.7,0.265700483 +07,15,92.7,0.028985507 +07,16,94,0 +07,17,100,0 +07,18,99.2,0 +07,19,95.8,0 +07,20,88.9,0 +07,21,83.3,0 +07,22,79.2,0 +07,23,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..f6b59233ce0 --- /dev/null +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -0,0 +1,275 @@ +# 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**: this value defines how many periods (for example, +# hours) we will optimize each time. For this example, we set the default value +# to 48 hours, meaning we will optimize two days each time. + +optimization_window = 48 + +# **Move Forward**: this value defines how many periods (for example, hours) we +# will move forward to optimize the next optimization window. For this example, +# we set the default value in 24 hours, meaning we will move one 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 + +# Let's explore the input data in file [rolling_horizon.csv](rolling_horizon.csv). +# We have a total time horizon of a week (that is, 168 hours), 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 (for example, 150 MW) to determine the solar +# production during the operation optimization step. + +solar_investment = 150 +time_series.solar_MW = solar_investment * time_series.solar_pu + +# 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. + +total_time_length = size(time_series, 1) + +# The total number of time windows we will solve for is: + +ceil(Int, total_time_length / move_forward) + +# 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. + +x_series = 1:total_time_length +y_series = [time_series.demand_MW, time_series.solar_MW] +plot_1 = Plots.plot(x_series, y_series; label = ["demand" "solar"]) +plot_2 = Plots.plot(x_series, y_series; label = false) +window = [0, optimization_window] +Plots.vspan!(plot_1, window; alpha = 0.25, label = false) +Plots.vspan!(plot_2, move_forward .+ window; alpha = 0.25, label = false) +text_1 = Plots.text("optimization\n window 1", :top, :left, 8) +Plots.annotate!(plot_1, 18, time_series.solar_MW[12], text_1) +text_2 = Plots.text("optimization\n window 2", :top, :left, 8) +Plots.annotate!(plot_2, 42, time_series.solar_MW[12], text_2) +Plots.plot( + plot_1, + plot_2; + layout = (2, 1), + linewidth = 3, + xticks = 0:12:total_time_length, + xlabel = "Hours", + ylabel = "MW", +) + +# ## JuMP model + +# We have all the information we need to create a JuMP model to solve a single +# window of our rolling horizon problem. + +model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) +set_silent(model) +@variables(model, begin + i == solar_investment + 0 <= r[1:optimization_window] + 0 <= p[1:optimization_window] <= 150 + 0 <= s[1:optimization_window] <= 40 + 0 <= c[1:optimization_window] <= 10 + 0 <= d[1:optimization_window] <= 10 + ## Initialize empty parameters. These values will get updated layer + D[t in 1:optimization_window] in Parameter(0) + A[t in 1:optimization_window] in Parameter(0) + So in Parameter(0) +end) +@objective(model, Min, 100 * i + 50 * sum(p)) +@constraints( + model, + begin + p .+ r .+ d .== D .+ c + s[1] == So + 0.9 * c[1] - d[1] / 0.9 + [t in 2:optimization_window], s[t] == s[t-1] + 0.9 * c[t] - d[t] / 0.9 + r .<= A .* i + end +) +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 (for example, 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 = Float64[] +renewable_production = Float64[] +storage_level = Float64[0.0] # Include an initial storage level + +# Now we can iterate across the windows of our rolling horizon problem, and at +# each window, we: + +# 1. update the parameters in the models +# 2. solve the model for that window +# 3. store the results for later analysis + +for offset in 0:move_forward:total_time_length-1 + ## Step 1: update the parameter values over the optimization_window + for t in 1:optimization_window + ## This row computation just let's us "wrap around" the `time_series` + ## DataFrame, so that the forecase for demand and solar PU in day 8 is + ## the same as day 1. In real models, you might choose to do something + ## different. + row = 1 + mod(offset + t, size(time_series, 1)) + set_parameter_value(model[:D][t], time_series[row, :demand_MW]) + set_parameter_value(model[:A][t], time_series[row, :solar_pu]) + end + set_parameter_value(model[:So], storage_level[end]) + ## Step 2: solve the model + optimize!(model) + ## Step 3: store the results of the move_forward values + push!(objective_function_per_window, objective_value(model)) + for t in 1:move_forward + push!(renewable_production, value(model[:r][t])) + push!(storage_level, value(model[:s][t])) + end +end + +# We can explore the outputs in the following graphs: + +Plots.plot( + objective_function_per_window ./ 10^3; + label = false, + linewidth = 3, + xlabel = "Window", + ylabel = "[000'] \$", +) + +#- + +Plots.plot( + [time_series.demand_MW, renewable_production, storage_level[2:end]]; + label = ["demand" "solar" "battery"], + linewidth = 3, + xlabel = "Hours", + ylabel = "MW", + xticks = 0:12:total_time_length, +) + +# ## Final remark + +# [ParametricOptInterface.jl](@ref) offers an easy way to update the parameters +# 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. From e4a1be32cfd1328884b894626d0e475808145360 Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 30 Aug 2024 13:02:45 +1200 Subject: [PATCH 02/11] Update --- docs/src/tutorials/algorithms/rolling_horizon.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/docs/src/tutorials/algorithms/rolling_horizon.jl b/docs/src/tutorials/algorithms/rolling_horizon.jl index f6b59233ce0..11b9e534b57 100644 --- a/docs/src/tutorials/algorithms/rolling_horizon.jl +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -38,18 +38,16 @@ 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 +# ### [Variables](@id rolling_horizon_variables) # # - Investment: $i \geq 0$ # - Renewable production: $r_t \geq 0$ @@ -64,7 +62,7 @@ import StatsPlots # - Renewable availability at time $t$: $A_t$ # - Initial storage: $S_0$ # -# ### Constraints +# ### [Constraints](@id rolling_horizon_constraints) # # 1. **Balance Constraint:** # From feec2c81b5388775e43bd4d92b93a609a7d9bea1 Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 30 Aug 2024 13:50:40 +1200 Subject: [PATCH 03/11] Fix accept.txt --- docs/styles/config/vocabularies/JuMP/accept.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/styles/config/vocabularies/JuMP/accept.txt b/docs/styles/config/vocabularies/JuMP/accept.txt index 5c73dbd1489..499869f1737 100644 --- a/docs/styles/config/vocabularies/JuMP/accept.txt +++ b/docs/styles/config/vocabularies/JuMP/accept.txt @@ -42,6 +42,7 @@ preprint README recurse reimplemented +renewables runtime(?s) [Ss]tacktrace subexpression(?s) @@ -256,6 +257,7 @@ Sungho Taccari Tanneau Teghem +Tejada Tillmann Ulungu Vandenberghe From f6501e971ca2d90c4d348dbd97f48279a2eb8634 Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 30 Aug 2024 13:54:29 +1200 Subject: [PATCH 04/11] Remove objective_function_per_window --- docs/src/tutorials/algorithms/rolling_horizon.jl | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/docs/src/tutorials/algorithms/rolling_horizon.jl b/docs/src/tutorials/algorithms/rolling_horizon.jl index 11b9e534b57..30bdeef8805 100644 --- a/docs/src/tutorials/algorithms/rolling_horizon.jl +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -211,7 +211,6 @@ model # 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 = Float64[] renewable_production = Float64[] storage_level = Float64[0.0] # Include an initial storage level @@ -237,7 +236,6 @@ for offset in 0:move_forward:total_time_length-1 ## Step 2: solve the model optimize!(model) ## Step 3: store the results of the move_forward values - push!(objective_function_per_window, objective_value(model)) for t in 1:move_forward push!(renewable_production, value(model[:r][t])) push!(storage_level, value(model[:s][t])) @@ -246,16 +244,6 @@ end # We can explore the outputs in the following graphs: -Plots.plot( - objective_function_per_window ./ 10^3; - label = false, - linewidth = 3, - xlabel = "Window", - ylabel = "[000'] \$", -) - -#- - Plots.plot( [time_series.demand_MW, renewable_production, storage_level[2:end]]; label = ["demand" "solar" "battery"], From dd2ae58d73951b640c6049ff1890af1aa5394b8c Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 30 Aug 2024 14:35:20 +1200 Subject: [PATCH 05/11] Update --- .../tutorials/algorithms/rolling_horizon.jl | 98 +++++++++---------- 1 file changed, 44 insertions(+), 54 deletions(-) diff --git a/docs/src/tutorials/algorithms/rolling_horizon.jl b/docs/src/tutorials/algorithms/rolling_horizon.jl index 30bdeef8805..49f53d789c4 100644 --- a/docs/src/tutorials/algorithms/rolling_horizon.jl +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -22,14 +22,15 @@ # # **This tutorial was originally contributed by Diego Tejada.** # +# The purpose of this tutorial is to demonstrate how to use [ParametricOptInterface.jl](@ref) +# to solve a rolling horizon optimization problem. +# # 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). +# As a motivating example, this tutorial models the operations of a power system +# with solar generation and a battery. # ## Required packages # @@ -44,67 +45,53 @@ import Plots # ## 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. +# The model is a simplified model of a power system's operations with battery +# storage. +# +# We model the system of a set of time-steps $t \in 1,\ldots,T$, where each time +# step is a period of one hour. # -# ### [Variables](@id rolling_horizon_variables) +# There are five types of decision variables in the model: # -# - 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 +# For the purpose of this tutorial, there are three parameters of interest: # # - Demand at time $t$: $D_t$ # - Renewable availability at time $t$: $A_t$ # - Initial storage: $S_0$ # -# ### [Constraints](@id rolling_horizon_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:** +# The objective function to minimize is the total cost of thermal generation: +# $$\min \sum_{t} O \cdot p_t$$ # -# $r_t \leq A_t \cdot i, \quad \forall t$ +# For the constraints, we must balance power generation and consumption in all +# time periods: +# $$p_t + r_t + d_t = D_t + c_t, \quad \forall t$$ # -# ### Objective Function +# We need to account for the dynamics of the battery storage: +# $$s_t = s_{t-1} + \eta^c \cdot c_t - \frac{d_t}{\eta^d}, \quad \forall t$$ +# with the boundary condition that $s_0 = S_0$. # -# The objective function to minimize is the total cost: -# -# $\min \left(I \cdot i + \sum_{t} O \cdot p_t\right)$ +# Finally, the level of renewable energy production is limited by the +# availability factor $A$ and the installed capacity $i$: +# $$r_t \leq A_t \cdot i, \quad \forall t$$ + +# Solving this problem with a large number of time steps is computationally +# challenging. 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 availability, and initial storage. By +# combining the solution of many smaller problems, we can recover a feasible +# solution to the full problem. However, because we don't optimize the full set +# of decisions in a single optimization problem, the recovered solution might be +# suboptimal. -# 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 +# There are two main parameters for a rolling horizon implementation: the # optimization window and the move forward. # **Optimization Window**: this value defines how many periods (for example, @@ -136,7 +123,11 @@ time_series = CSV.read(filename, DataFrames.DataFrame); # production during the operation optimization step. solar_investment = 150 -time_series.solar_MW = solar_investment * time_series.solar_pu + +# We multiple the level of solar investment by the time series of availability +# to get actual MW generated. + +time_series.solar_MW = solar_investment * time_series.solar_pu; # 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 @@ -181,7 +172,6 @@ Plots.plot( model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) set_silent(model) @variables(model, begin - i == solar_investment 0 <= r[1:optimization_window] 0 <= p[1:optimization_window] <= 150 0 <= s[1:optimization_window] <= 40 @@ -190,16 +180,16 @@ set_silent(model) ## Initialize empty parameters. These values will get updated layer D[t in 1:optimization_window] in Parameter(0) A[t in 1:optimization_window] in Parameter(0) - So in Parameter(0) + S_0 in Parameter(0) end) -@objective(model, Min, 100 * i + 50 * sum(p)) +@objective(model, Min, 50 * sum(p)) @constraints( model, begin p .+ r .+ d .== D .+ c - s[1] == So + 0.9 * c[1] - d[1] / 0.9 + s[1] == S_0 + 0.9 * c[1] - d[1] / 0.9 [t in 2:optimization_window], s[t] == s[t-1] + 0.9 * c[t] - d[t] / 0.9 - r .<= A .* i + r .<= A .* solar_investment end ) model @@ -232,7 +222,7 @@ for offset in 0:move_forward:total_time_length-1 set_parameter_value(model[:D][t], time_series[row, :demand_MW]) set_parameter_value(model[:A][t], time_series[row, :solar_pu]) end - set_parameter_value(model[:So], storage_level[end]) + set_parameter_value(model[:S_0], storage_level[end]) ## Step 2: solve the model optimize!(model) ## Step 3: store the results of the move_forward values @@ -242,7 +232,7 @@ for offset in 0:move_forward:total_time_length-1 end end -# We can explore the outputs in the following graphs: +# We can now plot the solution to the week-long problem: Plots.plot( [time_series.demand_MW, renewable_production, storage_level[2:end]]; From a49d1798a5c2457942ac49662d87531f9726bb08 Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 30 Aug 2024 16:15:56 +1200 Subject: [PATCH 06/11] Update --- .../tutorials/algorithms/rolling_horizon.jl | 32 +++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/docs/src/tutorials/algorithms/rolling_horizon.jl b/docs/src/tutorials/algorithms/rolling_horizon.jl index 49f53d789c4..006b5d0e062 100644 --- a/docs/src/tutorials/algorithms/rolling_horizon.jl +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -204,6 +204,28 @@ model renewable_production = Float64[] storage_level = Float64[0.0] # Include an initial storage level +# We'll also plot the solution at each of the time-steps to help visualize the +# solution to the rolling horizon problems. + +function plot_solution(model, offset) + plot = Plots.plot(; + ylabel = "MW", + xlims = (0, total_time_length), + xticks = 0:12:total_time_length, + ) + x = offset .+ (1:optimization_window) + y = hcat(value.(model[:p]), value.(model[:r]), value.(model[:d])) + if offset == 0 + Plots.areaplot!(x, y; label = ["thermal" "solar" "discharge"]) + Plots.areaplot!(x, -value.(model[:c]); label = "charge") + else + Plots.areaplot!(x, y; label = false) + Plots.areaplot!(x, -value.(model[:c]); label = false) + end + return plot +end +plots = Any[] + # Now we can iterate across the windows of our rolling horizon problem, and at # each window, we: @@ -230,6 +252,7 @@ for offset in 0:move_forward:total_time_length-1 push!(renewable_production, value(model[:r][t])) push!(storage_level, value(model[:s][t])) end + push!(plots, plot_solution(model, offset)) end # We can now plot the solution to the week-long problem: @@ -243,6 +266,15 @@ Plots.plot( xticks = 0:12:total_time_length, ) +# and visualize each of the rolling horizon subplots: + +Plots.plot( + plots...; + layout = (length(plots), 1), + size = (600, 800), + margin = 3Plots.mm, +) + # ## Final remark # [ParametricOptInterface.jl](@ref) offers an easy way to update the parameters From 4902d527ed97d558226a311b3a2bac339af52380 Mon Sep 17 00:00:00 2001 From: odow Date: Sat, 31 Aug 2024 15:34:51 +1200 Subject: [PATCH 07/11] update --- .../tutorials/algorithms/rolling_horizon.jl | 112 ++++++++++-------- 1 file changed, 60 insertions(+), 52 deletions(-) diff --git a/docs/src/tutorials/algorithms/rolling_horizon.jl b/docs/src/tutorials/algorithms/rolling_horizon.jl index 006b5d0e062..94a44f66220 100644 --- a/docs/src/tutorials/algorithms/rolling_horizon.jl +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -66,19 +66,24 @@ import Plots # - Initial storage: $S_0$ # # The objective function to minimize is the total cost of thermal generation: +# # $$\min \sum_{t} O \cdot p_t$$ # # For the constraints, we must balance power generation and consumption in all # time periods: -# $$p_t + r_t + d_t = D_t + c_t, \quad \forall t$$ +# +# $$p_t + r_t + d_t = D_t + c_t \forall t$$ # # We need to account for the dynamics of the battery storage: -# $$s_t = s_{t-1} + \eta^c \cdot c_t - \frac{d_t}{\eta^d}, \quad \forall t$$ +# +# $$s_t = s_{t-1} + \eta^c \cdot c_t - \frac{d_t}{\eta^d}, \forall t$$ +# # with the boundary condition that $s_0 = S_0$. # -# Finally, the level of renewable energy production is limited by the -# availability factor $A$ and the installed capacity $i$: -# $$r_t \leq A_t \cdot i, \quad \forall t$$ +# Finally, the level of renewable energy production is limited by the quantity +# of potential solar generation $A$: +# +# $$r_t \leq A_t, \quad \forall t$$ # Solving this problem with a large number of time steps is computationally # challenging. A common practice is to use the rolling horizon idea to solve @@ -98,14 +103,14 @@ import Plots # hours) we will optimize each time. For this example, we set the default value # to 48 hours, meaning we will optimize two days each time. -optimization_window = 48 +optimization_window = 48; # **Move Forward**: this value defines how many periods (for example, hours) we # will move forward to optimize the next optimization window. For this example, # we set the default value in 24 hours, meaning we will move one day ahead each # time. -move_forward = 24 +move_forward = 24; # Note that the move forward parameter must be lower or equal to the # optimization window parameter to work correctly. @@ -117,12 +122,12 @@ move_forward = 24 # demand, and a solar production profile. filename = joinpath(@__DIR__, "rolling_horizon.csv") -time_series = CSV.read(filename, DataFrames.DataFrame); +time_series = CSV.read(filename, DataFrames.DataFrame) # We define the solar investment (for example, 150 MW) to determine the solar # production during the operation optimization step. -solar_investment = 150 +solar_investment = 150; # We multiple the level of solar investment by the time series of availability # to get actual MW generated. @@ -189,7 +194,7 @@ end) p .+ r .+ d .== D .+ c s[1] == S_0 + 0.9 * c[1] - d[1] / 0.9 [t in 2:optimization_window], s[t] == s[t-1] + 0.9 * c[t] - d[t] / 0.9 - r .<= A .* solar_investment + r .<= A end ) model @@ -201,30 +206,15 @@ model # additional periods or hours beyond the "move forward" parameter to prevent the # storage from depleting entirely at the end of the specified hours. -renewable_production = Float64[] -storage_level = Float64[0.0] # Include an initial storage level - -# We'll also plot the solution at each of the time-steps to help visualize the -# solution to the rolling horizon problems. - -function plot_solution(model, offset) - plot = Plots.plot(; - ylabel = "MW", - xlims = (0, total_time_length), - xticks = 0:12:total_time_length, - ) - x = offset .+ (1:optimization_window) - y = hcat(value.(model[:p]), value.(model[:r]), value.(model[:d])) - if offset == 0 - Plots.areaplot!(x, y; label = ["thermal" "solar" "discharge"]) - Plots.areaplot!(x, -value.(model[:c]); label = "charge") - else - Plots.areaplot!(x, y; label = false) - Plots.areaplot!(x, -value.(model[:c]); label = false) - end - return plot -end -plots = Any[] +sol_complete = Dict( + :r => zeros(total_time_length), + :p => zeros(total_time_length), + :c => zeros(total_time_length), + :d => zeros(total_time_length), + ## The storage level is initialized with an initial value + :s => zeros(total_time_length + 1), +) +sol_windows = Pair{Int,Dict{Symbol,Vector{Float64}}}[] # Now we can iterate across the windows of our rolling horizon problem, and at # each window, we: @@ -237,44 +227,62 @@ for offset in 0:move_forward:total_time_length-1 ## Step 1: update the parameter values over the optimization_window for t in 1:optimization_window ## This row computation just let's us "wrap around" the `time_series` - ## DataFrame, so that the forecase for demand and solar PU in day 8 is + ## DataFrame, so that the forecast for demand and solar MW in day 8 is ## the same as day 1. In real models, you might choose to do something ## different. row = 1 + mod(offset + t, size(time_series, 1)) set_parameter_value(model[:D][t], time_series[row, :demand_MW]) - set_parameter_value(model[:A][t], time_series[row, :solar_pu]) + set_parameter_value(model[:A][t], time_series[row, :solar_MW]) end - set_parameter_value(model[:S_0], storage_level[end]) + set_parameter_value(model[:S_0], sol_complete[:s][end]) ## Step 2: solve the model optimize!(model) ## Step 3: store the results of the move_forward values for t in 1:move_forward - push!(renewable_production, value(model[:r][t])) - push!(storage_level, value(model[:s][t])) + for key in (:r, :p, :c, :d) + sol_complete[key][offset+t] = value(model[key][t]) + end + sol_complete[:s][offset+t+1] = value(model[:s][t]) end - push!(plots, plot_solution(model, offset)) + sol_window = Dict(key => value.(model[key]) for key in (:r, :p, :s, :c, :d)) + push!(sol_windows, offset => sol_window) end -# We can now plot the solution to the week-long problem: +# ## Solution -Plots.plot( - [time_series.demand_MW, renewable_production, storage_level[2:end]]; - label = ["demand" "solar" "battery"], - linewidth = 3, - xlabel = "Hours", - ylabel = "MW", - xticks = 0:12:total_time_length, -) +# Here is a function to plot the solution at each of the time-steps to help +# visualize the rolling horizon scheme: -# and visualize each of the rolling horizon subplots: +function plot_solution(sol; offset = 0, kwargs...) + plot = Plots.plot(; + ylabel = "MW", + xlims = (0, total_time_length), + xticks = 0:12:total_time_length, + kwargs..., + ) + y = hcat(sol[:p], sol[:r], sol[:d]) + x = offset .+ (1:size(y, 1)) + if offset == 0 + Plots.areaplot!(x, y; label = ["thermal" "solar" "discharge"]) + Plots.areaplot!(x, -sol[:c]; label = "charge") + else + Plots.areaplot!(x, y; label = false) + Plots.areaplot!(x, -sol[:c]; label = false) + end + return plot +end Plots.plot( - plots...; - layout = (length(plots), 1), + [plot_solution(sol; offset) for (offset, sol) in sol_windows]...; + layout = (length(sol_windows), 1), size = (600, 800), margin = 3Plots.mm, ) +# We can re-use the function to plot the recovered solution of the full problem: + +plot_solution(sol_complete; offset = 0, xlabel = "Hour") + # ## Final remark # [ParametricOptInterface.jl](@ref) offers an easy way to update the parameters From bb0fbb1dab91224033fb4d03010ae4c1d33f0111 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Sat, 31 Aug 2024 16:31:29 +1200 Subject: [PATCH 08/11] Update docs/src/tutorials/algorithms/rolling_horizon.jl --- docs/src/tutorials/algorithms/rolling_horizon.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/tutorials/algorithms/rolling_horizon.jl b/docs/src/tutorials/algorithms/rolling_horizon.jl index 94a44f66220..4ef9f1a7e24 100644 --- a/docs/src/tutorials/algorithms/rolling_horizon.jl +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -123,6 +123,7 @@ move_forward = 24; filename = joinpath(@__DIR__, "rolling_horizon.csv") time_series = CSV.read(filename, DataFrames.DataFrame) +time_series[1:21:end, :] # We define the solar investment (for example, 150 MW) to determine the solar # production during the operation optimization step. From 67098a5419977b1fbaae0f531c6faa25e6c40bde Mon Sep 17 00:00:00 2001 From: odow Date: Sat, 31 Aug 2024 19:53:58 +1200 Subject: [PATCH 09/11] Update --- .../tutorials/algorithms/rolling_horizon.jl | 37 +++++++++++-------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/docs/src/tutorials/algorithms/rolling_horizon.jl b/docs/src/tutorials/algorithms/rolling_horizon.jl index 4ef9f1a7e24..7c5a0c7146c 100644 --- a/docs/src/tutorials/algorithms/rolling_horizon.jl +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -72,7 +72,7 @@ import Plots # For the constraints, we must balance power generation and consumption in all # time periods: # -# $$p_t + r_t + d_t = D_t + c_t \forall t$$ +# $$p_t + r_t + d_t = D_t + c_t, \forall t$$ # # We need to account for the dynamics of the battery storage: # @@ -136,14 +136,14 @@ solar_investment = 150; time_series.solar_MW = solar_investment * time_series.solar_pu; # 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. +# horizon, such as the number of data points we have: total_time_length = size(time_series, 1) -# The total number of time windows we will solve for is: +# adn the number of windows that we are going to optimize given the problem's +# time horizon: -ceil(Int, total_time_length / move_forward) +(total_time_length + move_forward - optimization_window) / move_forward # 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 @@ -175,6 +175,12 @@ Plots.plot( # We have all the information we need to create a JuMP model to solve a single # window of our rolling horizon problem. +# As the optimizer, we use `POI.Optimizer`, which is part of +# [ParametricOptInterface.jl](@ref). `POI.Optimizer` converts the [`Parameter`](@ref) +# decision variables into constants in the underlying optimization model, and it +# efficiently updates the solver in-place when we call [`set_parameter_value`](@ref) +# which avoids having to rebuild the problem each time we call [`optimize!`](@ref). + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) set_silent(model) @variables(model, begin @@ -224,22 +230,21 @@ sol_windows = Pair{Int,Dict{Symbol,Vector{Float64}}}[] # 2. solve the model for that window # 3. store the results for later analysis -for offset in 0:move_forward:total_time_length-1 +offsets = 0:move_forward:total_time_length-optimization_window +for offset in offsets ## Step 1: update the parameter values over the optimization_window for t in 1:optimization_window - ## This row computation just let's us "wrap around" the `time_series` - ## DataFrame, so that the forecast for demand and solar MW in day 8 is - ## the same as day 1. In real models, you might choose to do something - ## different. - row = 1 + mod(offset + t, size(time_series, 1)) - set_parameter_value(model[:D][t], time_series[row, :demand_MW]) - set_parameter_value(model[:A][t], time_series[row, :solar_MW]) + set_parameter_value(model[:D][t], time_series[offset+t, :demand_MW]) + set_parameter_value(model[:A][t], time_series[offset+t, :solar_MW]) end - set_parameter_value(model[:S_0], sol_complete[:s][end]) + ## Set the starting storage level as the value from the end of the previous + ## solve. The `+1` accounts for the initial storage value in time step "t=0" + set_parameter_value(model[:S_0], sol_complete[:s][offset+1]) ## Step 2: solve the model optimize!(model) - ## Step 3: store the results of the move_forward values - for t in 1:move_forward + ## Step 3: store the results of the move_forward values, except in the last + ## horizon where we store the full `optimization_window`. + for t in 1:(offset == last(offsets) ? optimization_window : move_forward) for key in (:r, :p, :c, :d) sol_complete[key][offset+t] = value(model[key][t]) end From 0d330178cae8124fca4621397ac626d77d764b78 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Sat, 31 Aug 2024 20:49:13 +1200 Subject: [PATCH 10/11] Update docs/src/tutorials/algorithms/rolling_horizon.jl --- docs/src/tutorials/algorithms/rolling_horizon.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/tutorials/algorithms/rolling_horizon.jl b/docs/src/tutorials/algorithms/rolling_horizon.jl index 7c5a0c7146c..45de2fe509f 100644 --- a/docs/src/tutorials/algorithms/rolling_horizon.jl +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -140,7 +140,7 @@ time_series.solar_MW = solar_investment * time_series.solar_pu; total_time_length = size(time_series, 1) -# adn the number of windows that we are going to optimize given the problem's +# and the number of windows that we are going to optimize given the problem's # time horizon: (total_time_length + move_forward - optimization_window) / move_forward From 3772b30d3e8e31192b21afa717bc643b3f8cb49e Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Sun, 1 Sep 2024 08:16:18 +1200 Subject: [PATCH 11/11] Apply suggestions from code review --- docs/src/tutorials/algorithms/rolling_horizon.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorials/algorithms/rolling_horizon.jl b/docs/src/tutorials/algorithms/rolling_horizon.jl index 45de2fe509f..17943088cd2 100644 --- a/docs/src/tutorials/algorithms/rolling_horizon.jl +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -130,7 +130,7 @@ time_series[1:21:end, :] solar_investment = 150; -# We multiple the level of solar investment by the time series of availability +# We multiply the level of solar investment by the time series of availability # to get actual MW generated. time_series.solar_MW = solar_investment * time_series.solar_pu; @@ -189,7 +189,7 @@ set_silent(model) 0 <= s[1:optimization_window] <= 40 0 <= c[1:optimization_window] <= 10 0 <= d[1:optimization_window] <= 10 - ## Initialize empty parameters. These values will get updated layer + ## Initialize empty parameters. These values will get updated later D[t in 1:optimization_window] in Parameter(0) A[t in 1:optimization_window] in Parameter(0) S_0 in Parameter(0)