diff --git a/src/algorithms/admm_one_level.jl b/src/algorithms/admm_one_level.jl index a80fdb8..a040cf7 100644 --- a/src/algorithms/admm_one_level.jl +++ b/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 @@ -27,7 +27,7 @@ function admm_one_level( "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) + info.outer, info.objval, info.auglag, info.mismatch, info.primtol, info.dualres, info.dualtol) end info.status = :IterationLimit @@ -56,13 +56,13 @@ function admm_one_level( 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) + info.outer, info.objval, info.auglag, info.primres, info.primtol, info.dualres, info.dualtol) 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 diff --git a/src/algorithms/admm_two_level.jl b/src/algorithms/admm_two_level.jl index 559ffbc..611a466 100644 --- a/src/algorithms/admm_two_level.jl +++ b/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 diff --git a/src/models/acopf/acopf_admm_update_residual_cpu.jl b/src/models/acopf/acopf_admm_update_residual_cpu.jl index 9d8a29f..0eb2e11 100644 --- a/src/models/acopf/acopf_admm_update_residual_cpu.jl +++ b/src/models/acopf/acopf_admm_update_residual_cpu.jl @@ -13,7 +13,8 @@ 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 @@ -21,8 +22,18 @@ function admm_update_residual( 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 + diff --git a/src/models/acopf/acopf_admm_update_residual_gpu.jl b/src/models/acopf/acopf_admm_update_residual_gpu.jl index 6754247..04a0e17 100644 --- a/src/models/acopf/acopf_admm_update_residual_gpu.jl +++ b/src/models/acopf/acopf_admm_update_residual_gpu.jl @@ -11,7 +11,8 @@ 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 @@ -20,8 +21,18 @@ function admm_update_residual( 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) diff --git a/src/models/acopf/acopf_admm_update_residual_ka.jl b/src/models/acopf/acopf_admm_update_residual_ka.jl index 429c48f..fc09a35 100644 --- a/src/models/acopf/acopf_admm_update_residual_ka.jl +++ b/src/models/acopf/acopf_admm_update_residual_ka.jl @@ -13,7 +13,8 @@ end function admm_update_residual( env::AdmmEnv, mod::AbstractOPFModel, - device + device; + normalized=true ) sol, info = mod.solution, mod.info nblk_nvar = div(mod.nvar-1, 64)+1 @@ -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) diff --git a/src/models/mpacopf/mpacopf_admm_update_residual_cpu.jl b/src/models/mpacopf/mpacopf_admm_update_residual_cpu.jl index d1fde6c..fae503d 100644 --- a/src/models/mpacopf/mpacopf_admm_update_residual_cpu.jl +++ b/src/models/mpacopf/mpacopf_admm_update_residual_cpu.jl @@ -7,11 +7,15 @@ function admm_update_residual( 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) diff --git a/src/models/mpacopf/mpacopf_admm_update_residual_gpu.jl b/src/models/mpacopf/mpacopf_admm_update_residual_gpu.jl index 51fc4c5..74e2802 100644 --- a/src/models/mpacopf/mpacopf_admm_update_residual_gpu.jl +++ b/src/models/mpacopf/mpacopf_admm_update_residual_gpu.jl @@ -20,11 +20,15 @@ function admm_update_residual( 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(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) diff --git a/src/models/qpsub/qpsub_admm_update_residual_cpu.jl b/src/models/qpsub/qpsub_admm_update_residual_cpu.jl index 607fd8b..561bafa 100644 --- a/src/models/qpsub/qpsub_admm_update_residual_cpu.jl +++ b/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) diff --git a/src/models/qpsub/qpsub_admm_update_residual_gpu.jl b/src/models/qpsub/qpsub_admm_update_residual_gpu.jl index 98e7f33..8c27a1f 100644 --- a/src/models/qpsub/qpsub_admm_update_residual_gpu.jl +++ b/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 diff --git a/src/models/qpsub/qpsub_admm_update_residual_ka.jl b/src/models/qpsub/qpsub_admm_update_residual_ka.jl index 082bb1a..9500dff 100644 --- a/src/models/qpsub/qpsub_admm_update_residual_ka.jl +++ b/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 diff --git a/src/utils/environment.jl b/src/utils/environment.jl index 040fde6..6fa031a 100644 --- a/src/utils/environment.jl +++ b/src/utils/environment.jl @@ -333,6 +333,10 @@ mutable struct IterationInformation{U} objval::Float64 primres::Float64 dualres::Float64 + primtol::Float64 # normalized primal tolerance; this may be redundant to eps_pri below. + dualtol::Float64 # normalized dual tolerance + primsca::Float64 # primal normalization scalar + dualsca::Float64 # primal normalization scalar mismatch::Float64 auglag::Float64 eps_pri::Float64 @@ -365,6 +369,10 @@ function Base.fill!(info::IterationInformation, val) info.objval = val info.primres = val info.dualres = val + info.primtol = val + info.dualtol = val + info.primsca = val + info.dualsca = val info.mismatch = val info.auglag = val info.eps_pri = val @@ -389,6 +397,10 @@ function Base.copy(ref::IterationInformation{ComponentInformation}) info.objval = ref.objval info.primres = ref.primres info.dualres = ref.dualres + info.primtol = ref.primtol + info.dualtol = ref.dualtol + info.primsca = ref.primsca + info.dualsca = ref.dualsca info.mismatch = ref.mismatch info.auglag = ref.auglag info.eps_pri = ref.eps_pri