From 78aa8719e2910022420b209f9799c16b9004b035 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Sun, 17 Sep 2023 22:33:54 -0400 Subject: [PATCH 01/36] working model for feasible limits of load --- Project.toml | 8 ++- examples/powermodels/pglib_datagen.jl | 13 +--- src/worst_case.jl | 2 + test/worst_case.jl | 93 +++++++++++++++++++++++++++ 4 files changed, 102 insertions(+), 14 deletions(-) create mode 100644 src/worst_case.jl create mode 100644 test/worst_case.jl diff --git a/Project.toml b/Project.toml index 613929b..2910533 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "1.0.0-DEV" [deps] Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +Dualization = "191a621a-6537-11e9-281d-650236a99e60" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" @@ -18,14 +19,15 @@ ParametricOptInterface = "0.5" julia = "1.6" [extras] +Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" -Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +PGLib = "07a8691f-3d11-4330-951b-3c50f98338be" PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" [targets] -test = ["Test", "DelimitedFiles", "Downloads", "HiGHS", "PowerModels", "Flux", "DataFrames", "Clarabel"] +test = ["Test", "DelimitedFiles", "PGLib", "HiGHS", "PowerModels", "Flux", "DataFrames", "Clarabel", "Ipopt"] diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 019899e..8fffef3 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -1,4 +1,4 @@ -using Downloads +using PGLib using PowerModels using JuMP, HiGHS import ParametricOptInterface as POI @@ -56,15 +56,6 @@ function generate_dataset_pglib( ) # Download file matpower_case_name = case_name * ".m" - case_file_path = joinpath(data_dir, matpower_case_name) - mkpath(data_dir) - if download_files && !isfile(case_file_path) - Downloads.download( - "https://raw.githubusercontent.com/power-grid-lib/pglib-opf/01681386d084d8bd03b429abcd1ee6966f68b9a3/" * - matpower_case_name, - case_file_path, - ) - end # save folder data_sim_dir = joinpath(data_dir, string(network_formulation)) @@ -73,7 +64,7 @@ function generate_dataset_pglib( end # Read data - network_data = PowerModels.parse_file(case_file_path) + network_data = make_basic_network(pglib(matpower_case_name)) # The problem to iterate over model = Model(solver) diff --git a/src/worst_case.jl b/src/worst_case.jl new file mode 100644 index 0000000..2b7a243 --- /dev/null +++ b/src/worst_case.jl @@ -0,0 +1,2 @@ +using Dualization + diff --git a/test/worst_case.jl b/test/worst_case.jl new file mode 100644 index 0000000..68323d5 --- /dev/null +++ b/test/worst_case.jl @@ -0,0 +1,93 @@ +using PowerModels +using Ipopt +using PGLib +using JuMP +using Dualization +import JuMP.MOI as MOI + +function _moi_set_upper_bound( + moi_backend, + idx::MOI.VariableIndex, + upper::Number, +) + new_set = MOI.LessThan(upper) + + JuMP._moi_add_constraint(moi_backend, idx, new_set) + + return +end + +function _moi_set_lower_bound( + moi_backend, + idx::MOI.VariableIndex, + lower::Number, +) + new_set = MOI.GreaterThan(lower) + + JuMP._moi_add_constraint(moi_backend, idx, new_set) + + return +end + +# Load data +case_name = "pglib_opf_case5_pjm" +matpower_case_name = case_name * ".m" +network_data = pglib(matpower_case_name) + +# Case Parameters +network_formulation=DCPPowerModel +original_load = [l["pd"] for l in values(network_data["load"])] +max_total_load = sum(original_load) * 1.1 + +# Create model +solver = Ipopt.Optimizer +model = Model(solver) +load_var = @variable(model, load_var[i=1:length(original_load)]) +for (i, l) in enumerate(values(network_data["load"])) + l["pd"] = load_var[i] +end + +# Instantiate the model +pm = instantiate_model( + network_data, + network_formulation, + PowerModels.build_opf; + setting=Dict("output" => Dict("duals" => true)), + jump_model=model, +) + +# Solve the model +JuMP.optimize!(model) + +optimal_cost = JuMP.objective_value(model) + +# Dualize the model +dual_st = Dualization.dualize(JuMP.backend(model), + variable_parameters = MOI.VariableIndex[i.index for i in load_var] +) + +dual_model = dual_st.dual_model +primal_dual_map = dual_st.primal_dual_map + +for (i,l) in enumerate(load_var) + _moi_set_upper_bound(dual_model, primal_dual_map.primal_parameter[l.index], max_total_load) + _moi_set_lower_bound(dual_model, primal_dual_map.primal_parameter[l.index], 0.0) +end + +jump_dual_model = JuMP.Model() +map_moi_to_jump = MOI.copy_to(JuMP.backend(jump_dual_model), dual_model) +set_optimizer(jump_dual_model, solver) + +load_dual_indexes = [map_moi_to_jump[primal_dual_map.primal_parameter[l.index]] for l in load_var] + +JuMP.optimize!(jump_dual_model) + +optimal_dual_cost = JuMP.objective_value(jump_dual_model) + +optimal_loads = [MOI.get(JuMP.backend(jump_dual_model), MOI.VariablePrimal(), l) for l in load_dual_indexes] + +cost_diff = (optimal_dual_cost - optimal_cost) / optimal_cost + +load_diff = sum(abs.(optimal_loads - original_load)) ./ sum(original_load) + +aux_load = optimal_loads \ No newline at end of file From c30cb31cb52e8e02517240b580c31b54c13466d0 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Sun, 17 Sep 2023 23:10:49 -0400 Subject: [PATCH 02/36] update test --- test/worst_case.jl | 49 +++++++++++++++++++++++++++++++++------------- 1 file changed, 35 insertions(+), 14 deletions(-) diff --git a/test/worst_case.jl b/test/worst_case.jl index 68323d5..fc4eb9b 100644 --- a/test/worst_case.jl +++ b/test/worst_case.jl @@ -8,7 +8,7 @@ import JuMP.MOI as MOI function _moi_set_upper_bound( moi_backend, idx::MOI.VariableIndex, - upper::Number, + upper, ) new_set = MOI.LessThan(upper) @@ -20,7 +20,7 @@ end function _moi_set_lower_bound( moi_backend, idx::MOI.VariableIndex, - lower::Number, + lower, ) new_set = MOI.GreaterThan(lower) @@ -46,6 +46,7 @@ load_var = @variable(model, load_var[i=1:length(original_load)]) for (i, l) in enumerate(values(network_data["load"])) l["pd"] = load_var[i] end +load_moi_idx = MOI.VariableIndex[i.index for i in load_var] # Instantiate the model pm = instantiate_model( @@ -57,37 +58,57 @@ pm = instantiate_model( ) # Solve the model -JuMP.optimize!(model) +# JuMP.optimize!(model) -optimal_cost = JuMP.objective_value(model) +# optimal_cost = JuMP.objective_value(model) # Dualize the model dual_st = Dualization.dualize(JuMP.backend(model), - variable_parameters = MOI.VariableIndex[i.index for i in load_var] + variable_parameters = load_moi_idx ) dual_model = dual_st.dual_model primal_dual_map = dual_st.primal_dual_map -for (i,l) in enumerate(load_var) - _moi_set_upper_bound(dual_model, primal_dual_map.primal_parameter[l.index], max_total_load) - _moi_set_lower_bound(dual_model, primal_dual_map.primal_parameter[l.index], 0.0) -end +# for (i,l) in enumerate(load_var) +# _moi_set_upper_bound(dual_model, primal_dual_map.primal_parameter[l.index], max_total_load) +# _moi_set_lower_bound(dual_model, primal_dual_map.primal_parameter[l.index], 0.0) +# end jump_dual_model = JuMP.Model() map_moi_to_jump = MOI.copy_to(JuMP.backend(jump_dual_model), dual_model) set_optimizer(jump_dual_model, solver) -load_dual_indexes = [map_moi_to_jump[primal_dual_map.primal_parameter[l.index]] for l in load_var] +load_dual_idxs = [map_moi_to_jump[primal_dual_map.primal_parameter[l]].value for l in load_moi_idx] + +load_var_dual = JuMP.all_variables(jump_dual_model)[load_dual_idxs] + +for (i,l) in enumerate(load_var_dual) + @constraint(jump_dual_model, l <= original_load[i]*1.5) + @constraint(jump_dual_model, l >= 0.0) +end +for (i, l) in enumerate(values(network_data["load"])) + l["pd"] = load_var_dual[i] +end + +obj = objective_function(jump_dual_model) + +pm = instantiate_model( + network_data, + network_formulation, + PowerModels.build_opf; + setting=Dict("output" => Dict("duals" => true)), + jump_model=jump_dual_model, +) + +@objective(jump_dual_model, Max, obj) JuMP.optimize!(jump_dual_model) optimal_dual_cost = JuMP.objective_value(jump_dual_model) -optimal_loads = [MOI.get(JuMP.backend(jump_dual_model), MOI.VariablePrimal(), l) for l in load_dual_indexes] +optimal_loads = value.(load_var_dual) cost_diff = (optimal_dual_cost - optimal_cost) / optimal_cost -load_diff = sum(abs.(optimal_loads - original_load)) ./ sum(original_load) - -aux_load = optimal_loads \ No newline at end of file +load_diff = sum(abs.(optimal_loads - original_load)) ./ sum(original_load) \ No newline at end of file From 33059d6827056d1487112e7223dae2afe56cbd71 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Mon, 18 Sep 2023 09:27:58 -0400 Subject: [PATCH 03/36] update code --- test/worst_case.jl | 70 +++++++++++++++++++++------------------------- 1 file changed, 32 insertions(+), 38 deletions(-) diff --git a/test/worst_case.jl b/test/worst_case.jl index fc4eb9b..f82c48c 100644 --- a/test/worst_case.jl +++ b/test/worst_case.jl @@ -5,41 +5,18 @@ using JuMP using Dualization import JuMP.MOI as MOI -function _moi_set_upper_bound( - moi_backend, - idx::MOI.VariableIndex, - upper, -) - new_set = MOI.LessThan(upper) - - JuMP._moi_add_constraint(moi_backend, idx, new_set) - - return -end - -function _moi_set_lower_bound( - moi_backend, - idx::MOI.VariableIndex, - lower, -) - new_set = MOI.GreaterThan(lower) - - JuMP._moi_add_constraint(moi_backend, idx, new_set) - - return -end - # Load data case_name = "pglib_opf_case5_pjm" matpower_case_name = case_name * ".m" -network_data = pglib(matpower_case_name) +network_data_original = pglib(matpower_case_name) +network_data = deepcopy(network_data_original) # Case Parameters network_formulation=DCPPowerModel original_load = [l["pd"] for l in values(network_data["load"])] max_total_load = sum(original_load) * 1.1 -# Create model +# Create primal model solver = Ipopt.Optimizer model = Model(solver) load_var = @variable(model, load_var[i=1:length(original_load)]) @@ -57,11 +34,6 @@ pm = instantiate_model( jump_model=model, ) -# Solve the model -# JuMP.optimize!(model) - -# optimal_cost = JuMP.objective_value(model) - # Dualize the model dual_st = Dualization.dualize(JuMP.backend(model), variable_parameters = load_moi_idx @@ -70,11 +42,6 @@ dual_st = Dualization.dualize(JuMP.backend(model), dual_model = dual_st.dual_model primal_dual_map = dual_st.primal_dual_map -# for (i,l) in enumerate(load_var) -# _moi_set_upper_bound(dual_model, primal_dual_map.primal_parameter[l.index], max_total_load) -# _moi_set_lower_bound(dual_model, primal_dual_map.primal_parameter[l.index], 0.0) -# end - jump_dual_model = JuMP.Model() map_moi_to_jump = MOI.copy_to(JuMP.backend(jump_dual_model), dual_model) set_optimizer(jump_dual_model, solver) @@ -109,6 +76,33 @@ optimal_dual_cost = JuMP.objective_value(jump_dual_model) optimal_loads = value.(load_var_dual) -cost_diff = (optimal_dual_cost - optimal_cost) / optimal_cost +load_diff = sum(abs.(optimal_loads - original_load)) ./ sum(original_load) -load_diff = sum(abs.(optimal_loads - original_load)) ./ sum(original_load) \ No newline at end of file +# Create final primal model +network_data = deepcopy(network_data_original) +model = Model(solver) +for (i, l) in enumerate(values(network_data["load"])) + l["pd"] = optimal_loads[i] +end + +# Instantiate the model +pm = instantiate_model( + network_data, + network_formulation, + PowerModels.build_opf; + setting=Dict("output" => Dict("duals" => true)), + jump_model=model, +) + +# Solve the model +JuMP.optimize!(model) + +# Check the solution +optimal_final_cost = JuMP.objective_value(model) +termination_status = JuMP.termination_status(model) +solution_primal_status = JuMP.primal_status(model) +solution_dual_status = JuMP.dual_status(model) +termination_status == MOI.INFEASIBLE && @error("Optimal solution not found") +solution_primal_status != MOI.FEASIBLE_POINT && @error("Primal solution not found") +solution_dual_status != MOI.FEASIBLE_POINT && @error("Dual solution not found") +isapprox(optimal_final_cost, optimal_dual_cost; rtol=1e-4) || @warn("Final cost is not equal to dual cost") From 85bf5ee8505447352cbbdd7d07daea39cc673e1e Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Mon, 18 Sep 2023 12:46:16 -0400 Subject: [PATCH 04/36] add worstcaseiterator --- Project.toml | 3 +- examples/powermodels/pglib_datagen.jl | 4 +- src/L2O.jl | 1 + src/arrowrecorder.jl | 8 +++ src/csvrecorder.jl | 12 ++++ src/datasetgen.jl | 19 +++--- src/worst_case.jl | 88 ++++++++++++++++++++++++++- test/datasetgen.jl | 6 +- 8 files changed, 127 insertions(+), 14 deletions(-) diff --git a/Project.toml b/Project.toml index 2910533..e996c08 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "L2O" uuid = "e1d8bfa7-c465-446a-84b9-451470f6e76c" authors = ["andrewrosemberg and contributors"] -version = "1.0.0-DEV" +version = "1.2.0-DEV" [deps] Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45" @@ -14,6 +14,7 @@ UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [compat] Arrow = "2" CSV = "0.10" +Dualization = "0.5" JuMP = "1" ParametricOptInterface = "0.5" julia = "1.6" diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 8fffef3..2b25a59 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -113,7 +113,7 @@ function generate_dataset_pglib( end number_vars = length(variable_refs) recorder = Recorder{filetype}(file; primal_variables=variable_refs) - return solve_batch(model, problem_iterator, recorder), + return solve_batch(problem_iterator, recorder), number_vars, length(original_load), batch_id @@ -135,7 +135,7 @@ function test_pglib_datasetgen(path::AbstractString, case_name::AbstractString, # Check if the number of successfull solves is equal to the number of problems saved @test isfile(file_out) @test length(readdlm(file_out, ',')[:, 1]) == num_p * success_solves + 1 - @test length(readdlm(file_out, ',')[1, :]) == number_variables + 1 + @test length(readdlm(file_out, ',')[1, :]) == number_variables + 2 return file_in, file_out end diff --git a/src/L2O.jl b/src/L2O.jl index 65f46ec..8b7a669 100644 --- a/src/L2O.jl +++ b/src/L2O.jl @@ -2,6 +2,7 @@ module L2O using Arrow using CSV +using Dualization using JuMP using UUIDs import ParametricOptInterface as POI diff --git a/src/arrowrecorder.jl b/src/arrowrecorder.jl index 816d1b9..562314b 100644 --- a/src/arrowrecorder.jl +++ b/src/arrowrecorder.jl @@ -8,6 +8,13 @@ Base.string(::Type{ArrowFile}) = "arrow" Record optimization problem solution to an Arrow file. """ function record(recorder::Recorder{ArrowFile}, id::UUID) + model = if length(recorder.primal_variables) > 0 + owner_model(recorder.primal_variables[1]) + elseif length(recorder.dual_variables) > 0 + owner_model(recorder.dual_variables[1]) + else + @error("Recorder has no variables") + end return Arrow.append( recorder.filename, (; @@ -20,6 +27,7 @@ function record(recorder::Recorder{ArrowFile}, id::UUID) Symbol.("dual_" .* name.(recorder.dual_variables)), [[dual.(p)] for p in recorder.dual_variables], )..., + objective=[JuMP.objective_value(model)], ), ) end diff --git a/src/csvrecorder.jl b/src/csvrecorder.jl index ed98013..a526755 100644 --- a/src/csvrecorder.jl +++ b/src/csvrecorder.jl @@ -17,6 +17,7 @@ function record(recorder::Recorder{CSVFile}, id::UUID) for p in recorder.dual_variables write(f, ",dual_$(name(p))") end + write(f, ",objective") write(f, "\n") end end @@ -30,6 +31,17 @@ function record(recorder::Recorder{CSVFile}, id::UUID) val = dual.(p) write(f, ",$val") end + # save objective value + model = if length(recorder.primal_variables) > 0 + owner_model(recorder.primal_variables[1]) + elseif length(recorder.dual_variables) > 0 + owner_model(recorder.dual_variables[1]) + else + @error("Recorder has no variables") + end + obj = JuMP.objective_value(model) + write(f, ",$obj") + # end line write(f, "\n") end end diff --git a/src/datasetgen.jl b/src/datasetgen.jl index f1db811..8790182 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -21,21 +21,25 @@ mutable struct Recorder{T<:RecorderFile} end end +abstract type AbstractProblemIterator end + """ ProblemIterator(ids::Vector{UUID}, pairs::Dict{VariableRef, Vector{Real}}) Iterator for optimization problem instances. """ -struct ProblemIterator{T<:Real} +struct ProblemIterator{T<:Real} <: AbstractProblemIterator + model::JuMP.Model ids::Vector{UUID} pairs::Dict{VariableRef,Vector{T}} function ProblemIterator( ids::Vector{UUID}, pairs::Dict{VariableRef,Vector{T}} ) where {T<:Real} + model = JuMP.owner_model(first(keys(pairs))) for (p, val) in pairs @assert length(ids) == length(val) end - return new{T}(ids, pairs) + return new{T}(model, ids, pairs) end end @@ -85,13 +89,14 @@ function update_model!(model::JuMP.Model, pairs::Dict, idx::Integer) end """ - solve_and_record(model::JuMP.Model, problem_iterator::ProblemIterator, recorder::Recorder, idx::Integer) + solve_and_record(problem_iterator::ProblemIterator, recorder::Recorder, idx::Integer) Solve an optimization problem and record the solution. """ function solve_and_record( - model::JuMP.Model, problem_iterator::ProblemIterator, recorder::Recorder, idx::Integer + problem_iterator::ProblemIterator, recorder::Recorder, idx::Integer ) + model = problem_iterator.model update_model!(model, problem_iterator.pairs, idx) optimize!(model) if recorder.filterfn(model) @@ -102,16 +107,16 @@ function solve_and_record( end """ - solve_batch(model::JuMP.Model, problem_iterator::ProblemIterator, recorder::Recorder) + solve_batch(problem_iterator::ProblemIterator, recorder::Recorder) Solve a batch of optimization problems and record the solutions. """ function solve_batch( - model::JuMP.Model, problem_iterator::ProblemIterator, recorder::Recorder + problem_iterator::ProblemIterator, recorder::Recorder ) successfull_solves = sum( - solve_and_record(model, problem_iterator, recorder, idx) for + solve_and_record(problem_iterator, recorder, idx) for idx in 1:length(problem_iterator.ids) ) / length(problem_iterator.ids) diff --git a/src/worst_case.jl b/src/worst_case.jl index 2b7a243..1fe2f46 100644 --- a/src/worst_case.jl +++ b/src/worst_case.jl @@ -1,2 +1,88 @@ -using Dualization +struct WorstCaseProblemIterator <: AbstractProblemIterator + ids::Vector{UUID} + parameters::Function + primal_builder!::Function + set_iterator!::Function + optimizer + function WorstCaseProblemIterator( + ids::Vector{UUID}, + parameters::Function, + primal_builder!::Function, + set_iterator!::Function, + optimizer, + ) + return new(ids, parameters, primal_builder!, set_iterator!, optimizer) + end +end +function solve_and_record( + problem_iterator::WorstCaseProblemIterator, recorder::Recorder, idx::Integer +) + # Build Primal + model = JuMP.Model() + parameters = problem_iterator.parameters(model) + problem_iterator.primal_builder!(model, parameters) + + # Parameter indices + load_moi_idx = MOI.VariableIndex[i.index for i in parameters] + + # Dualize the model + dual_st = Dualization.dualize(JuMP.backend(model), + variable_parameters = load_moi_idx + ) + + dual_model = dual_st.dual_model + primal_dual_map = dual_st.primal_dual_map + + # Build Dual in JuMP + jump_dual_model = JuMP.Model() + map_moi_to_jump = MOI.copy_to(JuMP.backend(jump_dual_model), dual_model) + set_optimizer(jump_dual_model, problem_iterator.optimizer) + + # Get dual variables for the parameters + load_dual_idxs = [map_moi_to_jump[primal_dual_map.primal_parameter[l]].value for l in load_moi_idx] + load_var_dual = JuMP.all_variables(jump_dual_model)[load_dual_idxs] + + # Add constraints to the dual associated with the parameters + problem_iterator.set_iterator!(jump_dual_model, parameters, idx) + + # Get the objective function + obj = objective_function(jump_dual_model) + dual_sense = JuMP.objective_sense(jump_dual_model) + + # Inforce primal constraints + problem_iterator.primal_builder!(jump_dual_model, load_var_dual) + + # Re-set objective function in case primal_builder! overwrote it + @objective(jump_dual_model, dual_sense, obj) + + # Solve the dual + JuMP.optimize!(jump_dual_model) + optimal_loads = value.(load_var_dual) + optimal_dual_cost = JuMP.objective_value(jump_dual_model) + + # Create final primal model and solve + model = JuMP.Model(problem_iterator.optimizer) + problem_iterator.primal_builder!(jump_dual_model, optimal_loads) + JuMP.optimize!(model) + + # Check if method was effective + optimal_final_cost = JuMP.objective_value(model) + termination_status = JuMP.termination_status(model) + solution_primal_status = JuMP.primal_status(model) + solution_dual_status = JuMP.dual_status(model) + termination_status == MOI.INFEASIBLE && @error("Optimal solution not found") + solution_primal_status != MOI.FEASIBLE_POINT && @error("Primal solution not found") + solution_dual_status != MOI.FEASIBLE_POINT && @error("Dual solution not found") + + if !isapprox(optimal_final_cost, optimal_dual_cost; rtol=1e-4) + rtol = abs(optimal_final_cost - optimal_dual_cost) / optimal_final_cost * 100 + @warn "Final cost is not equal to dual cost by $(rtol) %" optimal_final_cost optimal_dual_cost + end + + if recorder.filterfn(model) + record(recorder, problem_iterator.ids[idx]) + return 1 + end + return 0 +end \ No newline at end of file diff --git a/test/datasetgen.jl b/test/datasetgen.jl index 27f764c..af56a3f 100644 --- a/test/datasetgen.jl +++ b/test/datasetgen.jl @@ -35,7 +35,7 @@ function testdataset_gen(path::AbstractString) ) # Solve all problems and record solutions - solve_batch(model, problem_iterator, recorder) + solve_batch(problem_iterator, recorder) # Check if file exists and has the correct number of rows and columns @test isfile(file_output) @@ -46,7 +46,7 @@ function testdataset_gen(path::AbstractString) rm(file_input) # test output file @test length(readdlm(file_output, ',')[:, 1]) == num_p + 1 - @test length(readdlm(file_output, ',')[1, :]) == 3 + @test length(readdlm(file_output, ',')[1, :]) == 4 rm(file_output) else # test input file @@ -56,7 +56,7 @@ function testdataset_gen(path::AbstractString) rm(file_input) # test output file df = Arrow.Table(file_output) - @test length(df) == 3 + @test length(df) == 4 @test length(df[1]) == num_p rm(file_output) end From 4ea4a9b85d1af943051e65b289c7f6d23f406c9e Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Mon, 18 Sep 2023 14:48:34 -0400 Subject: [PATCH 05/36] update --- examples/powermodels/pglib_datagen.jl | 2 +- src/L2O.jl | 3 +- src/arrowrecorder.jl | 8 +- src/csvrecorder.jl | 11 +- src/datasetgen.jl | 40 ++++-- src/worst_case.jl | 40 +++++- test/datasetgen.jl | 4 +- test/runtests.jl | 6 +- test/worst_case.jl | 177 +++++++++++--------------- 9 files changed, 161 insertions(+), 130 deletions(-) diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 2b25a59..01f92ab 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -38,7 +38,7 @@ function load_sampler( end """ - generate_dataset_pglib(data_dir::AbstractString, case_name::AbstractString; download_files::Bool=true, filetype::Type{RecorderFile}, + generate_dataset_pglib(data_dir::AbstractString, case_name::AbstractString; download_files::Bool=true, filetype::Type{FileType}, num_p::Int=10 ) diff --git a/src/L2O.jl b/src/L2O.jl index 8b7a669..94c3122 100644 --- a/src/L2O.jl +++ b/src/L2O.jl @@ -8,10 +8,11 @@ using UUIDs import ParametricOptInterface as POI import Base: string -export ArrowFile, CSVFile, ProblemIterator, Recorder, save, solve_batch +export ArrowFile, CSVFile, ProblemIterator, Recorder, save, solve_batch, WorstCaseProblemIterator include("datasetgen.jl") include("csvrecorder.jl") include("arrowrecorder.jl") +include("worst_case.jl") end diff --git a/src/arrowrecorder.jl b/src/arrowrecorder.jl index 562314b..0ed3420 100644 --- a/src/arrowrecorder.jl +++ b/src/arrowrecorder.jl @@ -1,4 +1,4 @@ -abstract type ArrowFile <: RecorderFile end +abstract type ArrowFile <: FileType end Base.string(::Type{ArrowFile}) = "arrow" @@ -7,7 +7,9 @@ Base.string(::Type{ArrowFile}) = "arrow" Record optimization problem solution to an Arrow file. """ -function record(recorder::Recorder{ArrowFile}, id::UUID) +function record(recorder::Recorder{ArrowFile}, id::UUID; input=false) + _filename = input ? filename_input(recorder) : filename(recorder) + model = if length(recorder.primal_variables) > 0 owner_model(recorder.primal_variables[1]) elseif length(recorder.dual_variables) > 0 @@ -16,7 +18,7 @@ function record(recorder::Recorder{ArrowFile}, id::UUID) @error("Recorder has no variables") end return Arrow.append( - recorder.filename, + _filename, (; id=[id], zip( diff --git a/src/csvrecorder.jl b/src/csvrecorder.jl index a526755..8b54a9f 100644 --- a/src/csvrecorder.jl +++ b/src/csvrecorder.jl @@ -1,4 +1,4 @@ -abstract type CSVFile <: RecorderFile end +abstract type CSVFile <: FileType end Base.string(::Type{CSVFile}) = "csv" @@ -7,9 +7,10 @@ Base.string(::Type{CSVFile}) = "csv" Record optimization problem solution to a CSV file. """ -function record(recorder::Recorder{CSVFile}, id::UUID) - if !isfile(recorder.filename) - open(recorder.filename, "w") do f +function record(recorder::Recorder{CSVFile}, id::UUID; input=false) + _filename = input ? filename_input(recorder) : filename(recorder) + if !isfile(_filename) + open(_filename, "w") do f write(f, "id") for p in recorder.primal_variables write(f, ",$(name(p))") @@ -21,7 +22,7 @@ function record(recorder::Recorder{CSVFile}, id::UUID) write(f, "\n") end end - open(recorder.filename, "a") do f + open(_filename, "a") do f write(f, "$id") for p in recorder.primal_variables val = value.(p) diff --git a/src/datasetgen.jl b/src/datasetgen.jl index 8790182..fedba31 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -1,26 +1,46 @@ -abstract type RecorderFile end +abstract type FileType end + +mutable struct RecorderFile{T<:FileType} + filename::String +end + +filename(recorder_file::RecorderFile) = recorder_file.filename """ Recorder(filename; primal_variables=[], dual_variables=[], filterfn=(model)-> termination_status(model) == MOI.OPTIMAL) Recorder of optimization problem solutions. """ -mutable struct Recorder{T<:RecorderFile} - filename::String - primal_variables::AbstractArray{VariableRef} - dual_variables::AbstractArray{ConstraintRef} +mutable struct Recorder{T<:FileType} + recorder_file::RecorderFile{T} + recorder_file_input::RecorderFile{T} + primal_variables::Vector{VariableRef} + dual_variables::Vector{ConstraintRef} filterfn::Function function Recorder{T}( filename::String; + filename_input::String=filename * "_input_", primal_variables=[], dual_variables=[], filterfn=(model) -> termination_status(model) == MOI.OPTIMAL, - ) where {T<:RecorderFile} - return new{T}(filename, primal_variables, dual_variables, filterfn) + ) where {T<:FileType} + return new{T}(RecorderFile{T}(filename), RecorderFile{T}(filename_input), primal_variables, dual_variables, filterfn) end end +filename(recorder::Recorder) = filename(recorder.recorder_file) + +filename_input(recorder::Recorder) = filename(recorder.recorder_file_input) + +function set_primal_variable!(recorder::Recorder, p::Vector{VariableRef}) + recorder.primal_variables = p +end + +function set_dual_variable!(recorder::Recorder, p::Vector{ConstraintRef}) + recorder.dual_variables = p +end + abstract type AbstractProblemIterator end """ @@ -55,7 +75,7 @@ Save optimization problem instances to a file. """ function save( problem_iterator::ProblemIterator, filename::String, file_type::Type{T} -) where {T<:RecorderFile} +) where {T<:FileType} return save( (; id=problem_iterator.ids, @@ -107,12 +127,12 @@ function solve_and_record( end """ - solve_batch(problem_iterator::ProblemIterator, recorder::Recorder) + solve_batch(problem_iterator::AbstractProblemIterator, recorder) Solve a batch of optimization problems and record the solutions. """ function solve_batch( - problem_iterator::ProblemIterator, recorder::Recorder + problem_iterator::AbstractProblemIterator, recorder ) successfull_solves = sum( diff --git a/src/worst_case.jl b/src/worst_case.jl index 1fe2f46..0236e76 100644 --- a/src/worst_case.jl +++ b/src/worst_case.jl @@ -1,3 +1,25 @@ +""" + WorstCaseProblemIterator + +An iterator that iterates over a set of problems and solves the worst case +problem for each problem in the set. The worst case problem is defined as +the problem that maximizes the objective function over the set of problems. + +The iterator is initialized with a set of problem ids, a function that +returns the parameters for a given problem, a function that builds the +primal problem, a function that builds the dual problem, and a function +that sets the iterator for the dual problem. + +The iterator returns the number of problems that were solved. + +# Arguments +- `ids::Vector{UUID}`: A vector of problem ids. +- `parameters::Function`: A function `(model) -> Vector{JuMP.Variable}` that + returns the parameters for a given problem. +- `primal_builder!::Function`: A function `(model, parameters) -> modified model` that builds the primal problem for a given set of parameters. +- `set_iterator!::Function`: A function `(model, parameters, idx) -> modified model` that modifies the dual problem for a given set of parameters and problem index. +- `optimizer`: The optimizer to use for the primal and dual problem. # TODO: Duplicate argument +""" struct WorstCaseProblemIterator <: AbstractProblemIterator ids::Vector{UUID} parameters::Function @@ -15,6 +37,17 @@ struct WorstCaseProblemIterator <: AbstractProblemIterator end end +""" + solve_and_record + +Solves the worst case problem for a given problem iterator and records the +solution if it passes the filter function. + +# Arguments +- `problem_iterator::WorstCaseProblemIterator`: The problem iterator. +- `recorder::Recorder`: The recorder. +- `idx::Integer`: The index of the problem to solve. +""" function solve_and_record( problem_iterator::WorstCaseProblemIterator, recorder::Recorder, idx::Integer ) @@ -61,9 +94,14 @@ function solve_and_record( optimal_loads = value.(load_var_dual) optimal_dual_cost = JuMP.objective_value(jump_dual_model) + # Save input + recorder.primal_variables = load_var_dual + recorder.dual_variables = [] + record(recorder, problem_iterator.ids[idx]; input=true) + # Create final primal model and solve model = JuMP.Model(problem_iterator.optimizer) - problem_iterator.primal_builder!(jump_dual_model, optimal_loads) + problem_iterator.primal_builder!(jump_dual_model, optimal_loads; recorder=recorder) JuMP.optimize!(model) # Check if method was effective diff --git a/test/datasetgen.jl b/test/datasetgen.jl index af56a3f..6ea1c52 100644 --- a/test/datasetgen.jl +++ b/test/datasetgen.jl @@ -1,10 +1,10 @@ """ - testdataset_gen(path::AbstractString) + test_problem_iterator(path::AbstractString) Test dataset generation for different filetypes """ -function testdataset_gen(path::AbstractString) +function test_problem_iterator(path::AbstractString) @testset "Dataset Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] # The problem to iterate over model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) diff --git a/test/runtests.jl b/test/runtests.jl index d742466..c869715 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,19 +7,23 @@ using L2O import ParametricOptInterface as POI using Test using UUIDs +using Ipopt const test_dir = dirname(@__FILE__) const examples_dir = joinpath(test_dir, "..", "examples") include(joinpath(test_dir, "datasetgen.jl")) +include(joinpath(test_dir, "worst_case.jl")) + include(joinpath(examples_dir, "powermodels", "pglib_datagen.jl")) include(joinpath(examples_dir, "flux", "test_flux_forecaster.jl")) @testset "L2O.jl" begin mktempdir() do path - testdataset_gen(path) + test_problem_iterator(path) + test_worst_case_problem_iterator(path) file_in, file_out = test_pglib_datasetgen(path, "pglib_opf_case5_pjm", 20) test_flux_forecaster(file_in, file_out) end diff --git a/test/worst_case.jl b/test/worst_case.jl index f82c48c..627a073 100644 --- a/test/worst_case.jl +++ b/test/worst_case.jl @@ -1,108 +1,73 @@ -using PowerModels -using Ipopt -using PGLib -using JuMP -using Dualization -import JuMP.MOI as MOI -# Load data -case_name = "pglib_opf_case5_pjm" -matpower_case_name = case_name * ".m" -network_data_original = pglib(matpower_case_name) -network_data = deepcopy(network_data_original) - -# Case Parameters -network_formulation=DCPPowerModel -original_load = [l["pd"] for l in values(network_data["load"])] -max_total_load = sum(original_load) * 1.1 - -# Create primal model -solver = Ipopt.Optimizer -model = Model(solver) -load_var = @variable(model, load_var[i=1:length(original_load)]) -for (i, l) in enumerate(values(network_data["load"])) - l["pd"] = load_var[i] -end -load_moi_idx = MOI.VariableIndex[i.index for i in load_var] - -# Instantiate the model -pm = instantiate_model( - network_data, - network_formulation, - PowerModels.build_opf; - setting=Dict("output" => Dict("duals" => true)), - jump_model=model, -) - -# Dualize the model -dual_st = Dualization.dualize(JuMP.backend(model), - variable_parameters = load_moi_idx -) - -dual_model = dual_st.dual_model -primal_dual_map = dual_st.primal_dual_map - -jump_dual_model = JuMP.Model() -map_moi_to_jump = MOI.copy_to(JuMP.backend(jump_dual_model), dual_model) -set_optimizer(jump_dual_model, solver) - -load_dual_idxs = [map_moi_to_jump[primal_dual_map.primal_parameter[l]].value for l in load_moi_idx] - -load_var_dual = JuMP.all_variables(jump_dual_model)[load_dual_idxs] - -for (i,l) in enumerate(load_var_dual) - @constraint(jump_dual_model, l <= original_load[i]*1.5) - @constraint(jump_dual_model, l >= 0.0) -end -for (i, l) in enumerate(values(network_data["load"])) - l["pd"] = load_var_dual[i] +""" + test_worst_case_problem_iterator(path::AbstractString) + +Test dataset generation using the worst case problem iterator for different filetypes. +""" +function test_worst_case_problem_iterator(path::AbstractString) + @testset "Dataset Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] + # The problem to iterate over + optimizer = () -> Ipopt.Optimizer() + parameter_factory = (model) -> [@variable(model, _p)] + function primal_builder!(model, parameters; recorder=nothing) + @variable(model, x) + @constraint(model, cons, x + parameters[1] >= 3) + @objective(model, Min, 2x) + + if !isnothing(recorder) + set_primal_variable!(recorder, [x]) + set_dual_variable!(recorder, [cons]) + end + end + function set_iterator!(model, parameters, idx) + @constraint(model, parameters[1] <= 1.0 * idx) + @constraint(model, parameters[1] >= 0.0) + end + num_p = 10 + problem_iterator = WorstCaseProblemIterator( + [uuid1() for _ in 1:num_p], + parameter_factory, + primal_builder!, + set_iterator!, + optimizer, + ) + + # file_names + file_input = joinpath(path, "test_input.$(string(filetype))") + file_output = joinpath(path, "test_output.$(string(filetype))") + + # The recorder + recorder = Recorder{filetype}( + file_output; filename_input=file_input, + primal_variables=[], dual_variables=[] + ) + + # Solve all problems and record solutions + solve_batch(problem_iterator, recorder) + + # Check if file exists and has the correct number of rows and columns + @test isfile(file_input) + @test isfile(file_output) + if filetype == CSVFile + # test input file + @test length(readdlm(file_input, ',')[:, 1]) == num_p + 1 + @test length(readdlm(file_input, ',')[1, :]) == 2 + rm(file_input) + # test output file + @test length(readdlm(file_output, ',')[:, 1]) == num_p + 1 + @test length(readdlm(file_output, ',')[1, :]) == 4 + rm(file_output) + else + # test input file + df = Arrow.Table(file_input) + @test length(df) == 2 + @test length(df[1]) == num_p + rm(file_input) + # test output file + df = Arrow.Table(file_output) + @test length(df) == 4 + @test length(df[1]) == num_p + rm(file_output) + end + end end - -obj = objective_function(jump_dual_model) - -pm = instantiate_model( - network_data, - network_formulation, - PowerModels.build_opf; - setting=Dict("output" => Dict("duals" => true)), - jump_model=jump_dual_model, -) - -@objective(jump_dual_model, Max, obj) - -JuMP.optimize!(jump_dual_model) - -optimal_dual_cost = JuMP.objective_value(jump_dual_model) - -optimal_loads = value.(load_var_dual) - -load_diff = sum(abs.(optimal_loads - original_load)) ./ sum(original_load) - -# Create final primal model -network_data = deepcopy(network_data_original) -model = Model(solver) -for (i, l) in enumerate(values(network_data["load"])) - l["pd"] = optimal_loads[i] -end - -# Instantiate the model -pm = instantiate_model( - network_data, - network_formulation, - PowerModels.build_opf; - setting=Dict("output" => Dict("duals" => true)), - jump_model=model, -) - -# Solve the model -JuMP.optimize!(model) - -# Check the solution -optimal_final_cost = JuMP.objective_value(model) -termination_status = JuMP.termination_status(model) -solution_primal_status = JuMP.primal_status(model) -solution_dual_status = JuMP.dual_status(model) -termination_status == MOI.INFEASIBLE && @error("Optimal solution not found") -solution_primal_status != MOI.FEASIBLE_POINT && @error("Primal solution not found") -solution_dual_status != MOI.FEASIBLE_POINT && @error("Dual solution not found") -isapprox(optimal_final_cost, optimal_dual_cost; rtol=1e-4) || @warn("Final cost is not equal to dual cost") From 94343c3d138cd320ba07a50e0d796b605bc6997f Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Mon, 18 Sep 2023 15:01:00 -0400 Subject: [PATCH 06/36] fix --- src/L2O.jl | 3 ++- src/datasetgen.jl | 2 +- src/worst_case.jl | 4 ++-- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/L2O.jl b/src/L2O.jl index 94c3122..2ea3dad 100644 --- a/src/L2O.jl +++ b/src/L2O.jl @@ -8,7 +8,8 @@ using UUIDs import ParametricOptInterface as POI import Base: string -export ArrowFile, CSVFile, ProblemIterator, Recorder, save, solve_batch, WorstCaseProblemIterator +export ArrowFile, CSVFile, ProblemIterator, Recorder, save, solve_batch, + WorstCaseProblemIterator, set_primal_variable!, set_dual_variable! include("datasetgen.jl") include("csvrecorder.jl") diff --git a/src/datasetgen.jl b/src/datasetgen.jl index fedba31..216fd9c 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -37,7 +37,7 @@ function set_primal_variable!(recorder::Recorder, p::Vector{VariableRef}) recorder.primal_variables = p end -function set_dual_variable!(recorder::Recorder, p::Vector{ConstraintRef}) +function set_dual_variable!(recorder::Recorder, p::Vector) recorder.dual_variables = p end diff --git a/src/worst_case.jl b/src/worst_case.jl index 0236e76..afd68e1 100644 --- a/src/worst_case.jl +++ b/src/worst_case.jl @@ -77,7 +77,7 @@ function solve_and_record( load_var_dual = JuMP.all_variables(jump_dual_model)[load_dual_idxs] # Add constraints to the dual associated with the parameters - problem_iterator.set_iterator!(jump_dual_model, parameters, idx) + problem_iterator.set_iterator!(jump_dual_model, load_var_dual, idx) # Get the objective function obj = objective_function(jump_dual_model) @@ -101,7 +101,7 @@ function solve_and_record( # Create final primal model and solve model = JuMP.Model(problem_iterator.optimizer) - problem_iterator.primal_builder!(jump_dual_model, optimal_loads; recorder=recorder) + problem_iterator.primal_builder!(model, optimal_loads; recorder=recorder) JuMP.optimize!(model) # Check if method was effective From 3a8152a84bdd85f50aaf83286260b17fe6c5899f Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Mon, 18 Sep 2023 17:38:12 -0400 Subject: [PATCH 07/36] add tests --- src/arrowrecorder.jl | 29 +++++++++++++++++------------ src/csvrecorder.jl | 10 +++++++--- src/datasetgen.jl | 8 +++++++- test/worst_case.jl | 2 +- 4 files changed, 32 insertions(+), 17 deletions(-) diff --git a/src/arrowrecorder.jl b/src/arrowrecorder.jl index 0ed3420..23061ea 100644 --- a/src/arrowrecorder.jl +++ b/src/arrowrecorder.jl @@ -17,20 +17,25 @@ function record(recorder::Recorder{ArrowFile}, id::UUID; input=false) else @error("Recorder has no variables") end + + df = (; + id=[id], + zip( + Symbol.(name.(recorder.primal_variables)), + [[value.(p)] for p in recorder.primal_variables], + )..., + zip( + Symbol.("dual_" .* name.(recorder.dual_variables)), + [[dual.(p)] for p in recorder.dual_variables], + )..., + ) + if !input + df=merge(df, (;objective=[JuMP.objective_value(model)])) + end + return Arrow.append( _filename, - (; - id=[id], - zip( - Symbol.(name.(recorder.primal_variables)), - [[value.(p)] for p in recorder.primal_variables], - )..., - zip( - Symbol.("dual_" .* name.(recorder.dual_variables)), - [[dual.(p)] for p in recorder.dual_variables], - )..., - objective=[JuMP.objective_value(model)], - ), + df, ) end diff --git a/src/csvrecorder.jl b/src/csvrecorder.jl index 8b54a9f..bb8fcef 100644 --- a/src/csvrecorder.jl +++ b/src/csvrecorder.jl @@ -18,7 +18,9 @@ function record(recorder::Recorder{CSVFile}, id::UUID; input=false) for p in recorder.dual_variables write(f, ",dual_$(name(p))") end - write(f, ",objective") + if !input + write(f, ",objective") + end write(f, "\n") end end @@ -40,8 +42,10 @@ function record(recorder::Recorder{CSVFile}, id::UUID; input=false) else @error("Recorder has no variables") end - obj = JuMP.objective_value(model) - write(f, ",$obj") + if !input + obj = JuMP.objective_value(model) + write(f, ",$obj") + end # end line write(f, "\n") end diff --git a/src/datasetgen.jl b/src/datasetgen.jl index 216fd9c..bf78e5c 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -6,6 +6,12 @@ end filename(recorder_file::RecorderFile) = recorder_file.filename +termination_status_filter(status) = status == MOI.OPTIMAL || status == MOI.SLOW_PROGRESS || status == MOI.LOCALLY_SOLVED || status == MOI.ITERATION_LIMIT +primal_status_filter(status) = status == MOI.FEASIBLE_POINT +dual_status_filter(status) = status == MOI.FEASIBLE_POINT + +filter_fn(model) = termination_status_filter(termination_status(model)) && primal_status_filter(primal_status(model)) && dual_status_filter(dual_status(model)) + """ Recorder(filename; primal_variables=[], dual_variables=[], filterfn=(model)-> termination_status(model) == MOI.OPTIMAL) @@ -23,7 +29,7 @@ mutable struct Recorder{T<:FileType} filename_input::String=filename * "_input_", primal_variables=[], dual_variables=[], - filterfn=(model) -> termination_status(model) == MOI.OPTIMAL, + filterfn=filter_fn, ) where {T<:FileType} return new{T}(RecorderFile{T}(filename), RecorderFile{T}(filename_input), primal_variables, dual_variables, filterfn) end diff --git a/test/worst_case.jl b/test/worst_case.jl index 627a073..231866e 100644 --- a/test/worst_case.jl +++ b/test/worst_case.jl @@ -5,7 +5,7 @@ Test dataset generation using the worst case problem iterator for different filetypes. """ function test_worst_case_problem_iterator(path::AbstractString) - @testset "Dataset Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] + @testset "Worst Case Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] # The problem to iterate over optimizer = () -> Ipopt.Optimizer() parameter_factory = (model) -> [@variable(model, _p)] From bb772ce92ef4ebda8f73bf4a4f27f5ccb6277755 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Mon, 18 Sep 2023 19:02:07 -0400 Subject: [PATCH 08/36] tests passing --- examples/powermodels/pglib_datagen.jl | 101 +++++++++++++++++--------- src/datasetgen.jl | 6 +- 2 files changed, 70 insertions(+), 37 deletions(-) diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 01f92ab..338a1c7 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -37,6 +37,52 @@ function load_sampler( return load_samples end +function load_parameter_factory(model, indices; load_set=nothing) + if isnothing(load_set) + return @variable( + model, _p[i in indices] + ) + end + return @variable( + model, _p[i in indices] in load_set + ) +end + +function pm_primal_builder!(model, parameters, network_data, network_formulation; recorder=nothing) + num_loads = length(network_data["load"]) + for (i, l) in enumerate(values(network_data["load"])) + l["pd"] = parameters[i] + l["qd"] = parameters[num_loads+i] + end + + # Instantiate the model + pm = instantiate_model( + network_data, + network_formulation, + PowerModels.build_opf; + setting=Dict("output" => Dict("duals" => true)), + jump_model=model, + ) + + if !isnothing(recorder) + variable_refs = return_variablerefs(pm) + for variableref in variable_refs + set_name(variableref, replace(name(variableref), "," => "_")) + end + set_primal_variable!(recorder, variable_refs) + # set_dual_variable!(recorder, [cons for cons in values(pm["constraint"])]) + end + + return nothing +end + +function load_set_iterator!(model, parameters, idx, original_load) + for (i, p) in enumerate(parameters) + @constraint(model, p <= original_load[i] * (1.0 + 0.01 * idx)) + @constraint(model, p >= 0.0) + end +end + """ generate_dataset_pglib(data_dir::AbstractString, case_name::AbstractString; download_files::Bool=true, filetype::Type{FileType}, num_p::Int=10 @@ -48,7 +94,6 @@ function generate_dataset_pglib( data_dir, case_name; filetype=CSVFile, - download_files=true, num_p=10, load_sampler=load_sampler, network_formulation=DCPPowerModel, @@ -71,34 +116,29 @@ function generate_dataset_pglib( MOI.set(model, MOI.Silent(), true) # Save original load value and Link POI - original_load = [l["pd"] for l in values(network_data["load"])] - p = @variable( - model, _p[i=1:length(network_data["load"])] in POI.Parameter.(original_load) - ) # vector of parameters - for (i, l) in enumerate(values(network_data["load"])) - l["pd"] = p[i] - end - - # Instantiate the model - pm = instantiate_model( - network_data, - network_formulation, - PowerModels.build_opf; - setting=Dict("output" => Dict("duals" => true)), - jump_model=model, + num_loads = length(network_data["load"]) + original_load = vcat( + [l["pd"] for l in values(network_data["load"])], + [l["qd"] for l in values(network_data["load"])], ) + p = load_parameter_factory(model, 1:(num_loads * 2); load_set=POI.Parameter.(original_load)) + + # Define batch id + batch_id = string(uuid1()) + @info "Batch ID: $batch_id" + + # Build model and Recorder + file = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_output_" * batch_id * "." * string(filetype)) + recorder = Recorder{filetype}(file) + pm_primal_builder!(model, p, network_data, network_formulation; recorder=recorder) # The problem iterator problem_iterator = ProblemIterator( Dict( - p .=> [ - load_sampler(original_load[i], num_p) for - i in 1:length(network_data["load"]) - ], + p .=> load_sampler.(original_load, num_p), ), ) - batch_id = string(uuid1()) - @info "Batch ID: $batch_id" + save( problem_iterator, joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_input_" * batch_id * "." * string(filetype)), @@ -106,23 +146,16 @@ function generate_dataset_pglib( ) # Solve the problem and return the number of successfull solves - file = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_output_" * batch_id * "." * string(filetype)) - variable_refs = return_variablerefs(pm) - for variableref in variable_refs - set_name(variableref, replace(name(variableref), "," => "_")) - end - number_vars = length(variable_refs) - recorder = Recorder{filetype}(file; primal_variables=variable_refs) return solve_batch(problem_iterator, recorder), - number_vars, - length(original_load), - batch_id + length(recorder.primal_variables), + length(original_load), + batch_id end function test_pglib_datasetgen(path::AbstractString, case_name::AbstractString, num_p::Int) @testset "Dataset Generation pglib case" begin network_formulation = DCPPowerModel - success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib( + success_solves, number_variables, number_parameters, batch_id = generate_dataset_pglib( path, case_name; num_p=num_p, network_formulation=network_formulation ) file_in = joinpath(path, string(network_formulation), case_name * "_" * string(network_formulation) * "_input_" * batch_id * ".csv") @@ -130,7 +163,7 @@ function test_pglib_datasetgen(path::AbstractString, case_name::AbstractString, # Check if problem iterator was saved @test isfile(file_in) @test length(readdlm(file_in, ',')[:, 1]) == num_p + 1 - @test length(readdlm(file_in, ',')[1, :]) == 1 + number_loads + @test length(readdlm(file_in, ',')[1, :]) == 1 + number_parameters # Check if the number of successfull solves is equal to the number of problems saved @test isfile(file_out) diff --git a/src/datasetgen.jl b/src/datasetgen.jl index bf78e5c..b936a9b 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -20,8 +20,8 @@ Recorder of optimization problem solutions. mutable struct Recorder{T<:FileType} recorder_file::RecorderFile{T} recorder_file_input::RecorderFile{T} - primal_variables::Vector{VariableRef} - dual_variables::Vector{ConstraintRef} + primal_variables::Vector + dual_variables::Vector filterfn::Function function Recorder{T}( @@ -39,7 +39,7 @@ filename(recorder::Recorder) = filename(recorder.recorder_file) filename_input(recorder::Recorder) = filename(recorder.recorder_file_input) -function set_primal_variable!(recorder::Recorder, p::Vector{VariableRef}) +function set_primal_variable!(recorder::Recorder, p::Vector) recorder.primal_variables = p end From 43e4839aed5fcba65d01f7a6943c50579b37bcd7 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Mon, 18 Sep 2023 19:45:48 -0400 Subject: [PATCH 09/36] fix test --- examples/powermodels/pglib_datagen.jl | 88 ++++++++++++++++++++++++++- src/worst_case.jl | 2 +- test/runtests.jl | 1 + 3 files changed, 87 insertions(+), 4 deletions(-) diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 338a1c7..17eb08a 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -99,9 +99,6 @@ function generate_dataset_pglib( network_formulation=DCPPowerModel, solver = () -> POI.Optimizer(HiGHS.Optimizer()), ) - # Download file - matpower_case_name = case_name * ".m" - # save folder data_sim_dir = joinpath(data_dir, string(network_formulation)) if !isdir(data_sim_dir) @@ -109,6 +106,7 @@ function generate_dataset_pglib( end # Read data + matpower_case_name = case_name * ".m" network_data = make_basic_network(pglib(matpower_case_name)) # The problem to iterate over @@ -152,6 +150,65 @@ function generate_dataset_pglib( batch_id end +function generate_worst_case_dataset(data_dir, + case_name; + filetype=CSVFile, + num_p=10, + network_formulation=DCPPowerModel, + optimizer = () -> Ipopt.Optimizer(), +) + # save folder + data_sim_dir = joinpath(data_dir, string(network_formulation)) + if !isdir(data_sim_dir) + mkdir(data_sim_dir) + end + + # Read data + matpower_case_name = case_name * ".m" + network_data = make_basic_network(pglib(matpower_case_name)) + + # Parameter factory + num_loads = length(network_data["load"]) + original_load = vcat( + [l["pd"] for l in values(network_data["load"])], + [l["qd"] for l in values(network_data["load"])], + ) + parameter_factory = (model) -> load_parameter_factory(model, 1:(num_loads * 2)) + + # Define batch id + batch_id = string(uuid1()) + @info "Batch ID: $batch_id" + + # Build model + primal_builder! = (model, parameters; recorder=nothing) -> pm_primal_builder!(model, parameters, network_data, network_formulation; recorder=recorder) + + # Set iterator + set_iterator! = (model, parameters, idx) -> load_set_iterator!(model, parameters, idx, original_load) + + # The problem iterator + problem_iterator = WorstCaseProblemIterator( + [uuid1() for _ in 1:num_p], + parameter_factory, + primal_builder!, + set_iterator!, + optimizer, + ) + + # File names + file_input = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_input_" * batch_id * "." * string(filetype)) + file_output = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_output_" * batch_id * "." * string(filetype)) + recorder = Recorder{filetype}( + file_output; filename_input=file_input, + primal_variables=[], dual_variables=[] + ) + + # Solve all problems and record solutions + return solve_batch(problem_iterator, recorder), + length(recorder.primal_variables), + length(original_load), + batch_id +end + function test_pglib_datasetgen(path::AbstractString, case_name::AbstractString, num_p::Int) @testset "Dataset Generation pglib case" begin network_formulation = DCPPowerModel @@ -173,3 +230,28 @@ function test_pglib_datasetgen(path::AbstractString, case_name::AbstractString, return file_in, file_out end end + +function test_generate_worst_case_dataset(path::AbstractString, case_name::AbstractString, num_p::Int) + @testset "Worst Case Dataset Generation pglib case" begin + network_formulation = DCPPowerModel + # Improve dataset + success_solves, number_variables, number_parameters, batch_id = generate_worst_case_dataset( + path, case_name; num_p=num_p, network_formulation=network_formulation + ) + + file_in = joinpath(path, string(network_formulation), case_name * "_" * string(network_formulation) * "_input_" * batch_id * ".csv") + file_out = joinpath(path, string(network_formulation), case_name * "_" * string(network_formulation) * "_output_" * batch_id * ".csv") + + # Check if problem iterator was saved + @test isfile(file_in) + @test length(readdlm(file_in, ',')[:, 1]) == num_p + 1 + @test length(readdlm(file_in, ',')[1, :]) == 1 + number_parameters + + # Check if the number of successfull solves is equal to the number of problems saved + @test isfile(file_out) + @test length(readdlm(file_out, ',')[:, 1]) == num_p * success_solves + 1 + @test length(readdlm(file_out, ',')[1, :]) == number_variables + 2 + + return file_in, file_out + end +end diff --git a/src/worst_case.jl b/src/worst_case.jl index afd68e1..01127ee 100644 --- a/src/worst_case.jl +++ b/src/worst_case.jl @@ -57,7 +57,7 @@ function solve_and_record( problem_iterator.primal_builder!(model, parameters) # Parameter indices - load_moi_idx = MOI.VariableIndex[i.index for i in parameters] + load_moi_idx = JuMP.index.(parameters) # Dualize the model dual_st = Dualization.dualize(JuMP.backend(model), diff --git a/test/runtests.jl b/test/runtests.jl index c869715..1d1ebd3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -25,6 +25,7 @@ include(joinpath(examples_dir, "flux", "test_flux_forecaster.jl")) test_problem_iterator(path) test_worst_case_problem_iterator(path) file_in, file_out = test_pglib_datasetgen(path, "pglib_opf_case5_pjm", 20) + file_in, file_out = test_generate_worst_case_dataset(path, "pglib_opf_case5_pjm", 20) test_flux_forecaster(file_in, file_out) end end From ec0a4c52f16c5a78a7f38e26f3f0aa59ff6447fd Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Mon, 18 Sep 2023 19:55:21 -0400 Subject: [PATCH 10/36] tests passing --- src/worst_case.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/worst_case.jl b/src/worst_case.jl index 01127ee..70f2e60 100644 --- a/src/worst_case.jl +++ b/src/worst_case.jl @@ -57,7 +57,7 @@ function solve_and_record( problem_iterator.primal_builder!(model, parameters) # Parameter indices - load_moi_idx = JuMP.index.(parameters) + load_moi_idx = Vector{MOI.VariableIndex}(JuMP.index.(parameters)) # Dualize the model dual_st = Dualization.dualize(JuMP.backend(model), From 664a202ed273635e6521aba962d06b620bf5e4b7 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Mon, 18 Sep 2023 23:40:44 -0400 Subject: [PATCH 11/36] update code --- examples/flux/flux_forecaster_script.jl | 4 +-- .../generate_full_datasets_script.jl | 35 +++++++++++++++---- examples/powermodels/pglib_datagen.jl | 23 ++++++++---- src/worst_case.jl | 15 ++++---- test/worst_case.jl | 6 ++-- 5 files changed, 60 insertions(+), 23 deletions(-) diff --git a/examples/flux/flux_forecaster_script.jl b/examples/flux/flux_forecaster_script.jl index 2d6d670..310f973 100644 --- a/examples/flux/flux_forecaster_script.jl +++ b/examples/flux/flux_forecaster_script.jl @@ -8,9 +8,9 @@ using PowerModels # Paths path_dataset = joinpath(pwd(), "examples", "powermodels", "data") -case_name = "pglib_opf_case5_pjm" +case_name = "pglib_opf_case300_ieee" filetype = ArrowFile -network_formulation = DCPPowerModel +network_formulation = SOCWRConicPowerModel case_file_path = joinpath(path, case_name, string(network_formulation)) # Load input and output data tables diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 4475844..5c229a0 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -1,3 +1,4 @@ +# run with: julia ./examples/powermodels/generate_full_datasets_script.jl using TestEnv TestEnv.activate() @@ -9,6 +10,7 @@ using PowerModels using Clarabel import JuMP.MOI as MOI import ParametricOptInterface as POI +using Ipopt cached = MOI.Bridges.full_bridge_optimizer( MOI.Utilities.CachingOptimizer( @@ -32,15 +34,36 @@ filetype = ArrowFile case_name = "pglib_opf_case300_ieee" network_formulation = SOCWRConicPowerModel case_file_path = joinpath(path, case_name) -solver = () -> POI.Optimizer(cached()) +solver = () -> POI.Optimizer(cached) # Generate dataset -success_solves = 0.0 +# global success_solves = 0.0 +# for i in 1:num_batches +# _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; +# num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=solver, +# load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01) +# ) +# global success_solves += _success_solves +# end +# success_solves /= num_batches + +# @info "Success solves Normal: $(success_solves)" + +# Generate worst case dataset + +function optimizer_factory() + IPO_OPT = Ipopt.Optimizer() + IPO = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(IPO_OPT) + return () -> IPO +end + +global success_solves = 0.0 for i in 1:num_batches - _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; - num_p=num_p, filetype=filetype, network_formulation=network_formulation, solver=solver, - load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01) + _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset(case_file_path, case_name; + num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer_factory=optimizer_factory, ) - success_solves += _success_solves + global success_solves += _success_solves end success_solves /= num_batches + +@info "Success solves Worst Case: $(success_solves) of $(num_batches * num_p)" diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 17eb08a..ec77d2c 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -19,6 +19,13 @@ function return_variablerefs(pm::AbstractPowerModel) ) end +function add_names_pm(pm::AbstractPowerModel) + variable_refs = return_variablerefs(pm) + for variableref in variable_refs + set_name(variableref, replace(name(variableref), "," => "_")) + end +end + """ load_sampler(original_load::T, num_p::Int, max_multiplier::T=3.0, min_multiplier::T=0.0, step_multiplier::T=0.1) @@ -97,20 +104,18 @@ function generate_dataset_pglib( num_p=10, load_sampler=load_sampler, network_formulation=DCPPowerModel, - solver = () -> POI.Optimizer(HiGHS.Optimizer()), + optimizer = () -> POI.Optimizer(HiGHS.Optimizer()), ) # save folder data_sim_dir = joinpath(data_dir, string(network_formulation)) - if !isdir(data_sim_dir) - mkdir(data_sim_dir) - end + mkpath(data_sim_dir) # Read data matpower_case_name = case_name * ".m" network_data = make_basic_network(pglib(matpower_case_name)) # The problem to iterate over - model = Model(solver) + model = Model(optimizer) MOI.set(model, MOI.Silent(), true) # Save original load value and Link POI @@ -150,12 +155,16 @@ function generate_dataset_pglib( batch_id end +function default_optimizer_factory() + return () -> Ipopt.Optimizer() +end + function generate_worst_case_dataset(data_dir, case_name; filetype=CSVFile, num_p=10, network_formulation=DCPPowerModel, - optimizer = () -> Ipopt.Optimizer(), + optimizer_factory = default_optimizer_factory, ) # save folder data_sim_dir = joinpath(data_dir, string(network_formulation)) @@ -191,7 +200,7 @@ function generate_worst_case_dataset(data_dir, parameter_factory, primal_builder!, set_iterator!, - optimizer, + optimizer_factory, ) # File names diff --git a/src/worst_case.jl b/src/worst_case.jl index 70f2e60..59884cf 100644 --- a/src/worst_case.jl +++ b/src/worst_case.jl @@ -70,7 +70,7 @@ function solve_and_record( # Build Dual in JuMP jump_dual_model = JuMP.Model() map_moi_to_jump = MOI.copy_to(JuMP.backend(jump_dual_model), dual_model) - set_optimizer(jump_dual_model, problem_iterator.optimizer) + set_optimizer(jump_dual_model, problem_iterator.optimizer()) # Get dual variables for the parameters load_dual_idxs = [map_moi_to_jump[primal_dual_map.primal_parameter[l]].value for l in load_moi_idx] @@ -95,12 +95,15 @@ function solve_and_record( optimal_dual_cost = JuMP.objective_value(jump_dual_model) # Save input - recorder.primal_variables = load_var_dual - recorder.dual_variables = [] - record(recorder, problem_iterator.ids[idx]; input=true) - + if recorder.filterfn(jump_dual_model) + recorder.primal_variables = load_var_dual + recorder.dual_variables = [] + record(recorder, problem_iterator.ids[idx]; input=true) + else + return 0 + end # Create final primal model and solve - model = JuMP.Model(problem_iterator.optimizer) + model = JuMP.Model(problem_iterator.optimizer()) problem_iterator.primal_builder!(model, optimal_loads; recorder=recorder) JuMP.optimize!(model) diff --git a/test/worst_case.jl b/test/worst_case.jl index 231866e..81f252e 100644 --- a/test/worst_case.jl +++ b/test/worst_case.jl @@ -7,7 +7,9 @@ Test dataset generation using the worst case problem iterator for different file function test_worst_case_problem_iterator(path::AbstractString) @testset "Worst Case Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] # The problem to iterate over - optimizer = () -> Ipopt.Optimizer() + function optimizer_factory() + return () -> Ipopt.Optimizer() + end parameter_factory = (model) -> [@variable(model, _p)] function primal_builder!(model, parameters; recorder=nothing) @variable(model, x) @@ -29,7 +31,7 @@ function test_worst_case_problem_iterator(path::AbstractString) parameter_factory, primal_builder!, set_iterator!, - optimizer, + optimizer_factory, ) # file_names From a494f6b4ca9a5c965738d6e9509859a12b566cfb Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Tue, 19 Sep 2023 01:19:30 -0400 Subject: [PATCH 12/36] update script --- .../powermodels/generate_full_datasets_script.jl | 15 ++++++++++----- examples/powermodels/pglib_datagen.jl | 8 +++++--- src/worst_case.jl | 13 ++++++++++--- 3 files changed, 25 insertions(+), 11 deletions(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 5c229a0..dba6a54 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -10,7 +10,8 @@ using PowerModels using Clarabel import JuMP.MOI as MOI import ParametricOptInterface as POI -using Ipopt +using Gurobi +# using QuadraticToBinary cached = MOI.Bridges.full_bridge_optimizer( MOI.Utilities.CachingOptimizer( @@ -31,9 +32,11 @@ num_p = 10 filetype = ArrowFile # Case name -case_name = "pglib_opf_case300_ieee" +case_name = "pglib_opf_case5_pjm" # "pglib_opf_case300_ieee" network_formulation = SOCWRConicPowerModel case_file_path = joinpath(path, case_name) +mkpath(case_file_path) + solver = () -> POI.Optimizer(cached) # Generate dataset @@ -52,15 +55,17 @@ solver = () -> POI.Optimizer(cached) # Generate worst case dataset function optimizer_factory() - IPO_OPT = Ipopt.Optimizer() - IPO = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(IPO_OPT) - return () -> IPO + IPO_OPT = Gurobi.Optimizer() # MadNLP.Optimizer(print_level=MadNLP.INFO, max_iter=100) + # IPO = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(IPO_OPT) + # MIP = QuadraticToBinary.Optimizer{Float64}(IPO) + return () -> IPO_OPT end global success_solves = 0.0 for i in 1:num_batches _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset(case_file_path, case_name; num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer_factory=optimizer_factory, + hook = (model) -> set_optimizer_attribute(model, "NonConvex", 2) ) global success_solves += _success_solves end diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index ec77d2c..683ba76 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -85,8 +85,8 @@ end function load_set_iterator!(model, parameters, idx, original_load) for (i, p) in enumerate(parameters) - @constraint(model, p <= original_load[i] * (1.0 + 0.01 * idx)) - @constraint(model, p >= 0.0) + @constraint(model, p <= original_load[i] * (1.0 + 0.1 * idx)) + @constraint(model, p >= original_load[i] * (1.0 - 0.1 * idx)) end end @@ -165,6 +165,7 @@ function generate_worst_case_dataset(data_dir, num_p=10, network_formulation=DCPPowerModel, optimizer_factory = default_optimizer_factory, + hook = nothing ) # save folder data_sim_dir = joinpath(data_dir, string(network_formulation)) @@ -201,6 +202,7 @@ function generate_worst_case_dataset(data_dir, primal_builder!, set_iterator!, optimizer_factory, + hook ) # File names @@ -253,7 +255,7 @@ function test_generate_worst_case_dataset(path::AbstractString, case_name::Abstr # Check if problem iterator was saved @test isfile(file_in) - @test length(readdlm(file_in, ',')[:, 1]) == num_p + 1 + @test length(readdlm(file_in, ',')[:, 1]) == num_p * success_solves + 1 @test length(readdlm(file_in, ',')[1, :]) == 1 + number_parameters # Check if the number of successfull solves is equal to the number of problems saved diff --git a/src/worst_case.jl b/src/worst_case.jl index 59884cf..55f2092 100644 --- a/src/worst_case.jl +++ b/src/worst_case.jl @@ -26,14 +26,16 @@ struct WorstCaseProblemIterator <: AbstractProblemIterator primal_builder!::Function set_iterator!::Function optimizer + hook function WorstCaseProblemIterator( ids::Vector{UUID}, parameters::Function, primal_builder!::Function, set_iterator!::Function, optimizer, + hook=nothing, ) - return new(ids, parameters, primal_builder!, set_iterator!, optimizer) + return new(ids, parameters, primal_builder!, set_iterator!, optimizer, hook) end end @@ -71,6 +73,9 @@ function solve_and_record( jump_dual_model = JuMP.Model() map_moi_to_jump = MOI.copy_to(JuMP.backend(jump_dual_model), dual_model) set_optimizer(jump_dual_model, problem_iterator.optimizer()) + if !isnothing(problem_iterator.hook) + problem_iterator.hook(jump_dual_model) + end # Get dual variables for the parameters load_dual_idxs = [map_moi_to_jump[primal_dual_map.primal_parameter[l]].value for l in load_moi_idx] @@ -91,8 +96,6 @@ function solve_and_record( # Solve the dual JuMP.optimize!(jump_dual_model) - optimal_loads = value.(load_var_dual) - optimal_dual_cost = JuMP.objective_value(jump_dual_model) # Save input if recorder.filterfn(jump_dual_model) @@ -102,6 +105,10 @@ function solve_and_record( else return 0 end + + optimal_loads = value.(load_var_dual) + optimal_dual_cost = JuMP.objective_value(jump_dual_model) + # Create final primal model and solve model = JuMP.Model(problem_iterator.optimizer()) problem_iterator.primal_builder!(model, optimal_loads; recorder=recorder) From edb522f16a66c91e35212adbda123f755afae048 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Tue, 19 Sep 2023 09:15:55 -0400 Subject: [PATCH 13/36] update generator --- examples/powermodels/generate_full_datasets_script.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index dba6a54..4761fe4 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -33,7 +33,7 @@ filetype = ArrowFile # Case name case_name = "pglib_opf_case5_pjm" # "pglib_opf_case300_ieee" -network_formulation = SOCWRConicPowerModel +network_formulation = DCPPowerModel # SOCWRConicPowerModel case_file_path = joinpath(path, case_name) mkpath(case_file_path) From 7b10f846374d83d886cb528c433326b7b5d5a5c2 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Tue, 19 Sep 2023 13:31:34 -0400 Subject: [PATCH 14/36] update code --- Project.toml | 8 +- examples/powermodels/pglib_datagen.jl | 6 +- src/L2O.jl | 4 + src/datasetgen.jl | 2 +- src/worst_case.jl | 18 +++-- src/worst_case_iter.jl | 87 +++++++++++++++++++++ test/runtests.jl | 13 +++- test/worst_case.jl | 87 ++++++++++++++++++++- test/worst_case_old.jl | 108 ++++++++++++++++++++++++++ 9 files changed, 315 insertions(+), 18 deletions(-) create mode 100644 src/worst_case_iter.jl create mode 100644 test/worst_case_old.jl diff --git a/Project.toml b/Project.toml index e996c08..92bb902 100644 --- a/Project.toml +++ b/Project.toml @@ -8,8 +8,10 @@ Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" Dualization = "191a621a-6537-11e9-281d-650236a99e60" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +Nonconvex = "01bcebdf-4d21-426d-b5c4-6132c1619978" ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] Arrow = "2" @@ -20,6 +22,7 @@ ParametricOptInterface = "0.5" julia = "1.6" [extras] +AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918" Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" @@ -29,6 +32,9 @@ Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" PGLib = "07a8691f-3d11-4330-951b-3c50f98338be" PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +NonconvexBayesian = "fb352abc-de7b-48de-9ebd-665b54b5d9b3" +NonconvexIpopt = "bf347577-a06d-49ad-a669-8c0e005493b8" +KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392" [targets] -test = ["Test", "DelimitedFiles", "PGLib", "HiGHS", "PowerModels", "Flux", "DataFrames", "Clarabel", "Ipopt"] +test = ["Test", "DelimitedFiles", "PGLib", "HiGHS", "PowerModels", "Flux", "DataFrames", "Clarabel", "Ipopt", "NonconvexBayesian", "NonconvexIpopt", "KernelFunctions", "AbstractGPs"] diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 683ba76..0e67d4b 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -115,7 +115,7 @@ function generate_dataset_pglib( network_data = make_basic_network(pglib(matpower_case_name)) # The problem to iterate over - model = Model(optimizer) + model = JuMP.Model(optimizer) MOI.set(model, MOI.Silent(), true) # Save original load value and Link POI @@ -201,8 +201,8 @@ function generate_worst_case_dataset(data_dir, parameter_factory, primal_builder!, set_iterator!, - optimizer_factory, - hook + optimizer_factory; + hook=hook ) # File names diff --git a/src/L2O.jl b/src/L2O.jl index 2ea3dad..615a6b6 100644 --- a/src/L2O.jl +++ b/src/L2O.jl @@ -8,6 +8,9 @@ using UUIDs import ParametricOptInterface as POI import Base: string +using Nonconvex +using Zygote + export ArrowFile, CSVFile, ProblemIterator, Recorder, save, solve_batch, WorstCaseProblemIterator, set_primal_variable!, set_dual_variable! @@ -15,5 +18,6 @@ include("datasetgen.jl") include("csvrecorder.jl") include("arrowrecorder.jl") include("worst_case.jl") +include("worst_case_iter.jl") end diff --git a/src/datasetgen.jl b/src/datasetgen.jl index b936a9b..8f2cc22 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -80,7 +80,7 @@ end Save optimization problem instances to a file. """ function save( - problem_iterator::ProblemIterator, filename::String, file_type::Type{T} + problem_iterator::AbstractProblemIterator, filename::String, file_type::Type{T} ) where {T<:FileType} return save( (; diff --git a/src/worst_case.jl b/src/worst_case.jl index 55f2092..bab4216 100644 --- a/src/worst_case.jl +++ b/src/worst_case.jl @@ -20,22 +20,24 @@ The iterator returns the number of problems that were solved. - `set_iterator!::Function`: A function `(model, parameters, idx) -> modified model` that modifies the dual problem for a given set of parameters and problem index. - `optimizer`: The optimizer to use for the primal and dual problem. # TODO: Duplicate argument """ -struct WorstCaseProblemIterator <: AbstractProblemIterator +struct WorstCaseProblemIterator{F} <: AbstractProblemIterator ids::Vector{UUID} parameters::Function primal_builder!::Function set_iterator!::Function - optimizer - hook + optimizer::F + hook::Union{Nothing, Function} + options::Any function WorstCaseProblemIterator( ids::Vector{UUID}, parameters::Function, primal_builder!::Function, set_iterator!::Function, - optimizer, - hook=nothing, - ) - return new(ids, parameters, primal_builder!, set_iterator!, optimizer, hook) + optimizer::F; + hook::Union{Nothing, Function}=nothing, + options::Any=nothing + ) where {F} + new{F}(ids, parameters, primal_builder!, set_iterator!, optimizer, hook, options) end end @@ -59,7 +61,7 @@ function solve_and_record( problem_iterator.primal_builder!(model, parameters) # Parameter indices - load_moi_idx = Vector{MOI.VariableIndex}(JuMP.index.(parameters)) + load_moi_idx = Vector{MOI.VariableIndex}(JuMP.index.(parameters)) # Dualize the model dual_st = Dualization.dualize(JuMP.backend(model), diff --git a/src/worst_case_iter.jl b/src/worst_case_iter.jl new file mode 100644 index 0000000..d6278fd --- /dev/null +++ b/src/worst_case_iter.jl @@ -0,0 +1,87 @@ +# Nonconvex needs a minimization objective function that only receives the decision vector. +function primal_objective(parameter_values, parameters, filter_fn; penalty=1e8) + model = owner_model(first(parameters)) + for (i, p) in enumerate(parameters) + update_model!(model, p, parameter_values[i]) + end + + # Solve the model + JuMP.optimize!(model) + + status = filter_fn(model) + + objective = if status + JuMP.objective_value(model) + else + - penalty + end + + return objective +end + +function save_input(parameter_values, _recorder, id) + recorder = deepcopy(_recorder) + recorder.primal_variables = parameter_values + recorder.dual_variables = [] + record(recorder, id; input=true) + return nothing +end + +mutable struct StorageCallbackObjective <: Function + success_solves::Int + fcalls::Int + parameters::Vector{VariableRef} + filter_fn::Function + recorder::Recorder +end +StorageCallbackObjective(parameters, filter_fn, recorder) = StorageCallbackObjective(0, 0, + parameters, + filter_fn, + recorder +) + +function (callback::StorageCallbackObjective)(parameter_values) + Zygote.@ignore callback.fcalls += 1 + + obj = primal_objective(parameter_values, callback.parameters, callback.filter_fn) + + Zygote.@ignore obj > 0 && Zygote.@ignore callback.success_solves += 1 + Zygote.@ignore begin + id = uuid1(); record(callback.recorder, id); save_input(parameter_values, recorder, id) + end + Zygote.@ignore @info "Iter: $(callback.fcalls):" obj + return - obj +end + +function solve_and_record( + problem_iterator::WorstCaseProblemIterator{T}, recorder::Recorder, idx::Integer; maxiter=1000, +) where {T<:NonconvexCore.AbstractOptimizer} + # Build Primal + model, parameters = problem_iterator.primal_builder!(;recorder=recorder) + (min_demands, max_demands, max_total_volume, starting_point) = problem_iterator.set_iterator!(idx) + + storage_objective_function = StorageCallbackObjective( + parameters, recorder.filterfn, + recorder + ) + + # Build Nonconvex optimization model: + model_non = Nonconvex.Model() + set_objective!(model_non, storage_objective_function, flags = [:expensive]) + addvar!(model_non, min_demands, max_demands) + add_ineq_constraint!(model_non, x -> sum(x) - max_total_volume) + + # Solution Method: Bayesian Optimization + + # Optimize model_non: + if !isnothing(problem_iterator.options) + optimize(model_non, problem_iterator.algorithm, starting_point; options = problem_iterator.options) + else + optimize(model_non, problem_iterator.algorithm, starting_point) + end + + # best_solution = r_bayes.minimizer + # best_profit = -r_bayes.minimum + + return storage_objective_function.success_solves / storage_objective_function.fcalls +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 1d1ebd3..cdca740 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,11 @@ using Test using UUIDs using Ipopt +using NonconvexIpopt # Nonconvex.@load Ipopt +using NonconvexBayesian # Nonconvex.@load BayesOpt +using AbstractGPs +using KernelFunctions + const test_dir = dirname(@__FILE__) const examples_dir = joinpath(test_dir, "..", "examples") @@ -22,10 +27,10 @@ include(joinpath(examples_dir, "flux", "test_flux_forecaster.jl")) @testset "L2O.jl" begin mktempdir() do path - test_problem_iterator(path) + # test_problem_iterator(path) test_worst_case_problem_iterator(path) - file_in, file_out = test_pglib_datasetgen(path, "pglib_opf_case5_pjm", 20) - file_in, file_out = test_generate_worst_case_dataset(path, "pglib_opf_case5_pjm", 20) - test_flux_forecaster(file_in, file_out) + # file_in, file_out = test_pglib_datasetgen(path, "pglib_opf_case5_pjm", 20) + # file_in, file_out = test_generate_worst_case_dataset(path, "pglib_opf_case5_pjm", 20) + # test_flux_forecaster(file_in, file_out) end end diff --git a/test/worst_case.jl b/test/worst_case.jl index 81f252e..2665c7a 100644 --- a/test/worst_case.jl +++ b/test/worst_case.jl @@ -1,11 +1,25 @@ +_BayesOptAlg = () -> BayesOptAlg(IpoptAlg()) + +function _bayes_options(maxiter) + BayesOptOptions( + sub_options = IpoptOptions(max_iter = 20, print_level = 0), + # ninit=Int(floor(maxiter / 5)), + maxiter = maxiter, ftol = 1e-4, ctol = 1e-5, initialize=true, postoptimize=false, + kernel= RationalKernel(α=2.27e8) ∘ ScaleTransform(0.01), + noise=0.001, + std_multiple=8.67e4, + fit_prior=false # not working with custom priors + ) +end + """ test_worst_case_problem_iterator(path::AbstractString) Test dataset generation using the worst case problem iterator for different filetypes. """ function test_worst_case_problem_iterator(path::AbstractString) - @testset "Worst Case Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] + @testset "Worst Case Dual Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] # The problem to iterate over function optimizer_factory() return () -> Ipopt.Optimizer() @@ -72,4 +86,75 @@ function test_worst_case_problem_iterator(path::AbstractString) rm(file_output) end end + + @testset "Worst Case Bayes Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] + function _primal_builder!(;recorder=nothing) + model = JuMP.Model(() -> POI.Optimizer(HiGHS.Optimizer())) + parameters = @variable(model, _p in POI.Parameter(1.0)) + @variable(model, x) + @constraint(model, cons, x + parameters >= 3) + @objective(model, Min, 2x) + + if !isnothing(recorder) + set_primal_variable!(recorder, [x]) + set_dual_variable!(recorder, [cons]) + end + + return model, [parameters] + end + function _set_iterator!(idx) + min_demands = [0.0] + max_demands = [1.0 * idx] + max_total_volume = 1.0 * idx + starting_point = [1.0] + return min_demands, max_demands, max_total_volume, starting_point + end + num_p = 10 + problem_iterator = WorstCaseProblemIterator( + [uuid1() for _ in 1:num_p], # will be ignored + () -> nothing, # will be ignored + _primal_builder!, + _set_iterator!, + _BayesOptAlg(); + options = _bayes_options(num_p) + ) + + # file_names + file_input = joinpath(path, "test_input.$(string(filetype))") + file_output = joinpath(path, "test_output.$(string(filetype))") + + # The recorder + recorder = Recorder{filetype}( + file_output; filename_input=file_input, + primal_variables=[], dual_variables=[] + ) + + # Solve all problems and record solutions + solve_batch(problem_iterator, recorder) + + # Check if file exists and has the correct number of rows and columns + @test isfile(file_input) + @test isfile(file_output) + if filetype == CSVFile + # test input file + @test length(readdlm(file_input, ',')[:, 1]) == num_p + 1 + @test length(readdlm(file_input, ',')[1, :]) == 2 + rm(file_input) + # test output file + @test length(readdlm(file_output, ',')[:, 1]) == num_p + 1 + @test length(readdlm(file_output, ',')[1, :]) == 4 + rm(file_output) + else + # test input file + df = Arrow.Table(file_input) + @test length(df) == 2 + @test length(df[1]) == num_p + rm(file_input) + # test output file + df = Arrow.Table(file_output) + @test length(df) == 4 + @test length(df[1]) == num_p + rm(file_output) + end + end end diff --git a/test/worst_case_old.jl b/test/worst_case_old.jl new file mode 100644 index 0000000..42af394 --- /dev/null +++ b/test/worst_case_old.jl @@ -0,0 +1,108 @@ +using PowerModels +using Ipopt +using PGLib +using JuMP +using Dualization +import JuMP.MOI as MOI + +# Load data +case_name = "pglib_opf_case5_pjm" +matpower_case_name = case_name * ".m" +network_data_original = pglib(matpower_case_name) +network_data = deepcopy(network_data_original) + +# Case Parameters +network_formulation=DCPPowerModel +original_load = [l["pd"] for l in values(network_data["load"])] +max_total_load = sum(original_load) * 1.1 + +# Create primal model +solver = Ipopt.Optimizer +model = JuMP.Model(solver) +load_var = @variable(model, load_var[i=1:length(original_load)]) +for (i, l) in enumerate(values(network_data["load"])) + l["pd"] = load_var[i] +end +load_moi_idx = MOI.VariableIndex[i.index for i in load_var] + +# Instantiate the model +pm = instantiate_model( + network_data, + network_formulation, + PowerModels.build_opf; + setting=Dict("output" => Dict("duals" => true)), + jump_model=model, +) + +# Dualize the model +dual_st = Dualization.dualize(JuMP.backend(model), + variable_parameters = load_moi_idx +) + +dual_model = dual_st.dual_model +primal_dual_map = dual_st.primal_dual_map + +jump_dual_model = JuMP.Model() +map_moi_to_jump = MOI.copy_to(JuMP.backend(jump_dual_model), dual_model) +set_optimizer(jump_dual_model, solver) + +load_dual_idxs = [map_moi_to_jump[primal_dual_map.primal_parameter[l]].value for l in load_moi_idx] + +load_var_dual = JuMP.all_variables(jump_dual_model)[load_dual_idxs] + +for (i,l) in enumerate(load_var_dual) + @constraint(jump_dual_model, l <= original_load[i]*1.5) + @constraint(jump_dual_model, l >= 0.0) +end +for (i, l) in enumerate(values(network_data["load"])) + l["pd"] = load_var_dual[i] +end + +obj = objective_function(jump_dual_model) + +pm = instantiate_model( + network_data, + network_formulation, + PowerModels.build_opf; + setting=Dict("output" => Dict("duals" => true)), + jump_model=jump_dual_model, +) + +@objective(jump_dual_model, Max, obj) + +JuMP.optimize!(jump_dual_model) + +optimal_dual_cost = JuMP.objective_value(jump_dual_model) + +optimal_loads = value.(load_var_dual) + +load_diff = sum(abs.(optimal_loads - original_load)) ./ sum(original_load) + +# Create final primal model +network_data = deepcopy(network_data_original) +model = JuMP.Model(solver) +for (i, l) in enumerate(values(network_data["load"])) + l["pd"] = optimal_loads[i] +end + +# Instantiate the model +pm = instantiate_model( + network_data, + network_formulation, + PowerModels.build_opf; + setting=Dict("output" => Dict("duals" => true)), + jump_model=model, +) + +# Solve the model +JuMP.optimize!(model) + +# Check the solution +optimal_final_cost = JuMP.objective_value(model) +termination_status = JuMP.termination_status(model) +solution_primal_status = JuMP.primal_status(model) +solution_dual_status = JuMP.dual_status(model) +termination_status == MOI.INFEASIBLE && @error("Optimal solution not found") +solution_primal_status != MOI.FEASIBLE_POINT && @error("Primal solution not found") +solution_dual_status != MOI.FEASIBLE_POINT && @error("Dual solution not found") +isapprox(optimal_final_cost, optimal_dual_cost; rtol=1e-4) || @warn("Final cost is not equal to dual cost") From b07d392fb5e650017a1ff6cfeaa5e52871df4997 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Tue, 19 Sep 2023 13:47:48 -0400 Subject: [PATCH 15/36] tests passing bayes --- src/worst_case_iter.jl | 10 +++++----- test/worst_case.jl | 8 ++++---- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/worst_case_iter.jl b/src/worst_case_iter.jl index d6278fd..09be80c 100644 --- a/src/worst_case_iter.jl +++ b/src/worst_case_iter.jl @@ -19,9 +19,9 @@ function primal_objective(parameter_values, parameters, filter_fn; penalty=1e8) return objective end -function save_input(parameter_values, _recorder, id) +function save_input(parameters, _recorder, id) recorder = deepcopy(_recorder) - recorder.primal_variables = parameter_values + recorder.primal_variables = parameters recorder.dual_variables = [] record(recorder, id; input=true) return nothing @@ -47,7 +47,7 @@ function (callback::StorageCallbackObjective)(parameter_values) Zygote.@ignore obj > 0 && Zygote.@ignore callback.success_solves += 1 Zygote.@ignore begin - id = uuid1(); record(callback.recorder, id); save_input(parameter_values, recorder, id) + id = uuid1(); record(callback.recorder, id); save_input(callback.parameters, callback.recorder, id) end Zygote.@ignore @info "Iter: $(callback.fcalls):" obj return - obj @@ -75,9 +75,9 @@ function solve_and_record( # Optimize model_non: if !isnothing(problem_iterator.options) - optimize(model_non, problem_iterator.algorithm, starting_point; options = problem_iterator.options) + optimize(model_non, problem_iterator.optimizer, starting_point; options = problem_iterator.options) else - optimize(model_non, problem_iterator.algorithm, starting_point) + optimize(model_non, problem_iterator.optimizer, starting_point) end # best_solution = r_bayes.minimizer diff --git a/test/worst_case.jl b/test/worst_case.jl index 2665c7a..e0c6aaf 100644 --- a/test/worst_case.jl +++ b/test/worst_case.jl @@ -137,23 +137,23 @@ function test_worst_case_problem_iterator(path::AbstractString) @test isfile(file_output) if filetype == CSVFile # test input file - @test length(readdlm(file_input, ',')[:, 1]) == num_p + 1 + @test length(readdlm(file_input, ',')[:, 1]) >= num_p + 1 @test length(readdlm(file_input, ',')[1, :]) == 2 rm(file_input) # test output file - @test length(readdlm(file_output, ',')[:, 1]) == num_p + 1 + @test length(readdlm(file_output, ',')[:, 1]) >= num_p + 1 @test length(readdlm(file_output, ',')[1, :]) == 4 rm(file_output) else # test input file df = Arrow.Table(file_input) @test length(df) == 2 - @test length(df[1]) == num_p + @test length(df[1]) >= num_p rm(file_input) # test output file df = Arrow.Table(file_output) @test length(df) == 4 - @test length(df[1]) == num_p + @test length(df[1]) >= num_p rm(file_output) end end From d9bd22fe7cbfbf91e720622e5dbb2d85c2087e31 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Tue, 19 Sep 2023 14:35:11 -0400 Subject: [PATCH 16/36] add pglib_bayes --- Project.toml | 7 +- examples/powermodels/pglib_datagen.jl | 104 +++++++++++++++++++++++++- src/worst_case_iter.jl | 2 +- test/datasetgen.jl | 2 +- test/runtests.jl | 23 +++++- test/worst_case.jl | 43 ++++------- 6 files changed, 140 insertions(+), 41 deletions(-) diff --git a/Project.toml b/Project.toml index 92bb902..d7258dc 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" Dualization = "191a621a-6537-11e9-281d-650236a99e60" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Nonconvex = "01bcebdf-4d21-426d-b5c4-6132c1619978" ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" @@ -29,12 +30,12 @@ DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392" +NonconvexBayesian = "fb352abc-de7b-48de-9ebd-665b54b5d9b3" +NonconvexIpopt = "bf347577-a06d-49ad-a669-8c0e005493b8" PGLib = "07a8691f-3d11-4330-951b-3c50f98338be" PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -NonconvexBayesian = "fb352abc-de7b-48de-9ebd-665b54b5d9b3" -NonconvexIpopt = "bf347577-a06d-49ad-a669-8c0e005493b8" -KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392" [targets] test = ["Test", "DelimitedFiles", "PGLib", "HiGHS", "PowerModels", "Flux", "DataFrames", "Clarabel", "Ipopt", "NonconvexBayesian", "NonconvexIpopt", "KernelFunctions", "AbstractGPs"] diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 0e67d4b..1e39f8b 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -2,6 +2,7 @@ using PGLib using PowerModels using JuMP, HiGHS import ParametricOptInterface as POI +using LinearAlgebra """ return_variablerefs(pm::AbstractPowerModel) @@ -80,7 +81,7 @@ function pm_primal_builder!(model, parameters, network_data, network_formulation # set_dual_variable!(recorder, [cons for cons in values(pm["constraint"])]) end - return nothing + return model, parameters end function load_set_iterator!(model, parameters, idx, original_load) @@ -159,6 +160,78 @@ function default_optimizer_factory() return () -> Ipopt.Optimizer() end +function generate_worst_case_dataset_bayes(data_dir, + case_name; + filetype=CSVFile, + num_p=10, + network_formulation=DCPPowerModel, + optimizer = () -> POI.Optimizer(HiGHS.Optimizer()), +) + # save folder + data_sim_dir = joinpath(data_dir, string(network_formulation)) + if !isdir(data_sim_dir) + mkdir(data_sim_dir) + end + + # Read data + matpower_case_name = case_name * ".m" + network_data = make_basic_network(pglib(matpower_case_name)) + + # The problem to iterate over + model = JuMP.Model(optimizer) + MOI.set(model, MOI.Silent(), true) + + # Save original load value and Link POI + num_loads = length(network_data["load"]) + original_load = vcat( + [l["pd"] for l in values(network_data["load"])], + [l["qd"] for l in values(network_data["load"])], + ) + p = load_parameter_factory(model, 1:(num_loads * 2); load_set=POI.Parameter.(original_load)) + + # Define batch id + batch_id = string(uuid1()) + @info "Batch ID: $batch_id" + + # Build model + function _primal_builder!(;recorder=nothing) + pm_primal_builder!(model, p, network_data, network_formulation; recorder=recorder) + end + + # Set iterator + function _set_iterator!(idx) + min_demands = zeros(num_loads * 2) + max_demands = ones(num_loads * 2) .+ 0.1 * idx + max_total_volume = norm(original_load, 2) ^ 2 + starting_point = original_load + return min_demands, max_demands, max_total_volume, starting_point + end + + # The problem iterator + problem_iterator = WorstCaseProblemIterator( + [uuid1() for _ in 1:num_p], # will be ignored + () -> nothing, # will be ignored + _primal_builder!, + _set_iterator!, + _BayesOptAlg(); + options = _bayes_options(num_p) + ) + + # File names + file_input = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_input_" * batch_id * "." * string(filetype)) + file_output = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_output_" * batch_id * "." * string(filetype)) + recorder = Recorder{filetype}( + file_output; filename_input=file_input, + primal_variables=[], dual_variables=[] + ) + + # Solve all problems and record solutions + return solve_batch(problem_iterator, recorder), + length(recorder.primal_variables), + length(original_load), + batch_id +end + function generate_worst_case_dataset(data_dir, case_name; filetype=CSVFile, @@ -255,12 +328,37 @@ function test_generate_worst_case_dataset(path::AbstractString, case_name::Abstr # Check if problem iterator was saved @test isfile(file_in) - @test length(readdlm(file_in, ',')[:, 1]) == num_p * success_solves + 1 + @test length(readdlm(file_in, ',')[:, 1]) >= num_p * success_solves + 1 @test length(readdlm(file_in, ',')[1, :]) == 1 + number_parameters # Check if the number of successfull solves is equal to the number of problems saved @test isfile(file_out) - @test length(readdlm(file_out, ',')[:, 1]) == num_p * success_solves + 1 + @test length(readdlm(file_out, ',')[:, 1]) >= num_p * success_solves + 1 + @test length(readdlm(file_out, ',')[1, :]) == number_variables + 2 + + return file_in, file_out + end +end + +function test_generate_worst_case_dataset_bayes(path::AbstractString, case_name::AbstractString, num_p::Int) + @testset "Worst Case Dataset Generation pglib case" begin + network_formulation = DCPPowerModel + # Improve dataset + success_solves, number_variables, number_parameters, batch_id = generate_worst_case_dataset_bayes( + path, case_name; num_p=num_p, network_formulation=network_formulation + ) + + file_in = joinpath(path, string(network_formulation), case_name * "_" * string(network_formulation) * "_input_" * batch_id * ".csv") + file_out = joinpath(path, string(network_formulation), case_name * "_" * string(network_formulation) * "_output_" * batch_id * ".csv") + + # Check if problem iterator was saved + @test isfile(file_in) + @test length(readdlm(file_in, ',')[:, 1]) >= num_p * success_solves + 1 + @test length(readdlm(file_in, ',')[1, :]) == 1 + number_parameters + + # Check if the number of successfull solves is equal to the number of problems saved + @test isfile(file_out) + @test length(readdlm(file_out, ',')[:, 1]) >= num_p * success_solves + 1 @test length(readdlm(file_out, ',')[1, :]) == number_variables + 2 return file_in, file_out diff --git a/src/worst_case_iter.jl b/src/worst_case_iter.jl index 09be80c..e7906b3 100644 --- a/src/worst_case_iter.jl +++ b/src/worst_case_iter.jl @@ -69,7 +69,7 @@ function solve_and_record( model_non = Nonconvex.Model() set_objective!(model_non, storage_objective_function, flags = [:expensive]) addvar!(model_non, min_demands, max_demands) - add_ineq_constraint!(model_non, x -> sum(x) - max_total_volume) + add_ineq_constraint!(model_non, x -> sum(x .^ 2) - max_total_volume ^ 2) # Solution Method: Bayesian Optimization diff --git a/test/datasetgen.jl b/test/datasetgen.jl index 6ea1c52..e62070c 100644 --- a/test/datasetgen.jl +++ b/test/datasetgen.jl @@ -7,7 +7,7 @@ Test dataset generation for different filetypes function test_problem_iterator(path::AbstractString) @testset "Dataset Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] # The problem to iterate over - model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + model = JuMP.Model(() -> POI.Optimizer(HiGHS.Optimizer())) @variable(model, x) p = @variable(model, _p in POI.Parameter(1.0)) @constraint(model, cons, x + _p >= 3) diff --git a/test/runtests.jl b/test/runtests.jl index cdca740..e647e26 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,20 @@ using KernelFunctions const test_dir = dirname(@__FILE__) const examples_dir = joinpath(test_dir, "..", "examples") +_BayesOptAlg = () -> BayesOptAlg(IpoptAlg()) + +function _bayes_options(maxiter) + BayesOptOptions( + sub_options = IpoptOptions(max_iter = 20, print_level = 0), + # ninit=Int(floor(maxiter / 5)), + maxiter = maxiter, ftol = 1e-4, ctol = 1e-5, initialize=true, postoptimize=false, + kernel= RationalKernel(α=2.27e8) ∘ ScaleTransform(0.01), + noise=0.001, + std_multiple=8.67e4, + fit_prior=false # not working with custom priors + ) +end + include(joinpath(test_dir, "datasetgen.jl")) include(joinpath(test_dir, "worst_case.jl")) @@ -27,10 +41,11 @@ include(joinpath(examples_dir, "flux", "test_flux_forecaster.jl")) @testset "L2O.jl" begin mktempdir() do path - # test_problem_iterator(path) + test_problem_iterator(path) test_worst_case_problem_iterator(path) - # file_in, file_out = test_pglib_datasetgen(path, "pglib_opf_case5_pjm", 20) - # file_in, file_out = test_generate_worst_case_dataset(path, "pglib_opf_case5_pjm", 20) - # test_flux_forecaster(file_in, file_out) + file_in, file_out = test_pglib_datasetgen(path, "pglib_opf_case5_pjm", 20) + file_in, file_out = test_generate_worst_case_dataset(path, "pglib_opf_case5_pjm", 20) + file_in, file_out = test_generate_worst_case_dataset_bayes(path, "pglib_opf_case5_pjm", 20) + test_flux_forecaster(file_in, file_out) end end diff --git a/test/worst_case.jl b/test/worst_case.jl index e0c6aaf..5b42614 100644 --- a/test/worst_case.jl +++ b/test/worst_case.jl @@ -1,24 +1,9 @@ - -_BayesOptAlg = () -> BayesOptAlg(IpoptAlg()) - -function _bayes_options(maxiter) - BayesOptOptions( - sub_options = IpoptOptions(max_iter = 20, print_level = 0), - # ninit=Int(floor(maxiter / 5)), - maxiter = maxiter, ftol = 1e-4, ctol = 1e-5, initialize=true, postoptimize=false, - kernel= RationalKernel(α=2.27e8) ∘ ScaleTransform(0.01), - noise=0.001, - std_multiple=8.67e4, - fit_prior=false # not working with custom priors - ) -end - """ test_worst_case_problem_iterator(path::AbstractString) Test dataset generation using the worst case problem iterator for different filetypes. """ -function test_worst_case_problem_iterator(path::AbstractString) +function test_worst_case_problem_iterator(path::AbstractString, num_p=10) @testset "Worst Case Dual Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] # The problem to iterate over function optimizer_factory() @@ -39,7 +24,7 @@ function test_worst_case_problem_iterator(path::AbstractString) @constraint(model, parameters[1] <= 1.0 * idx) @constraint(model, parameters[1] >= 0.0) end - num_p = 10 + problem_iterator = WorstCaseProblemIterator( [uuid1() for _ in 1:num_p], parameter_factory, @@ -59,30 +44,30 @@ function test_worst_case_problem_iterator(path::AbstractString) ) # Solve all problems and record solutions - solve_batch(problem_iterator, recorder) + successfull_solves = solve_batch(problem_iterator, recorder) # Check if file exists and has the correct number of rows and columns @test isfile(file_input) @test isfile(file_output) if filetype == CSVFile # test input file - @test length(readdlm(file_input, ',')[:, 1]) == num_p + 1 + @test length(readdlm(file_input, ',')[:, 1]) >= num_p *successfull_solves + 1 @test length(readdlm(file_input, ',')[1, :]) == 2 rm(file_input) # test output file - @test length(readdlm(file_output, ',')[:, 1]) == num_p + 1 + @test length(readdlm(file_output, ',')[:, 1]) == num_p * successfull_solves + 1 @test length(readdlm(file_output, ',')[1, :]) == 4 rm(file_output) else # test input file df = Arrow.Table(file_input) @test length(df) == 2 - @test length(df[1]) == num_p + @test length(df[1]) >= num_p * successfull_solves rm(file_input) # test output file df = Arrow.Table(file_output) @test length(df) == 4 - @test length(df[1]) == num_p + @test length(df[1]) >= num_p rm(file_output) end end @@ -109,14 +94,14 @@ function test_worst_case_problem_iterator(path::AbstractString) starting_point = [1.0] return min_demands, max_demands, max_total_volume, starting_point end - num_p = 10 + problem_iterator = WorstCaseProblemIterator( [uuid1() for _ in 1:num_p], # will be ignored () -> nothing, # will be ignored _primal_builder!, _set_iterator!, _BayesOptAlg(); - options = _bayes_options(num_p) + options = _bayes_options(floor(Int, num_p / 5)) ) # file_names @@ -130,30 +115,30 @@ function test_worst_case_problem_iterator(path::AbstractString) ) # Solve all problems and record solutions - solve_batch(problem_iterator, recorder) + successfull_solves = solve_batch(problem_iterator, recorder) # Check if file exists and has the correct number of rows and columns @test isfile(file_input) @test isfile(file_output) if filetype == CSVFile # test input file - @test length(readdlm(file_input, ',')[:, 1]) >= num_p + 1 + @test length(readdlm(file_input, ',')[:, 1]) >= num_p * successfull_solves + 1 @test length(readdlm(file_input, ',')[1, :]) == 2 rm(file_input) # test output file - @test length(readdlm(file_output, ',')[:, 1]) >= num_p + 1 + @test length(readdlm(file_output, ',')[:, 1]) >= num_p * successfull_solves + 1 @test length(readdlm(file_output, ',')[1, :]) == 4 rm(file_output) else # test input file df = Arrow.Table(file_input) @test length(df) == 2 - @test length(df[1]) >= num_p + @test length(df[1]) >= num_p * successfull_solves rm(file_input) # test output file df = Arrow.Table(file_output) @test length(df) == 4 - @test length(df[1]) >= num_p + @test length(df[1]) >= num_p * successfull_solves rm(file_output) end end From be9efe61f64191f3a6057fc82dca4b2ddf6a8cc8 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Tue, 19 Sep 2023 14:41:57 -0400 Subject: [PATCH 17/36] update dataset generation --- .../generate_full_datasets_script.jl | 31 +++++++++++++------ src/worst_case_iter.jl | 6 ++-- 2 files changed, 25 insertions(+), 12 deletions(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 4761fe4..71188be 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -33,7 +33,7 @@ filetype = ArrowFile # Case name case_name = "pglib_opf_case5_pjm" # "pglib_opf_case300_ieee" -network_formulation = DCPPowerModel # SOCWRConicPowerModel +network_formulation = SOCWRConicPowerModel # SOCWRConicPowerModel # DCPPowerModel case_file_path = joinpath(path, case_name) mkpath(case_file_path) @@ -54,18 +54,29 @@ solver = () -> POI.Optimizer(cached) # Generate worst case dataset -function optimizer_factory() - IPO_OPT = Gurobi.Optimizer() # MadNLP.Optimizer(print_level=MadNLP.INFO, max_iter=100) - # IPO = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(IPO_OPT) - # MIP = QuadraticToBinary.Optimizer{Float64}(IPO) - return () -> IPO_OPT -end +# function optimizer_factory() +# IPO_OPT = Gurobi.Optimizer() # MadNLP.Optimizer(print_level=MadNLP.INFO, max_iter=100) +# # IPO = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(IPO_OPT) +# # MIP = QuadraticToBinary.Optimizer{Float64}(IPO) +# return () -> IPO_OPT +# end + +# global success_solves = 0.0 +# for i in 1:num_batches +# _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset(case_file_path, case_name; +# num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer_factory=optimizer_factory, +# hook = (model) -> set_optimizer_attribute(model, "NonConvex", 2) +# ) +# global success_solves += _success_solves +# end +# success_solves /= num_batches + +# @info "Success solves Worst Case: $(success_solves) of $(num_batches * num_p)" global success_solves = 0.0 for i in 1:num_batches - _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset(case_file_path, case_name; - num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer_factory=optimizer_factory, - hook = (model) -> set_optimizer_attribute(model, "NonConvex", 2) + _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset_bayes(case_file_path, case_name; + num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=solver, ) global success_solves += _success_solves end diff --git a/src/worst_case_iter.jl b/src/worst_case_iter.jl index e7906b3..23afa3b 100644 --- a/src/worst_case_iter.jl +++ b/src/worst_case_iter.jl @@ -45,9 +45,11 @@ function (callback::StorageCallbackObjective)(parameter_values) obj = primal_objective(parameter_values, callback.parameters, callback.filter_fn) - Zygote.@ignore obj > 0 && Zygote.@ignore callback.success_solves += 1 Zygote.@ignore begin - id = uuid1(); record(callback.recorder, id); save_input(callback.parameters, callback.recorder, id) + if obj > 0 + callback.success_solves += 1 + id = uuid1(); record(callback.recorder, id); save_input(callback.parameters, callback.recorder, id) + end end Zygote.@ignore @info "Iter: $(callback.fcalls):" obj return - obj From ee1c4c019c219b6d665bc8f487f067583f9c3b09 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Tue, 19 Sep 2023 16:29:38 -0400 Subject: [PATCH 18/36] working --- .../generate_full_datasets_script.jl | 36 +++++++++++++++---- examples/powermodels/pglib_datagen.jl | 28 +++++++++------ src/L2O.jl | 1 + src/datasetgen.jl | 14 +++++++- src/worst_case_iter.jl | 15 +++++--- 5 files changed, 71 insertions(+), 23 deletions(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 71188be..2913312 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -11,9 +11,17 @@ using Clarabel import JuMP.MOI as MOI import ParametricOptInterface as POI using Gurobi + +using NonconvexIpopt # Nonconvex.@load Ipopt +using NonconvexBayesian # Nonconvex.@load BayesOpt +using AbstractGPs +using KernelFunctions + # using QuadraticToBinary -cached = MOI.Bridges.full_bridge_optimizer( +########## SOLVERS ########## + +cached = () -> MOI.Bridges.full_bridge_optimizer( MOI.Utilities.CachingOptimizer( MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), Clarabel.Optimizer(), @@ -21,6 +29,24 @@ cached = MOI.Bridges.full_bridge_optimizer( Float64, ) +POI_cached_optimizer() = POI.Optimizer(cached()) + +_BayesOptAlg = () -> BayesOptAlg(IpoptAlg()) + +function _bayes_options(maxiter) + BayesOptOptions( + sub_options = IpoptOptions(max_iter = 20, print_level = 0), + # ninit=Int(floor(maxiter / 5)), + maxiter = maxiter, ftol = 1e-4, ctol = 1e-5, initialize=true, postoptimize=false, + kernel= RationalKernel(α=2.27e8) ∘ ScaleTransform(0.01), + noise=0.001, + std_multiple=8.67e4, + fit_prior=false # not working with custom priors + ) +end + +########## DATASET GENERATION ########## + # Paths path_powermodels = joinpath(pwd(), "examples", "powermodels") path = joinpath(path_powermodels, "data") @@ -37,13 +63,11 @@ network_formulation = SOCWRConicPowerModel # SOCWRConicPowerModel # DCPPowerMode case_file_path = joinpath(path, case_name) mkpath(case_file_path) -solver = () -> POI.Optimizer(cached) - # Generate dataset # global success_solves = 0.0 # for i in 1:num_batches # _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; -# num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=solver, +# num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, # load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01) # ) # global success_solves += _success_solves @@ -76,10 +100,10 @@ solver = () -> POI.Optimizer(cached) global success_solves = 0.0 for i in 1:num_batches _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset_bayes(case_file_path, case_name; - num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=solver, + num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, ) global success_solves += _success_solves end success_solves /= num_batches -@info "Success solves Worst Case: $(success_solves) of $(num_batches * num_p)" +@info "Success solves Worst Case: $(success_solves * 100) of $(num_batches * num_p)" diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 1e39f8b..ea42edd 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -79,9 +79,10 @@ function pm_primal_builder!(model, parameters, network_data, network_formulation end set_primal_variable!(recorder, variable_refs) # set_dual_variable!(recorder, [cons for cons in values(pm["constraint"])]) + return model, parameters, variable_refs end - return model, parameters + return nothing end function load_set_iterator!(model, parameters, idx, original_load) @@ -178,7 +179,7 @@ function generate_worst_case_dataset_bayes(data_dir, network_data = make_basic_network(pglib(matpower_case_name)) # The problem to iterate over - model = JuMP.Model(optimizer) + model = JuMP.Model(() -> optimizer()) MOI.set(model, MOI.Silent(), true) # Save original load value and Link POI @@ -193,9 +194,22 @@ function generate_worst_case_dataset_bayes(data_dir, batch_id = string(uuid1()) @info "Batch ID: $batch_id" + # File names + file_input = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_input_" * batch_id * "." * string(filetype)) + file_output = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_output_" * batch_id * "." * string(filetype)) + recorder = Recorder{filetype}( + file_output; filename_input=file_input, + primal_variables=[], dual_variables=[] + ) + # Build model + model, parameters, variable_refs = pm_primal_builder!(model, p, network_data, network_formulation; recorder=recorder) function _primal_builder!(;recorder=nothing) - pm_primal_builder!(model, p, network_data, network_formulation; recorder=recorder) + if !isnothing(recorder) + set_primal_variable!(recorder, variable_refs) + end + + return model, parameters end # Set iterator @@ -217,14 +231,6 @@ function generate_worst_case_dataset_bayes(data_dir, options = _bayes_options(num_p) ) - # File names - file_input = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_input_" * batch_id * "." * string(filetype)) - file_output = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_output_" * batch_id * "." * string(filetype)) - recorder = Recorder{filetype}( - file_output; filename_input=file_input, - primal_variables=[], dual_variables=[] - ) - # Solve all problems and record solutions return solve_batch(problem_iterator, recorder), length(recorder.primal_variables), diff --git a/src/L2O.jl b/src/L2O.jl index 615a6b6..a0a0b2f 100644 --- a/src/L2O.jl +++ b/src/L2O.jl @@ -6,6 +6,7 @@ using Dualization using JuMP using UUIDs import ParametricOptInterface as POI +import JuMP.MOI as MOI import Base: string using Nonconvex diff --git a/src/datasetgen.jl b/src/datasetgen.jl index 8f2cc22..b72b810 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -36,8 +36,20 @@ mutable struct Recorder{T<:FileType} end filename(recorder::Recorder) = filename(recorder.recorder_file) - filename_input(recorder::Recorder) = filename(recorder.recorder_file_input) +get_primal_variables(recorder::Recorder) = recorder.primal_variables +get_dual_variables(recorder::Recorder) = recorder.dual_variables +get_filterfn(recorder::Recorder) = recorder.filterfn + +function similar(recorder::Recorder{T}) where {T<:FileType} + return Recorder{T}( + filename(recorder); + filename_input=filename_input(recorder), + primal_variables=get_primal_variables(recorder), + dual_variables=get_dual_variables(recorder), + filterfn=get_filterfn(recorder), + ) +end function set_primal_variable!(recorder::Recorder, p::Vector) recorder.primal_variables = p diff --git a/src/worst_case_iter.jl b/src/worst_case_iter.jl index 23afa3b..1a0d2f1 100644 --- a/src/worst_case_iter.jl +++ b/src/worst_case_iter.jl @@ -6,9 +6,14 @@ function primal_objective(parameter_values, parameters, filter_fn; penalty=1e8) end # Solve the model - JuMP.optimize!(model) + status = false + try + JuMP.optimize!(model) - status = filter_fn(model) + status = filter_fn(model) + catch e + @warn "Error in solve_and_record: $e" + end objective = if status JuMP.objective_value(model) @@ -20,9 +25,9 @@ function primal_objective(parameter_values, parameters, filter_fn; penalty=1e8) end function save_input(parameters, _recorder, id) - recorder = deepcopy(_recorder) - recorder.primal_variables = parameters - recorder.dual_variables = [] + recorder = similar(_recorder) + set_primal_variable!(recorder, parameters) + set_dual_variable!(recorder, []) record(recorder, id; input=true) return nothing end From 3649c019d9e7e3aa6658a28cd97b96e3ff6a06a5 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Tue, 19 Sep 2023 17:36:35 -0400 Subject: [PATCH 19/36] end --- examples/powermodels/generate_full_datasets_script.jl | 6 +++--- examples/powermodels/pglib_datagen.jl | 6 +++--- src/worst_case.jl | 6 ++++-- src/worst_case_iter.jl | 8 ++++++-- 4 files changed, 16 insertions(+), 10 deletions(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 2913312..b48a6a0 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -37,7 +37,7 @@ function _bayes_options(maxiter) BayesOptOptions( sub_options = IpoptOptions(max_iter = 20, print_level = 0), # ninit=Int(floor(maxiter / 5)), - maxiter = maxiter, ftol = 1e-4, ctol = 1e-5, initialize=true, postoptimize=false, + maxiter = maxiter * 100, ftol = 1e-4, ctol = 1e-5, initialize=true, postoptimize=false, kernel= RationalKernel(α=2.27e8) ∘ ScaleTransform(0.01), noise=0.001, std_multiple=8.67e4, @@ -54,11 +54,11 @@ include(joinpath(path_powermodels, "pglib_datagen.jl")) # Parameters num_batches = 1 -num_p = 10 +num_p = 5 filetype = ArrowFile # Case name -case_name = "pglib_opf_case5_pjm" # "pglib_opf_case300_ieee" +case_name = "pglib_opf_case300_ieee" # "pglib_opf_case300_ieee" network_formulation = SOCWRConicPowerModel # SOCWRConicPowerModel # DCPPowerModel case_file_path = joinpath(path, case_name) mkpath(case_file_path) diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index ea42edd..63d6b5c 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -215,9 +215,9 @@ function generate_worst_case_dataset_bayes(data_dir, # Set iterator function _set_iterator!(idx) min_demands = zeros(num_loads * 2) - max_demands = ones(num_loads * 2) .+ 0.1 * idx - max_total_volume = norm(original_load, 2) ^ 2 - starting_point = original_load + max_demands = original_load .+ ones(num_loads * 2) .* 0.1 * idx + max_total_volume = norm(max_demands, 2) ^ 2 + starting_point = original_load ./ 10 return min_demands, max_demands, max_total_volume, starting_point end diff --git a/src/worst_case.jl b/src/worst_case.jl index bab4216..9f97b07 100644 --- a/src/worst_case.jl +++ b/src/worst_case.jl @@ -28,6 +28,7 @@ struct WorstCaseProblemIterator{F} <: AbstractProblemIterator optimizer::F hook::Union{Nothing, Function} options::Any + ext::Dict function WorstCaseProblemIterator( ids::Vector{UUID}, parameters::Function, @@ -35,9 +36,10 @@ struct WorstCaseProblemIterator{F} <: AbstractProblemIterator set_iterator!::Function, optimizer::F; hook::Union{Nothing, Function}=nothing, - options::Any=nothing + options::Any=nothing, + ext::Dict=Dict() ) where {F} - new{F}(ids, parameters, primal_builder!, set_iterator!, optimizer, hook, options) + new{F}(ids, parameters, primal_builder!, set_iterator!, optimizer, hook, options, ext) end end diff --git a/src/worst_case_iter.jl b/src/worst_case_iter.jl index 1a0d2f1..a7fc9f7 100644 --- a/src/worst_case_iter.jl +++ b/src/worst_case_iter.jl @@ -72,6 +72,10 @@ function solve_and_record( recorder ) + if haskey(problem_iterator.ext, :best_solution) + starting_point = problem_iterator.ext[:best_solution] + end + # Build Nonconvex optimization model: model_non = Nonconvex.Model() set_objective!(model_non, storage_objective_function, flags = [:expensive]) @@ -81,13 +85,13 @@ function solve_and_record( # Solution Method: Bayesian Optimization # Optimize model_non: - if !isnothing(problem_iterator.options) + r_bayes = if !isnothing(problem_iterator.options) optimize(model_non, problem_iterator.optimizer, starting_point; options = problem_iterator.options) else optimize(model_non, problem_iterator.optimizer, starting_point) end - # best_solution = r_bayes.minimizer + problem_iterator.ext[:best_solution] = r_bayes.minimizer # best_profit = -r_bayes.minimum return storage_objective_function.success_solves / storage_objective_function.fcalls From 69a68f1765b4a809fe74b7d2ad1b98251d4086ab Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Tue, 19 Sep 2023 19:40:36 -0400 Subject: [PATCH 20/36] update code --- .../powermodels/generate_full_datasets_script.jl | 6 ++++++ examples/powermodels/pglib_datagen.jl | 16 ++++++++++------ src/worst_case_iter.jl | 4 +--- 3 files changed, 17 insertions(+), 9 deletions(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index b48a6a0..436b822 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -17,6 +17,8 @@ using NonconvexBayesian # Nonconvex.@load BayesOpt using AbstractGPs using KernelFunctions +using NonconvexNLopt + # using QuadraticToBinary ########## SOLVERS ########## @@ -45,6 +47,9 @@ function _bayes_options(maxiter) ) end +algorithm = NLoptAlg(:LN_BOBYQA) +options = NLoptOptions(maxeval=30) + ########## DATASET GENERATION ########## # Paths @@ -101,6 +106,7 @@ global success_solves = 0.0 for i in 1:num_batches _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset_bayes(case_file_path, case_name; num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, + algorithm=algorithm, options=options, ) global success_solves += _success_solves end diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 63d6b5c..58e9c66 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -167,6 +167,8 @@ function generate_worst_case_dataset_bayes(data_dir, num_p=10, network_formulation=DCPPowerModel, optimizer = () -> POI.Optimizer(HiGHS.Optimizer()), + algorithm = _BayesOptAlg(), + options = _bayes_options(num_p), ) # save folder data_sim_dir = joinpath(data_dir, string(network_formulation)) @@ -214,10 +216,12 @@ function generate_worst_case_dataset_bayes(data_dir, # Set iterator function _set_iterator!(idx) - min_demands = zeros(num_loads * 2) - max_demands = original_load .+ ones(num_loads * 2) .* 0.1 * idx - max_total_volume = norm(max_demands, 2) ^ 2 - starting_point = original_load ./ 10 + _min_demands = original_load .- ones(num_loads * 2) .* 0.1 * idx + _max_demands = original_load .+ ones(num_loads * 2) .* 0.1 * idx + min_demands = min.(_min_demands, _max_demands) + max_demands = max.(_min_demands, _max_demands) + max_total_volume = (norm(max_demands, 2) + norm(min_demands, 2)) ^ 2 + starting_point = original_load return min_demands, max_demands, max_total_volume, starting_point end @@ -227,8 +231,8 @@ function generate_worst_case_dataset_bayes(data_dir, () -> nothing, # will be ignored _primal_builder!, _set_iterator!, - _BayesOptAlg(); - options = _bayes_options(num_p) + algorithm; + options = options ) # Solve all problems and record solutions diff --git a/src/worst_case_iter.jl b/src/worst_case_iter.jl index a7fc9f7..e3bf75c 100644 --- a/src/worst_case_iter.jl +++ b/src/worst_case_iter.jl @@ -80,9 +80,7 @@ function solve_and_record( model_non = Nonconvex.Model() set_objective!(model_non, storage_objective_function, flags = [:expensive]) addvar!(model_non, min_demands, max_demands) - add_ineq_constraint!(model_non, x -> sum(x .^ 2) - max_total_volume ^ 2) - - # Solution Method: Bayesian Optimization + # add_ineq_constraint!(model_non, x -> sum(x .^ 2) - max_total_volume ^ 2) # Optimize model_non: r_bayes = if !isnothing(problem_iterator.options) From d875f532aa88249e705260a051fded71a82b2c7b Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Wed, 20 Sep 2023 09:43:30 -0400 Subject: [PATCH 21/36] update code --- examples/flux/flux_forecaster_script.jl | 3 ++- examples/powermodels/generate_full_datasets_script.jl | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/flux/flux_forecaster_script.jl b/examples/flux/flux_forecaster_script.jl index 310f973..884e43a 100644 --- a/examples/flux/flux_forecaster_script.jl +++ b/examples/flux/flux_forecaster_script.jl @@ -5,13 +5,14 @@ using Arrow using Flux using DataFrames using PowerModels +using L2O # Paths path_dataset = joinpath(pwd(), "examples", "powermodels", "data") case_name = "pglib_opf_case300_ieee" filetype = ArrowFile network_formulation = SOCWRConicPowerModel -case_file_path = joinpath(path, case_name, string(network_formulation)) +case_file_path = joinpath(path_dataset, case_name, string(network_formulation)) # Load input and output data tables iter_files = readdir(joinpath(case_file_path)) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 436b822..83f598d 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -59,7 +59,7 @@ include(joinpath(path_powermodels, "pglib_datagen.jl")) # Parameters num_batches = 1 -num_p = 5 +num_p = 200 filetype = ArrowFile # Case name From dbef0b0d96999538facf1427831b51b419493267 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Fri, 22 Sep 2023 16:21:54 -0400 Subject: [PATCH 22/36] update names --- Project.toml | 6 ++--- .../generate_full_datasets_script.jl | 25 +------------------ examples/powermodels/pglib_datagen.jl | 12 ++++----- src/worst_case_iter.jl | 6 ++--- test/runtests.jl | 21 +++------------- test/worst_case.jl | 6 ++--- 6 files changed, 18 insertions(+), 58 deletions(-) diff --git a/Project.toml b/Project.toml index d7258dc..875ef7f 100644 --- a/Project.toml +++ b/Project.toml @@ -30,12 +30,10 @@ DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" -KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392" -NonconvexBayesian = "fb352abc-de7b-48de-9ebd-665b54b5d9b3" -NonconvexIpopt = "bf347577-a06d-49ad-a669-8c0e005493b8" +NonconvexNLopt = "b43a31b8-ff9b-442d-8e31-c163daa8ab75" PGLib = "07a8691f-3d11-4330-951b-3c50f98338be" PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "DelimitedFiles", "PGLib", "HiGHS", "PowerModels", "Flux", "DataFrames", "Clarabel", "Ipopt", "NonconvexBayesian", "NonconvexIpopt", "KernelFunctions", "AbstractGPs"] +test = ["Test", "DelimitedFiles", "PGLib", "HiGHS", "PowerModels", "Flux", "DataFrames", "Clarabel", "Ipopt", "NonconvexNLopt"] diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 83f598d..6917fcb 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -12,11 +12,6 @@ import JuMP.MOI as MOI import ParametricOptInterface as POI using Gurobi -using NonconvexIpopt # Nonconvex.@load Ipopt -using NonconvexBayesian # Nonconvex.@load BayesOpt -using AbstractGPs -using KernelFunctions - using NonconvexNLopt # using QuadraticToBinary @@ -33,23 +28,6 @@ cached = () -> MOI.Bridges.full_bridge_optimizer( POI_cached_optimizer() = POI.Optimizer(cached()) -_BayesOptAlg = () -> BayesOptAlg(IpoptAlg()) - -function _bayes_options(maxiter) - BayesOptOptions( - sub_options = IpoptOptions(max_iter = 20, print_level = 0), - # ninit=Int(floor(maxiter / 5)), - maxiter = maxiter * 100, ftol = 1e-4, ctol = 1e-5, initialize=true, postoptimize=false, - kernel= RationalKernel(α=2.27e8) ∘ ScaleTransform(0.01), - noise=0.001, - std_multiple=8.67e4, - fit_prior=false # not working with custom priors - ) -end - -algorithm = NLoptAlg(:LN_BOBYQA) -options = NLoptOptions(maxeval=30) - ########## DATASET GENERATION ########## # Paths @@ -104,9 +82,8 @@ mkpath(case_file_path) global success_solves = 0.0 for i in 1:num_batches - _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset_bayes(case_file_path, case_name; + _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset_Nonconvex(case_file_path, case_name; num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, - algorithm=algorithm, options=options, ) global success_solves += _success_solves end diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 58e9c66..c3cf026 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -161,14 +161,14 @@ function default_optimizer_factory() return () -> Ipopt.Optimizer() end -function generate_worst_case_dataset_bayes(data_dir, +function generate_worst_case_dataset_Nonconvex(data_dir, case_name; filetype=CSVFile, num_p=10, network_formulation=DCPPowerModel, optimizer = () -> POI.Optimizer(HiGHS.Optimizer()), - algorithm = _BayesOptAlg(), - options = _bayes_options(num_p), + algorithm = NLoptAlg(:LN_BOBYQA), + options = NLoptOptions(maxeval=10), ) # save folder data_sim_dir = joinpath(data_dir, string(network_formulation)) @@ -350,11 +350,11 @@ function test_generate_worst_case_dataset(path::AbstractString, case_name::Abstr end end -function test_generate_worst_case_dataset_bayes(path::AbstractString, case_name::AbstractString, num_p::Int) - @testset "Worst Case Dataset Generation pglib case" begin +function test_generate_worst_case_dataset_Nonconvex(path::AbstractString, case_name::AbstractString, num_p::Int) + @testset "WC Nonconvex Dataset Generation pglib case" begin network_formulation = DCPPowerModel # Improve dataset - success_solves, number_variables, number_parameters, batch_id = generate_worst_case_dataset_bayes( + success_solves, number_variables, number_parameters, batch_id = generate_worst_case_dataset_Nonconvex( path, case_name; num_p=num_p, network_formulation=network_formulation ) diff --git a/src/worst_case_iter.jl b/src/worst_case_iter.jl index e3bf75c..0579154 100644 --- a/src/worst_case_iter.jl +++ b/src/worst_case_iter.jl @@ -83,14 +83,14 @@ function solve_and_record( # add_ineq_constraint!(model_non, x -> sum(x .^ 2) - max_total_volume ^ 2) # Optimize model_non: - r_bayes = if !isnothing(problem_iterator.options) + r_Nonconvex = if !isnothing(problem_iterator.options) optimize(model_non, problem_iterator.optimizer, starting_point; options = problem_iterator.options) else optimize(model_non, problem_iterator.optimizer, starting_point) end - problem_iterator.ext[:best_solution] = r_bayes.minimizer - # best_profit = -r_bayes.minimum + problem_iterator.ext[:best_solution] = r_Nonconvex.minimizer + # best_profit = -r_Nonconvex.minimum return storage_objective_function.success_solves / storage_objective_function.fcalls end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index e647e26..015075c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,27 +9,12 @@ using Test using UUIDs using Ipopt -using NonconvexIpopt # Nonconvex.@load Ipopt -using NonconvexBayesian # Nonconvex.@load BayesOpt -using AbstractGPs -using KernelFunctions +using NonconvexNLopt const test_dir = dirname(@__FILE__) const examples_dir = joinpath(test_dir, "..", "examples") -_BayesOptAlg = () -> BayesOptAlg(IpoptAlg()) - -function _bayes_options(maxiter) - BayesOptOptions( - sub_options = IpoptOptions(max_iter = 20, print_level = 0), - # ninit=Int(floor(maxiter / 5)), - maxiter = maxiter, ftol = 1e-4, ctol = 1e-5, initialize=true, postoptimize=false, - kernel= RationalKernel(α=2.27e8) ∘ ScaleTransform(0.01), - noise=0.001, - std_multiple=8.67e4, - fit_prior=false # not working with custom priors - ) -end + include(joinpath(test_dir, "datasetgen.jl")) @@ -45,7 +30,7 @@ include(joinpath(examples_dir, "flux", "test_flux_forecaster.jl")) test_worst_case_problem_iterator(path) file_in, file_out = test_pglib_datasetgen(path, "pglib_opf_case5_pjm", 20) file_in, file_out = test_generate_worst_case_dataset(path, "pglib_opf_case5_pjm", 20) - file_in, file_out = test_generate_worst_case_dataset_bayes(path, "pglib_opf_case5_pjm", 20) + file_in, file_out = test_generate_worst_case_dataset_Nonconvex(path, "pglib_opf_case5_pjm", 20) test_flux_forecaster(file_in, file_out) end end diff --git a/test/worst_case.jl b/test/worst_case.jl index 5b42614..39b487a 100644 --- a/test/worst_case.jl +++ b/test/worst_case.jl @@ -72,7 +72,7 @@ function test_worst_case_problem_iterator(path::AbstractString, num_p=10) end end - @testset "Worst Case Bayes Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] + @testset "Worst Case NonConvex Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] function _primal_builder!(;recorder=nothing) model = JuMP.Model(() -> POI.Optimizer(HiGHS.Optimizer())) parameters = @variable(model, _p in POI.Parameter(1.0)) @@ -100,8 +100,8 @@ function test_worst_case_problem_iterator(path::AbstractString, num_p=10) () -> nothing, # will be ignored _primal_builder!, _set_iterator!, - _BayesOptAlg(); - options = _bayes_options(floor(Int, num_p / 5)) + NLoptAlg(:LN_BOBYQA); + options = NLoptOptions(maxeval=10) ) # file_names From 1a41033f08802d6c38406cd5b9ad119f3dc90787 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Fri, 22 Sep 2023 17:22:11 -0400 Subject: [PATCH 23/36] update code --- .../generate_full_datasets_script.jl | 37 ++++++++++++++----- examples/powermodels/pglib_datagen.jl | 7 +++- src/datasetgen.jl | 33 +++++++++++------ src/worst_case.jl | 27 +++++++++----- src/worst_case_iter.jl | 2 +- 5 files changed, 72 insertions(+), 34 deletions(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 6917fcb..11b41b1 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -47,17 +47,34 @@ case_file_path = joinpath(path, case_name) mkpath(case_file_path) # Generate dataset -# global success_solves = 0.0 -# for i in 1:num_batches -# _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; -# num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, -# load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01) -# ) -# global success_solves += _success_solves -# end -# success_solves /= num_batches +global success_solves = 0.0 +for i in 1:num_batches + _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; + num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, + load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01) + ) + global success_solves += _success_solves +end +success_solves /= num_batches + +@info "Success solves Normal: $(success_solves)" + +# Generate Line-search dataset + +early_stop_fn = (model, status, recorder) -> !status + +global success_solves = 0.0 +for i in 1:num_batches + _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; + num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, + load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01), + early_stop_fn=early_stop_fn + ) + global success_solves += _success_solves +end +success_solves /= num_batches -# @info "Success solves Normal: $(success_solves)" +@info "Success solves Normal: $(success_solves)" # Generate worst case dataset diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index c3cf026..ea8140e 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -107,6 +107,8 @@ function generate_dataset_pglib( load_sampler=load_sampler, network_formulation=DCPPowerModel, optimizer = () -> POI.Optimizer(HiGHS.Optimizer()), + filterfn=L2O.filter_fn, + early_stop_fn = (model, status, recorder) -> false, ) # save folder data_sim_dir = joinpath(data_dir, string(network_formulation)) @@ -134,14 +136,15 @@ function generate_dataset_pglib( # Build model and Recorder file = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_output_" * batch_id * "." * string(filetype)) - recorder = Recorder{filetype}(file) + recorder = Recorder{filetype}(file; filterfn=filterfn) pm_primal_builder!(model, p, network_data, network_formulation; recorder=recorder) # The problem iterator problem_iterator = ProblemIterator( Dict( p .=> load_sampler.(original_load, num_p), - ), + ); + early_stop=early_stop_fn ) save( diff --git a/src/datasetgen.jl b/src/datasetgen.jl index b72b810..246ad19 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -70,20 +70,21 @@ struct ProblemIterator{T<:Real} <: AbstractProblemIterator model::JuMP.Model ids::Vector{UUID} pairs::Dict{VariableRef,Vector{T}} + early_stop::Function function ProblemIterator( - ids::Vector{UUID}, pairs::Dict{VariableRef,Vector{T}} + ids::Vector{UUID}, pairs::Dict{VariableRef,Vector{T}}, early_stop::Function=(args...) -> false ) where {T<:Real} model = JuMP.owner_model(first(keys(pairs))) for (p, val) in pairs @assert length(ids) == length(val) end - return new{T}(model, ids, pairs) + return new{T}(model, ids, pairs, early_stop) end end -function ProblemIterator(pairs::Dict{VariableRef,Vector{T}}) where {T<:Real} +function ProblemIterator(pairs::Dict{VariableRef,Vector{T}}; early_stop::Function=(args...) -> false) where {T<:Real} ids = [uuid1() for _ in 1:length(first(values(pairs)))] - return ProblemIterator(ids, pairs) + return ProblemIterator(ids, pairs, early_stop) end """ @@ -137,11 +138,13 @@ function solve_and_record( model = problem_iterator.model update_model!(model, problem_iterator.pairs, idx) optimize!(model) - if recorder.filterfn(model) + status = recorder.filterfn(model) + early_stop_bool = problem_iterator.early_stop(model, status, recorder) + if status record(recorder, problem_iterator.ids[idx]) - return 1 + return 1, early_stop_bool end - return 0 + return 0, early_stop_bool end """ @@ -152,11 +155,17 @@ Solve a batch of optimization problems and record the solutions. function solve_batch( problem_iterator::AbstractProblemIterator, recorder ) - successfull_solves = - sum( - solve_and_record(problem_iterator, recorder, idx) for - idx in 1:length(problem_iterator.ids) - ) / length(problem_iterator.ids) + successfull_solves = 0.0 + for idx in 1:length(problem_iterator.ids) + _success_bool, early_stop_bool = solve_and_record(problem_iterator, recorder, idx) + if _success_bool == 1 + successfull_solves += 1 + end + if early_stop_bool + break + end + end + successfull_solves = successfull_solves / length(problem_iterator.ids) @info "Recorded $(successfull_solves * 100) % of $(length(problem_iterator.ids)) problems" return successfull_solves diff --git a/src/worst_case.jl b/src/worst_case.jl index 9f97b07..698d09f 100644 --- a/src/worst_case.jl +++ b/src/worst_case.jl @@ -29,6 +29,7 @@ struct WorstCaseProblemIterator{F} <: AbstractProblemIterator hook::Union{Nothing, Function} options::Any ext::Dict + early_stop::Function function WorstCaseProblemIterator( ids::Vector{UUID}, parameters::Function, @@ -37,9 +38,10 @@ struct WorstCaseProblemIterator{F} <: AbstractProblemIterator optimizer::F; hook::Union{Nothing, Function}=nothing, options::Any=nothing, - ext::Dict=Dict() + ext::Dict=Dict(), + early_stop::Function=(args...) -> false ) where {F} - new{F}(ids, parameters, primal_builder!, set_iterator!, optimizer, hook, options, ext) + new{F}(ids, parameters, primal_builder!, set_iterator!, optimizer, hook, options, ext, early_stop) end end @@ -102,12 +104,17 @@ function solve_and_record( JuMP.optimize!(jump_dual_model) # Save input + status = recorder.filterfn(jump_dual_model) + early_stop_bool = problem_iterator.early_stop(jump_dual_model, status, recorder) + if early_stop_bool + return 0, early_stop_bool + end if recorder.filterfn(jump_dual_model) recorder.primal_variables = load_var_dual recorder.dual_variables = [] record(recorder, problem_iterator.ids[idx]; input=true) else - return 0 + return 0, early_stop_bool end optimal_loads = value.(load_var_dual) @@ -118,14 +125,16 @@ function solve_and_record( problem_iterator.primal_builder!(model, optimal_loads; recorder=recorder) JuMP.optimize!(model) + termination_status = recorder.filterfn(model) + early_stop_bool = problem_iterator.early_stop(model, termination_status, recorder) + # Check if method was effective optimal_final_cost = JuMP.objective_value(model) - termination_status = JuMP.termination_status(model) solution_primal_status = JuMP.primal_status(model) solution_dual_status = JuMP.dual_status(model) - termination_status == MOI.INFEASIBLE && @error("Optimal solution not found") - solution_primal_status != MOI.FEASIBLE_POINT && @error("Primal solution not found") - solution_dual_status != MOI.FEASIBLE_POINT && @error("Dual solution not found") + termination_status == MOI.INFEASIBLE && @warn("Optimal solution not found") + solution_primal_status != MOI.FEASIBLE_POINT && @warn("Primal solution not found") + solution_dual_status != MOI.FEASIBLE_POINT && @warn("Dual solution not found") if !isapprox(optimal_final_cost, optimal_dual_cost; rtol=1e-4) rtol = abs(optimal_final_cost - optimal_dual_cost) / optimal_final_cost * 100 @@ -134,7 +143,7 @@ function solve_and_record( if recorder.filterfn(model) record(recorder, problem_iterator.ids[idx]) - return 1 + return 1, early_stop_bool end - return 0 + return 0, early_stop_bool end \ No newline at end of file diff --git a/src/worst_case_iter.jl b/src/worst_case_iter.jl index 0579154..6a590fa 100644 --- a/src/worst_case_iter.jl +++ b/src/worst_case_iter.jl @@ -92,5 +92,5 @@ function solve_and_record( problem_iterator.ext[:best_solution] = r_Nonconvex.minimizer # best_profit = -r_Nonconvex.minimum - return storage_objective_function.success_solves / storage_objective_function.fcalls + return storage_objective_function.success_solves / storage_objective_function.fcalls, false end \ No newline at end of file From 585371179b33c4b3cce67bbd989ec6c28471c4c9 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Fri, 22 Sep 2023 19:02:40 -0400 Subject: [PATCH 24/36] update --- .../generate_full_datasets_script.jl | 74 ++++++++++------- examples/powermodels/pglib_datagen.jl | 24 +++--- src/arrowrecorder.jl | 5 +- src/csvrecorder.jl | 3 +- src/datasetgen.jl | 3 +- test/datasetgen.jl | 82 +++++++++++-------- 6 files changed, 117 insertions(+), 74 deletions(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 11b41b1..b8aa8b5 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -36,45 +36,63 @@ path = joinpath(path_powermodels, "data") include(joinpath(path_powermodels, "pglib_datagen.jl")) # Parameters -num_batches = 1 -num_p = 200 -filetype = ArrowFile +filetype = CSVFile # CSVFile # ArrowFile # Case name -case_name = "pglib_opf_case300_ieee" # "pglib_opf_case300_ieee" +case_name = "pglib_opf_case5_pjm" # pglib_opf_case300_ieee # pglib_opf_case5_pjm network_formulation = SOCWRConicPowerModel # SOCWRConicPowerModel # DCPPowerModel case_file_path = joinpath(path, case_name) mkpath(case_file_path) # Generate dataset -global success_solves = 0.0 -for i in 1:num_batches - _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; - num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, - load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01) - ) - global success_solves += _success_solves -end -success_solves /= num_batches +# num_batches = 1 +# num_p = 200 +# global success_solves = 0.0 +# for i in 1:num_batches +# _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; +# num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, +# load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01) +# ) +# global success_solves += _success_solves +# end +# success_solves /= num_batches -@info "Success solves Normal: $(success_solves)" +# @info "Success solves Normal: $(success_solves)" # Generate Line-search dataset +matpower_case_name = case_name * ".m" +network_data = make_basic_network(pglib(matpower_case_name)) + early_stop_fn = (model, status, recorder) -> !status +step_multiplier = 1.01 +num_loads = length(network_data["load"]) +num_batches = num_loads * 2 + 1 +num_p = 10 + +function line_sampler(_o, n, idx, num_inputs, ibatc) + @info "ibatc: $(ibatc)" n idx num_inputs + if (idx == ibatc) || (idx == num_inputs + 1) + return [_o * step_multiplier ^ j for j in 1:n] + else + return ones(n) + end +end global success_solves = 0.0 -for i in 1:num_batches - _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; +global batch_id = string(uuid1()) +for ibatc in 1:num_batches + _success_solves, number_variables, number_loads, b_id = generate_dataset_pglib(case_file_path, case_name; num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, - load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01), - early_stop_fn=early_stop_fn + load_sampler= (_o, n, idx, num_inputs) -> line_sampler(_o, n, idx, num_inputs, ibatc), + early_stop_fn=early_stop_fn, + batch_id=batch_id, ) global success_solves += _success_solves end success_solves /= num_batches -@info "Success solves Normal: $(success_solves)" +@info "Success solves: $(success_solves * 100) % of $(num_batches * num_p)" # Generate worst case dataset @@ -97,13 +115,13 @@ success_solves /= num_batches # @info "Success solves Worst Case: $(success_solves) of $(num_batches * num_p)" -global success_solves = 0.0 -for i in 1:num_batches - _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset_Nonconvex(case_file_path, case_name; - num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, - ) - global success_solves += _success_solves -end -success_solves /= num_batches +# global success_solves = 0.0 +# for i in 1:num_batches +# _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset_Nonconvex(case_file_path, case_name; +# num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, +# ) +# global success_solves += _success_solves +# end +# success_solves /= num_batches -@info "Success solves Worst Case: $(success_solves * 100) of $(num_batches * num_p)" +# @info "Success solves Worst Case: $(success_solves * 100) of $(num_batches * num_p)" diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index ea8140e..dd24028 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -34,11 +34,13 @@ Load sampling """ function load_sampler( original_load::T, - num_p::Int; + num_p::F, + ::F, + ::F; max_multiplier::T=2.5, min_multiplier::T=0.0, step_multiplier::T=0.1, -) where {T<:Real} +) where {T<:Real, F<:Integer} # Load sampling load_samples = original_load * rand(min_multiplier:step_multiplier:max_multiplier, num_p) @@ -109,7 +111,10 @@ function generate_dataset_pglib( optimizer = () -> POI.Optimizer(HiGHS.Optimizer()), filterfn=L2O.filter_fn, early_stop_fn = (model, status, recorder) -> false, + batch_id = string(uuid1()) ) + @info "Batch ID: $batch_id" + # save folder data_sim_dir = joinpath(data_dir, string(network_formulation)) mkpath(data_sim_dir) @@ -124,15 +129,12 @@ function generate_dataset_pglib( # Save original load value and Link POI num_loads = length(network_data["load"]) + num_inputs = num_loads * 2 original_load = vcat( [l["pd"] for l in values(network_data["load"])], [l["qd"] for l in values(network_data["load"])], ) - p = load_parameter_factory(model, 1:(num_loads * 2); load_set=POI.Parameter.(original_load)) - - # Define batch id - batch_id = string(uuid1()) - @info "Batch ID: $batch_id" + p = load_parameter_factory(model, 1:(num_inputs); load_set=POI.Parameter.(original_load)) # Build model and Recorder file = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_output_" * batch_id * "." * string(filetype)) @@ -140,10 +142,12 @@ function generate_dataset_pglib( pm_primal_builder!(model, p, network_data, network_formulation; recorder=recorder) # The problem iterator + pairs = Dict{VariableRef, Vector{Float64}}() + for i in 1:num_inputs + pairs[p[i]] = load_sampler(original_load[i], num_p, i, num_inputs) + end problem_iterator = ProblemIterator( - Dict( - p .=> load_sampler.(original_load, num_p), - ); + pairs; early_stop=early_stop_fn ) diff --git a/src/arrowrecorder.jl b/src/arrowrecorder.jl index 23061ea..a7cefbd 100644 --- a/src/arrowrecorder.jl +++ b/src/arrowrecorder.jl @@ -40,5 +40,8 @@ function record(recorder::Recorder{ArrowFile}, id::UUID; input=false) end function save(table::NamedTuple, filename::String, ::Type{ArrowFile}) - return Arrow.write(filename, table) + Arrow.append( + filename, + table, + ) end diff --git a/src/csvrecorder.jl b/src/csvrecorder.jl index bb8fcef..4006e0b 100644 --- a/src/csvrecorder.jl +++ b/src/csvrecorder.jl @@ -52,5 +52,6 @@ function record(recorder::Recorder{CSVFile}, id::UUID; input=false) end function save(table::NamedTuple, filename::String, ::Type{CSVFile}) - return CSV.write(filename, table) + CSV.write(filename, table; append=isfile(filename)) + return nothing end diff --git a/src/datasetgen.jl b/src/datasetgen.jl index 246ad19..68700eb 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -95,7 +95,7 @@ Save optimization problem instances to a file. function save( problem_iterator::AbstractProblemIterator, filename::String, file_type::Type{T} ) where {T<:FileType} - return save( + save( (; id=problem_iterator.ids, zip( @@ -105,6 +105,7 @@ function save( filename, file_type, ) + return nothing end """ diff --git a/test/datasetgen.jl b/test/datasetgen.jl index e62070c..fd87be7 100644 --- a/test/datasetgen.jl +++ b/test/datasetgen.jl @@ -15,50 +15,66 @@ function test_problem_iterator(path::AbstractString) # The problem iterator num_p = 10 - @test_throws AssertionError ProblemIterator( - [uuid1() for _ in 1:num_p], Dict(p => collect(1.0:3.0)) - ) - @test_throws MethodError ProblemIterator( - collect(1.0:3.0), Dict(p => collect(1.0:3.0)) - ) + @testset "Problem Iterator Builder" begin + @test_throws AssertionError ProblemIterator( + [uuid1() for _ in 1:num_p], Dict(p => collect(1.0:3.0)) + ) + @test_throws MethodError ProblemIterator( + collect(1.0:3.0), Dict(p => collect(1.0:3.0)) + ) + problem_iterator = ProblemIterator(Dict(p => collect(1.0:num_p))) + file_input = joinpath(path, "test_input.$(string(filetype))") # file path + save(problem_iterator, file_input, filetype) + @test isfile(file_input) + end problem_iterator = ProblemIterator(Dict(p => collect(1.0:num_p))) file_input = joinpath(path, "test_input.$(string(filetype))") # file path - save(problem_iterator, file_input, filetype) - @test isfile(file_input) # The recorder file_output = joinpath(path, "test_output.$(string(filetype))") # file path - @test Recorder{filetype}(file_output; primal_variables=[x]) isa Recorder{filetype} - @test Recorder{filetype}(file_output; dual_variables=[cons]) isa Recorder{filetype} + @testset "Recorder Builder" begin + @test Recorder{filetype}(file_output; primal_variables=[x]) isa Recorder{filetype} + @test Recorder{filetype}(file_output; dual_variables=[cons]) isa Recorder{filetype} + end recorder = Recorder{filetype}( file_output; primal_variables=[x], dual_variables=[cons] ) # Solve all problems and record solutions - solve_batch(problem_iterator, recorder) + @testset "early_stop" begin + problem_iterator = ProblemIterator(Dict(p => collect(1.0:num_p)); + early_stop=(args...) -> true + ) + successfull_solves = solve_batch(problem_iterator, recorder) + @test num_p * successfull_solves == 1 + end + + @testset "solve_batch" begin + successfull_solves = solve_batch(problem_iterator, recorder) - # Check if file exists and has the correct number of rows and columns - @test isfile(file_output) - if filetype == CSVFile - # test input file - @test length(readdlm(file_input, ',')[:, 1]) == num_p + 1 - @test length(readdlm(file_input, ',')[1, :]) == 2 - rm(file_input) - # test output file - @test length(readdlm(file_output, ',')[:, 1]) == num_p + 1 - @test length(readdlm(file_output, ',')[1, :]) == 4 - rm(file_output) - else - # test input file - df = Arrow.Table(file_input) - @test length(df) == 2 - @test length(df[1]) == num_p - rm(file_input) - # test output file - df = Arrow.Table(file_output) - @test length(df) == 4 - @test length(df[1]) == num_p - rm(file_output) + # Check if file exists and has the correct number of rows and columns + @test isfile(file_output) + if filetype == CSVFile + # test input file + @test length(readdlm(file_input, ',')[:, 1]) == num_p + 1 # 1 from header + @test length(readdlm(file_input, ',')[1, :]) == 2 # 2 parameter + rm(file_input) + # test output file + @test length(readdlm(file_output, ',')[:, 1]) == num_p * successfull_solves + 2 # 1 from early_stop and 1 from header + @test length(readdlm(file_output, ',')[1, :]) == 4 + rm(file_output) + else + # test input file + df = Arrow.Table(file_input) + @test length(df) == 2 # 2 parameter + @test length(df[1]) == num_p + rm(file_input) + # test output file + df = Arrow.Table(file_output) + @test length(df) == 4 + @test length(df[1]) == num_p * successfull_solves + 1 # 1 from early_stop + rm(file_output) + end end end end From f68c4932a254b993fecf7b16784456b065feab5c Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Fri, 22 Sep 2023 19:13:45 -0400 Subject: [PATCH 25/36] update --- examples/powermodels/generate_full_datasets_script.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index b8aa8b5..d215d26 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -72,7 +72,8 @@ num_p = 10 function line_sampler(_o, n, idx, num_inputs, ibatc) @info "ibatc: $(ibatc)" n idx num_inputs - if (idx == ibatc) || (idx == num_inputs + 1) + if (idx == ibatc) || (ibatc == num_inputs + 1) + @info "Minha vez" [_o * step_multiplier ^ j for j in 1:n] return [_o * step_multiplier ^ j for j in 1:n] else return ones(n) From f3b34f1204c40b0ac725d5ca5f958674c904bd9c Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Fri, 22 Sep 2023 21:45:25 -0400 Subject: [PATCH 26/36] update code --- examples/flux/flux_forecaster_script.jl | 5 ++--- .../generate_full_datasets_script.jl | 10 ++++------ examples/powermodels/pglib_datagen.jl | 17 +++++++++-------- src/arrowrecorder.jl | 5 +++-- src/csvrecorder.jl | 11 +++++++++-- src/datasetgen.jl | 17 ++++++++++------- 6 files changed, 37 insertions(+), 28 deletions(-) diff --git a/examples/flux/flux_forecaster_script.jl b/examples/flux/flux_forecaster_script.jl index 884e43a..8e2cd14 100644 --- a/examples/flux/flux_forecaster_script.jl +++ b/examples/flux/flux_forecaster_script.jl @@ -8,10 +8,9 @@ using PowerModels using L2O # Paths -path_dataset = joinpath(pwd(), "examples", "powermodels", "data") -case_name = "pglib_opf_case300_ieee" +case_name = "pglib_opf_case5_pjm" # pglib_opf_case300_ieee # pglib_opf_case5_pjm +network_formulation = SOCWRConicPowerModel # SOCWRConicPowerModel # DCPPowerModel filetype = ArrowFile -network_formulation = SOCWRConicPowerModel case_file_path = joinpath(path_dataset, case_name, string(network_formulation)) # Load input and output data tables diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index d215d26..607e7e6 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -36,7 +36,7 @@ path = joinpath(path_powermodels, "data") include(joinpath(path_powermodels, "pglib_datagen.jl")) # Parameters -filetype = CSVFile # CSVFile # ArrowFile +filetype = ArrowFile # CSVFile # ArrowFile # Case name case_name = "pglib_opf_case5_pjm" # pglib_opf_case300_ieee # pglib_opf_case5_pjm @@ -51,7 +51,7 @@ mkpath(case_file_path) # for i in 1:num_batches # _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; # num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, -# load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01) +# internal_load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01) # ) # global success_solves += _success_solves # end @@ -71,10 +71,8 @@ num_batches = num_loads * 2 + 1 num_p = 10 function line_sampler(_o, n, idx, num_inputs, ibatc) - @info "ibatc: $(ibatc)" n idx num_inputs if (idx == ibatc) || (ibatc == num_inputs + 1) - @info "Minha vez" [_o * step_multiplier ^ j for j in 1:n] - return [_o * step_multiplier ^ j for j in 1:n] + return [_o * step_multiplier ^ (j-1) for j in 1:n] else return ones(n) end @@ -85,7 +83,7 @@ global batch_id = string(uuid1()) for ibatc in 1:num_batches _success_solves, number_variables, number_loads, b_id = generate_dataset_pglib(case_file_path, case_name; num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, - load_sampler= (_o, n, idx, num_inputs) -> line_sampler(_o, n, idx, num_inputs, ibatc), + internal_load_sampler= (_o, n, idx, num_inputs) -> line_sampler(_o, n, idx, num_inputs, ibatc), early_stop_fn=early_stop_fn, batch_id=batch_id, ) diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index dd24028..10c096e 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -50,17 +50,18 @@ end function load_parameter_factory(model, indices; load_set=nothing) if isnothing(load_set) return @variable( - model, _p[i in indices] + model, _p[i=indices] ) end return @variable( - model, _p[i in indices] in load_set + model, _p[i=indices] in load_set ) end function pm_primal_builder!(model, parameters, network_data, network_formulation; recorder=nothing) num_loads = length(network_data["load"]) - for (i, l) in enumerate(values(network_data["load"])) + for (str_i, l) in network_data["load"] + i = parse(Int, str_i) l["pd"] = parameters[i] l["qd"] = parameters[num_loads+i] end @@ -106,7 +107,7 @@ function generate_dataset_pglib( case_name; filetype=CSVFile, num_p=10, - load_sampler=load_sampler, + internal_load_sampler=load_sampler, network_formulation=DCPPowerModel, optimizer = () -> POI.Optimizer(HiGHS.Optimizer()), filterfn=L2O.filter_fn, @@ -131,10 +132,10 @@ function generate_dataset_pglib( num_loads = length(network_data["load"]) num_inputs = num_loads * 2 original_load = vcat( - [l["pd"] for l in values(network_data["load"])], - [l["qd"] for l in values(network_data["load"])], + [network_data["load"]["$l"]["pd"] for l=1:num_loads], + [network_data["load"]["$l"]["qd"] for l=1:num_loads], ) - p = load_parameter_factory(model, 1:(num_inputs); load_set=POI.Parameter.(original_load)) + p = load_parameter_factory(model, 1:num_inputs; load_set=POI.Parameter.(original_load)) # Build model and Recorder file = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_output_" * batch_id * "." * string(filetype)) @@ -144,7 +145,7 @@ function generate_dataset_pglib( # The problem iterator pairs = Dict{VariableRef, Vector{Float64}}() for i in 1:num_inputs - pairs[p[i]] = load_sampler(original_load[i], num_p, i, num_inputs) + pairs[p[i]] = internal_load_sampler(original_load[i], num_p, i, num_inputs) end problem_iterator = ProblemIterator( pairs; diff --git a/src/arrowrecorder.jl b/src/arrowrecorder.jl index a7cefbd..8374d17 100644 --- a/src/arrowrecorder.jl +++ b/src/arrowrecorder.jl @@ -39,9 +39,10 @@ function record(recorder::Recorder{ArrowFile}, id::UUID; input=false) ) end -function save(table::NamedTuple, filename::String, ::Type{ArrowFile}) +function save(table::NamedTuple, filename::String, ::Type{ArrowFile}; kwargs...) Arrow.append( filename, - table, + table; + kwargs... ) end diff --git a/src/csvrecorder.jl b/src/csvrecorder.jl index 4006e0b..35b6374 100644 --- a/src/csvrecorder.jl +++ b/src/csvrecorder.jl @@ -51,7 +51,14 @@ function record(recorder::Recorder{CSVFile}, id::UUID; input=false) end end -function save(table::NamedTuple, filename::String, ::Type{CSVFile}) - CSV.write(filename, table; append=isfile(filename)) +function save(table::NamedTuple, filename::String, ::Type{CSVFile}; kwargs...) + isappend = isfile(filename) + mode = isappend ? "append" : "write" + @info "Saving CSV file to $filename - Mode: $mode" + if isappend + @warn "Appending to existing CSV file is NOT RECOMMENDED" + # CSV does not take into account the column names when appending + end + CSV.write(filename, table; append=isappend) return nothing end diff --git a/src/datasetgen.jl b/src/datasetgen.jl index 68700eb..2e08058 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -95,15 +95,18 @@ Save optimization problem instances to a file. function save( problem_iterator::AbstractProblemIterator, filename::String, file_type::Type{T} ) where {T<:FileType} + kys = keys(problem_iterator.pairs) + df = (; + id=problem_iterator.ids, + ) + for ky in kys + df = merge(df, (; Symbol(ky) => problem_iterator.pairs[ky])) + end save( - (; - id=problem_iterator.ids, - zip( - Symbol.(name.(keys(problem_iterator.pairs))), values(problem_iterator.pairs) - )..., - ), + df, filename, - file_type, + file_type; + dictencode=true, ) return nothing end From 760247e391fb9932bf390a4f703480153a2885f5 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Fri, 22 Sep 2023 22:08:55 -0400 Subject: [PATCH 27/36] update --- examples/flux/flux_forecaster_script.jl | 1 + src/datasetgen.jl | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/flux/flux_forecaster_script.jl b/examples/flux/flux_forecaster_script.jl index 8e2cd14..1565596 100644 --- a/examples/flux/flux_forecaster_script.jl +++ b/examples/flux/flux_forecaster_script.jl @@ -11,6 +11,7 @@ using L2O case_name = "pglib_opf_case5_pjm" # pglib_opf_case300_ieee # pglib_opf_case5_pjm network_formulation = SOCWRConicPowerModel # SOCWRConicPowerModel # DCPPowerModel filetype = ArrowFile +path_dataset = joinpath(pwd(), "examples", "powermodels", "data") case_file_path = joinpath(path_dataset, case_name, string(network_formulation)) # Load input and output data tables diff --git a/src/datasetgen.jl b/src/datasetgen.jl index 2e08058..02fb8ca 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -95,7 +95,7 @@ Save optimization problem instances to a file. function save( problem_iterator::AbstractProblemIterator, filename::String, file_type::Type{T} ) where {T<:FileType} - kys = keys(problem_iterator.pairs) + kys = sort(collect(keys(problem_iterator.pairs)); by=(v) -> index(v).value) df = (; id=problem_iterator.ids, ) From 9b8adf70bdcb00b57f8a32e08acad9cf6f164240 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Fri, 22 Sep 2023 22:16:27 -0400 Subject: [PATCH 28/36] update --- examples/powermodels/generate_full_datasets_script.jl | 4 ++-- src/csvrecorder.jl | 4 ---- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 607e7e6..53111dd 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -36,10 +36,10 @@ path = joinpath(path_powermodels, "data") include(joinpath(path_powermodels, "pglib_datagen.jl")) # Parameters -filetype = ArrowFile # CSVFile # ArrowFile +filetype = CSVFile # CSVFile # ArrowFile # Case name -case_name = "pglib_opf_case5_pjm" # pglib_opf_case300_ieee # pglib_opf_case5_pjm +case_name = "pglib_opf_case300_ieee" # pglib_opf_case300_ieee # pglib_opf_case5_pjm network_formulation = SOCWRConicPowerModel # SOCWRConicPowerModel # DCPPowerModel case_file_path = joinpath(path, case_name) mkpath(case_file_path) diff --git a/src/csvrecorder.jl b/src/csvrecorder.jl index 35b6374..c51385d 100644 --- a/src/csvrecorder.jl +++ b/src/csvrecorder.jl @@ -55,10 +55,6 @@ function save(table::NamedTuple, filename::String, ::Type{CSVFile}; kwargs...) isappend = isfile(filename) mode = isappend ? "append" : "write" @info "Saving CSV file to $filename - Mode: $mode" - if isappend - @warn "Appending to existing CSV file is NOT RECOMMENDED" - # CSV does not take into account the column names when appending - end CSV.write(filename, table; append=isappend) return nothing end From f920399e7fa79821d639f63c3a2790046bb73731 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Fri, 22 Sep 2023 22:52:04 -0400 Subject: [PATCH 29/36] update code --- .../generate_full_datasets_script.jl | 8 ++++---- src/datasetgen.jl | 19 ++++++++++++++++--- 2 files changed, 20 insertions(+), 7 deletions(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 53111dd..02bbae7 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -67,14 +67,14 @@ network_data = make_basic_network(pglib(matpower_case_name)) early_stop_fn = (model, status, recorder) -> !status step_multiplier = 1.01 num_loads = length(network_data["load"]) -num_batches = num_loads * 2 + 1 -num_p = 10 +num_batches = num_loads + 1 +num_p = 2 function line_sampler(_o, n, idx, num_inputs, ibatc) - if (idx == ibatc) || (ibatc == num_inputs + 1) + if (idx == ibatc) || (idx - num_loads == ibatc) || (ibatc == num_inputs + 1) return [_o * step_multiplier ^ (j-1) for j in 1:n] else - return ones(n) + return ones(n) * _o end end diff --git a/src/datasetgen.jl b/src/datasetgen.jl index 02fb8ca..c898d41 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -6,9 +6,22 @@ end filename(recorder_file::RecorderFile) = recorder_file.filename -termination_status_filter(status) = status == MOI.OPTIMAL || status == MOI.SLOW_PROGRESS || status == MOI.LOCALLY_SOLVED || status == MOI.ITERATION_LIMIT -primal_status_filter(status) = status == MOI.FEASIBLE_POINT -dual_status_filter(status) = status == MOI.FEASIBLE_POINT +const ACCEPTED_TERMINATION_STATUSES = [ + MOI.OPTIMAL, + MOI.SLOW_PROGRESS, + MOI.LOCALLY_SOLVED, + MOI.ITERATION_LIMIT, + MOI.ALMOST_OPTIMAL, +] + +DECISION_STATUS = [ + MOI.FEASIBLE_POINT, + MOI.NEARLY_FEASIBLE_POINT +] + +termination_status_filter(status) = in(status, ACCEPTED_TERMINATION_STATUSES) +primal_status_filter(status) = in(status, DECISION_STATUS) +dual_status_filter(status) = in(status, DECISION_STATUS) filter_fn(model) = termination_status_filter(termination_status(model)) && primal_status_filter(primal_status(model)) && dual_status_filter(dual_status(model)) From 8f1163071f09954eccdd3c7d2cd4a24af2fb42ec Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Fri, 22 Sep 2023 23:00:59 -0400 Subject: [PATCH 30/36] update tests --- examples/flux/flux_forecaster_script.jl | 2 +- examples/powermodels/generate_full_datasets_script.jl | 2 +- test/datasetgen.jl | 5 ++++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/examples/flux/flux_forecaster_script.jl b/examples/flux/flux_forecaster_script.jl index 1565596..659a352 100644 --- a/examples/flux/flux_forecaster_script.jl +++ b/examples/flux/flux_forecaster_script.jl @@ -8,7 +8,7 @@ using PowerModels using L2O # Paths -case_name = "pglib_opf_case5_pjm" # pglib_opf_case300_ieee # pglib_opf_case5_pjm +case_name = "pglib_opf_case300_ieee" # pglib_opf_case300_ieee # pglib_opf_case5_pjm network_formulation = SOCWRConicPowerModel # SOCWRConicPowerModel # DCPPowerModel filetype = ArrowFile path_dataset = joinpath(pwd(), "examples", "powermodels", "data") diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 02bbae7..0451b63 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -36,7 +36,7 @@ path = joinpath(path_powermodels, "data") include(joinpath(path_powermodels, "pglib_datagen.jl")) # Parameters -filetype = CSVFile # CSVFile # ArrowFile +filetype = ArrowFile # CSVFile # ArrowFile # Case name case_name = "pglib_opf_case300_ieee" # pglib_opf_case300_ieee # pglib_opf_case5_pjm diff --git a/test/datasetgen.jl b/test/datasetgen.jl index fd87be7..dbe067f 100644 --- a/test/datasetgen.jl +++ b/test/datasetgen.jl @@ -42,10 +42,13 @@ function test_problem_iterator(path::AbstractString) # Solve all problems and record solutions @testset "early_stop" begin + recorder_dual = Recorder{filetype}( + file_output; dual_variables=[cons] + ) problem_iterator = ProblemIterator(Dict(p => collect(1.0:num_p)); early_stop=(args...) -> true ) - successfull_solves = solve_batch(problem_iterator, recorder) + successfull_solves = solve_batch(problem_iterator, recorder_dual) @test num_p * successfull_solves == 1 end From 9539e2c7a0d74d15ab4e9d1285a8e03bffeac069 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Fri, 22 Sep 2023 23:29:23 -0400 Subject: [PATCH 31/36] update tests --- test/datasetgen.jl | 7 ++++--- test/worst_case.jl | 10 ++++++---- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/test/datasetgen.jl b/test/datasetgen.jl index dbe067f..15d6959 100644 --- a/test/datasetgen.jl +++ b/test/datasetgen.jl @@ -42,8 +42,9 @@ function test_problem_iterator(path::AbstractString) # Solve all problems and record solutions @testset "early_stop" begin + file_dual_output = joinpath(path, "test_$(string(uuid1()))_output.$(string(filetype))") # file path recorder_dual = Recorder{filetype}( - file_output; dual_variables=[cons] + file_dual_output; dual_variables=[cons] ) problem_iterator = ProblemIterator(Dict(p => collect(1.0:num_p)); early_stop=(args...) -> true @@ -63,7 +64,7 @@ function test_problem_iterator(path::AbstractString) @test length(readdlm(file_input, ',')[1, :]) == 2 # 2 parameter rm(file_input) # test output file - @test length(readdlm(file_output, ',')[:, 1]) == num_p * successfull_solves + 2 # 1 from early_stop and 1 from header + @test length(readdlm(file_output, ',')[:, 1]) == num_p * successfull_solves + 1 # 1 from header @test length(readdlm(file_output, ',')[1, :]) == 4 rm(file_output) else @@ -75,7 +76,7 @@ function test_problem_iterator(path::AbstractString) # test output file df = Arrow.Table(file_output) @test length(df) == 4 - @test length(df[1]) == num_p * successfull_solves + 1 # 1 from early_stop + @test length(df[1]) == num_p * successfull_solves rm(file_output) end end diff --git a/test/worst_case.jl b/test/worst_case.jl index 39b487a..31b1ada 100644 --- a/test/worst_case.jl +++ b/test/worst_case.jl @@ -34,8 +34,9 @@ function test_worst_case_problem_iterator(path::AbstractString, num_p=10) ) # file_names - file_input = joinpath(path, "test_input.$(string(filetype))") - file_output = joinpath(path, "test_output.$(string(filetype))") + batch_id = string(uuid1()) + file_input = joinpath(path, "test_$(batch_id)_input.$(string(filetype))") + file_output = joinpath(path, "test_$(batch_id)_output.$(string(filetype))") # The recorder recorder = Recorder{filetype}( @@ -105,8 +106,9 @@ function test_worst_case_problem_iterator(path::AbstractString, num_p=10) ) # file_names - file_input = joinpath(path, "test_input.$(string(filetype))") - file_output = joinpath(path, "test_output.$(string(filetype))") + batch_id = string(uuid1()) + file_input = joinpath(path, "test_$(batch_id)_input.$(string(filetype))") + file_output = joinpath(path, "test_$(batch_id)_output.$(string(filetype))") # The recorder recorder = Recorder{filetype}( From 10af4f8b7b951583776a68a03734190fc88a3738 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Fri, 22 Sep 2023 23:40:43 -0400 Subject: [PATCH 32/36] rm legacy --- test/worst_case_old.jl | 108 ----------------------------------------- 1 file changed, 108 deletions(-) delete mode 100644 test/worst_case_old.jl diff --git a/test/worst_case_old.jl b/test/worst_case_old.jl deleted file mode 100644 index 42af394..0000000 --- a/test/worst_case_old.jl +++ /dev/null @@ -1,108 +0,0 @@ -using PowerModels -using Ipopt -using PGLib -using JuMP -using Dualization -import JuMP.MOI as MOI - -# Load data -case_name = "pglib_opf_case5_pjm" -matpower_case_name = case_name * ".m" -network_data_original = pglib(matpower_case_name) -network_data = deepcopy(network_data_original) - -# Case Parameters -network_formulation=DCPPowerModel -original_load = [l["pd"] for l in values(network_data["load"])] -max_total_load = sum(original_load) * 1.1 - -# Create primal model -solver = Ipopt.Optimizer -model = JuMP.Model(solver) -load_var = @variable(model, load_var[i=1:length(original_load)]) -for (i, l) in enumerate(values(network_data["load"])) - l["pd"] = load_var[i] -end -load_moi_idx = MOI.VariableIndex[i.index for i in load_var] - -# Instantiate the model -pm = instantiate_model( - network_data, - network_formulation, - PowerModels.build_opf; - setting=Dict("output" => Dict("duals" => true)), - jump_model=model, -) - -# Dualize the model -dual_st = Dualization.dualize(JuMP.backend(model), - variable_parameters = load_moi_idx -) - -dual_model = dual_st.dual_model -primal_dual_map = dual_st.primal_dual_map - -jump_dual_model = JuMP.Model() -map_moi_to_jump = MOI.copy_to(JuMP.backend(jump_dual_model), dual_model) -set_optimizer(jump_dual_model, solver) - -load_dual_idxs = [map_moi_to_jump[primal_dual_map.primal_parameter[l]].value for l in load_moi_idx] - -load_var_dual = JuMP.all_variables(jump_dual_model)[load_dual_idxs] - -for (i,l) in enumerate(load_var_dual) - @constraint(jump_dual_model, l <= original_load[i]*1.5) - @constraint(jump_dual_model, l >= 0.0) -end -for (i, l) in enumerate(values(network_data["load"])) - l["pd"] = load_var_dual[i] -end - -obj = objective_function(jump_dual_model) - -pm = instantiate_model( - network_data, - network_formulation, - PowerModels.build_opf; - setting=Dict("output" => Dict("duals" => true)), - jump_model=jump_dual_model, -) - -@objective(jump_dual_model, Max, obj) - -JuMP.optimize!(jump_dual_model) - -optimal_dual_cost = JuMP.objective_value(jump_dual_model) - -optimal_loads = value.(load_var_dual) - -load_diff = sum(abs.(optimal_loads - original_load)) ./ sum(original_load) - -# Create final primal model -network_data = deepcopy(network_data_original) -model = JuMP.Model(solver) -for (i, l) in enumerate(values(network_data["load"])) - l["pd"] = optimal_loads[i] -end - -# Instantiate the model -pm = instantiate_model( - network_data, - network_formulation, - PowerModels.build_opf; - setting=Dict("output" => Dict("duals" => true)), - jump_model=model, -) - -# Solve the model -JuMP.optimize!(model) - -# Check the solution -optimal_final_cost = JuMP.objective_value(model) -termination_status = JuMP.termination_status(model) -solution_primal_status = JuMP.primal_status(model) -solution_dual_status = JuMP.dual_status(model) -termination_status == MOI.INFEASIBLE && @error("Optimal solution not found") -solution_primal_status != MOI.FEASIBLE_POINT && @error("Primal solution not found") -solution_dual_status != MOI.FEASIBLE_POINT && @error("Dual solution not found") -isapprox(optimal_final_cost, optimal_dual_cost; rtol=1e-4) || @warn("Final cost is not equal to dual cost") From 85d86f4f351ddec9fec7de5b495a98d06f3ea48f Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Sat, 23 Sep 2023 13:26:08 -0400 Subject: [PATCH 33/36] update code --- .gitignore | 1 + examples/flux/flux_forecaster_script.jl | 1 + .../generate_full_datasets_script.jl | 179 +++++++++--------- examples/powermodels/pglib_datagen.jl | 49 ++++- 4 files changed, 135 insertions(+), 95 deletions(-) diff --git a/.gitignore b/.gitignore index 6cc3760..25db443 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,4 @@ Manifest.toml *.arrow *.m *.csv +*.config.toml diff --git a/examples/flux/flux_forecaster_script.jl b/examples/flux/flux_forecaster_script.jl index 659a352..8cf8b34 100644 --- a/examples/flux/flux_forecaster_script.jl +++ b/examples/flux/flux_forecaster_script.jl @@ -16,6 +16,7 @@ case_file_path = joinpath(path_dataset, case_name, string(network_formulation)) # Load input and output data tables iter_files = readdir(joinpath(case_file_path)) +iter_files = filter(x -> occursin(string(ArrowFile), x), iter_files) file_ins = [joinpath(case_file_path, file) for file in iter_files if occursin("input", file)] file_outs = [joinpath(case_file_path, file) for file in iter_files if occursin("output", file)] batch_ids = [split(split(file, "_")[end], ".")[1] for file in file_ins] diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 0451b63..763fc57 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -1,22 +1,27 @@ -# run with: julia ./examples/powermodels/generate_full_datasets_script.jl +# run with: julia ./examples/powermodels/generate_full_datasets_script.jl "./examples/powermodels/data/pglib_opf_case300_ieee/case300.config.toml" +config_path = ARGS[1] + using TestEnv -TestEnv.activate() +TestEnv.activate(dirname(config_path)) + +########## SCRIPT REQUIRED PACKAGES ########## -using Arrow using L2O +using Arrow using Test using UUIDs using PowerModels -using Clarabel import JuMP.MOI as MOI import ParametricOptInterface as POI -using Gurobi +using TOML -using NonconvexNLopt +## SOLVER PACKAGES ## -# using QuadraticToBinary +using Clarabel +using Gurobi +using NonconvexNLopt -########## SOLVERS ########## +########## POI SOLVER ########## cached = () -> MOI.Bridges.full_bridge_optimizer( MOI.Utilities.CachingOptimizer( @@ -28,99 +33,91 @@ cached = () -> MOI.Bridges.full_bridge_optimizer( POI_cached_optimizer() = POI.Optimizer(cached()) -########## DATASET GENERATION ########## +########## PARAMETERS ########## + +config = TOML.parsefile(config_path) +path = config["export_dir"] -# Paths -path_powermodels = joinpath(pwd(), "examples", "powermodels") -path = joinpath(path_powermodels, "data") +path_powermodels = joinpath(@__FILE__, "data") # TODO: Make it a submodule include(joinpath(path_powermodels, "pglib_datagen.jl")) -# Parameters -filetype = ArrowFile # CSVFile # ArrowFile +filetype = ArrowFile -# Case name -case_name = "pglib_opf_case300_ieee" # pglib_opf_case300_ieee # pglib_opf_case5_pjm -network_formulation = SOCWRConicPowerModel # SOCWRConicPowerModel # DCPPowerModel +case_name = config["case_name"] case_file_path = joinpath(path, case_name) mkpath(case_file_path) +network_formulation=ARGS[2] + +########## SAMPLER DATASET GENERATION ########## + +if haskey(config, "sampler") + num_batches = config["sampler"]["num_batches"] + num_p = config["sampler"]["num_samples"] + global success_solves = 0.0 + for i in 1:num_batches + _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; + num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, + internal_load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01) + ) + global success_solves += _success_solves + end + success_solves /= num_batches + + @info "Success solves Normal: $(success_solves)" +end -# Generate dataset -# num_batches = 1 -# num_p = 200 -# global success_solves = 0.0 -# for i in 1:num_batches -# _success_solves, number_variables, number_loads, batch_id = generate_dataset_pglib(case_file_path, case_name; -# num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, -# internal_load_sampler= (_o, n) -> load_sampler(_o, n, max_multiplier=1.25, min_multiplier=0.8, step_multiplier=0.01) -# ) -# global success_solves += _success_solves -# end -# success_solves /= num_batches - -# @info "Success solves Normal: $(success_solves)" - -# Generate Line-search dataset - -matpower_case_name = case_name * ".m" -network_data = make_basic_network(pglib(matpower_case_name)) - -early_stop_fn = (model, status, recorder) -> !status -step_multiplier = 1.01 -num_loads = length(network_data["load"]) -num_batches = num_loads + 1 -num_p = 2 - -function line_sampler(_o, n, idx, num_inputs, ibatc) - if (idx == ibatc) || (idx - num_loads == ibatc) || (ibatc == num_inputs + 1) - return [_o * step_multiplier ^ (j-1) for j in 1:n] - else - return ones(n) * _o +########## LINE SEARCH DATASET GENERATION ########## + +if haskey(config, "line_search") + network_data = make_basic_network(pglib(case_name * ".m")) + step_multiplier = 1.01 + num_loads = length(network_data["load"]) + num_batches = num_loads + 1 + num_p = config["line_search"]["num_samples"] + early_stop_fn = (model, status, recorder) -> !status + + global success_solves = 0.0 + global batch_id = string(uuid1()) + for ibatc in 1:num_batches + _success_solves, number_variables, number_loads, b_id = generate_dataset_pglib(case_file_path, case_name; + num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, + internal_load_sampler= (_o, n, idx, num_inputs) -> line_sampler(_o, n, idx, num_inputs, ibatc; step_multiplier=step_multiplier), + early_stop_fn=early_stop_fn, + batch_id=batch_id, + ) + global success_solves += _success_solves end + success_solves /= num_batches + + @info "Success solves: $(success_solves * 100) % of $(num_batches * num_p)" end -global success_solves = 0.0 -global batch_id = string(uuid1()) -for ibatc in 1:num_batches - _success_solves, number_variables, number_loads, b_id = generate_dataset_pglib(case_file_path, case_name; +########## WORST CASE DUAL DATASET GENERATION ########## +if haskey(config, "worst_case_dual") + num_p = config["worst_case_dual"]["num_samples"] + function optimizer_factory() + IPO_OPT = Gurobi.Optimizer() + # IPO_OPT = MadNLP.Optimizer(print_level=MadNLP.INFO, max_iter=100) + # IPO = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(IPO_OPT) + # MIP = QuadraticToBinary.Optimizer{Float64}(IPO) + return () -> IPO_OPT + end + + success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset(case_file_path, case_name; + num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer_factory=optimizer_factory, + hook = (model) -> set_optimizer_attribute(model, "NonConvex", 2) + ) + + @info "Success solves Worst Case: $(success_solves) of $(num_p)" +end + +########## WORST CASE NONCONVEX DATASET GENERATION ########## +if haskey(config, "worst_case_nonconvex") + num_p = config["worst_case_nonconvex"]["num_samples"] + + success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset_Nonconvex(case_file_path, case_name; num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, - internal_load_sampler= (_o, n, idx, num_inputs) -> line_sampler(_o, n, idx, num_inputs, ibatc), - early_stop_fn=early_stop_fn, - batch_id=batch_id, ) - global success_solves += _success_solves + + @info "Success solves Worst Case: $(success_solves * 100) of $(num_p)" end -success_solves /= num_batches - -@info "Success solves: $(success_solves * 100) % of $(num_batches * num_p)" - -# Generate worst case dataset - -# function optimizer_factory() -# IPO_OPT = Gurobi.Optimizer() # MadNLP.Optimizer(print_level=MadNLP.INFO, max_iter=100) -# # IPO = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(IPO_OPT) -# # MIP = QuadraticToBinary.Optimizer{Float64}(IPO) -# return () -> IPO_OPT -# end - -# global success_solves = 0.0 -# for i in 1:num_batches -# _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset(case_file_path, case_name; -# num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer_factory=optimizer_factory, -# hook = (model) -> set_optimizer_attribute(model, "NonConvex", 2) -# ) -# global success_solves += _success_solves -# end -# success_solves /= num_batches - -# @info "Success solves Worst Case: $(success_solves) of $(num_batches * num_p)" - -# global success_solves = 0.0 -# for i in 1:num_batches -# _success_solves, number_variables, number_loads, batch_id = generate_worst_case_dataset_Nonconvex(case_file_path, case_name; -# num_p=num_p, filetype=filetype, network_formulation=network_formulation, optimizer=POI_cached_optimizer, -# ) -# global success_solves += _success_solves -# end -# success_solves /= num_batches - -# @info "Success solves Worst Case: $(success_solves * 100) of $(num_batches * num_p)" diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 10c096e..fe447ba 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -47,6 +47,37 @@ function load_sampler( return load_samples end +""" + line_sampler(original_parameter::T, num_p::Int, parameter_index::F, num_inputs::F, line_index::F; step_multiplier::T=0.1) + +line_sampler is a function to help generate a dataset for varying parameter values. It has two modes: + - If line_index is not outside the parameter index range: + Return an incremental vector for the parameter at parameter_index and an unchanged parameter for the rest; + - If line_index is outside the parameter index range: + Return an incremental vector for all parameters. +""" +function line_sampler( + original_parameter::T, + num_p::F, + parameter_index::F, + num_inputs::F, + line_index::F; + step_multiplier::T=1.01, +) where {T<:Real, F<:Integer} + # parameter sampling + num_parameters = floor(Int, num_inputs / 2) + if (parameter_index == line_index) || (parameter_index - num_parameters == line_index) || (line_index == num_inputs + 1) + return [original_parameter * step_multiplier ^ (j) for j in 1:num_p] + else + return ones(num_p) * original_parameter + end +end + +""" + load_parameter_factory(model, indices; load_set=nothing) + +load_parameter_factory is a function to help generate a parameter vector for the problem iterator. +""" function load_parameter_factory(model, indices; load_set=nothing) if isnothing(load_set) return @variable( @@ -58,6 +89,11 @@ function load_parameter_factory(model, indices; load_set=nothing) ) end +""" + pm_primal_builder!(model, parameters, network_data, network_formulation; recorder=nothing) + +pm_primal_builder! is a function to help build a PowerModels model and update recorder. +""" function pm_primal_builder!(model, parameters, network_data, network_formulation; recorder=nothing) num_loads = length(network_data["load"]) for (str_i, l) in network_data["load"] @@ -88,6 +124,11 @@ function pm_primal_builder!(model, parameters, network_data, network_formulation return nothing end +""" + load_set_iterator!(model, parameters, idx, original_load) + +load_set_iterator! is a function to help iterate over the load set. Used in the worst case generator. +""" function load_set_iterator!(model, parameters, idx, original_load) for (i, p) in enumerate(parameters) @constraint(model, p <= original_load[i] * (1.0 + 0.1 * idx)) @@ -165,10 +206,6 @@ function generate_dataset_pglib( batch_id end -function default_optimizer_factory() - return () -> Ipopt.Optimizer() -end - function generate_worst_case_dataset_Nonconvex(data_dir, case_name; filetype=CSVFile, @@ -250,6 +287,10 @@ function generate_worst_case_dataset_Nonconvex(data_dir, batch_id end +function default_optimizer_factory() + return () -> Ipopt.Optimizer() +end + function generate_worst_case_dataset(data_dir, case_name; filetype=CSVFile, From 337efe9940e412e4220410812d6ace2748a8c78d Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Sat, 23 Sep 2023 13:27:32 -0400 Subject: [PATCH 34/36] add toml example --- .gitignore | 1 - .../case300.config.toml | 19 +++++++++++++++++++ 2 files changed, 19 insertions(+), 1 deletion(-) create mode 100644 examples/powermodels/data/pglib_opf_case300_ieee/case300.config.toml diff --git a/.gitignore b/.gitignore index 25db443..6cc3760 100644 --- a/.gitignore +++ b/.gitignore @@ -8,4 +8,3 @@ Manifest.toml *.arrow *.m *.csv -*.config.toml diff --git a/examples/powermodels/data/pglib_opf_case300_ieee/case300.config.toml b/examples/powermodels/data/pglib_opf_case300_ieee/case300.config.toml new file mode 100644 index 0000000..ae5dcc6 --- /dev/null +++ b/examples/powermodels/data/pglib_opf_case300_ieee/case300.config.toml @@ -0,0 +1,19 @@ +# Name of the reference PGLib case. Must be a valid PGLib case name. +case_name = "pglib_opf_case300_ieee" + +# Directory where instance/solution files are exported +# must be a valid directory +export_dir = "EXPORT_DIR" + +# [sampler] +# num_batches = 1 +# num_samples = 10 + +[line_search] +num_samples = 10 + +# [worst_case_dual] +# num_samples = 10 + +# [worst_case_nonconvex] +# num_samples = 10 From d9d08899236e85539eed82e2bc53011bcee77192 Mon Sep 17 00:00:00 2001 From: Andrew Rosemberg Date: Sat, 23 Sep 2023 16:31:50 -0400 Subject: [PATCH 35/36] add line example --- examples/flux/flux_forecaster_script.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/flux/flux_forecaster_script.jl b/examples/flux/flux_forecaster_script.jl index 8cf8b34..cd5ec4c 100644 --- a/examples/flux/flux_forecaster_script.jl +++ b/examples/flux/flux_forecaster_script.jl @@ -42,6 +42,9 @@ output_data_test = DataFrame(output_table_test) output_variables_train = output_data_train[!, Not(:id)] input_features_train = innerjoin(input_data_train, output_data_train[!, [:id]], on = :id)[!, Not(:id)] # just use success solves +num_loads = floor(Int,size(input_features_train,2)/2) +total_volume=[sum(sqrt(input_features_train[i,l]^2 + input_features_train[i,l+num_loads]^2) for l in 1:num_loads) for i in 1:size(input_features_train,1) ] + output_variables_test = output_data_test[!, Not(:id)] input_features_test = innerjoin(input_data_test, output_data_test[!, [:id]], on = :id)[!, Not(:id)] # just use success solves From 2cac65bb45a80dc80cf2e1ca5925f5f7aa6bc272 Mon Sep 17 00:00:00 2001 From: Andrew David Werner Rosemberg Date: Sat, 23 Sep 2023 18:39:07 -0400 Subject: [PATCH 36/36] update code --- examples/powermodels/generate_full_datasets_script.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 763fc57..503240a 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -1,8 +1,10 @@ # run with: julia ./examples/powermodels/generate_full_datasets_script.jl "./examples/powermodels/data/pglib_opf_case300_ieee/case300.config.toml" config_path = ARGS[1] +import Pkg; Pkg.activate(".") + using TestEnv -TestEnv.activate(dirname(config_path)) +TestEnv.activate() ########## SCRIPT REQUIRED PACKAGES ########## @@ -15,6 +17,8 @@ import JuMP.MOI as MOI import ParametricOptInterface as POI using TOML +PowerModels.silence() + ## SOLVER PACKAGES ## using Clarabel @@ -38,7 +42,7 @@ POI_cached_optimizer() = POI.Optimizer(cached()) config = TOML.parsefile(config_path) path = config["export_dir"] -path_powermodels = joinpath(@__FILE__, "data") # TODO: Make it a submodule +path_powermodels = joinpath(dirname(@__FILE__)) # TODO: Make it a submodule include(joinpath(path_powermodels, "pglib_datagen.jl")) filetype = ArrowFile @@ -46,7 +50,7 @@ filetype = ArrowFile case_name = config["case_name"] case_file_path = joinpath(path, case_name) mkpath(case_file_path) -network_formulation=ARGS[2] +network_formulation= eval(Symbol(ARGS[2])) ########## SAMPLER DATASET GENERATION ##########