From 02fe75116b03fcb3e40c73b031ea90c1101b2975 Mon Sep 17 00:00:00 2001 From: Sungho Shin Date: Fri, 26 Jan 2024 20:26:40 -0600 Subject: [PATCH 1/4] Added ExaModels.jl (#57) * Added ExaModels.jl example * Update project files * added examodels to testing * [hotfix] "ma27" option removed * "typo fix README" --------- Co-authored-by: Oscar Dowson --- Manifest.toml | 22 +++- Project.toml | 3 + README.md | 2 + examodels.jl | 285 +++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 5 files changed, 312 insertions(+), 1 deletion(-) create mode 100755 examodels.jl diff --git a/Manifest.toml b/Manifest.toml index 4ccaad9..7efda40 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.0" manifest_format = "2.0" -project_hash = "f82d4f25fcfe31a45d52dae8a3020a2e2d59cb05" +project_hash = "4dde0c1ed630e05fd3041df844c0a01316b0a2cd" [[deps.ADNLPModels]] deps = ["ColPack", "ForwardDiff", "LinearAlgebra", "NLPModels", "Requires", "ReverseDiff", "SparseArrays"] @@ -449,6 +449,26 @@ git-tree-sha1 = "d6863c556f1142a061532e79f611aa46be201686" uuid = "90fa49ef-747e-5e6f-a989-263ba693cf1a" version = "0.5.2" +[[deps.ExaModels]] +deps = ["NLPModels", "SolverCore"] +git-tree-sha1 = "054b99d8f8d19a81a1fbbad74223af16f8012f2a" +uuid = "1037b233-b668-4ce9-9b63-f9f681f55dd2" +version = "0.4.2" + + [deps.ExaModels.extensions] + ExaModelsAMDGPU = "AMDGPU" + ExaModelsCUDA = "CUDA" + ExaModelsKernelAbstractions = "KernelAbstractions" + ExaModelsOneAPI = "oneAPI" + ExaModelsSpecialFunctions = "SpecialFunctions" + + [deps.ExaModels.weakdeps] + AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" + CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" + SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" + oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b" + [[deps.ExponentialUtilities]] deps = ["Adapt", "ArrayInterface", "GPUArraysCore", "GenericSchur", "LinearAlgebra", "PrecompileTools", "Printf", "SparseArrays", "libblastrampoline_jll"] git-tree-sha1 = "602e4585bcbd5a25bc06f514724593d13ff9e862" diff --git a/Project.toml b/Project.toml index c374a3f..4bcf335 100644 --- a/Project.toml +++ b/Project.toml @@ -1,5 +1,6 @@ [deps] ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a" +ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" @@ -15,8 +16,10 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] ADNLPModels = "0.7" +ExaModels = "0.4" Ipopt = "1" JuMP = "1.15" +NLPModelsIpopt = "0.10" Nonconvex = "2" NonconvexIpopt = "0.4" Optim = "1" diff --git a/README.md b/README.md index 58759bf..b0a7b3c 100644 --- a/README.md +++ b/README.md @@ -50,12 +50,14 @@ In this formulation the complex voltage terms expand into the following expressi ### AC-OPF Implementations Each of these files is designed to be _stand-alone_ and can be tested with minimal package dependencies. Consequently, there is some code replication between implementations. +- `examodels.jl`: implementation using [ExaModels](https://github.com/exanauts/ExaModels.jl) - `jump.jl`: implementation using [JuMP](https://github.com/jump-dev/JuMP.jl) - `nlpmodels.jl`: implementation using [ADNLPModels](https://github.com/JuliaSmoothOptimizers/ADNLPModels.jl) - `nonconvex.jl`: implementation using [Nonconvex](https://github.com/JuliaNonconvex/Nonconvex.jl) - `optim.jl`: implementation using [Optim](https://github.com/JuliaNLSolvers/Optim.jl) - `optimization.jl`: implementation using [Optimization](https://github.com/SciML/Optimization.jl) + ### Other Files - `data/*`: small example datasets - `test/*`: basic tests for the primary AC-OPF implementations diff --git a/examodels.jl b/examodels.jl new file mode 100755 index 0000000..7d12146 --- /dev/null +++ b/examodels.jl @@ -0,0 +1,285 @@ +#!/usr/bin/env julia +###### AC-OPF using ExaModels ###### +# +# implementation reference: https://github.com/lanl-ansi/PowerModelsAnnex.jl/blob/master/src/model/ac-opf.jl +# only the built-in AD library is supported +# + +import PowerModels +import ExaModels +import NLPModelsIpopt +import LinearAlgebra + +const constraint_tol = 1e-6 + +function solve_opf(file_name) + time_data_start = time() + + data = PowerModels.parse_file(file_name) + PowerModels.standardize_cost_terms!(data, order=2) + PowerModels.calc_thermal_limits!(data) + ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] + + arcdict = Dict(a => k for (k, a) in enumerate(ref[:arcs])) + busdict = Dict(k => i for (i, (k, v)) in enumerate(ref[:bus])) + gendict = Dict(k => i for (i, (k, v)) in enumerate(ref[:gen])) + branchdict = Dict(k => i for (i, (k, v)) in enumerate(ref[:branch])) + + data = ( + bus = [ + begin + bus_loads = [ref[:load][l] for l in ref[:bus_loads][k]] + bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][k]] + pd = sum(load["pd"] for load in bus_loads; init = 0.0) + gs = sum(shunt["gs"] for shunt in bus_shunts; init = 0.0) + qd = sum(load["qd"] for load in bus_loads; init = 0.0) + bs = sum(shunt["bs"] for shunt in bus_shunts; init = 0.0) + (i = busdict[k], pd = pd, gs = gs, qd = qd, bs = bs) + end for (k, v) in ref[:bus] + ], + gen = [ + ( + i = gendict[k], + cost1 = v["cost"][1], + cost2 = v["cost"][2], + cost3 = v["cost"][3], + bus = busdict[v["gen_bus"]], + ) for (k, v) in ref[:gen] + ], + arc = [ + (i = k, rate_a = ref[:branch][l]["rate_a"], bus = busdict[i]) for + (k, (l, i, j)) in enumerate(ref[:arcs]) + ], + branch = [ + begin + f_idx = arcdict[i, branch["f_bus"], branch["t_bus"]] + t_idx = arcdict[i, branch["t_bus"], branch["f_bus"]] + g, b = PowerModels.calc_branch_y(branch) + tr, ti = PowerModels.calc_branch_t(branch) + ttm = tr^2 + ti^2 + g_fr = branch["g_fr"] + b_fr = branch["b_fr"] + g_to = branch["g_to"] + b_to = branch["b_to"] + c1 = (-g * tr - b * ti) / ttm + c2 = (-b * tr + g * ti) / ttm + c3 = (-g * tr + b * ti) / ttm + c4 = (-b * tr - g * ti) / ttm + c5 = (g + g_fr) / ttm + c6 = (b + b_fr) / ttm + c7 = (g + g_to) + c8 = (b + b_to) + ( + i = branchdict[i], + j = 1, + f_idx = f_idx, + t_idx = t_idx, + f_bus = busdict[branch["f_bus"]], + t_bus = busdict[branch["t_bus"]], + c1 = c1, + c2 = c2, + c3 = c3, + c4 = c4, + c5 = c5, + c6 = c6, + c7 = c7, + c8 = c8, + rate_a_sq = branch["rate_a"]^2, + ) + end for (i, branch) in ref[:branch] + ], + ref_buses = [busdict[i] for (i, k) in ref[:ref_buses]], + vmax = [v["vmax"] for (k, v) in ref[:bus]], + vmin = [v["vmin"] for (k, v) in ref[:bus]], + pmax = [v["pmax"] for (k, v) in ref[:gen]], + pmin = [v["pmin"] for (k, v) in ref[:gen]], + qmax = [v["qmax"] for (k, v) in ref[:gen]], + qmin = [v["qmin"] for (k, v) in ref[:gen]], + rate_a = [ref[:branch][l]["rate_a"] for (k, (l, i, j)) in enumerate(ref[:arcs])], + angmax = [b["angmax"] for (i, b) in ref[:branch]], + angmin = [b["angmin"] for (i, b) in ref[:branch]], + ) + + data_load_time = time() - time_data_start + + + time_model_start = time() + + w = ExaModels.ExaCore() + + va = ExaModels.variable(w, length(data.bus);) + + vm = ExaModels.variable( + w, + length(data.bus); + start = fill!(similar(data.bus, Float64), 1.0), + lvar = data.vmin, + uvar = data.vmax, + ) + pg = ExaModels.variable(w, length(data.gen); lvar = data.pmin, uvar = data.pmax) + + qg = ExaModels.variable(w, length(data.gen); lvar = data.qmin, uvar = data.qmax) + + p = ExaModels.variable(w, length(data.arc); lvar = -data.rate_a, uvar = data.rate_a) + + q = ExaModels.variable(w, length(data.arc); lvar = -data.rate_a, uvar = data.rate_a) + + o = ExaModels.objective( + w, + g.cost1 * pg[g.i]^2 + g.cost2 * pg[g.i] + g.cost3 for g in data.gen + ) + + c1 = ExaModels.constraint(w, va[i] for i in data.ref_buses) + + c2 = ExaModels.constraint( + w, + p[b.f_idx] - b.c5 * vm[b.f_bus]^2 - + b.c3 * (vm[b.f_bus] * vm[b.t_bus] * cos(va[b.f_bus] - va[b.t_bus])) - + b.c4 * (vm[b.f_bus] * vm[b.t_bus] * sin(va[b.f_bus] - va[b.t_bus])) for + b in data.branch + ) + + c3 = ExaModels.constraint( + w, + q[b.f_idx] + + b.c6 * vm[b.f_bus]^2 + + b.c4 * (vm[b.f_bus] * vm[b.t_bus] * cos(va[b.f_bus] - va[b.t_bus])) - + b.c3 * (vm[b.f_bus] * vm[b.t_bus] * sin(va[b.f_bus] - va[b.t_bus])) for + b in data.branch + ) + + c4 = ExaModels.constraint( + w, + p[b.t_idx] - b.c7 * vm[b.t_bus]^2 - + b.c1 * (vm[b.t_bus] * vm[b.f_bus] * cos(va[b.t_bus] - va[b.f_bus])) - + b.c2 * (vm[b.t_bus] * vm[b.f_bus] * sin(va[b.t_bus] - va[b.f_bus])) for + b in data.branch + ) + + c5 = ExaModels.constraint( + w, + q[b.t_idx] + + b.c8 * vm[b.t_bus]^2 + + b.c2 * (vm[b.t_bus] * vm[b.f_bus] * cos(va[b.t_bus] - va[b.f_bus])) - + b.c1 * (vm[b.t_bus] * vm[b.f_bus] * sin(va[b.t_bus] - va[b.f_bus])) for + b in data.branch + ) + + c6 = ExaModels.constraint( + w, + va[b.f_bus] - va[b.t_bus] for b in data.branch; + lcon = data.angmin, + ucon = data.angmax, + ) + c7 = ExaModels.constraint( + w, + p[b.f_idx]^2 + q[b.f_idx]^2 - b.rate_a_sq for b in data.branch; + lcon = fill!(similar(data.branch, Float64, length(data.branch)), -Inf), + ) + c8 = ExaModels.constraint( + w, + p[b.t_idx]^2 + q[b.t_idx]^2 - b.rate_a_sq for b in data.branch; + lcon = fill!(similar(data.branch, Float64, length(data.branch)), -Inf), + ) + + c9 = ExaModels.constraint(w, b.pd + b.gs * vm[b.i]^2 for b in data.bus) + + c10 = ExaModels.constraint(w, b.qd - b.bs * vm[b.i]^2 for b in data.bus) + + c11 = ExaModels.constraint!(w, c9, a.bus => p[a.i] for a in data.arc) + c12 = ExaModels.constraint!(w, c10, a.bus => q[a.i] for a in data.arc) + + c13 = ExaModels.constraint!(w, c9, g.bus => -pg[g.i] for g in data.gen) + c14 = ExaModels.constraint!(w, c10, g.bus => -qg[g.i] for g in data.gen) + + model = ExaModels.ExaModel(w) + + model_build_time = time() - time_model_start + + + time_solve_start = time() + + result = NLPModelsIpopt.ipopt(model) + + cost = result.objective + + feasible = result.status == :first_order + + solve_time = time() - time_solve_start + total_time = time() - time_data_start + + # TODO: once a built-in timer is implemented in ExaModels, report callback timing + + # total_callback_time = + # nlp_block.evaluator.eval_objective_timer + + # nlp_block.evaluator.eval_objective_gradient_timer + + # nlp_block.evaluator.eval_constraint_timer + + # nlp_block.evaluator.eval_constraint_jacobian_timer + + # nlp_block.evaluator.eval_hessian_lagrangian_timer + + println("") + println("\033[1mSummary\033[0m") + println(" case........: $(file_name)") + println(" variables...: $(model.meta.nvar)") + println(" constraints.: $(model.meta.ncon)") + println(" feasible....: $(feasible)") + println(" cost........: $(round(Int, cost))") + println(" total time..: $(total_time)") + println(" data time.: $(data_load_time)") + println(" build time: $(model_build_time)") + println(" solve time: $(solve_time)") + # println(" callbacks: $(total_callback_time)") + println("") + # println(" callbacks time:") + # println(" * obj.....: $(nlp_block.evaluator.eval_objective_timer)") + # println(" * grad....: $(nlp_block.evaluator.eval_objective_gradient_timer)") + # println(" * cons....: $(nlp_block.evaluator.eval_constraint_timer)") + # println(" * jac.....: $(nlp_block.evaluator.eval_constraint_jacobian_timer)") + # println(" * hesslag.: $(nlp_block.evaluator.eval_hessian_lagrangian_timer)") + # println("") + + va_sol = ExaModels.solution(result, va) + va_dict = Dict("va_$(i)" => va_sol[busdict[i]] for (i,b) in enumerate(data.bus)) + + vm_sol = ExaModels.solution(result, vm) + vm_dict = Dict("vm_$(i)" => vm_sol[busdict[i]] for (i,b) in enumerate(data.bus)) + + pg_sol = ExaModels.solution(result, pg) + pg_dict = Dict("pg_$(i)" => pg_sol[gendict[i]] for (i,b) in enumerate(data.gen)) + + qg_sol = ExaModels.solution(result, qg) + qg_dict = Dict("qg_$(i)" => qg_sol[gendict[i]] for (i,b) in enumerate(data.gen)) + + p_sol = ExaModels.solution(result, p) + p_dict = Dict("p_$(write_out_tuple(ref[:arcs][i]))" => p_sol[i] for (i,b) in enumerate(data.arc)) + + q_sol = ExaModels.solution(result, q) + q_dict = Dict("q_$(write_out_tuple(ref[:arcs][i]))" => q_sol[i] for (i,b) in enumerate(data.arc)) + + println(ref[:arcs]) + return Dict( + "case" => file_name, + "variables" => model.meta.nvar, + "constraints" => model.meta.ncon, + "feasible" => feasible, + "cost" => cost, + "time_total" => total_time, + "time_data" => data_load_time, + "time_build" => model_build_time, + "time_solve" => solve_time, + # "time_callbacks" => total_callback_time, + "solution" => Dict( + va_dict..., + vm_dict..., + pg_dict..., + qg_dict..., + p_dict..., + q_dict...), + ) +end + +write_out_tuple((i,j,k)) = "$(i)_$(j)_$(k)" + +if isinteractive() == false + solve_opf("$(@__DIR__)/data/pglib_opf_case5_pjm.m") +end diff --git a/test/runtests.jl b/test/runtests.jl index d2ffaf8..2d6e800 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,7 @@ include("validator.jl") @testset "Rosetta OPF" begin @testset "$framework" for framework in [ + "examodels", "jump", "nlpmodels", "nonconvex", From c685c2cf60c1ee49fd1a686ef5a3f7b7edeaa794 Mon Sep 17 00:00:00 2001 From: Carleton Coffrin Date: Fri, 26 Jan 2024 20:01:19 -0700 Subject: [PATCH 2/4] remove extra println, update to opf_warmup case --- examodels.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/examodels.jl b/examodels.jl index 7d12146..163b596 100755 --- a/examodels.jl +++ b/examodels.jl @@ -1,7 +1,7 @@ #!/usr/bin/env julia ###### AC-OPF using ExaModels ###### # -# implementation reference: https://github.com/lanl-ansi/PowerModelsAnnex.jl/blob/master/src/model/ac-opf.jl +# implementation reference: https://exanauts.github.io/ExaModels.jl/stable/guide/ # only the built-in AD library is supported # @@ -10,8 +10,6 @@ import ExaModels import NLPModelsIpopt import LinearAlgebra -const constraint_tol = 1e-6 - function solve_opf(file_name) time_data_start = time() @@ -200,7 +198,7 @@ function solve_opf(file_name) time_solve_start = time() result = NLPModelsIpopt.ipopt(model) - + cost = result.objective feasible = result.status == :first_order @@ -256,7 +254,6 @@ function solve_opf(file_name) q_sol = ExaModels.solution(result, q) q_dict = Dict("q_$(write_out_tuple(ref[:arcs][i]))" => q_sol[i] for (i,b) in enumerate(data.arc)) - println(ref[:arcs]) return Dict( "case" => file_name, "variables" => model.meta.nvar, @@ -281,5 +278,5 @@ end write_out_tuple((i,j,k)) = "$(i)_$(j)_$(k)" if isinteractive() == false - solve_opf("$(@__DIR__)/data/pglib_opf_case5_pjm.m") + solve_opf("$(@__DIR__)/data/opf_warmup.m") end From 13a41828aff757167e9c9c032da458faa81ae71c Mon Sep 17 00:00:00 2001 From: Carleton Coffrin Date: Fri, 26 Jan 2024 21:17:42 -0700 Subject: [PATCH 3/4] fix solution building --- examodels.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examodels.jl b/examodels.jl index 163b596..85747db 100755 --- a/examodels.jl +++ b/examodels.jl @@ -237,16 +237,16 @@ function solve_opf(file_name) # println("") va_sol = ExaModels.solution(result, va) - va_dict = Dict("va_$(i)" => va_sol[busdict[i]] for (i,b) in enumerate(data.bus)) + va_dict = Dict("va_$(i)" => va_sol[b.i] for (i,b) in enumerate(data.bus)) vm_sol = ExaModels.solution(result, vm) - vm_dict = Dict("vm_$(i)" => vm_sol[busdict[i]] for (i,b) in enumerate(data.bus)) + vm_dict = Dict("vm_$(i)" => vm_sol[b.i] for (i,b) in enumerate(data.bus)) pg_sol = ExaModels.solution(result, pg) - pg_dict = Dict("pg_$(i)" => pg_sol[gendict[i]] for (i,b) in enumerate(data.gen)) + pg_dict = Dict("pg_$(i)" => pg_sol[b.i] for (i,b) in enumerate(data.gen)) qg_sol = ExaModels.solution(result, qg) - qg_dict = Dict("qg_$(i)" => qg_sol[gendict[i]] for (i,b) in enumerate(data.gen)) + qg_dict = Dict("qg_$(i)" => qg_sol[b.i] for (i,b) in enumerate(data.gen)) p_sol = ExaModels.solution(result, p) p_dict = Dict("p_$(write_out_tuple(ref[:arcs][i]))" => p_sol[i] for (i,b) in enumerate(data.arc)) From 846db98550322967b1adac91f7a35a1f8a0b7390 Mon Sep 17 00:00:00 2001 From: Sungho Shin Date: Fri, 26 Jan 2024 23:07:41 -0600 Subject: [PATCH 4/4] ci failure fix (#59) --- examodels.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examodels.jl b/examodels.jl index 85747db..42ca209 100755 --- a/examodels.jl +++ b/examodels.jl @@ -32,12 +32,12 @@ function solve_opf(file_name) gs = sum(shunt["gs"] for shunt in bus_shunts; init = 0.0) qd = sum(load["qd"] for load in bus_loads; init = 0.0) bs = sum(shunt["bs"] for shunt in bus_shunts; init = 0.0) - (i = busdict[k], pd = pd, gs = gs, qd = qd, bs = bs) + (i = busdict[k], j = k, pd = pd, gs = gs, qd = qd, bs = bs) end for (k, v) in ref[:bus] ], gen = [ ( - i = gendict[k], + i = gendict[k], j = k, cost1 = v["cost"][1], cost2 = v["cost"][2], cost3 = v["cost"][3], @@ -237,16 +237,16 @@ function solve_opf(file_name) # println("") va_sol = ExaModels.solution(result, va) - va_dict = Dict("va_$(i)" => va_sol[b.i] for (i,b) in enumerate(data.bus)) + va_dict = Dict("va_$(b.j)" => va_sol[b.i] for (i,b) in enumerate(data.bus)) vm_sol = ExaModels.solution(result, vm) - vm_dict = Dict("vm_$(i)" => vm_sol[b.i] for (i,b) in enumerate(data.bus)) + vm_dict = Dict("vm_$(b.j)" => vm_sol[b.i] for (i,b) in enumerate(data.bus)) pg_sol = ExaModels.solution(result, pg) - pg_dict = Dict("pg_$(i)" => pg_sol[b.i] for (i,b) in enumerate(data.gen)) + pg_dict = Dict("pg_$(b.j)" => pg_sol[b.i] for (i,b) in enumerate(data.gen)) qg_sol = ExaModels.solution(result, qg) - qg_dict = Dict("qg_$(i)" => qg_sol[b.i] for (i,b) in enumerate(data.gen)) + qg_dict = Dict("qg_$(b.j)" => qg_sol[b.i] for (i,b) in enumerate(data.gen)) p_sol = ExaModels.solution(result, p) p_dict = Dict("p_$(write_out_tuple(ref[:arcs][i]))" => p_sol[i] for (i,b) in enumerate(data.arc))