From 4783914c1f0a4b8b1c9d40567c33f3d2655d93c4 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Sun, 1 Sep 2024 09:10:43 +1200 Subject: [PATCH] [docs] add Rolling horizon tutorial (#3815) Co-authored-by: datejada --- docs/Project.toml | 5 + docs/make.jl | 1 + .../tutorials/algorithms/rolling_horizon.csv | 169 ++++++++++ .../tutorials/algorithms/rolling_horizon.jl | 297 ++++++++++++++++++ .../config/vocabularies/JuMP/accept.txt | 2 + 5 files changed, 474 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..17943088cd2 --- /dev/null +++ b/docs/src/tutorials/algorithms/rolling_horizon.jl @@ -0,0 +1,297 @@ +# 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 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. +# +# As a motivating example, this tutorial models the operations of a power system +# with solar generation and a battery. + +# ## Required packages +# +# This tutorial uses the following packages + +using JuMP +import CSV +import DataFrames +import HiGHS +import ParametricOptInterface as POI +import Plots + +# ## The optimization model +# +# 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. +# +# There are five types of decision variables in the model: +# +# - 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}$ +# +# 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$ +# +# 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, \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}, \forall t$$ +# +# with the boundary condition that $s_0 = S_0$. +# +# 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 +# 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. + +# ## Parameter definition and input data + +# 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, +# 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) +time_series[1:21:end, :] + +# We define the solar investment (for example, 150 MW) to determine the solar +# production during the operation optimization step. + +solar_investment = 150; + +# 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; + +# In addition, we can determine some basic information about the rolling +# horizon, such as the number of data points we have: + +total_time_length = size(time_series, 1) + +# 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 + +# 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. + +# 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 + 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 later + D[t in 1:optimization_window] in Parameter(0) + A[t in 1:optimization_window] in Parameter(0) + S_0 in Parameter(0) +end) +@objective(model, Min, 50 * sum(p)) +@constraints( + model, + begin + 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 + 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. + +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: + +# 1. update the parameters in the models +# 2. solve the model for that window +# 3. store the results for later analysis + +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 + 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 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, 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 + sol_complete[:s][offset+t+1] = value(model[:s][t]) + end + sol_window = Dict(key => value.(model[key]) for key in (:r, :p, :s, :c, :d)) + push!(sol_windows, offset => sol_window) +end + +# ## Solution + +# Here is a function to plot the solution at each of the time-steps to help +# visualize the rolling horizon scheme: + +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( + [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 +# 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. 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