Normalized residuals and tolerances

src/algorithms/admm_one_level.jl

@@ -6,8 +6,8 @@ function admm_one_level(
     info = mod.info
     sol = mod.solution
 
-    sqrt_d = sqrt(mod.nvar)
-    OUTER_TOL = sqrt_d*(par.outer_eps) #adjusted outer loop tolerance
+    info.primtol = par.RELTOL
+    info.dualtol = par.RELTOL
 
     fill!(info, 0)
     info.mismatch = Inf
@@ -23,11 +23,11 @@ function admm_one_level(
 
     if par.verbose > 0
         admm_update_residual(env, mod, device)
-        @printf("%8s  %10s  %10s  %10s  %10s %10s %10s\n",
-                "Iter", "Objval", "Auglag", "PrimRes", "PrimTol", "DualRes", "DualTol")
+        @printf("%8s  %10s  %10s  %10s  %10s  %10s  %10s %10s %10s %10s\n",
+                "Iter", "Objval", "Auglag", "PrimScal", "PrimRes", "PrimTol", "DualScal", "DualRes", "DualTol", "rho")
 
-        @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)
+        @printf("%8d  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e %10.3e  %10.3e  %10.3e\n",
+                info.outer, info.objval, info.auglag, info.primsca, info.primres, info.primtol, info.dualsca, info.dualres, info.dualtol, sol.rho[1])
     end
 
     info.status = :IterationLimit
@@ -51,22 +51,31 @@ function admm_one_level(
 
             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")
+                    @printf("%8s  %10s  %10s  %10s  %10s  %10s  %10s %10s %10s %10s\n",
+                            "Iter", "Objval", "Auglag", "PrimScal", "PrimRes", "PrimTol", "DualScal", "DualRes", "DualTol", "rho")
                 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)
+                @printf("%8d  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e %10.3e  %10.3e  %10.3e\n",
+                info.outer, info.objval, info.auglag, info.primsca, info.primres, info.primtol, info.dualsca, info.dualres, info.dualtol, sol.rho[1])
             end
 
         end # while inner
 
         # mismatch: x-xbar
-        if info.mismatch <= OUTER_TOL && info.dualres <= OUTER_TOL*norm(sol.rho, device)/sqrt_d
+        if info.primres <= info.primtol # && info.dualres <= info.dualtol
             info.status = :Solved
             break
         end
 
+        # residual balancing
+        if par.rb_switch
+            if info.primres > par.rb_beta1 * info.dualres
+                sol.rho .= min.(sol.rho * par.rb_tau, 1e+16)
+            elseif par.rb_beta2 * info.primres < info.dualres
+                sol.rho .= max.(sol.rho / par.rb_tau, 1e-16)
+            end
+        end
+
     end # while outer
     end # @timed
 

src/algorithms/admm_two_level.jl

@@ -42,7 +42,7 @@ function admm_two_level(
             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)
+            info.eps_pri = min(sqrt_d / (2500*info.outer) / info.primsca, info.primres)
 
             if par.verbose > 0
                 if (info.cumul % 50) == 0

src/interface/solve_qpsub.jl

@@ -34,7 +34,8 @@ function solve_qpsub(
     use_linelimit = true,
     use_projection = false,
     tight_factor = 1.0,
-    outer_eps = 2 * 1e-4,
+    ABSTOL = 2 * 1e-4,
+    RELTOL = 1e-4,
     gpu_no = 0,
     verbose = 1,
     onelevel = true,
@@ -102,7 +103,8 @@ function solve_qpsub(
 
     env.params.scale = scale
     env.params.obj_scale = obj_scale
-    env.params.outer_eps = outer_eps
+    env.params.ABSTOL = ABSTOL
+    env.params.RELTOL = RELTOL
     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)

src/models/acopf/acopf_admm_update_residual_cpu.jl

@@ -13,16 +13,27 @@ The primal and dual residuals are stored in `mod.solution.rp` and `mod.solution.
 function admm_update_residual(
     env::AdmmEnv{Float64,Array{Float64,1},Array{Int,1},Array{Float64,2}},
     mod::AbstractOPFModel{Float64,Array{Float64,1},Array{Int,1},Array{Float64,2}},
-    device::Nothing=nothing
+    device::Nothing=nothing;
+    normalized=true
 )
     sol, info, data, par, grid_data = mod.solution, mod.info, env.data, env.params, mod.grid_data
 
     sol.rp .= sol.u_curr .- sol.v_curr .+ sol.z_curr
     sol.rd .= sol.z_curr .- sol.z_prev
     sol.Ax_plus_By .= sol.rp .- sol.z_curr
 
+    info.primsca = max(norm(sol.u_curr), norm(sol.v_curr), norm(sol.z_curr))
+    info.dualsca = norm(sol.l_curr)
     info.primres = norm(sol.rp)
     info.dualres = norm(sol.rd)
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.primsca
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.dualsca
+    if normalized
+        info.primres /= info.primsca
+        info.dualres /= info.dualsca
+        info.primtol /= info.primsca
+        info.dualtol /= info.dualsca
+    end
     info.norm_z_curr = norm(sol.z_curr)
     info.mismatch = norm(sol.Ax_plus_By)
     info.objval = sum(data.generators[g].coeff[data.generators[g].n-2]*(grid_data.baseMVA*sol.u_curr[mod.gen_start+2*(g-1)])^2 +

src/models/acopf/acopf_admm_update_residual_gpu.jl

@@ -11,17 +11,28 @@ end
 function admm_update_residual(
     env::AdmmEnv{Float64,CuArray{Float64,1},CuArray{Int,1},CuArray{Float64,2}},
     mod::AbstractOPFModel{Float64,CuArray{Float64,1},CuArray{Int,1},CuArray{Float64,2}},
-    device::Nothing=nothing
+    device::Nothing=nothing;
+    normalized=true
 )
-    sol, info = mod.solution, mod.info
+    sol, info, par = mod.solution, mod.info, env.params
 
     @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) compute_primal_residual_kernel(mod.nvar, sol.rp, sol.u_curr, sol.v_curr, sol.z_curr)
     @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) vector_difference(mod.nvar, sol.rd, sol.z_curr, sol.z_prev)
     CUDA.synchronize()
     CUDA.@sync @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) vector_difference(mod.nvar, sol.Ax_plus_By, sol.rp, sol.z_curr)
 
+    info.primsca = max(CUDA.norm(sol.u_curr), CUDA.norm(sol.v_curr), CUDA.norm(sol.z_curr))
+    info.dualsca = CUDA.norm(sol.l_curr)
     info.primres = CUDA.norm(sol.rp)
     info.dualres = CUDA.norm(sol.rd)
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.primsca
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.dualsca
+    if normalized
+        info.primres /= info.primsca
+        info.dualres /= info.dualsca
+        info.primtol /= info.primsca
+        info.dualtol /= info.dualsca
+    end
     info.norm_z_curr = CUDA.norm(sol.z_curr)
     info.mismatch = CUDA.norm(sol.Ax_plus_By)
 

src/models/acopf/acopf_admm_update_residual_ka.jl

@@ -13,9 +13,10 @@ end
 function admm_update_residual(
     env::AdmmEnv,
     mod::AbstractOPFModel,
-    device
+    device;
+    normalized=true
 )
-    sol, info = mod.solution, mod.info
+    sol, info, par = mod.solution, mod.info, env.params
     nblk_nvar = div(mod.nvar-1, 64)+1
     compute_primal_residual_kernel_ka(device,64,64*nblk_nvar)(
         mod.nvar, sol.rp, sol.u_curr, sol.v_curr, sol.z_curr
@@ -30,8 +31,18 @@ function admm_update_residual(
     )
     KA.synchronize(device)
 
+    info.primsca = max(norm(sol.u_curr, device), norm(sol.v_curr, device), norm(sol.z_curr, device))
+    info.dualsca = norm(sol.l_curr, device)
     info.primres = norm(sol.rp, device)
     info.dualres = norm(sol.rd, device)
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.primsca
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL + par.RELTOL * info.dualsca
+    if normalized
+        info.primres /= info.primsca
+        info.dualres /= info.dualsca
+        info.primtol /= info.primsca
+        info.dualtol /= info.dualsca
+    end
     info.norm_z_curr = norm(sol.z_curr, device)
     info.mismatch = norm(sol.Ax_plus_By, device)
 

src/models/mpacopf/mpacopf_admm_update_residual_cpu.jl

@@ -3,15 +3,19 @@ function admm_update_residual(
     mod::ModelMpacopf{Float64,Array{Float64,1},Array{Int,1},Array{Float64,2}},
     device::Nothing=nothing
 )
-    info = mod.info
+    info, par = mod.info, env.params
 
     info.primres = 0.0
     info.dualres = 0.0
+    info.primsca = 0.0
+    info.dualsca = 0.0
     info.norm_z_curr = 0.0
     info.mismatch = 0.0
 
     for i=1:mod.len_horizon
-        admm_update_residual(env, mod.models[i])
+        admm_update_residual(env, mod.models[i], normalized = false)
+        info.primsca = max(info.primsca, mod.models[i].info.primsca)
+        info.dualsca = max(info.dualsca, mod.models[i].info.dualsca)
     #=
         info.primres += mod.models[i].info.primres^2
         info.dualres += mod.models[i].info.dualres^2
@@ -31,6 +35,8 @@ function admm_update_residual(
         mod.models[i].info.dualres = sqrt(mod.models[i].info.dualres^2 + norm(sol_ramp.rd)^2)
         mod.models[i].info.norm_z_curr = sqrt(mod.models[i].info.norm_z_curr^2 + norm(sol_ramp.z_curr)^2)
         mod.models[i].info.mismatch = sqrt(mod.models[i].info.mismatch^2 + norm(sol_ramp.Ax_plus_By)^2)
+        info.primsca = max(info.primsca, norm(sol_ramp.u_curr), norm(v_curr), norm(sol_ramp.z_curr))
+        info.dualsca = max(info.dualsca, norm(sol_ramp.l_curr))
     #=
         info.primres += norm(sol_ramp.rp)^2
         info.dualres += norm(sol_ramp.rd)^2
@@ -45,6 +51,10 @@ function admm_update_residual(
         info.norm_z_curr = max(info.norm_z_curr, mod.models[i].info.norm_z_curr)
         info.mismatch = max(info.mismatch, mod.models[i].info.mismatch)
     end
+    info.primres /= info.primsca
+    info.dualres /= info.dualsca
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL / info.primsca + par.RELTOL
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL / info.dualsca + par.RELTOL
 
     #=
     info.primres = sqrt(info.primres)

src/models/mpacopf/mpacopf_admm_update_residual_gpu.jl

@@ -16,15 +16,19 @@ function admm_update_residual(
     mod::ModelMpacopf{Float64,CuArray{Float64,1},CuArray{Int,1},CuArray{Float64,2}},
     device::Nothing=nothing
 )
-    info = mod.info
+    info, par = mod.info, env.params
 
     info.primres = 0.0
     info.dualres = 0.0
+    info.primsca = 0.0
+    info.dualsca = 0.0
     info.norm_z_curr = 0.0
     info.mismatch = 0.0
 
     for i=1:mod.len_horizon
-        admm_update_residual(env, mod.models[i])
+        admm_update_residual(env, mod.models[i], normalized = false)
+        info.primsca = max(info.primsca, mod.models[i].info.primsca)
+        info.dualsca = max(info.dualsca, mod.models[i].info.dualsca)
     #=
         info.primres += mod.models[i].info.primres^2
         info.dualres += mod.models[i].info.dualres^2
@@ -44,6 +48,8 @@ function admm_update_residual(
         mod.models[i].info.dualres = sqrt(mod.models[i].info.dualres^2 + CUDA.norm(sol_ramp.rd)^2)
         mod.models[i].info.norm_z_curr = sqrt(mod.models[i].info.norm_z_curr^2 + CUDA.norm(sol_ramp.z_curr)^2)
         mod.models[i].info.mismatch = sqrt(mod.models[i].info.mismatch^2 + CUDA.norm(sol_ramp.Ax_plus_By)^2)
+        info.primsca = max(info.primsca, CUDA.norm(sol_ramp.u_curr), CUDA.norm(sol_ramp.v_curr), CUDA.norm(sol_ramp.z_curr))
+        info.dualsca = max(info.dualsca, CUDA.norm(sol_ramp.l_curr))
     #=
         info.primres += norm(sol_ramp.rp)^2
         info.dualres += norm(sol_ramp.rd)^2
@@ -58,6 +64,10 @@ function admm_update_residual(
         info.norm_z_curr = max(info.norm_z_curr, mod.models[i].info.norm_z_curr)
         info.mismatch = max(info.mismatch, mod.models[i].info.mismatch)
     end
+    info.primres /= info.primsca
+    info.dualres /= info.dualsca
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL / info.primsca + par.RELTOL
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL / info.dualsca + par.RELTOL
 
     #=
     info.primres = sqrt(info.primres)

src/models/qpsub/qpsub_admm_update_residual_cpu.jl

@@ -19,8 +19,12 @@ function admm_update_residual(
     sol.rd .= sol.rho .* (sol.v_curr - mod.v_prev)#from Boyd's single-level admm
     sol.Ax_plus_By .= sol.rp #x-xbar
 
-    info.primres = norm(sol.rp)
-    info.dualres = norm(sol.rd)
+    info.primsca = max(norm(sol.u_curr), norm(sol.v_curr))
+    info.dualsca = norm(sol.l_curr)
+    info.primres = norm(sol.rp) / info.primsca
+    info.dualres = norm(sol.rd) / info.dualsca
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL / info.primsca + par.RELTOL
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL / info.dualsca + par.RELTOL
     info.mismatch = norm(sol.Ax_plus_By)
 
 

src/models/qpsub/qpsub_admm_update_residual_gpu.jl

@@ -42,11 +42,12 @@ function admm_update_residual(
 
     @cuda threads=64 blocks=(div(mod.nvar-1, 64)+1) copy_data_kernel(mod.nvar, sol.Ax_plus_By, sol.rp) # from gpu utility
 
-
-    info.primres = CUDA.norm(sol.rp)
-
-    info.dualres = CUDA.norm(sol.rd)
-
+    info.primsca = max(CUDA.norm(sol.u_curr), CUDA.norm(sol.v_curr))
+    info.dualsca = CUDA.norm(sol.l_curr)
+    info.primres = CUDA.norm(sol.rp) / info.primsca
+    info.dualres = CUDA.norm(sol.rd) / info.dualsca
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL / info.primsca + par.RELTOL
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL / info.dualsca + par.RELTOL
     info.mismatch = CUDA.norm(sol.Ax_plus_By)
 
     return

src/models/qpsub/qpsub_admm_update_residual_ka.jl

@@ -51,11 +51,12 @@ function admm_update_residual(
     ) # from gpu utility
     KA.synchronize(device)
 
-
-    info.primres = norm(sol.rp, device)
-
-    info.dualres = norm(sol.rd, device)
-
+    info.primsca = max(norm(sol.u_curr, device), norm(sol.v_curr, device))
+    info.dualsca = norm(sol.l_curr, device)
+    info.primres = norm(sol.rp, device) / info.primsca
+    info.dualres = norm(sol.rd, device) / info.dualsca
+    info.primtol = sqrt(mod.nvar) * par.ABSTOL / info.primsca + par.RELTOL
+    info.dualtol = sqrt(mod.nvar) * par.ABSTOL / info.dualsca + par.RELTOL
     info.mismatch = norm(sol.Ax_plus_By, device)
 
     return

src/models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_gpu.jl

@@ -51,6 +51,8 @@ function auglag_linelimit_qpsub(Hs, l_curr, rho, u_curr, v_curr, z_curr, YffR, Y
     Atron = CuDynamicSharedArray(Float64, (6,6), (5*n + n^2 + 136)*sizeof(Float64))
     btron = CuDynamicSharedArray(Float64, 6, (5*n + n^2 + 172)*sizeof(Float64))
 
+    fill!(Ctron,0.0)
+    CUDA.sync_threads()
 
     #initialization: variable wrt branch structure wrt Exatron
     x[1] = 0.0
@@ -108,7 +110,6 @@ function auglag_linelimit_qpsub(Hs, l_curr, rho, u_curr, v_curr, z_curr, YffR, Y
     inv12 = -LH_1h[lineidx,2]/prod
     inv21 = -LH_1i[lineidx,1]/prod
     inv22 = LH_1h[lineidx,1]/prod
-    fill!(Ctron,0.0)
     Ctron[1,1] = 1.0
     Ctron[2,2] = 1.0
     Ctron[3,1] = 0.0
@@ -162,7 +163,7 @@ function auglag_linelimit_qpsub(Hs, l_curr, rho, u_curr, v_curr, z_curr, YffR, Y
 
 
     while !terminate
-         it += 1
+        it += 1
 
         eval_A_b_auglag_branch_kernel_gpu_qpsub_red(Hbr, bbr, A_aug, Atron, btron, scale,vec_1j,vec_1k,membuf,lineidx, Ctron,dtron,RH_1j,RH_1k)
 
@@ -215,7 +216,7 @@ function auglag_linelimit_qpsub(Hs, l_curr, rho, u_curr, v_curr, z_curr, YffR, Y
 
         if it >= max_auglag #maximum iteration for auglag
             if tx == 1
-                @cuprintln("max iteration reach at block = ",lineidx, "and threads = ",tx)
+                @cuprintln("max_auglag reached for line with cnorm = ", cnorm, " for line = ", lineidx, " with mu ", mu)
             end
             terminate = true
         end

src/models/qpsub/qpsub_auglag_Ab_linelimit_kernel_red_ka.jl

@@ -53,6 +53,10 @@
     Atron = @localmem Float64 (6,6)
     btron = @localmem Float64 (6,)
 
+    fill!(Ctron,0.0)
+    fill!(A_aug, 0.0)
+    fill!(Atron, 0.0)
+    CUDA.sync_threads()
 
     #initialization: variable wrt branch structure wrt Exatron
     x[1] = 0.0
@@ -112,7 +116,6 @@
     inv12 = -LH_1h[lineidx,2]/prod
     inv21 = -LH_1i[lineidx,1]/prod
     inv22 = LH_1h[lineidx,1]/prod
-    fill!(Ctron,0.0)
     Ctron[1,1] = 1.0
     Ctron[2,2] = 1.0
     Ctron[3,1] = 0.0
@@ -166,7 +169,7 @@
 
 
     while !terminate
-         it += 1
+        it += 1
 
         eval_A_b_auglag_branch_kernel_gpu_qpsub_red(
             Hbr, bbr, A_aug, Atron, btron, scale,vec_1j,vec_1k,membuf,lineidx, Ctron,dtron,RH_1j,RH_1k,

src/models/qpsub/qpsub_eval_Ab_linelimit_kernel_gpu.jl

@@ -72,9 +72,12 @@ end
 
 function eval_A_b_auglag_branch_kernel_gpu_qpsub_red(Hbr, bbr, A_aug, Atron, btron, scale,vec_1j,vec_1k,membuf,lineidx, Ctron,dtron,RH_1j,RH_1k)
     tx = threadIdx().x
+
+    fill!(A_aug, 0.0)
+    fill!(Atron, 0.0)
+    CUDA.sync_threads()
+
     if tx == 1
-        fill!(A_aug, 0.0)
-        fill!(Atron, 0.0)
 
         @inbounds begin
                     for i = 3:8