Various updates

odow · Sep 19, 2021 · 5c87dea · 5c87dea
1 parent 62b2b84
commit 5c87dea
Show file tree

Hide file tree

Showing 5 changed files with 57 additions and 80 deletions.
diff --git a/papers/biobjective_sddp/BiObjectiveSDDP.jl b/papers/biobjective_sddp/BiObjectiveSDDP.jl
@@ -6,10 +6,7 @@
 module BiObjectiveSDDP
 
 # Include the Gurobi-specific versions of get_BinvA and get_basis.
-function include_gurobi_specific_functions()
-    include(joinpath(@__DIR__, "gurobi.jl"))
-    return
-end
+include(joinpath(@__DIR__, "gurobi.jl"))
 
 import Printf
 import SDDP
@@ -395,6 +392,7 @@ function get_next_lambda(
     # Quickly optimize `dest` to obtain a basis. Note: for sanity sake, the
     # ObjectiveValue of `dest` after this should be identical to the objective
     # value of node.subproblem (modulo numerical tolerances).
+    MOI.set(dest, MOI.Silent(), true)
     MOI.optimize!(dest)
     # Get the next lambda using MOI calls defined above.
     return get_next_lambda(
@@ -467,7 +465,7 @@ function print_iteration_header(io)
     )
     println(
         io,
-        "          BI-OBJECTIVE SDDP.jl (c) Oscar Dowson, 2019-19          ",
+        "          BI-OBJECTIVE SDDP.jl (c) Oscar Dowson, 2019-21          ",
     )
     println(io)
     println(io, "      Iterations")
@@ -602,7 +600,12 @@ function bi_objective_sddp(
                 major_iterations += 1
             end
             SDDP.set_trade_off_weight(model, λ)
-            SDDP.train(model; run_numerical_stability_report = false, kwargs...)
+            SDDP.train(
+                model;
+                run_numerical_stability_report = false,
+                add_to_existing_cuts = true,
+                kwargs...,
+            )
             minor_iterations += 1
             sddp_iterations += length(model.most_recent_training_results.log)
             if bi_objective_post_train_callback !== nothing

diff --git a/papers/biobjective_sddp/Project.toml b/papers/biobjective_sddp/Project.toml
@@ -0,0 +1,4 @@
+[deps]
+Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SDDP = "f4570300-c277-11e8-125c-4912f86ce65d"
diff --git a/papers/biobjective_sddp/brazilian_example.jl b/papers/biobjective_sddp/brazilian_example.jl
@@ -3,12 +3,12 @@
 #  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/.
 
-using SDDP
-using Gurobi
-using Plots
-using BiObjectiveSDDP
+include(joinpath(@__DIR__, "BiObjectiveSDDP.jl"))
+using .BiObjectiveSDDP
 
-BiObjectiveSDDP.include_gurobi_specific_functions()
+using SDDP
+import Gurobi
+import Plots
 
 const OBJ_1_SCALING = 0.01
 const OBJ_2_SCALING = 0.1
@@ -26,10 +26,8 @@ function create_model()
         @variables(
             sp,
             begin
-                # State variables.
                 0 <= storedEnergy[i = 1:4] <= storedEnergy_ub[i],
                 (SDDP.State, initial_value = storedEnergy_initial[i])
-                # Control variables.
                 0 <= spillEnergy[i = 1:4]
                 0 <= hydroGeneration[i = 1:4] <= hydro_ub[i]
                 thermal_lb[i][j] <=
@@ -39,14 +37,12 @@ function create_model()
                 0 <=
                 deficit[i = 1:4, j = 1:4] <=
                 demand[month][i] * deficit_ub[j]
-                # Dummy variables for helpers.
                 inflow[i = 1:4]
             end
         )
         @constraints(
             sp,
             begin
-                # Model constraints.
                 [i = 1:4],
                 sum(deficit[i, :]) +
                 hydroGeneration[i] +
@@ -84,7 +80,8 @@ function create_model()
         SDDP.initialize_biobjective_subproblem(sp)
         SDDP.parameterize(sp, Ω) do ω
             JuMP.fix.(inflow, ω)
-            return SDDP.set_biobjective_functions(sp, objective_1, objective_2)
+            SDDP.set_biobjective_functions(sp, objective_1, objective_2)
+            return
         end
     end
     return model
@@ -110,60 +107,43 @@ function extract_objectives(simulation)
 end
 
 function plot_objective_space(simulations, simulation_weights)
-    plot(title = "Objective Space")
+    Plots.plot(title = "Objective Space")
     for λ in simulation_weights
-        scatter!(
+        Plots.scatter!(
             extract_objectives(simulations[λ])...,
             label = "\\lambda = $(λ)",
             alpha = 0.4,
         )
     end
-    p = plot!(xlabel = "Deficit cost", ylabel = "Thermal cost")
-    savefig("objective_space.pdf")
+    p = Plots.plot!(xlabel = "Deficit cost", ylabel = "Thermal cost")
+    Plots.savefig("objective_space.pdf")
     return p
 end
 
 function plot_weight_space(weights, bounds, simulations)
-    plot(weights, bounds, label = "")
+    Plots.plot(weights, bounds, label = "")
     for (λ, sim) in simulations
         obj_1, obj_2 = extract_objectives(simulations[λ])
         weighted_sum = λ .* obj_1 .+ (1 - λ) .* obj_2
-        scatter!(
+        Plots.scatter!(
             fill(λ, length(weighted_sum)),
             weighted_sum,
             alpha = 0.4,
             label = "",
         )
     end
-    p = plot!(xlabel = "Weight \\lambda", ylabel = "Weighted-sum")
-    savefig("weight_space.pdf")
+    p = Plots.plot!(xlabel = "Weight \\lambda", ylabel = "Weighted-sum")
+    Plots.savefig("weight_space.pdf")
     return p
 end
 
-# function plot_publication(simulations, simulation_weights)
-#     plts = [
-#         SDDP.publication_plot(simulations[λ]; title = "λ = $(λ)") do data
-#             return sum(data[:storedEnergy][i].out for i in 1:4)
-#         end
-#         for λ in simulation_weights
-#     ]
-#     plot(plts...,
-#         xlabel = "Stage",
-#         layout = (1, length(keys)),
-#         margin = 5Plots.mm,
-#         size = (1000, 300),
-#         ylim = (0, 2.5e5),
-#         ylabel = "Stored energy"
-#     )
-# end
-
-import Logging
-Logging.global_logger(Logging.ConsoleLogger(stdout, Logging.Debug))
+# import Logging
+# Logging.global_logger(Logging.ConsoleLogger(stdout, Logging.Debug))
 
 model = create_model()
 lower_bound, weights, bounds = BiObjectiveSDDP.bi_objective_sddp(
     model,
-    with_optimizer(Gurobi.Optimizer, GUROBI_ENV, OutputFlag = 0);
+    () -> Gurobi.Optimizer(GUROBI_ENV);
     # BiObjectiveSDDP kwargs ...
     bi_objective_minor_iteration_limit = 60,
     bi_objective_lambda_atol = 1e-3,

diff --git a/papers/biobjective_sddp/gurobi.jl b/papers/biobjective_sddp/gurobi.jl
@@ -8,36 +8,41 @@
 #
 # See `BiObjectiveSDDP.get_BinvA` and `BiObjectiveSDDP.get_basis`.
 
-import BiObjectiveSDDP
 import Gurobi
 import SparseArrays
 
 ####
 ####    New functions for the Gurobi C API.
 ####
 
-function BiObjectiveSDDP.get_basis(model::Gurobi.Optimizer)
+function get_basis(model::Gurobi.Optimizer)
     p = Ref{Cint}()
-    @assert GRBgetintattr(model, "NumConstrs", p) == 0
+    @assert Gurobi.GRBgetintattr(model, "NumConstrs", p) == 0
     bhead = zeros(Cint, p[])
     ret = Gurobi.GRBgetBasisHead(model, bhead)
     @assert ret == 0
     bhead .+= 1
     return bhead
 end
 
-function BiObjectiveSDDP.get_BinvA(model::Gurobi.Optimizer)
+mutable struct GRBsvec
+    len::Cint
+    ind::Ptr{Cint}
+    val::Ptr{Cdouble}
+end
+
+function get_BinvA(model::Gurobi.Optimizer)
     p = Ref{Cint}()
-    @assert GRBgetintattr(model, "NumConstrs", p) == 0
+    @assert Gurobi.GRBgetintattr(model, "NumConstrs", p) == 0
     ncon = p[]
-    @assert GRBgetintattr(model, "NumVars", p) == 0
+    @assert Gurobi.GRBgetintattr(model, "NumVars", p) == 0
     nvar = p[]
     function _GRBBinvRowi(model::Gurobi.Optimizer, i::Int)
         ind = zeros(Cint, ncon + nvar)
         val = zeros(Cdouble, ncon + nvar)
-        x = Gurobi.GRBsvec(0, pointer(ind), pointer(val))
-        GC.@preserve ind val begin
-            @assert Gurobi.GRBBinvRowi(model, i, x) == 0
+        x = GRBsvec(0, pointer(ind), pointer(val))
+        GC.@preserve ind val x begin
+            @assert Gurobi.GRBBinvRowi(model, i, pointer_from_objref(x)) == 0
         end
         return ind[1:x.len], val[1:x.len]
     end

diff --git a/papers/biobjective_sddp/simple_example.jl b/papers/biobjective_sddp/simple_example.jl
@@ -3,25 +3,22 @@
 #  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/.
 
-using SDDP
-using Gurobi
-using BiObjectiveSDDP
-using Random
+include(joinpath(@__DIR__, "BiObjectiveSDDP.jl"))
+using .BiObjectiveSDDP
 
-BiObjectiveSDDP.include_gurobi_specific_functions()
+using SDDP
+import Gurobi
+import Random
 
 const GUROBI_ENV = Gurobi.Env()
 
 function create_model()
     model = SDDP.LinearPolicyGraph(
         stages = 2,
         lower_bound = 0.0,
-        optimizer = with_optimizer(
-            Gurobi.Optimizer,
-            GUROBI_ENV,
-            OutputFlag = 0,
-        ),
+        optimizer = () -> Gurobi.Optimizer(GUROBI_ENV),
     ) do sp, t
+        set_silent(sp)
         @variable(sp, x >= 0, SDDP.State, initial_value = 0.0)
         if t == 1
             @expression(sp, objective_1, 2 * x.out)
@@ -36,34 +33,22 @@ function create_model()
             @expression(sp, objective_1, y)
             @expression(sp, objective_2, 3 * y)
         end
-        SDDP.add_objective_state(
-            sp,
-            initial_value = 1.0,
-            lower_bound = 0.0,
-            upper_bound = 1.0,
-            lipschitz = 10.0,
-        ) do y, ω
-            return y
-        end
+        SDDP.initialize_biobjective_subproblem(sp)
         SDDP.parameterize(sp, [nothing]) do ω
-            λ = SDDP.objective_state(sp)
-            @stageobjective(sp, λ * objective_1 + (1 - λ) * objective_2)
+            SDDP.set_biobjective_functions(sp, objective_1, objective_2)
         end
     end
     return model
 end
 
-import Logging
-Logging.global_logger(Logging.ConsoleLogger(stdout, Logging.Debug))
-
 Random.seed!(1)
 
 bounds_for_reporting = Tuple{Float64,Float64,Float64}[]
 
 model = create_model()
 lower_bound, weights, bounds = BiObjectiveSDDP.bi_objective_sddp(
     model,
-    with_optimizer(Gurobi.Optimizer, GUROBI_ENV, OutputFlag = 0);
+    () -> Gurobi.Optimizer(GUROBI_ENV);
     # BiObjectiveSDDP kwargs ...
     bi_objective_sddp_iteration_limit = 20,
     bi_objective_lower_bound = 0.0,
@@ -74,14 +59,14 @@ lower_bound, weights, bounds = BiObjectiveSDDP.bi_objective_sddp(
         begin
             upper_bound = BiObjectiveSDDP.surrogate_upper_bound(
                 model,
-                with_optimizer(Gurobi.Optimizer, GUROBI_ENV, OutputFlag = 0);
+                () -> Gurobi.Optimizer(GUROBI_ENV);
                 global_lower_bound = 0.0,
                 lambda_minimum_step = 1e-4,
                 lambda_atol = 1e-4,
             )
             lower_bound, _, _ = BiObjectiveSDDP.surrogate_lower_bound(
                 model,
-                with_optimizer(Gurobi.Optimizer, GUROBI_ENV, OutputFlag = 0);
+                () -> Gurobi.Optimizer(GUROBI_ENV);
                 global_lower_bound = 0.0,
                 lambda_minimum_step = 1e-4,
                 lambda_atol = 1e-4,