Add DiffOptt

Project.toml

@@ -5,13 +5,17 @@ version = "1.0.0-DEV"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+DiffOpt = "930fe3bc-9c6b-11ea-2d94-6184641e85e7"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 julia = "1"
@@ -21,6 +25,7 @@ CUDA_Runtime_jll = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
 
 [targets]
-test = ["Test", "HiGHS", "Zygote"]
+test = ["Test", "SCS", "Zygote"]

examples/HydroPowerModels/gen_inputs_l2O_hydropowermodels.jl

@@ -21,7 +21,7 @@ subproblems, state_params_in, state_params_out, uncertainty_samples, initial_sta
     joinpath(HydroPowerModels_dir, case_name), formulation_file; num_stages=num_stages
 )
 
-det_equivalent, uncertainty_samples = DecisionRules.deterministic_equivalent(subproblems, state_params_in, state_params_out, initial_state, uncertainty_samples)
+det_equivalent, uncertainty_samples = DecisionRules.deterministic_equivalent!(JuMP.Model(), subproblems, state_params_in, state_params_out, initial_state, uncertainty_samples)
 
 # Remove state imposing constraints
 

examples/HydroPowerModels/load_hydropowermodels.jl

@@ -57,7 +57,7 @@ function read_inflow(file::String, nHyd::Int; num_stages=nothing)
     return vector_inflows, nCen, num_stages
 end
 
-function build_hydropowermodels(case_folder::AbstractString, subproblem_file::AbstractString; num_stages=nothing, param_type=:Var, penalty=nothing) # :Param, :Cons, :Var
+function build_hydropowermodels(case_folder::AbstractString, subproblem_file::AbstractString; num_stages=nothing, penalty=nothing)
     hydro_file = JSON.parsefile(joinpath(case_folder, "hydro.json"))["Hydrogenerators"]
     nHyd = length(hydro_file)
     vector_inflows, nCen, num_stages = read_inflow(joinpath(case_folder, "inflows.csv"), nHyd; num_stages=num_stages)
@@ -77,9 +77,9 @@ function build_hydropowermodels(case_folder::AbstractString, subproblem_file::Ab
             delete(subproblems[t], con)
         end
         state_params_in[t], state_param_out, inflow = find_reservoirs_and_inflow(subproblems[t])
-        state_params_in[t] = variable_to_parameter.(subproblems[t], state_params_in[t], param_type=param_type)
-        state_params_out[t] = [variable_to_parameter(subproblems[t], state_param_out[i]; deficit=_deficit[i], param_type=param_type) for i in 1:nHyd]
-        inflow = [(variable_to_parameter(subproblems[t], inflow[i]; param_type=param_type),  vector_inflows[i][t, :] .+ 0.0) for i in 1:nHyd]
+        state_params_in[t] = variable_to_parameter.(subproblems[t], state_params_in[t])
+        state_params_out[t] = [variable_to_parameter(subproblems[t], state_param_out[i]; deficit=_deficit[i]) for i in 1:nHyd]
+        inflow = [(variable_to_parameter(subproblems[t], inflow[i]),  vector_inflows[i][t, :] .+ 0.0) for i in 1:nHyd]
         uncertainty_samples[t] = inflow
     end
 

examples/HydroPowerModels/test_dr_hydropowermodels.jl

@@ -6,10 +6,9 @@ using Gurobi
 using MosekTools
 using Ipopt, HSL_jll
 import MathOptSymbolicAD
-import ParametricOptInterface as POI
 using JLD2
 using HydroPowerModels
-
+using DiffOpt
 
 HydroPowerModels_dir = dirname(@__FILE__)
 include(joinpath(HydroPowerModels_dir, "load_hydropowermodels.jl"))
@@ -40,26 +39,16 @@ HydroPowerModels.gather_useful_info!(data)
 # Build MSP
 
 subproblems, state_params_in, state_params_out, uncertainty_samples, initial_state, max_volume = build_hydropowermodels(    
-    joinpath(HydroPowerModels_dir, case_name), formulation_file; num_stages=num_stages, param_type=:Var
+    joinpath(HydroPowerModels_dir, case_name), formulation_file; num_stages=num_stages
 )
 
-det_equivalent, uncertainty_samples = DecisionRules.deterministic_equivalent(subproblems, state_params_in, state_params_out, initial_state, uncertainty_samples)
-
-set_optimizer(det_equivalent, Gurobi.Optimizer)
-
-# set_optimizer(det_equivalent, Mosek.Optimizer)
-
-# set_optimizer(det_equivalent, optimizer_with_attributes(Ipopt.Optimizer, 
-#     "print_level" => 0,
-#     "hsllib" => HSL_jll.libhsl_path,
-#     "linear_solver" => "ma27"
-# ))
+det_equivalent = DiffOpt.diff_model(optimizer_with_attributes(Ipopt.Optimizer, 
+    "print_level" => 0,
+    "hsllib" => HSL_jll.libhsl_path,
+    "linear_solver" => "ma27"
+))
 
-# set_attribute(
-#     det_equivalent,
-#     MOI.AutomaticDifferentiationBackend(),
-#     MathOptSymbolicAD.DefaultBackend(),
-# )
+det_equivalent, uncertainty_samples = DecisionRules.deterministic_equivalent!(det_equivalent, subproblems, state_params_in, state_params_out, initial_state, uncertainty_samples)
 
 num_hydro = length(initial_state)
 
@@ -103,13 +92,13 @@ for i in 1:num_samples
         inflow_var = inflow_var[findfirst(x -> occursin("_inflow[$j]", JuMP.name(x)), inflow_var)]
         inflows[i, j, t] = uncertainty_s[t][inflow_var]
     end
-    objective_values[i] = simulate_multistage(
+    simulate_multistage(
         det_equivalent, state_params_in, state_params_out, 
         initial_state, uncertainty_s, 
         models;
-        ensure_feasibility=(x_out, x_in, _sa) -> ensure_feasibility(x_out, x_in, _sa, max_volume),
-        _objective_value = DecisionRules.get_objective_no_target_deficit,
+        ensure_feasibility=(x_out, x_in, _sa) -> ensure_feasibility(x_out, x_in, _sa, max_volume)
     )
+    objective_values[i] = DecisionRules.get_objective_no_target_deficit(det_equivalent)
     for _var in record_variables_names
         num_vars = length(find_variables(det_equivalent, [_var; r"#1$"]))
         for j in 1:num_vars, t in 1:num_stages

examples/HydroPowerModels/train_dr_hydropowermodels.jl

@@ -20,6 +20,7 @@ using Ipopt, HSL_jll # Gurobi, MosekTools, Ipopt, MadNLP
 # import ParametricOptInterface as POI
 using Wandb, Dates, Logging
 using JLD2
+using DiffOpt
 
 HydroPowerModels_dir = dirname(@__FILE__)
 include(joinpath(HydroPowerModels_dir, "load_hydropowermodels.jl"))
@@ -52,41 +53,16 @@ pre_trained_model = nothing #joinpath(HydroPowerModels_dir, case_name, formulati
 # Build MSP
 
 subproblems, state_params_in, state_params_out, uncertainty_samples, initial_state, max_volume = build_hydropowermodels(    
-    joinpath(HydroPowerModels_dir, case_name), formulation_file; num_stages=num_stages, param_type=:Var
+    joinpath(HydroPowerModels_dir, case_name), formulation_file; num_stages=num_stages
 )
 
-det_equivalent, uncertainty_samples = DecisionRules.deterministic_equivalent(subproblems, state_params_in, state_params_out, initial_state, uncertainty_samples) #; model = JuMP.Model(() -> POI_cached_optimizer()))
-
-set_optimizer(det_equivalent, optimizer_with_attributes(Ipopt.Optimizer, 
+det_equivalent = DiffOpt.diff_model(optimizer_with_attributes(Ipopt.Optimizer, 
     "print_level" => 0,
     "hsllib" => HSL_jll.libhsl_path,
     "linear_solver" => "ma27"
 ))
 
-# set_optimizer(det_equivalent, Gurobi.Optimizer)
-
-# set_optimizer(det_equivalent, Mosek.Optimizer)
-
-# ipopt = Ipopt.Optimizer()
-# MOI.set(ipopt, MOI.RawOptimizerAttribute("print_level"), 0)
-# MOI.set(ipopt, MOI.RawOptimizerAttribute("hsllib"), HSL_jll.libhsl_path)
-# MOI.set(ipopt, MOI.RawOptimizerAttribute("linear_solver"), "ma97")
-# cached =
-#     () -> MOI.Bridges.full_bridge_optimizer(
-#         MOI.Utilities.CachingOptimizer(
-#             MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
-#             ipopt,
-#         ),
-#         Float64,
-# )
-# POI_cached_optimizer() = POI.Optimizer(cached())
-
-# set_optimizer(det_equivalent, () -> POI.Optimizer(Ipopt.Optimizer()))
-
-# set_optimizer(det_equivalent, () -> POI_cached_optimizer())
-# set_optimizer(det_equivalent, () -> Mosek.Optimizer())
-# set_attribute(det_equivalent, "QUIET", true)
-# set_attributes(det_equivalent, "OutputFlag" => 0)
+det_equivalent, uncertainty_samples = DecisionRules.deterministic_equivalent!(det_equivalent, subproblems, state_params_in, state_params_out, initial_state, uncertainty_samples)
 
 num_hydro = length(initial_state)
 # for subproblem in subproblems

examples/RL/hydro_pre_trained.jl

@@ -56,7 +56,7 @@ end
     formulation_file = formulation * ".mof.json"
     Random.seed!(1234)
     subproblems, state_params_in, state_params_out, uncertainty_samples, initial_state, max_volume = build_hydropowermodels(    
-        joinpath(HydroPowerModels_dir, case_name), formulation_file; num_stages=num_stages, param_type=:Var, penalty=penalty
+        joinpath(HydroPowerModels_dir, case_name), formulation_file; num_stages=num_stages, penalty=penalty
     )
 
     for subproblem in subproblems
@@ -88,7 +88,7 @@ end
             uncertainty_sample = [(uncertainty_samples[j][i][1], rain_state[i]) for i in 1:length(rain)]
             simulate_stage(subproblems[j], state_params_in[j], state_params_out[j], uncertainty_sample, state_in, state_out)
             r = _objective_value(subproblems[j])
-            next_volume = DecisionRules.get_next_state(subproblems[j], state_params_out[j], state_in, state_out)
+            next_volume = DecisionRules.get_next_state(subproblems[j], state_params_in[j], state_params_out[j], state_in, state_out)
             @info "Stage t=$j" sum(state_in) sum(rain_state) sum(state_out) sum(next_volume) r
             sp = [next_volume; rain; j+1]
             o = rain #+ next_volume - state_out

examples/RL/hydropowermodels_rl.jl

@@ -58,7 +58,7 @@ end
     formulation_file = formulation * ".mof.json"
     Random.seed!(1234)
     subproblems, state_params_in, state_params_out, uncertainty_samples, initial_state, max_volume = build_hydropowermodels(    
-        joinpath(HydroPowerModels_dir, case_name), formulation_file; num_stages=num_stages, param_type=:Var, penalty=penalty
+        joinpath(HydroPowerModels_dir, case_name), formulation_file; num_stages=num_stages, penalty=penalty
     )
 
     for subproblem in subproblems
@@ -95,7 +95,7 @@ end
             uncertainty_sample = [(uncertainty_samples[j][i][1], rain_state[i]) for i in 1:length(rain)]
             simulate_stage(subproblems[j], state_params_in[j], state_params_out[j], uncertainty_sample, state_in, state_out)
             r = _objective_value(subproblems[j])
-            sp = DecisionRules.get_next_state(subproblems[j], state_params_out[j], state_in, state_out)
+            sp = DecisionRules.get_next_state(subproblems[j], state_params_in[j], state_params_out[j], state_in, state_out)
             @info "Stage t=$j" sum(state_in) sum(rain_state) sum(state_out) sum(sp) r
             sp = [sp; rain; j+1]
             o = sp # no hidden state

examples/rocket_control/build_rocket_problem.jl

@@ -5,6 +5,7 @@ using Statistics
 
 using JuMP
 import Ipopt, HSL_jll
+using DiffOpt
 # import Plots
 
 
@@ -32,8 +33,7 @@ function build_rocket_problem(;
     # program, we need to use a nonlinear solver like Ipopt. We cannot use a linear
     # solver like HiGHS.
 
-    det_equivalent = Model()
-    set_optimizer(det_equivalent, optimizer_with_attributes(Ipopt.Optimizer, 
+    det_equivalent = DiffOpt.diff_model(optimizer_with_attributes(Ipopt.Optimizer, 
         "print_level" => 0,
         "hsllib" => HSL_jll.libhsl_path,
         "linear_solver" => "ma27"

src/DecisionRules.jl

@@ -4,8 +4,10 @@ using JuMP
 import ParametricOptInterface as POI
 using Flux
 using JLD2
+using ChainRules: @ignore_derivatives
 using ChainRulesCore
 import ChainRulesCore.rrule
+using DiffOpt
 
 export simulate_multistage, sample, train_multistage, simulate_states, simulate_stage, dense_multilayer_nn, variable_to_parameter, create_deficit!, 
     SaveBest, find_variables, identity