Update

docs/Project.toml

@@ -49,6 +49,7 @@ Documenter = "=1.4.1"
 DocumenterCitations = "1"
 Dualization = "0.5"
 Enzyme = "0.12.14"
+ForwardDiff = "0.10"
 GLPK = "=1.2.1"
 HTTP = "1.5.4"
 HiGHS = "=1.9.2"
@@ -57,11 +58,14 @@ Ipopt = "=1.6.5"
 JSON = "0.21"
 JSONSchema = "1"
 Literate = "2.8"
+MarkdownAST = "0.1"
 MathOptInterface = "=1.31.0"
 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"

docs/src/tutorials/algorithms/rolling_horizon.jl

@@ -105,19 +105,20 @@ import StatsPlots
 # 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
-# (for example, 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**: 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** (move_forward): It defines how many periods (for example, hours) we
+# **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 24h, meaning we will move 1 day ahead each time.
+# we set the default value in 24 hours, meaning we will move one day ahead each
+# time.
 
 move_forward = 24
 
@@ -127,32 +128,32 @@ move_forward = 24
 @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 (i.e., 168h), an electricity demand,
-# and a solar production profile.
+# 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.
 
-## 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)
-println("number of windows:", ceil(Int, total_time_length / move_forward))
-
-#-
 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"])
@@ -174,10 +175,10 @@ Plots.plot(
     ylabel = "MW",
 )
 
-# ## Rolling horizon first window
-#
-# We have all the information we need to create and optimize the first window in
-# the model.
+# ## 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)
@@ -193,6 +194,7 @@ set_silent(model)
     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
@@ -202,7 +204,6 @@ end)
         r .<= A .* i
     end
 )
-@objective(model, Min, 100 * i + 50 * sum(p))
 model
 
 # After the optimization, we can store the results in vectors. It's important to
@@ -213,32 +214,35 @@ model
 # storage from depleting entirely at the end of the specified hours.
 
 objective_function_per_window = Float64[]
-renewable_production = zeros(total_time_length)
-storage_level = zeros(total_time_length)
+renewable_production = Float64[]
+storage_level = Float64[0.0]  # Include an initial storage level
 
-# ### Rolling horizon for the following windows
-#
-# For the following windows on the horizon, we:
-#
-# 1. Update the parameters in the models using the ParametricOptInterface.jl
-# 2. Solve the model for that window
-# 3. Store the results
+# 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
-        row = mod1(offset + t, size(time_series, 1))
+        ## 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], get(storage_level, offset, 0))
+    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
-        renewable_production[offset+t] = value(model[:r][t])
-        storage_level[offset+t] = value(model[:s][t])
+        push!(renewable_production, value(model[:r][t]))
+        push!(storage_level, value(model[:s][t]))
     end
 end
 
@@ -255,16 +259,17 @@ Plots.plot(
 #-
 
 Plots.plot(
-    [time_series.demand_MW, renewable_production, storage_level];
+    [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 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.
+# ## 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.