diff --git a/Project.toml b/Project.toml index 613929b..875ef7f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,31 +1,39 @@ 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" 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" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] Arrow = "2" CSV = "0.10" +Dualization = "0.5" JuMP = "1" 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" -Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +NonconvexNLopt = "b43a31b8-ff9b-442d-8e31-c163daa8ab75" +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", "NonconvexNLopt"] diff --git a/examples/flux/flux_forecaster_script.jl b/examples/flux/flux_forecaster_script.jl index 2d6d670..cd5ec4c 100644 --- a/examples/flux/flux_forecaster_script.jl +++ b/examples/flux/flux_forecaster_script.jl @@ -5,16 +5,18 @@ using Arrow using Flux using DataFrames using PowerModels +using L2O # Paths -path_dataset = joinpath(pwd(), "examples", "powermodels", "data") -case_name = "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 -network_formulation = DCPPowerModel -case_file_path = joinpath(path, case_name, string(network_formulation)) +path_dataset = joinpath(pwd(), "examples", "powermodels", "data") +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] @@ -40,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 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 diff --git a/examples/powermodels/generate_full_datasets_script.jl b/examples/powermodels/generate_full_datasets_script.jl index 4475844..503240a 100644 --- a/examples/powermodels/generate_full_datasets_script.jl +++ b/examples/powermodels/generate_full_datasets_script.jl @@ -1,16 +1,33 @@ +# 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() -using Arrow +########## SCRIPT REQUIRED PACKAGES ########## + using L2O +using Arrow using Test using UUIDs using PowerModels -using Clarabel import JuMP.MOI as MOI import ParametricOptInterface as POI +using TOML + +PowerModels.silence() + +## SOLVER PACKAGES ## -cached = MOI.Bridges.full_bridge_optimizer( +using Clarabel +using Gurobi +using NonconvexNLopt + +########## POI SOLVER ########## + +cached = () -> MOI.Bridges.full_bridge_optimizer( MOI.Utilities.CachingOptimizer( MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), Clarabel.Optimizer(), @@ -18,29 +35,93 @@ cached = MOI.Bridges.full_bridge_optimizer( Float64, ) -# Paths -path_powermodels = joinpath(pwd(), "examples", "powermodels") -path = joinpath(path_powermodels, "data") +POI_cached_optimizer() = POI.Optimizer(cached()) + +########## PARAMETERS ########## + +config = TOML.parsefile(config_path) +path = config["export_dir"] + +path_powermodels = joinpath(dirname(@__FILE__)) # TODO: Make it a submodule include(joinpath(path_powermodels, "pglib_datagen.jl")) -# Parameters -num_batches = 1 -num_p = 10 filetype = ArrowFile -# Case name -case_name = "pglib_opf_case300_ieee" -network_formulation = SOCWRConicPowerModel +case_name = config["case_name"] case_file_path = joinpath(path, case_name) -solver = () -> POI.Optimizer(cached()) - -# Generate dataset -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) +mkpath(case_file_path) +network_formulation= eval(Symbol(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 + +########## 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 + +########## 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) ) - success_solves += _success_solves + + @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, + ) + + @info "Success solves Worst Case: $(success_solves * 100) of $(num_p)" end -success_solves /= num_batches diff --git a/examples/powermodels/pglib_datagen.jl b/examples/powermodels/pglib_datagen.jl index 019899e..fe447ba 100644 --- a/examples/powermodels/pglib_datagen.jl +++ b/examples/powermodels/pglib_datagen.jl @@ -1,7 +1,8 @@ -using Downloads +using PGLib using PowerModels using JuMP, HiGHS import ParametricOptInterface as POI +using LinearAlgebra """ return_variablerefs(pm::AbstractPowerModel) @@ -19,6 +20,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) @@ -26,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) @@ -38,7 +48,96 @@ function load_sampler( end """ - generate_dataset_pglib(data_dir::AbstractString, case_name::AbstractString; download_files::Bool=true, filetype::Type{RecorderFile}, + 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( + model, _p[i=indices] + ) + end + return @variable( + model, _p[i=indices] in load_set + ) +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"] + i = parse(Int, str_i) + 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"])]) + return model, parameters, variable_refs + end + + 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)) + @constraint(model, p >= original_load[i] * (1.0 - 0.1 * idx)) + end +end + +""" + generate_dataset_pglib(data_dir::AbstractString, case_name::AbstractString; download_files::Bool=true, filetype::Type{FileType}, num_p::Int=10 ) @@ -48,24 +147,74 @@ function generate_dataset_pglib( data_dir, case_name; filetype=CSVFile, - download_files=true, num_p=10, - load_sampler=load_sampler, + internal_load_sampler=load_sampler, network_formulation=DCPPowerModel, - solver = () -> POI.Optimizer(HiGHS.Optimizer()), + optimizer = () -> POI.Optimizer(HiGHS.Optimizer()), + filterfn=L2O.filter_fn, + early_stop_fn = (model, status, recorder) -> false, + batch_id = string(uuid1()) ) - # Download file + @info "Batch ID: $batch_id" + + # save folder + data_sim_dir = joinpath(data_dir, string(network_formulation)) + mkpath(data_sim_dir) + + # Read data 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, - ) + 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"]) + num_inputs = num_loads * 2 + original_load = vcat( + [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)) + + # Build model and Recorder + file = joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_output_" * batch_id * "." * string(filetype)) + recorder = Recorder{filetype}(file; filterfn=filterfn) + 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]] = internal_load_sampler(original_load[i], num_p, i, num_inputs) end + problem_iterator = ProblemIterator( + pairs; + early_stop=early_stop_fn + ) + save( + problem_iterator, + joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_input_" * batch_id * "." * string(filetype)), + filetype, + ) + + # Solve the problem and return the number of successfull solves + return solve_batch(problem_iterator, recorder), + length(recorder.primal_variables), + length(original_load), + batch_id +end + +function generate_worst_case_dataset_Nonconvex(data_dir, + case_name; + filetype=CSVFile, + num_p=10, + network_formulation=DCPPowerModel, + optimizer = () -> POI.Optimizer(HiGHS.Optimizer()), + algorithm = NLoptAlg(:LN_BOBYQA), + options = NLoptOptions(maxeval=10), +) # save folder data_sim_dir = joinpath(data_dir, string(network_formulation)) if !isdir(data_sim_dir) @@ -73,65 +222,140 @@ function generate_dataset_pglib( end # Read data - network_data = PowerModels.parse_file(case_file_path) + matpower_case_name = case_name * ".m" + network_data = make_basic_network(pglib(matpower_case_name)) # The problem to iterate over - model = Model(solver) + model = JuMP.Model(() -> optimizer()) 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 + 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" - # Instantiate the model - pm = instantiate_model( - network_data, - network_formulation, - PowerModels.build_opf; - setting=Dict("output" => Dict("duals" => true)), - jump_model=model, + # 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) + if !isnothing(recorder) + set_primal_variable!(recorder, variable_refs) + end + + return model, parameters + end + + # Set iterator + function _set_iterator!(idx) + _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 + # The problem iterator - problem_iterator = ProblemIterator( - Dict( - p .=> [ - load_sampler(original_load[i], num_p) for - i in 1:length(network_data["load"]) - ], - ), + problem_iterator = WorstCaseProblemIterator( + [uuid1() for _ in 1:num_p], # will be ignored + () -> nothing, # will be ignored + _primal_builder!, + _set_iterator!, + algorithm; + options = options + ) + + # Solve all problems and record solutions + return solve_batch(problem_iterator, recorder), + length(recorder.primal_variables), + length(original_load), + 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_factory = default_optimizer_factory, + hook = nothing +) + # 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" - save( - problem_iterator, - joinpath(data_sim_dir, case_name * "_" * string(network_formulation) * "_input_" * batch_id * "." * string(filetype)), - filetype, + + # 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_factory; + hook=hook ) - # 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(model, problem_iterator, recorder), - number_vars, - length(original_load), - 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=[] + ) + + # 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 - 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") @@ -139,12 +363,62 @@ 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) @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 +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 * 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 + end +end + +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_Nonconvex( + 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 end diff --git a/src/L2O.jl b/src/L2O.jl index 65f46ec..a0a0b2f 100644 --- a/src/L2O.jl +++ b/src/L2O.jl @@ -2,15 +2,23 @@ module L2O using Arrow using CSV +using Dualization using JuMP using UUIDs import ParametricOptInterface as POI +import JuMP.MOI as MOI import Base: string -export ArrowFile, CSVFile, ProblemIterator, Recorder, save, solve_batch +using Nonconvex +using Zygote + +export ArrowFile, CSVFile, ProblemIterator, Recorder, save, solve_batch, + WorstCaseProblemIterator, set_primal_variable!, set_dual_variable! include("datasetgen.jl") include("csvrecorder.jl") include("arrowrecorder.jl") +include("worst_case.jl") +include("worst_case_iter.jl") end diff --git a/src/arrowrecorder.jl b/src/arrowrecorder.jl index 816d1b9..8374d17 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,23 +7,42 @@ 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 + owner_model(recorder.dual_variables[1]) + 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( - recorder.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], - )..., - ), + _filename, + df, ) end -function save(table::NamedTuple, filename::String, ::Type{ArrowFile}) - return Arrow.write(filename, table) +function save(table::NamedTuple, filename::String, ::Type{ArrowFile}; kwargs...) + Arrow.append( + filename, + table; + kwargs... + ) end diff --git a/src/csvrecorder.jl b/src/csvrecorder.jl index ed98013..c51385d 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))") @@ -17,10 +18,13 @@ function record(recorder::Recorder{CSVFile}, id::UUID) for p in recorder.dual_variables write(f, ",dual_$(name(p))") end + if !input + write(f, ",objective") + end 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) @@ -30,10 +34,27 @@ 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 + if !input + obj = JuMP.objective_value(model) + write(f, ",$obj") + end + # end line write(f, "\n") end end -function save(table::NamedTuple, filename::String, ::Type{CSVFile}) - return CSV.write(filename, table) +function save(table::NamedTuple, filename::String, ::Type{CSVFile}; kwargs...) + isappend = isfile(filename) + mode = isappend ? "append" : "write" + @info "Saving CSV file to $filename - Mode: $mode" + CSV.write(filename, table; append=isappend) + return nothing end diff --git a/src/datasetgen.jl b/src/datasetgen.jl index f1db811..c898d41 100644 --- a/src/datasetgen.jl +++ b/src/datasetgen.jl @@ -1,47 +1,103 @@ -abstract type RecorderFile end +abstract type FileType end + +mutable struct RecorderFile{T<:FileType} + filename::String +end + +filename(recorder_file::RecorderFile) = recorder_file.filename + +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)) """ 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 + dual_variables::Vector 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) + filterfn=filter_fn, + ) 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) +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 +end + +function set_dual_variable!(recorder::Recorder, p::Vector) + recorder.dual_variables = p +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}} + 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}(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 """ @@ -50,18 +106,22 @@ end Save optimization problem instances to a file. """ function save( - problem_iterator::ProblemIterator, filename::String, file_type::Type{T} -) where {T<:RecorderFile} - return save( - (; - id=problem_iterator.ids, - zip( - Symbol.(name.(keys(problem_iterator.pairs))), values(problem_iterator.pairs) - )..., - ), + problem_iterator::AbstractProblemIterator, filename::String, file_type::Type{T} +) where {T<:FileType} + kys = sort(collect(keys(problem_iterator.pairs)); by=(v) -> index(v).value) + df = (; + id=problem_iterator.ids, + ) + for ky in kys + df = merge(df, (; Symbol(ky) => problem_iterator.pairs[ky])) + end + save( + df, filename, - file_type, + file_type; + dictencode=true, ) + return nothing end """ @@ -85,35 +145,44 @@ 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) + 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 """ - solve_batch(model::JuMP.Model, problem_iterator::ProblemIterator, recorder::Recorder) + solve_batch(problem_iterator::AbstractProblemIterator, recorder) Solve a batch of optimization problems and record the solutions. """ function solve_batch( - model::JuMP.Model, problem_iterator::ProblemIterator, recorder::Recorder + problem_iterator::AbstractProblemIterator, recorder ) - successfull_solves = - sum( - solve_and_record(model, 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 new file mode 100644 index 0000000..698d09f --- /dev/null +++ b/src/worst_case.jl @@ -0,0 +1,149 @@ +""" + 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{F} <: AbstractProblemIterator + ids::Vector{UUID} + parameters::Function + primal_builder!::Function + set_iterator!::Function + optimizer::F + hook::Union{Nothing, Function} + options::Any + ext::Dict + early_stop::Function + function WorstCaseProblemIterator( + ids::Vector{UUID}, + parameters::Function, + primal_builder!::Function, + set_iterator!::Function, + optimizer::F; + hook::Union{Nothing, Function}=nothing, + options::Any=nothing, + ext::Dict=Dict(), + early_stop::Function=(args...) -> false + ) where {F} + new{F}(ids, parameters, primal_builder!, set_iterator!, optimizer, hook, options, ext, early_stop) + 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 +) + # Build Primal + model = JuMP.Model() + parameters = problem_iterator.parameters(model) + problem_iterator.primal_builder!(model, parameters) + + # Parameter indices + load_moi_idx = Vector{MOI.VariableIndex}(JuMP.index.(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()) + 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] + 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, load_var_dual, 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) + + # 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, early_stop_bool + 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) + 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) + solution_primal_status = JuMP.primal_status(model) + solution_dual_status = JuMP.dual_status(model) + 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 + @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, early_stop_bool + end + 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 new file mode 100644 index 0000000..6a590fa --- /dev/null +++ b/src/worst_case_iter.jl @@ -0,0 +1,96 @@ +# 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 + status = false + try + JuMP.optimize!(model) + + status = filter_fn(model) + catch e + @warn "Error in solve_and_record: $e" + end + + objective = if status + JuMP.objective_value(model) + else + - penalty + end + + return objective +end + +function save_input(parameters, _recorder, id) + recorder = similar(_recorder) + set_primal_variable!(recorder, parameters) + set_dual_variable!(recorder, []) + 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 begin + 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 +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 + ) + + 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]) + addvar!(model_non, min_demands, max_demands) + # add_ineq_constraint!(model_non, x -> sum(x .^ 2) - max_total_volume ^ 2) + + # Optimize model_non: + 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_Nonconvex.minimizer + # best_profit = -r_Nonconvex.minimum + + return storage_objective_function.success_solves / storage_objective_function.fcalls, false +end \ No newline at end of file diff --git a/test/datasetgen.jl b/test/datasetgen.jl index 27f764c..15d6959 100644 --- a/test/datasetgen.jl +++ b/test/datasetgen.jl @@ -1,13 +1,13 @@ """ - 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())) + 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) @@ -15,50 +15,70 @@ function testdataset_gen(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(model, problem_iterator, recorder) + @testset "early_stop" begin + file_dual_output = joinpath(path, "test_$(string(uuid1()))_output.$(string(filetype))") # file path + recorder_dual = Recorder{filetype}( + file_dual_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_dual) + @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, :]) == 3 - 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) == 3 - @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 + 1 # 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 + rm(file_output) + end end end end diff --git a/test/runtests.jl b/test/runtests.jl index d742466..015075c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,20 +7,30 @@ using L2O import ParametricOptInterface as POI using Test using UUIDs +using Ipopt + +using NonconvexNLopt 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) + file_in, file_out = test_generate_worst_case_dataset(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 new file mode 100644 index 0000000..31b1ada --- /dev/null +++ b/test/worst_case.jl @@ -0,0 +1,147 @@ +""" + 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, num_p=10) + @testset "Worst Case Dual Generation Type: $filetype" for filetype in [CSVFile, ArrowFile] + # The problem to iterate over + function optimizer_factory() + return () -> Ipopt.Optimizer() + end + 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 + + problem_iterator = WorstCaseProblemIterator( + [uuid1() for _ in 1:num_p], + parameter_factory, + primal_builder!, + set_iterator!, + optimizer_factory, + ) + + # file_names + 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}( + file_output; filename_input=file_input, + primal_variables=[], dual_variables=[] + ) + + # Solve all problems and record solutions + 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 *successfull_solves + 1 + @test length(readdlm(file_input, ',')[1, :]) == 2 + rm(file_input) + # test output file + @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 * successfull_solves + 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 + + @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)) + @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 + + problem_iterator = WorstCaseProblemIterator( + [uuid1() for _ in 1:num_p], # will be ignored + () -> nothing, # will be ignored + _primal_builder!, + _set_iterator!, + NLoptAlg(:LN_BOBYQA); + options = NLoptOptions(maxeval=10) + ) + + # file_names + 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}( + file_output; filename_input=file_input, + primal_variables=[], dual_variables=[] + ) + + # Solve all problems and record solutions + 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 * successfull_solves + 1 + @test length(readdlm(file_input, ',')[1, :]) == 2 + rm(file_input) + # test output file + @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 * successfull_solves + rm(file_input) + # test output file + df = Arrow.Table(file_output) + @test length(df) == 4 + @test length(df[1]) >= num_p * successfull_solves + rm(file_output) + end + end +end