Add AbstractKKTSystem structure

benchmark/.Project.toml

@@ -4,4 +4,5 @@ Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655"
+PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"

benchmark/README.md

@@ -1,5 +1,9 @@
-To run benchmarks, run the following script.
+Running benchmark requires pacakge ArgParse. First install ArgParse with:
 ```
-# to compare current version and master using 10 cores and CUTEst and PowerModels.jl problems,
-julia runbenchmarks.jl 10 power cutest current master 
+julia -e 'import Pkg; Pkg.add("ArgParse")'
+```
+
+To get a help for running benchmark, run:
+```
+julia runbenchmark.jl --help
 ```

benchmark/benchmark-cutest.jl

@@ -1,6 +1,6 @@
 include("config.jl")
 
-function get_status(code::Symbol) 
+function get_status(code::Symbol)
     if code == :first_order
         return 1
     elseif code == :acceptable
@@ -12,14 +12,11 @@ end
 
 @everywhere using CUTEst
 
-if SOLVER == "master"
-    @everywhere solver = nlp -> madnlp(nlp,linear_solver=MadNLPMa57,max_wall_time=900.,tol=1e-6)
-    @everywhere using MadNLP, MadNLPHSL
-elseif SOLVER == "current"
-    @everywhere solver = nlp -> madnlp(nlp,linear_solver=MadNLPMa57,max_wall_time=900.,tol=1e-6)
+if SOLVER == "master" || SOLVER == "current"
     @everywhere using MadNLP, MadNLPHSL
+    @everywhere solver = nlp -> madnlp(nlp,linear_solver=MadNLPMa57,max_wall_time=900.,tol=1e-6, print_level=PRINT_LEVEL)
 elseif SOLVER == "ipopt"
-    @everywhere solver = nlp -> ipopt(nlp,linear_solver="ma57",max_cpu_time=900.,tol=1e-6)
+    @everywhere solver = nlp -> ipopt(nlp,linear_solver="ma57",max_cpu_time=900.,tol=1e-6, print_level=PRINT_LEVEL)
     @everywhere using NLPModelsIpopt
 elseif SOLVER == "knitro"
     # TODO
@@ -33,47 +30,64 @@ end
     finalize(CUTEstModel(name))
 end
 
-@everywhere function evalmodel(name,solver)
+@everywhere function evalmodel(name,solver;gcoff=false)
     println("Solving $name")
     nlp = CUTEstModel(name; decode=false)
-    t = @elapsed begin
-        retval = solver(nlp)
+    try
+        gcoff && GC.enable(false);
+        mem = @allocated begin
+            t = @elapsed begin
+                retval = solver(nlp)
+            end
+        end
+        gcoff && GC.enable(true);
+        retval.elapsed_time = t
+        retval.solver_specific[:mem] = mem
+        finalize(nlp)
+        return retval
+    catch e
+        finalize(nlp)
+        throw(e)
     end
-    retval.elapsed_time = t
-    finalize(nlp)
-    return retval
 end
 
 function benchmark(solver,probs;warm_up_probs = [])
     println("Warming up (forcing JIT compile)")
     broadcast(decodemodel,warm_up_probs)
-    [remotecall_fetch.(prob->evalmodel(prob,solver),i,warm_up_probs) for i in procs() if i!= 1]
+    [remotecall_fetch.(prob->evalmodel(prob,solver;gcoff=GCOFF),i,warm_up_probs) for i in procs() if i!= 1]
 
     println("Decoding problems")
     broadcast(decodemodel,probs)
-    
+
     println("Solving problems")
     retvals = pmap(prob->evalmodel(prob,solver),probs)
-    time   = [retval.elapsed_time for retval in retvals]
     status = [get_status(retval.status) for retval in retvals]
-    time,status
+    time   = [retval.elapsed_time for retval in retvals]
+    mem    = [retval.solver_specific[:mem] for retval in retvals]
+    iter   = [retval.iter for retval in retvals]
+    status,time,mem,iter
 end
 
 exclude = [
     "PFIT1","PFIT2","PFIT4","DENSCHNE","SPECANNE","DJTL", "EG3","OET7",
     "PRIMAL3","TAX213322","TAXR213322","TAX53322","TAXR53322","HIMMELP2","MOSARQP2","LUKVLE11",
     "CYCLOOCT","CYCLOOCF","LIPPERT1","GAUSSELM","A2NSSSSL",
-    "YATP1LS","YATP2LS","YATP1CLS","YATP2CLS","BA-L52LS","BA-L73LS","BA-L21LS","CRESC132"    
+    "YATP1LS","YATP2LS","YATP1CLS","YATP2CLS","BA-L52LS","BA-L73LS","BA-L21LS","CRESC132"
 ]
 
 
-probs = CUTEst.select()
+if QUICK
+    probs = readdlm("cutest-quick-names.csv")[:]
+else
+    probs = CUTEst.select()
+end
 
-println(probs)
 filter!(e->!(e in exclude),probs)
 
-time,status = benchmark(solver,probs;warm_up_probs = ["EIGMINA"])
+status,time,mem,iter = benchmark(solver,probs;warm_up_probs = ["EIGMINA"])
 
 writedlm("name-cutest.csv",probs,',')
+writedlm("status-cutest-$(SOLVER).csv",status),','
 writedlm("time-cutest-$(SOLVER).csv",time,',')
-writedlm("status-cutest-$(SOLVER).csv",status),',' 
+writedlm("mem-cutest-$(SOLVER).csv",mem,',')
+writedlm("iter-cutest-$(SOLVER).csv",iter,',')

benchmark/benchmark-power.jl

@@ -1,76 +1,122 @@
-const PGLIB_PATH = ENV["PGLIB_PATH"]
+if haskey(ENV, "PGLIB_PATH")
+    const PGLIB_PATH = ENV["PGLIB_PATH"]
+else
+    error("Unable to find path to PGLIB benchmark.\n"*
+          "Please set environment variable `PGLIB_PATH` to run benchmark with PowerModels.jl")
+end
 
 include("config.jl")
 
-using MathOptInterface; const MOI=MathOptInterface
-function get_status(code::MOI.TerminationStatusCode) 
-    if code == MOI.LOCALLY_SOLVED
-        return 1
-    elseif code == MOI.ALMOST_OPTIMAL
-        return 2
-    else
-        return 3
+@everywhere begin
+    using PowerModels, MathOptInterface, JuMP
+    const MOI = MathOptInterface
+
+    PowerModels.silence()
+
+    function evalmodel(prob,solver;gcoff=false)
+        println("Solving $(get_name(prob))")
+        gcoff && GC.enable(false);
+        return solver(prob)
     end
-end
 
-@everywhere using PowerModels
-@everywhere PowerModels.silence()
+    function get_status(code::MOI.TerminationStatusCode)
+        if code == MOI.LOCALLY_SOLVED
+            return 1
+        elseif code == MOI.ALMOST_OPTIMAL
+            return 2
+        else
+            return 3
+        end
+    end
+
+    get_name(pm) = "$(pm.data["name"])-$(typeof(pm))"
+end
 
-if SOLVER == "master"
-    @everywhere solver = prob ->
-        run_opf(joinpath(PGLIB_PATH,prob[1]), prob[2],
-                ()->MadNLP.Optimizer(linear_solver=MadNLPMa57,max_wall_time=900.,tol=1e-6))
-    @everywhere using MadNLP, MadNLPHSL
-elseif SOLVER == "current"
-    @everywhere solver = prob ->
-        run_opf(joinpath(PGLIB_PATH,prob[1]), prob[2],
-                ()->MadNLP.Optimizer(linear_solver=MadNLPMa57,max_wall_time=900.,tol=1e-6))
-    @everywhere using MadNLP, MadNLPHSL
+if SOLVER == "master" || SOLVER == "current"
+    @everywhere begin
+        using MadNLP, MadNLPHSL
+        solver = pm -> begin
+            set_optimizer(pm.model,()->
+                MadNLP.Optimizer(linear_solver=MadNLPMa57,max_wall_time=900.,tol=1e-6, print_level=PRINT_LEVEL))
+            mem=@allocated begin
+                t=@elapsed begin
+                    optimize_model!(pm)
+                end
+            end
+            return get_status(termination_status(pm.model)),t,mem,barrier_iterations(pm.model)
+        end
+    end
 elseif SOLVER == "ipopt"
-    @everywhere solver = prob ->
-        run_opf(joinpath(PGLIB_PATH,prob[1]), prob[2],
-                ()->Ipopt.Optimizer(linear_solver="ma57",max_cpu_time=900.,tol=1e-6))
-    @everywhere using Ipopt
+    @everywhere begin
+        using Ipopt
+
+        const ITER = [-1]
+        function ipopt_callback(
+            prob::IpoptProblem,alg_mod::Cint,iter_count::Cint,obj_value::Float64,
+            inf_pr::Float64,inf_du::Float64,mu::Float64,d_norm::Float64,
+            regularization_size::Float64,alpha_du::Float64,alpha_pr::Float64,ls_trials::Cint)
+
+            ITER[] += 1
+            return true
+        end
+
+        solver = pm -> begin
+            ITER[] = 0
+            set_optimizer(pm.model,()->
+                Ipopt.Optimizer(linear_solver="ma57",max_cpu_time=900.,tol=1e-6, print_level=PRINT_LEVEL))
+            MOI.set(pm.model, Ipopt.CallbackFunction(), ipopt_callback)
+            mem=@allocated begin
+                t=@elapsed begin
+                    optimize_model!(pm)
+                end
+            end
+            return get_status(termination_status(pm.model)),t,mem,ITER[]
+        end
+    end
 elseif SOLVER == "knitro"
     # TODO
 else
     error("Proper SOLVER should be given")
 end
 
 
-@everywhere function evalmodel(prob,solver)
-    println("Solving $prob")
-    t = @elapsed begin
-        retval = solver(prob)
-    end
-    retval["solve_time"] = t
-    return retval
-end
-
 function benchmark(solver,probs;warm_up_probs = [])
     println("Warming up (forcing JIT compile)")
-    println(warm_up_probs)
+    println(get_name.(warm_up_probs))
     rs = [remotecall.(solver,i,warm_up_probs) for i in procs() if i!= 1]
     ws = [wait.(r) for r in rs]
     fs= [fetch.(r) for r in rs]
-    
+
     println("Solving problems")
-    retvals = pmap(prob->evalmodel(prob,solver),probs)
-    time   = [retval["solve_time"] for retval in retvals]
-    status = [get_status(retval["termination_status"]) for retval in retvals]
-    time,status
+    retvals = pmap(prob->evalmodel(prob,solver;gcoff=GCOFF),probs)
+
+    status = [status for (status,time,mem,iter) in retvals]
+    time   = [time for (status,time,mem,iter) in retvals]
+    mem    = [mem for (status,time,mem,iter) in retvals]
+    iter   = [iter for (status,time,mem,iter) in retvals]
+
+    return status,time,mem,iter
 end
 
-cases = filter!(e->occursin("pglib_opf_case",e),readdir(PGLIB_PATH))
-types = [ACPPowerModel, ACRPowerModel, ACTPowerModel,
-         DCPPowerModel, DCMPPowerModel, NFAPowerModel,
-         DCPLLPowerModel,LPACCPowerModel, SOCWRPowerModel,
-         QCRMPowerModel,QCLSPowerModel]
-probs = [(case,type) for case in cases for type in types]
+if QUICK
+    cases = filter!(e->(occursin("pglib_opf_case",e) && occursin("pegase",e)),readdir(PGLIB_PATH))
+    types = [ACPPowerModel, ACRPowerModel]
+else
+    cases = filter!(e->occursin("pglib_opf_case",e),readdir(PGLIB_PATH))
+    types = [ACPPowerModel, ACRPowerModel, ACTPowerModel,
+             DCPPowerModel, DCMPPowerModel, NFAPowerModel,
+             DCPLLPowerModel,LPACCPowerModel, SOCWRPowerModel,
+             QCRMPowerModel,QCLSPowerModel]
+end
+probs = [instantiate_model(joinpath(PGLIB_PATH,case),type,PowerModels.build_opf) for case in cases for type in types]
 name =  ["$case-$type" for case in cases for type in types]
 
-time,status = benchmark(solver,probs;warm_up_probs = [("pglib_opf_case1888_rte.m", ACPPowerModel)])
+status,time,mem,iter = benchmark(solver,probs;warm_up_probs = [
+    instantiate_model(joinpath(PGLIB_PATH,"pglib_opf_case1888_rte.m"), ACPPowerModel,PowerModels.build_opf)
+])
 
 writedlm("name-power.csv",name,',')
+writedlm("status-power-$(SOLVER).csv",status)
 writedlm("time-power-$(SOLVER).csv",time)
-writedlm("status-power-$(SOLVER).csv",status) 
+writedlm("mem-power-$(SOLVER).csv",mem)
+writedlm("iter-power-$(SOLVER).csv",iter)

benchmark/config.jl

@@ -2,6 +2,9 @@ using Pkg, Distributed, DelimitedFiles
 
 const NP = ARGS[1]
 const SOLVER = ARGS[2]
+const VERBOSE = ARGS[3] == "true"
+const QUICK = ARGS[4] == "true"
+const GCOFF = ARGS[5] == "true"
 
 addprocs(parse(Int,NP),exeflags="--project=.")
 Pkg.instantiate()
@@ -19,3 +22,16 @@ elseif SOLVER == "knitro"
 else
     error("Proper ARGS should be given")
 end
+
+# Set verbose option
+if SOLVER == "ipopt"
+    const PRINT_LEVEL = VERBOSE ? 5 : 0
+elseif SOLVER == "knitro"
+    const PRINT_LEVEL = VERBOSE ? 3 : 0
+else
+    using MadNLP
+    const PRINT_LEVEL = VERBOSE ? MadNLP.INFO : MadNLP.ERROR
+end
+
+# Set quick option
+