From 0b83293b734ff94dd45e8414ae8f0b46551d7b0e Mon Sep 17 00:00:00 2001 From: odow Date: Sun, 14 May 2023 11:25:58 +1200 Subject: [PATCH] Add PowerModels example --- Project.toml | 4 +- docs/Project.toml | 1 + docs/src/apireference.md | 2 + docs/src/examples/air_conditioning_forward.jl | 7 +- docs/src/examples/pglib_opf.jl | 89 ++++++++++++++ docs/src/examples/pglib_opf_case5_pjm.m | 116 ++++++++++++++++++ src/algorithm.jl | 10 +- src/alternative_forward.jl | 81 ++++++------ src/plugins/parallel_schemes.jl | 3 + 9 files changed, 269 insertions(+), 44 deletions(-) create mode 100644 docs/src/examples/pglib_opf.jl create mode 100644 docs/src/examples/pglib_opf_case5_pjm.m diff --git a/Project.toml b/Project.toml index e44c6cde2..09a1e820c 100644 --- a/Project.toml +++ b/Project.toml @@ -21,6 +21,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] HTTP = "1" HiGHS = "1" +Ipopt = "1" JSON = "0.21" JSONSchema = "1" JuMP = "1" @@ -33,9 +34,10 @@ julia = "1.6" [extras] Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Downloads", "HiGHS", "JSONSchema", "Pkg", "Test"] +test = ["Downloads", "HiGHS", "Ipopt", "JSONSchema", "Pkg", "Test"] diff --git a/docs/Project.toml b/docs/Project.toml index 46c49a74d..6fddee964 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/src/apireference.md b/docs/src/apireference.md index 4c3bf90d8..26ed0823e 100644 --- a/docs/src/apireference.md +++ b/docs/src/apireference.md @@ -77,6 +77,8 @@ SDDP.AbstractForwardPass SDDP.DefaultForwardPass SDDP.RevisitingForwardPass SDDP.RiskAdjustedForwardPass +SDDP.AlternativeForwardPass +SDDP.AlternativePostIterationCallback ``` ### Risk Measures diff --git a/docs/src/examples/air_conditioning_forward.jl b/docs/src/examples/air_conditioning_forward.jl index 2105ccfd7..bc1be57a5 100644 --- a/docs/src/examples/air_conditioning_forward.jl +++ b/docs/src/examples/air_conditioning_forward.jl @@ -32,6 +32,11 @@ end convex = create_air_conditioning_model(; convex = true) non_convex = create_air_conditioning_model(; convex = false) -SDDP.train_with_forward_model(non_convex, convex; iteration_limit = 10) +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) diff --git a/docs/src/examples/pglib_opf.jl b/docs/src/examples/pglib_opf.jl new file mode 100644 index 000000000..633bcc920 --- /dev/null +++ b/docs/src/examples/pglib_opf.jl @@ -0,0 +1,89 @@ +# 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 teh 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) + +# ## Training a non-convex model + +# We can also build and train the true 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) + +# ## Combining convex and non-convex models + +# As a compromise, we can combine the convex and non-convex models. +convex = build_model(PowerModels.DCPPowerModel) +non_convex = build_model(filename, PowerModels.ACPPowerModel) +SDDP.train( + convex; + forward_pass = SDDP.AlternativeForwardPass(non_convex), + post_iteration_callback = SDDP.AlternativePostIterationCallback(non_convex), + iteration_limit = 10, +) diff --git a/docs/src/examples/pglib_opf_case5_pjm.m b/docs/src/examples/pglib_opf_case5_pjm.m new file mode 100644 index 000000000..7fa81a5a5 --- /dev/null +++ b/docs/src/examples/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 (me.brennan@ieee.org) 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 === diff --git a/src/algorithm.jl b/src/algorithm.jl index cf166bb1f..58502d8b3 100644 --- a/src/algorithm.jl +++ b/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 diff --git a/src/alternative_forward.jl b/src/alternative_forward.jl index 9cbcf2d98..9697389ac 100644 --- a/src/alternative_forward.jl +++ b/src/alternative_forward.jl @@ -3,21 +3,41 @@ # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. """ - AlternativeForwardPass(; - forward_model::SDDP.PolicyGraph{T}, + AlternativeForwardPass( + forward_model::SDDP.PolicyGraph{T}; + forward_pass::AbstractForwardPass = DefaultForwardPass(), ) -A forward pass for holding an alternative foward model. This is useful for -calculating cuts and simulating the policy under different models. -For example, simulating the forward passes with the AC-OPF to obtain forward trajectories, -and then refine the value functions on the backward pass using some convex approximation -(e.g., DC with line losses) as in the following paper: - - Rosemberg, Andrew Rosemberg, Alexandre Street, Joaquim Dias Garcia, Davi M. Valladão, Thuener Silva, and Oscar Dowson. -"Assessing the cost of network simplifications in long-term hydrothermal dispatch planning models." -IEEE Transactions on Sustainable Energy 13, no. 1 (2021): 196-206. +A forward pass that simulates using `forward_model`, which may be different to +the model used in the backwards pass. + +When using this forward pass, you should almost always pass +[`SDDP.AlternativePostIterationCallback`](@ref) to the `post_iteration_callback` +argument of [`SDDP.train`](@ref). + +This forward pass is most useful when the `forward_model` is non-convex and we +use a convex approximation of the model in the backward pass. + +For example, in optimal power flow models, we can use an AC-OPF formulation as +the `forward_model` and a DC-OPF formulation as the backward model. + +For more details see the paper: + +Rosemberg, A., and Street, A., and Garcia, J.D., and Valladão, D.M., and Silva, +T., and Dowson, O. (2021). Assessing the cost of network simplifications in +long-term hydrothermal dispatch planning models. IEEE Transactions on +Sustainable Energy. 13(1), 196-206. """ struct AlternativeForwardPass{T} <: AbstractForwardPass model::PolicyGraph{T} + forward_pass::AbstractForwardPass + + function AlternativeForwardPass( + model::PolicyGraph{T}; + forward_pass::AbstractForwardPass = DefaultForwardPass(), + ) where {T} + return new{T}(model, forward_pass) + end end function forward_pass( @@ -25,41 +45,20 @@ function forward_pass( options::Options, pass::AlternativeForwardPass{T}, ) where {T} - return forward_pass(pass.model, options, DefaultForwardPass()) + return forward_pass(pass.model, options, pass.forward_pass) end -struct AlternativeParallelScheme{T} <: AbstractParallelScheme +""" + AlternativePostIterationCallback(forward_model::PolicyGraph) + +A post-iteration callback that should be used whenever [`SDDP.AlternativeForwardPass`](@ref) +is used. +""" +struct AlternativePostIterationCallback{T} model::PolicyGraph{T} end -Base.show(io::IO, ::AlternativeParallelScheme) = print(io, "alternative") - -interrupt(::AlternativeParallelScheme) = nothing - -function master_loop( - scheme::AlternativeParallelScheme{T}, - model::PolicyGraph{T}, - options::Options, -) where {T} - _initialize_solver(model; throw_error = false) - while true - result = iteration(model, options) - slave_update(scheme.model, result) - log_iteration(options) - if result.has_converged - return result.status - end - end +function (callback::AlternativePostIterationCallback)(result::IterationResult) + slave_update(callback.model, result) return end - -function train_with_forward_model(nonconvex, convex; kwargs...) - @assert isempty(setdiff(keys(nonconvex.initial_root_state), keys(convex.initial_root_state))) - @assert isempty(setdiff(keys(convex.initial_root_state), keys(nonconvex.initial_root_state))) - return SDDP.train( - convex; - forward_pass = AlternativeForwardPass(nonconvex), - parallel_scheme = AlternativeParallelScheme(nonconvex), - kwargs..., - ) -end diff --git a/src/plugins/parallel_schemes.jl b/src/plugins/parallel_schemes.jl index 631bd3008..951e9298d 100644 --- a/src/plugins/parallel_schemes.jl +++ b/src/plugins/parallel_schemes.jl @@ -40,6 +40,7 @@ function master_loop( _initialize_solver(model; throw_error = false) while true result = iteration(model, options) + options.post_iteration_callback(result) log_iteration(options) if result.has_converged return result.status @@ -167,6 +168,7 @@ function slave_loop( results_to_add = IterationResult{T}[] while true result = iteration(model, options) + options.post_iteration_callback(result) # The next four lines are subject to a race condition: if the master closes # `results` _after_ the call to `isopen` and _before_` the call to `put!` has # executed, we get an `InvalidStateException`. This gets trapped in the outer @@ -243,6 +245,7 @@ function master_loop( # implementation anyway. while async.use_master && !isready(results) result = iteration(model, options) + options.post_iteration_callback(result) for (_, ch) in updates put!(ch, result) end