implement the QP solver (#43)

* Setup TronQP

* implement qpsub

* update qpsub model and functions

* one level admm

* update solve_qpsub to Hessian and linearization input

* update models for SQP integration

* add second order correction

* indicator for SOC, FR, LF

* enable poststep

* fix poststep for SQP integration

* implement gpu branch kernel with first iteration test passed

* implement full ExaAdmm GPU tests

* implement multiplier extraction on GPU with tests

* add line_res for new SQOPF initialization

* Update Project.toml

---------

Co-authored-by: lukeli1990 <[email protected]>
Co-authored-by: Kibaek Kim <[email protected]>

.github/workflows/ci.yml

@@ -35,7 +35,6 @@ jobs:
     env:
       CUDA_VISIBLE_DEVICES: 1
       JULIA_DEPOT_PATH: /scratch/github-actions/julia_depot_exaadmmm
-      JULIA_CUDA_USE_BINARYBUILDER: true
     runs-on: self-hosted
     strategy:
       matrix:

.gitignore

@@ -1,3 +1,5 @@
 Manifest.toml
 docs/build
 docs/Manifest.toml
+.vscode/settings.json
+Manifest.toml

Artifacts.toml

@@ -0,0 +1,6 @@
+[ExaData]
+git-tree-sha1 = "df547b1e5d707fdb88721eee9bf84f9211e1fb50"
+lazy = true
+    [[ExaData.download]]
+    url = "https://web.cels.anl.gov/~mschanen/ExaData-e1b6123.tar.gz"
+    sha256 = "f6e22f26fbf273e97625a25b1b55f07d0eb0e15283483c545dda6eda84c48fdd"

Project.toml

@@ -1,7 +1,7 @@
 name = "ExaAdmm"
 uuid = "4d6a948c-1075-4240-a564-361a5d4e22a2"
-authors = ["Youngdae Kim <[email protected]>", "Kibaek Kim <[email protected]>", "Weiqi Zhang <[email protected]>", "François Pacaud <[email protected]>", "Michel Schanen <[email protected]>"]
-version = "0.4.0"
+authors = ["Youngdae Kim <[email protected]>", "Kibaek Kim <[email protected]>", "Weiqi Zhang <[email protected]>", "Bowen Li <[email protected]>", "François Pacaud <[email protected]>", "Michel Schanen <[email protected]>"]
+version = "0.4.1"
 
 [deps]
 AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
@@ -19,7 +19,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 AMDGPU = "0.4"
-CUDA = "4"
+CUDA = "4.1"
 ExaTron = "3"
 FileIO = "1.14"
 KernelAbstractions = "0.9"

src/ExaAdmm.jl

@@ -27,6 +27,7 @@ include("utils/utilities_gpu.jl")
 include("utils/utilities_ka.jl")
 
 include("algorithms/admm_two_level.jl")
+include("algorithms/admm_one_level.jl")
 
 # ----------------------------------------
 # A single period ACOPF implementation
@@ -158,4 +159,37 @@ include("models/mpec/mpec_admm_update_residual_gpu.jl")
 include("models/mpec/mpec_admm_update_lz_gpu.jl")
 include("models/mpec/mpec_admm_prepoststep_gpu.jl")
 =#
+
+# Interface to use ADMM solving QPsub.
+include("interface/solve_qpsub.jl")
+
+# Define "struct ModelAcopf" for encapsulating an ACOPF model.
+include("models/qpsub/qpsub_model.jl")
+
+# CPU
+include("models/qpsub/qpsub_init_solution_cpu.jl")
+include("models/qpsub/qpsub_generator_kernel_cpu.jl")
+include("models/qpsub/qpsub_admm_update_x_cpu.jl")
+include("models/qpsub/qpsub_admm_update_xbar_cpu.jl")
+include("models/qpsub/qpsub_admm_update_l_single_cpu.jl")
+include("models/qpsub/qpsub_admm_update_residual_cpu.jl")
+include("models/qpsub/qpsub_admm_prepoststep_cpu.jl")
+include("models/qpsub/qpsub_eval_Ab_linelimit_kernel_cpu.jl")
+include("models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_cpu.jl")
+include("models/qpsub/qpsub_auglag_tron_linelimit_kernel_cpu.jl")
+
+# GPU
+include("models/qpsub/qpsub_init_solution_gpu.jl")
+include("models/qpsub/qpsub_generator_kernel_gpu.jl")
+include("models/qpsub/qpsub_admm_update_x_gpu.jl")
+include("models/qpsub/qpsub_admm_update_xbar_gpu.jl")
+include("models/qpsub/qpsub_admm_update_l_single_gpu.jl")
+include("models/qpsub/qpsub_admm_update_residual_gpu.jl")
+include("models/qpsub/qpsub_eval_Ab_linelimit_kernel_gpu.jl")
+include("models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_gpu.jl")
+include("models/qpsub/qpsub_tron_linelimit_kernel.jl")
+include("models/qpsub/qpsub_admm_prepoststep_gpu.jl")
+
+
+
 end # module

src/algorithms/admm_one_level.jl

@@ -0,0 +1,80 @@
+function admm_one_level(
+    env::AdmmEnv, mod::AbstractOPFModel
+)
+    par = env.params
+    info = mod.info
+    sol = mod.solution
+
+    sqrt_d = sqrt(mod.nvar)
+    OUTER_TOL = sqrt_d*(par.outer_eps) #adjusted outer loop tolerance
+
+    fill!(info, 0)
+    info.mismatch = Inf
+
+    #eliminate second level
+    info.norm_z_prev = info.norm_z_curr = 0
+    par.initial_beta = 0
+    par.beta = 0
+    sol.lz .= 0
+    sol.z_curr .= 0
+    sol.z_prev .= 0
+    par.inner_iterlim = 1
+
+    if par.verbose > 0
+        admm_update_residual(env, mod)
+        @printf("%8s  %10s  %10s  %10s  %10s %10s %10s\n",
+                "Iter", "Objval", "Auglag", "PrimRes", "PrimTol", "DualRes", "DualTol")
+
+        @printf("%8d  %10.3e  %10.3e  %10.3e  %10.3e %10.3e  %10.3e\n",
+                info.outer, info.objval, info.auglag, info.mismatch, OUTER_TOL, info.dualres, OUTER_TOL*norm(sol.rho)/sqrt_d)
+    end
+
+    info.status = :IterationLimit
+
+    overall_time = @timed begin
+    while info.outer < par.outer_iterlim
+        admm_increment_outer(env, mod)
+
+
+        admm_increment_reset_inner(env, mod)
+        while info.inner < par.inner_iterlim
+            admm_increment_inner(env, mod)
+
+            admm_update_x(env, mod)
+
+            admm_update_xbar(env, mod)
+
+            admm_update_l_single(env, mod)
+
+            admm_update_residual(env, mod)
+
+            if par.verbose > 0
+                if (info.cumul % 50) == 0
+                    @printf("%8s  %10s  %10s  %10s  %10s  %10s  %10s\n",
+                            "Iter", "Objval", "Auglag", "PrimRes", "PrimTol", "DualRes", "DualTol")
+                end
+
+                @printf("%8d  %10.3e  %10.3e  %10.3e  %10.3e %10.3e  %10.3e\n",
+                        info.outer, info.objval, info.auglag, info.mismatch, OUTER_TOL, info.dualres, OUTER_TOL*norm(sol.rho)/sqrt_d)
+            end
+
+        end # while inner
+
+        # mismatch: x-xbar
+        if info.mismatch <= OUTER_TOL && info.dualres <= OUTER_TOL*norm(sol.rho)/sqrt_d
+            info.status = :Solved
+            break
+        end
+
+    end # while outer
+    end # @timed
+
+    info.time_overall = overall_time.time
+    admm_poststep(env, mod)
+
+    if par.verbose > 0
+        print_statistics(env, mod)
+    end
+
+    return
+end

src/algorithms/admm_two_level.jl

@@ -5,7 +5,7 @@ function admm_two_level(
     info = mod.info
 
     sqrt_d = sqrt(mod.nvar)
-    OUTER_TOL = sqrt_d*(par.outer_eps)
+    OUTER_TOL = sqrt_d*(par.outer_eps) #adjusted outer loop tolerance
 
     fill!(info, 0)
     info.mismatch = Inf
@@ -41,6 +41,7 @@ function admm_two_level(
             admm_update_l(env, mod, device)
             admm_update_residual(env, mod, device)
 
+            # an adjusting termination criteria for inner loop (i.e., inner loop is not solved to exact)
             info.eps_pri = sqrt_d / (2500*info.outer)
 
             if par.verbose > 0
@@ -55,18 +56,21 @@ function admm_two_level(
                         info.dualres, info.norm_z_curr, info.mismatch, OUTER_TOL, par.beta)
             end
 
+            # primres: x-xbar+z_curr
             if info.primres <= info.eps_pri #|| info.dualres <= par.DUAL_TOL
                 break
             end
         end # while inner
 
+        # mismatch: x-xbar
         if info.mismatch <= OUTER_TOL
             info.status = :Solved
             break
         end
 
         admm_update_lz(env, mod, device)
 
+        # if z_curr too large vs z_prev, increase penalty
         if info.norm_z_curr > par.theta*info.norm_z_prev
             par.beta = min(par.inc_c*par.beta, 1e24)
         end
@@ -81,4 +85,4 @@ function admm_two_level(
     end
 
     return
-end
+end

src/interface/solve_acopf.jl

@@ -37,6 +37,8 @@ function solve_acopf(case::String;
     env.params.outer_eps = outer_eps
     env.params.outer_iterlim = outer_iterlim
     env.params.inner_iterlim = inner_iterlim
+    #memory allocation for CuDynamicSharedArray
+    env.params.shmem_size = sizeof(Float64)*(14*mod.n+3*mod.n^2) + sizeof(Int)*(4*mod.n)
 
     admm_two_level(env, mod, isa(ka_device, KA.CPU) ? nothing : ka_device)
 

src/interface/solve_pf.jl

@@ -70,4 +70,4 @@ function solve_pf(pf::PowerFlow; tol=1e-6, max_iter=50)
     @printf("Elapsed time (secs). . . %8.2f\n", timed.time)
 
     return pf
-end
+end

src/interface/solve_qpsub.jl

@@ -0,0 +1,109 @@
+function solve_qpsub(
+    case::String,
+    Hs,
+    LH_1h,
+    RH_1h,
+    LH_1i,
+    RH_1i,
+    LH_1j,
+    RH_1j,
+    LH_1k,
+    RH_1k,
+    ls,
+    us,
+    pgmax,
+    pgmin,
+    qgmax,
+    qgmin,
+    c1,
+    c2,
+    Pd,
+    Qd,
+    initial_beta;
+    case_format = "matpower",
+    outer_iterlim = 20,
+    inner_iterlim = 1000,
+    rho_pq = 400.0,
+    rho_va = 40000.0,
+    obj_scale = 1.0,
+    scale = 1e-4,
+    storage_ratio = 0.0,
+    storage_charge_max = 1.0,
+    use_gpu = false,
+    use_linelimit = true,
+    use_projection = false,
+    tight_factor = 1.0,
+    outer_eps = 2 * 1e-4,
+    gpu_no = 0,
+    verbose = 1,
+    onelevel = true,
+)
+    T = Float64
+    TD = Array{Float64,1}
+    TI = Array{Int,1}
+    TM = Array{Float64,2}
+    if use_gpu
+        CUDA.device!(gpu_no)
+        TD = CuArray{Float64,1}
+        TI = CuArray{Int,1}
+        TM = CuArray{Float64,2}
+    end
+
+    env = AdmmEnv{T,TD,TI,TM}(
+        case,
+        rho_pq,
+        rho_va;
+        case_format = case_format,
+        use_gpu = use_gpu,
+        use_linelimit = use_linelimit,
+        use_projection = use_projection,
+        tight_factor = tight_factor,
+        gpu_no = gpu_no,
+        storage_ratio = storage_ratio,
+        storage_charge_max = storage_charge_max,
+        verbose = verbose,
+    )
+    mod = ModelQpsub{T,TD,TI,TM}(env)
+
+    data = mod.grid_data
+
+    mod.Hs = copy(Hs)
+    mod.LH_1h = copy(LH_1h)
+    mod.RH_1h = copy(RH_1h)
+    mod.LH_1i = copy(LH_1i)
+    mod.RH_1i = copy(RH_1i)
+    mod.LH_1j = copy(LH_1j)
+    mod.RH_1j = copy(RH_1j)
+    mod.LH_1k = copy(LH_1k)
+    mod.RH_1k = copy(RH_1k)
+    mod.ls = copy(ls)
+    mod.us = copy(us)
+
+    mod.qpsub_pgmax = copy(pgmax)
+    mod.qpsub_pgmin = copy(pgmin)
+    mod.qpsub_qgmax = copy(qgmax)
+    mod.qpsub_qgmin = copy(qgmin)
+
+    mod.qpsub_c1 = copy(c1)
+    mod.qpsub_c2 = copy(c2)
+    mod.qpsub_Pd = copy(Pd)
+    mod.qpsub_Qd = copy(Qd)
+
+    env.params.scale = scale
+    env.params.obj_scale = obj_scale
+    env.params.outer_eps = outer_eps
+    env.params.outer_iterlim = outer_iterlim
+    env.params.inner_iterlim = inner_iterlim
+    env.params.shmem_size = sizeof(Float64) * (16 * mod.n + 4 * mod.n^2 + 178) + sizeof(Int) * (4 * mod.n)
+    env.params.initial_beta = initial_beta #use my initial beta 
+
+    init_solution!(mod, mod.solution, env.initial_rho_pq, env.initial_rho_va)
+
+    #one-level or two-level admm 
+    if onelevel
+        admm_one_level(env, mod)
+    else
+        @warn "two-level ADMM is not implemented in QPsub"
+    end
+    return env, mod
+end

src/models/qpsub/qpsub_admm_increment.jl

@@ -0,0 +1,44 @@
+"""
+    admm_increment_outer()
+
+increment outer iteration
+"""
+
+function admm_increment_outer(
+    env::AdmmEnv,
+    mod::ModelQpsub
+)
+    mod.info.outer += 1
+    return
+end
+
+
+"""
+    admm_increment_reset_inner()
+
+reset inner iteration when outer iteration just incremented
+"""
+
+function admm_increment_reset_inner(
+    env::AdmmEnv,
+    mod::ModelQpsub
+)
+    mod.info.inner = 0
+    return
+end
+
+
+"""
+    admm_increment_inner()
+
+increment inner iteration
+"""
+
+function admm_increment_inner(
+    env::AdmmEnv,
+    mod::ModelQpsub
+)
+    mod.info.inner += 1
+    mod.info.cumul += 1
+    return
+end