ExaModels

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"

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"

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

examodels.jl

@@ -0,0 +1,282 @@
+#!/usr/bin/env julia
+###### AC-OPF using ExaModels ######
+#
+# implementation reference: https://exanauts.github.io/ExaModels.jl/stable/guide/
+# only the built-in AD library is supported
+#
+
+import PowerModels
+import ExaModels
+import NLPModelsIpopt
+import LinearAlgebra
+
+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], j = k, pd = pd, gs = gs, qd = qd, bs = bs)
+            end for (k, v) in ref[:bus]
+        ],
+        gen = [
+            (
+                i = gendict[k], j = 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_$(b.j)" => va_sol[b.i] for (i,b) in enumerate(data.bus))
+
+    vm_sol = ExaModels.solution(result, vm)
+    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_$(b.j)" => pg_sol[b.i] for (i,b) in enumerate(data.gen))
+
+    qg_sol = ExaModels.solution(result, qg)
+    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))
+
+    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))
+
+    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/opf_warmup.m")
+end

test/runtests.jl

@@ -4,6 +4,7 @@ include("validator.jl")
 
 @testset "Rosetta OPF" begin
     @testset "$framework" for framework in [
+        "examodels",
         "jump",
         "nlpmodels",
         "nonconvex",