Add AlternativeForwardPass

.gitignore

@@ -11,3 +11,4 @@ Manifest.toml
 docs/src/tutorial/*.json
 docs/src/*/*.ipynb
 docs/src/tutorial/spaghetti_plot.html
+convex.cuts.json

docs/Project.toml

@@ -11,6 +11,7 @@ JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655"
 SDDP = "f4570300-c277-11e8-125c-4912f86ce65d"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -24,3 +25,4 @@ JSON = "0.21"
 JuMP = "1"
 Literate = "2"
 Plots = "1"
+PowerModels = "0.19"

docs/make.jl

@@ -132,6 +132,7 @@ Documenter.makedocs(
             "tutorial/warnings.md",
             "tutorial/arma.md",
             "tutorial/objective_states.md",
+            "tutorial/pglib_opf.md",
             "tutorial/mdps.md",
             "tutorial/example_newsvendor.md",
             "tutorial/example_reservoir.md",

docs/src/apireference.md

@@ -77,6 +77,8 @@ SDDP.AbstractForwardPass
 SDDP.DefaultForwardPass
 SDDP.RevisitingForwardPass
 SDDP.RiskAdjustedForwardPass
+SDDP.AlternativeForwardPass
+SDDP.AlternativePostIterationCallback
 ```
 
 ### Risk Measures

docs/src/examples/air_conditioning_forward.jl

@@ -0,0 +1,42 @@
+#  Copyright (c) 2017-23, Oscar Dowson and SDDP.jl contributors.        #src
+#  This Source Code Form is subject to the terms of the Mozilla Public  #src
+#  License, v. 2.0. If a copy of the MPL was not distributed with this  #src
+#  file, You can obtain one at http://mozilla.org/MPL/2.0/.             #src
+
+# # Training with a different forward model
+
+using SDDP
+import HiGHS
+import Test
+
+function create_air_conditioning_model(; convex::Bool)
+    return SDDP.LinearPolicyGraph(
+        stages = 3,
+        lower_bound = 0.0,
+        optimizer = HiGHS.Optimizer,
+    ) do sp, t
+        @variable(sp, 0 <= x <= 100, SDDP.State, initial_value = 0)
+        @variable(sp, 0 <= u_production <= 200)
+        @variable(sp, u_overtime >= 0)
+        if !convex
+            set_integer(x.out)
+            set_integer(u_production)
+            set_integer(u_overtime)
+        end
+        @constraint(sp, demand, x.in - x.out + u_production + u_overtime == 0)
+        Ω = [[100.0], [100.0, 300.0], [100.0, 300.0]]
+        SDDP.parameterize(ω -> JuMP.set_normalized_rhs(demand, ω), sp, Ω[t])
+        @stageobjective(sp, 100 * u_production + 300 * u_overtime + 50 * x.out)
+    end
+end
+
+convex = create_air_conditioning_model(; convex = true)
+non_convex = create_air_conditioning_model(; convex = false)
+SDDP.train(
+    convex;
+    forward_pass = SDDP.AlternativeForwardPass(non_convex),
+    post_iteration_callback = SDDP.AlternativePostIterationCallback(non_convex),
+    iteration_limit = 10,
+)
+Test.@test isapprox(SDDP.calculate_bound(non_convex), 62_500.0, atol = 0.1)
+Test.@test isapprox(SDDP.calculate_bound(convex), 62_500.0, atol = 0.1)

docs/src/tutorial/pglib_opf.jl

@@ -0,0 +1,127 @@
+#  Copyright (c) 2017-23, Oscar Dowson and SDDP.jl contributors.        #src
+#  This Source Code Form is subject to the terms of the Mozilla Public  #src
+#  License, v. 2.0. If a copy of the MPL was not distributed with this  #src
+#  file, You can obtain one at http://mozilla.org/MPL/2.0/.             #src
+
+# # Alternative forward models
+
+# This example demonstrates how to train convex and non-convex models.
+
+# This example uses the following packages:
+
+using SDDP
+import Ipopt
+import PowerModels
+import Test
+
+# ## Formulation
+
+# For our model, we build a simple optimal power flow model with a single
+# hydro-electric generator.
+
+# The formulation of our optimal power flow problem depends on `model_type`,
+# which must be on of the `PowerModels` formulations.
+
+function build_model(model_type)
+    filename = joinpath(@__DIR__, "pglib_opf_case5_pjm.m")
+    data = PowerModels.parse_file(filename)
+    return SDDP.PolicyGraph(
+        SDDP.UnicyclicGraph(0.95);
+        sense = :Min,
+        lower_bound = 0.0,
+        optimizer = Ipopt.Optimizer,
+    ) do sp, t
+        power_model = PowerModels.instantiate_model(
+            data,
+            model_type,
+            PowerModels.build_opf;
+            jump_model = sp,
+        )
+        ## Now add hydro power models. Assume that generator 5 is hydro, and the
+        ## rest are thermal.
+        pg = power_model.var[:it][:pm][:nw][0][:pg][5]
+        sp[:pg] = pg
+        @variable(sp, x >= 0, SDDP.State, initial_value = 10.0)
+        @variable(sp, deficit >= 0)
+        @constraint(sp, balance, x.out == x.in - pg + deficit)
+        @stageobjective(sp, objective_function(sp) + 1e6 * deficit)
+        SDDP.parameterize(sp, [0, 2, 5]) do ω
+            return SDDP.set_normalized_rhs(balance, ω)
+        end
+        return
+    end
+end
+
+# ## Training a convex model
+
+# We can build and train a convex approximation of the optimal power flow
+# problem.
+
+# The problem with the convex model is that it does not accurately simulate the
+# true dynamics of the problem. Therefore, it under-estimates the true cost of
+# operation.
+
+convex = build_model(PowerModels.DCPPowerModel)
+SDDP.train(convex; iteration_limit = 10)
+
+# To more accurately simulate the dynamics of the problem, a common approach is
+# to write the cuts representing the policy to a file, and then read them into
+# a non-convex model:
+
+SDDP.write_cuts_to_file(convex, "convex.cuts.json")
+non_convex = build_model(PowerModels.ACPPowerModel)
+SDDP.read_cuts_from_file(non_convex, "convex.cuts.json")
+
+# Now we can simulate `non_convex` to evaluate the policy.
+
+result = SDDP.simulate(non_convex, 1)
+
+# A problem with reading and writing the cuts to file is that the cuts have been
+# generated from trial points of the convex model. Therefore, the policy may be
+# arbitrarily bad at points visited by the non-convex model.
+
+# ## Training a non-convex model
+
+# We can also build and train a non-convex formulation of the optimal power flow
+# problem.
+
+# The problem with the non-convex model is that because it is non-convex,
+# SDDP.jl may find a sub-optimal policy. Therefore, it may over-estimate the
+# true cost of operation.
+
+non_convex = build_model(PowerModels.ACPPowerModel)
+SDDP.train(non_convex; iteration_limit = 10)
+result = SDDP.simulate(non_convex, 1)
+
+# ## Combining convex and non-convex models
+
+# To summarize, training with the convex model constructs cuts at points that
+# may never be visited by the non-convex model, and training with the non-convex
+# model may construct arbitrarily poor cuts because a key assumption of SDDP is
+# convexity.
+
+# As a compromise, we can train a policy using a combination of the convex and
+# non-convex models; we'll use the non-convex model to generate trial points on
+# the forward pass, and we'll use the convex model to build cuts on the backward
+# pass.
+
+convex = build_model(PowerModels.DCPPowerModel)
+
+#-
+
+non_convex = build_model(PowerModels.ACPPowerModel)
+
+# To do so, we train `convex` using the [`SDDP.AlternativeForwardPass`](@ref)
+# forward pass, which simulates the model using `non_convex`, and we use
+# [`SDDP.AlternativePostIterationCallback`](@ref) as a post-iteration callback,
+# which copies cuts from the `convex` model back into the `non_convex` model.
+
+SDDP.train(
+    convex;
+    forward_pass = SDDP.AlternativeForwardPass(non_convex),
+    post_iteration_callback = SDDP.AlternativePostIterationCallback(non_convex),
+    iteration_limit = 10,
+)
+
+# In practice, if we were to simulate `non_convex` now, we should obtain a
+# better policy than either of the two previous approaches.

docs/src/tutorial/pglib_opf_case5_pjm.m

@@ -0,0 +1,116 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%                                                                  %%%%%
+%%%%    IEEE PES Power Grid Library - Optimal Power Flow - v21.07     %%%%%
+%%%%          (https://github.com/power-grid-lib/pglib-opf)           %%%%%
+%%%%               Benchmark Group - Typical Operations               %%%%%
+%%%%                         29 - July - 2021                         %%%%%
+%%%%                                                                  %%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%   CASE5  Power flow data for modified 5 bus, 5 gen case based on PJM 5-bus system
+%   Please see CASEFORMAT for details on the case file format.
+%
+%   Based on data from ...
+%     F.Li and R.Bo, "Small Test Systems for Power System Economic Studies",
+%     Proceedings of the 2010 IEEE Power & Energy Society General Meeting
+%
+%   Created by Rui Bo in 2006, modified in 2010, 2014.
+%
+%   Copyright (c) 2010 by The Institute of Electrical and Electronics Engineers (IEEE)
+%   Licensed under the Creative Commons Attribution 4.0
+%   International license, http://creativecommons.org/licenses/by/4.0/
+%
+%   Contact M.E. Brennan ([email protected]) for inquries on further reuse of
+%   this dataset.
+%
+function mpc = pglib_opf_case5_pjm
+mpc.version = '2';
+mpc.baseMVA = 100.0;
+
+%% area data
+%	area	refbus
+mpc.areas = [
+	1	 4;
+];
+
+%% bus data
+%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
+mpc.bus = [
+	1	 2	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 230.0	 1	    1.10000	    0.90000;
+	2	 1	 300.0	 98.61	 0.0	 0.0	 1	    1.00000	    0.00000	 230.0	 1	    1.10000	    0.90000;
+	3	 2	 300.0	 98.61	 0.0	 0.0	 1	    1.00000	    0.00000	 230.0	 1	    1.10000	    0.90000;
+	4	 3	 400.0	 131.47	 0.0	 0.0	 1	    1.00000	    0.00000	 230.0	 1	    1.10000	    0.90000;
+	5	 2	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 230.0	 1	    1.10000	    0.90000;
+];
+
+%% generator data
+%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin
+mpc.gen = [
+	1	 20.0	 0.0	 30.0	 -30.0	 1.0	 100.0	 1	 40.0	 0.0;
+	1	 85.0	 0.0	 127.5	 -127.5	 1.0	 100.0	 1	 170.0	 0.0;
+	3	 260.0	 0.0	 390.0	 -390.0	 1.0	 100.0	 1	 520.0	 0.0;
+	4	 100.0	 0.0	 150.0	 -150.0	 1.0	 100.0	 1	 200.0	 0.0;
+	5	 300.0	 0.0	 450.0	 -450.0	 1.0	 100.0	 1	 600.0	 0.0;
+];
+
+%% generator cost data
+%	2	startup	shutdown	n	c(n-1)	...	c0
+mpc.gencost = [
+	2	 0.0	 0.0	 3	   0.000000	  14.000000	   0.000000;
+	2	 0.0	 0.0	 3	   0.000000	  15.000000	   0.000000;
+	2	 0.0	 0.0	 3	   0.000000	  30.000000	   0.000000;
+	2	 0.0	 0.0	 3	   0.000000	  40.000000	   0.000000;
+	2	 0.0	 0.0	 3	   0.000000	  10.000000	   0.000000;
+];
+
+%% branch data
+%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
+mpc.branch = [
+	1	 2	 0.00281	 0.0281	 0.00712	 400.0	 400.0	 400.0	 0.0	 0.0	 1	 -30.0	 30.0;
+	1	 4	 0.00304	 0.0304	 0.00658	 426	 426	 426	 0.0	 0.0	 1	 -30.0	 30.0;
+	1	 5	 0.00064	 0.0064	 0.03126	 426	 426	 426	 0.0	 0.0	 1	 -30.0	 30.0;
+	2	 3	 0.00108	 0.0108	 0.01852	 426	 426	 426	 0.0	 0.0	 1	 -30.0	 30.0;
+	3	 4	 0.00297	 0.0297	 0.00674	 426	 426	 426	 0.0	 0.0	 1	 -30.0	 30.0;
+	4	 5	 0.00297	 0.0297	 0.00674	 240.0	 240.0	 240.0	 0.0	 0.0	 1	 -30.0	 30.0;
+];
+
+% INFO    : === Translation Options ===
+% INFO    : Phase Angle Bound:           30.0 (deg.)
+% INFO    : Line Capacity Model:         stat
+% INFO    : Setting Flat Start
+% INFO    : Line Capacity PAB:           15.0 (deg.)
+% INFO    :
+% INFO    : === Generator Bounds Update Notes ===
+% INFO    :
+% INFO    : === Base KV Replacement Notes ===
+% INFO    :
+% INFO    : === Transformer Setting Replacement Notes ===
+% INFO    :
+% INFO    : === Line Capacity Stat Model Notes ===
+% INFO    : Updated Thermal Rating: on line 1-4 : Rate A, Rate B, Rate C , 9900.0, 0.0, 0.0 -> 426
+% INFO    : Updated Thermal Rating: on line 1-5 : Rate A, Rate B, Rate C , 9900.0, 0.0, 0.0 -> 426
+% INFO    : Updated Thermal Rating: on line 2-3 : Rate A, Rate B, Rate C , 9900.0, 0.0, 0.0 -> 426
+% INFO    : Updated Thermal Rating: on line 3-4 : Rate A, Rate B, Rate C , 9900.0, 0.0, 0.0 -> 426
+% INFO    :
+% INFO    : === Line Capacity Monotonicity Notes ===
+% INFO    :
+% INFO    : === Voltage Setpoint Replacement Notes ===
+% INFO    : Bus 1	: V=1.0, theta=0.0 -> V=1.0, theta=0.0
+% INFO    : Bus 2	: V=1.0, theta=0.0 -> V=1.0, theta=0.0
+% INFO    : Bus 3	: V=1.0, theta=0.0 -> V=1.0, theta=0.0
+% INFO    : Bus 4	: V=1.0, theta=0.0 -> V=1.0, theta=0.0
+% INFO    : Bus 5	: V=1.0, theta=0.0 -> V=1.0, theta=0.0
+% INFO    :
+% INFO    : === Generator Setpoint Replacement Notes ===
+% INFO    : Gen at bus 1	: Pg=40.0, Qg=0.0 -> Pg=20.0, Qg=0.0
+% INFO    : Gen at bus 1	: Vg=1.0 -> Vg=1.0
+% INFO    : Gen at bus 1	: Pg=170.0, Qg=0.0 -> Pg=85.0, Qg=0.0
+% INFO    : Gen at bus 1	: Vg=1.0 -> Vg=1.0
+% INFO    : Gen at bus 3	: Pg=323.49, Qg=0.0 -> Pg=260.0, Qg=0.0
+% INFO    : Gen at bus 3	: Vg=1.0 -> Vg=1.0
+% INFO    : Gen at bus 4	: Pg=0.0, Qg=0.0 -> Pg=100.0, Qg=0.0
+% INFO    : Gen at bus 4	: Vg=1.0 -> Vg=1.0
+% INFO    : Gen at bus 5	: Pg=466.51, Qg=0.0 -> Pg=300.0, Qg=0.0
+% INFO    : Gen at bus 5	: Vg=1.0 -> Vg=1.0
+% INFO    :
+% INFO    : === Writing Matpower Case File Notes ===

src/SDDP.jl

@@ -62,6 +62,7 @@ include("visualization/value_functions.jl")
 # Other solvers.
 include("deterministic_equivalent.jl")
 include("biobjective.jl")
+include("alternative_forward.jl")
 
 include("Experimental.jl")
 include("MSPFormat.jl")

src/algorithm.jl

@@ -99,7 +99,7 @@ struct Options{T}
     duality_handler::AbstractDualityHandler
     # A callback called after the forward pass.
     forward_pass_callback::Any
-
+    post_iteration_callback::Any
     # Internal function: users should never construct this themselves.
     function Options(
         model::PolicyGraph{T},
@@ -119,6 +119,7 @@ struct Options{T}
         forward_pass::AbstractForwardPass = DefaultForwardPass(),
         duality_handler::AbstractDualityHandler = ContinuousConicDuality(),
         forward_pass_callback = x -> nothing,
+        post_iteration_callback = result -> nothing,
     ) where {T}
         return new{T}(
             initial_state,
@@ -140,6 +141,7 @@ struct Options{T}
             forward_pass,
             duality_handler,
             forward_pass_callback,
+            post_iteration_callback,
         )
     end
 end
@@ -913,6 +915,10 @@ Train the policy for `model`.
  - `duality_handler::AbstractDualityHandler`: specify a duality handler to use
    when creating cuts.
 
+ - `post_iteration_callback::Function`: a callback with the signature
+   `post_iteration_callback(::IterationResult)` that is evaluated after each
+   iteration of the algorithm.
+
 There is also a special option for infinite horizon problems
 
  - `cycle_discretization_delta`: the maximum distance between states allowed on
@@ -943,6 +949,7 @@ function train(
     add_to_existing_cuts::Bool = false,
     duality_handler::AbstractDualityHandler = SDDP.ContinuousConicDuality(),
     forward_pass_callback::Function = (x) -> nothing,
+    post_iteration_callback = result -> nothing,
 )
     function log_frequency_f(log::Vector{Log})
         if mod(length(log), log_frequency) != 0
@@ -1063,6 +1070,7 @@ function train(
         forward_pass,
         duality_handler,
         forward_pass_callback,
+        post_iteration_callback,
     )
     status = :not_solved
     try