From 2870d904176ea94cbb8096d47d2fc4374bffb1cb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Pacaud?= Date: Wed, 27 Sep 2023 15:42:51 +0200 Subject: [PATCH] Add preventive and corrective SCOPF (#67) - migrate to MadNLP 0.6 - migrate to CUDA 4.* - migrate to ExaPF 0.9 - add preventive and corrective SCOPF - add script for real-time SCOPF --------- Co-authored-by: Michel Schanen --- .github/workflows/action.yml | 1 - Project.toml | 14 +- docs/make.jl | 1 + lib/ArgosCUDA.jl/Project.toml | 1 - lib/ArgosCUDA.jl/src/ArgosCUDA.jl | 4 +- lib/ArgosCUDA.jl/src/kernels.jl | 42 +- lib/ArgosCUDA.jl/src/reduction.jl | 16 +- scripts/jump/scopf.jl | 239 ++++++++ scripts/scopf/Artifacts.toml | 6 + scripts/scopf/LocalPreferences.toml | 2 + scripts/scopf/Manifest.toml | 793 +++++++++++++++++++++++++ scripts/scopf/Project.toml | 17 + scripts/scopf/benchmark.jl | 102 ++++ scripts/scopf/benchmark_biegler.jl | 93 +++ scripts/scopf/benchmark_full_space.jl | 91 +++ scripts/scopf/common.jl | 119 ++++ scripts/scopf/config.jl | 53 ++ scripts/scopf/corrective_scopf.jl | 119 ++++ scripts/scopf/ipopt.jl | 27 + scripts/scopf/kkt.jl | 173 ++++++ scripts/scopf/preventive_scopf.jl | 81 +++ scripts/scopf/screening.jl | 74 +++ scripts/scopf/tracking.jl | 47 ++ src/Evaluators/Evaluators.jl | 1 + src/Evaluators/corrective_evaluator.jl | 280 +++++++++ src/Evaluators/stoch_evaluator.jl | 57 +- src/KKT/auglag_kkt.jl | 2 +- src/KKT/reduced_newton.jl | 38 +- test/Algorithms/MadNLP_wrapper.jl | 9 +- test/Evaluators/TestEvaluators.jl | 14 + test/Evaluators/api.jl | 7 +- test/runtests.jl | 3 +- 32 files changed, 2454 insertions(+), 72 deletions(-) create mode 100644 scripts/jump/scopf.jl create mode 100644 scripts/scopf/Artifacts.toml create mode 100644 scripts/scopf/LocalPreferences.toml create mode 100644 scripts/scopf/Manifest.toml create mode 100644 scripts/scopf/Project.toml create mode 100644 scripts/scopf/benchmark.jl create mode 100644 scripts/scopf/benchmark_biegler.jl create mode 100644 scripts/scopf/benchmark_full_space.jl create mode 100644 scripts/scopf/common.jl create mode 100644 scripts/scopf/config.jl create mode 100644 scripts/scopf/corrective_scopf.jl create mode 100644 scripts/scopf/ipopt.jl create mode 100644 scripts/scopf/kkt.jl create mode 100644 scripts/scopf/preventive_scopf.jl create mode 100644 scripts/scopf/screening.jl create mode 100644 scripts/scopf/tracking.jl create mode 100644 src/Evaluators/corrective_evaluator.jl diff --git a/.github/workflows/action.yml b/.github/workflows/action.yml index d18ba185..ac6f2fef 100644 --- a/.github/workflows/action.yml +++ b/.github/workflows/action.yml @@ -39,7 +39,6 @@ jobs: env: CUDA_VISIBLE_DEVICES: 1 JULIA_DEPOT_PATH: /scratch/github-actions/julia_depot_argos - JULIA_CUDA_USE_BINARYBUILDER: false runs-on: self-hosted strategy: matrix: diff --git a/Project.toml b/Project.toml index 515cfc4a..c1059fe4 100644 --- a/Project.toml +++ b/Project.toml @@ -14,20 +14,18 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] -CUDA = "^3.8" -CUDAKernels = "0.4" -ExaPF = "~0.8" +CUDA = "4.1" +ExaPF = "~0.9.3" FiniteDiff = "2.7" Ipopt = "1" -KernelAbstractions = "0.8" -MadNLP = "~0.5.2" +KernelAbstractions = "0.9" +MadNLP = "0.7" MathOptInterface = "1" -NLPModels = "0.18" +NLPModels = ">= 0.18" julia = "1.6" [extras] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" -CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" @@ -36,4 +34,4 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "CUDA", "CUDAKernels", "DelimitedFiles", "FiniteDiff", "Ipopt", "LazyArtifacts", "Random"] +test = ["Test", "CUDA", "DelimitedFiles", "FiniteDiff", "Ipopt", "LazyArtifacts", "Random"] diff --git a/docs/make.jl b/docs/make.jl index 00a44c75..57f285a3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -19,6 +19,7 @@ makedocs( ), modules = [Argos], repo = "https://github.com/exanauts/Argos.jl/blob/{commit}{path}#{line}", + doctest=false, checkdocs = :exports, pages = [ "Home" => "index.md", diff --git a/lib/ArgosCUDA.jl/Project.toml b/lib/ArgosCUDA.jl/Project.toml index ed9c3a32..f9647636 100644 --- a/lib/ArgosCUDA.jl/Project.toml +++ b/lib/ArgosCUDA.jl/Project.toml @@ -6,7 +6,6 @@ version = "0.1.0" [deps] Argos = "ef244971-cf80-42b0-9762-2c2c832df5d5" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" -CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57" CUSOLVERRF = "a8cc9031-bad2-4722-94f5-40deabb4245c" ExaPF = "0cf0e50c-a82e-488f-ac7e-41ffdff1b8aa" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" diff --git a/lib/ArgosCUDA.jl/src/ArgosCUDA.jl b/lib/ArgosCUDA.jl/src/ArgosCUDA.jl index a1276306..430b27b6 100644 --- a/lib/ArgosCUDA.jl/src/ArgosCUDA.jl +++ b/lib/ArgosCUDA.jl/src/ArgosCUDA.jl @@ -7,8 +7,8 @@ using CUDA using CUDA.CUSPARSE using CUSOLVERRF -using KernelAbstractions -using CUDAKernels +import KernelAbstractions as KA +import KernelAbstractions: @kernel, @index using ExaPF const LS = ExaPF.LinearSolvers diff --git a/lib/ArgosCUDA.jl/src/kernels.jl b/lib/ArgosCUDA.jl/src/kernels.jl index 17906077..207807b9 100644 --- a/lib/ArgosCUDA.jl/src/kernels.jl +++ b/lib/ArgosCUDA.jl/src/kernels.jl @@ -5,8 +5,8 @@ end function Argos.transfer2tril!(hessvals::AbstractVector, H::CuSparseMatrixCSR, csc2tril) Hz = nonzeros(H) ndrange = (length(hessvals),) - ev = _map2vec_kernel!(CUDADevice())(hessvals, Hz, csc2tril; ndrange=ndrange) - wait(ev) + _map2vec_kernel!(CUDABackend())(hessvals, Hz, csc2tril; ndrange=ndrange) + KA.synchronize(CUDABackend()) end @@ -17,8 +17,8 @@ end function Argos.fixed!(dest::CuVector, ind_fixed, val::Number) length(ind_fixed) == 0 && return g_ind_fixed = ind_fixed |> CuArray - ev = _fixed_kernel!(CUDADevice())(dest, g_ind_fixed, val; ndrange=length(ind_fixed)) - wait(ev) + _fixed_kernel!(CUDABackend())(dest, g_ind_fixed, val; ndrange=length(ind_fixed)) + KA.synchronize(CUDABackend()) end @@ -30,8 +30,8 @@ function Argos.copy_index!(dest::CuVector{T}, src::CuVector{T}, idx) where T @assert length(dest) == length(idx) ndrange = (length(dest),) idx_d = CuVector(idx) - ev = _copy_index_kernel!(CUDADevice())(dest, src, idx_d; ndrange=ndrange) - wait(ev) + _copy_index_kernel!(CUDABackend())(dest, src, idx_d; ndrange=ndrange) + KA.synchronize(CUDABackend()) end @@ -43,8 +43,8 @@ end function Argos.fixed_diag!(dest::CuMatrix, ind_fixed, val::Number) length(ind_fixed) == 0 && return g_ind_fixed = ind_fixed |> CuArray - ev = _fixed_diag_kernel!(CUDADevice())(dest, g_ind_fixed, val; ndrange=length(ind_fixed)) - wait(ev) + _fixed_diag_kernel!(CUDABackend())(dest, g_ind_fixed, val; ndrange=length(ind_fixed)) + KA.synchronize(CUDABackend()) end @kernel function _transfer_auglag_hessian!(dest, H, J, ρ, n, m) @@ -70,8 +70,8 @@ function Argos.transfer_auglag_hessian!( @assert size(ρ, 1) == m ndrange = (n+m, n) - ev = _transfer_auglag_hessian!(CUDADevice())(dest, H, J, ρ, n, m, ndrange=ndrange, dependencies=Event(CUDADevice())) - wait(ev) + _transfer_auglag_hessian!(CUDABackend())(dest, H, J, ρ, n, m, ndrange=ndrange) + KA.synchronize(CUDABackend()) return end @@ -84,11 +84,11 @@ function Argos.set_batch_tangents!(seeds::CuMatrix, offset, n, n_batches) @assert offset + n_batches <= n ndrange = (n_batches) fill!(seeds, 0.0) - ev = _batch_tangents_kernel!(CUDADevice())( + _batch_tangents_kernel!(CUDABackend())( seeds, offset, n_batches; - ndrange=ndrange, dependencies=Event(CUDADevice()), + ndrange=ndrange, ) - wait(ev) + KA.synchronize(CUDABackend()) return end @@ -115,11 +115,11 @@ function Argos.tgtmul!( k = size(z, 2) ndrange = (n, k) y .*= beta - ev = _tgtmul_1_kernel!(CUDADevice())( + _tgtmul_1_kernel!(CUDABackend())( y, A.rowPtr, A.colVal, A.nzVal, z, w, alpha, nz, nw; - ndrange=ndrange, dependencies=Event(CUDADevice()), + ndrange=ndrange, ) - wait(ev) + KA.synchronize(CUDABackend()) end @@ -150,11 +150,11 @@ function Argos.tgtmul!( ndrange = (n, k) yx .*= beta yu .*= beta - ev = _tgtmul_2_kernel!(CUDADevice())( + _tgtmul_2_kernel!(CUDABackend())( yx, yu, A.rowPtr, A.colVal, A.nzVal, z, w, alpha, nz, nw; - ndrange=ndrange, dependencies=Event(CUDADevice()), + ndrange=ndrange, ) - wait(ev) + KA.synchronize(CUDABackend()) end @@ -171,11 +171,11 @@ end function Argos.update!(K::Argos.HJDJ, A, D, Σ) m = size(A, 1) - ev = _scale_transpose_kernel!(CUDADevice())( + ev = _scale_transpose_kernel!(CUDABackend())( K.Jt.nzVal, A.rowPtr, A.colVal, A.nzVal, D, K.transperm, ndrange=(m, 1), ) - wait(ev) + KA.synchronize(ev) spgemm!('N', 'N', 1.0, K.Jt, A, 0.0, K.JtJ, 'O') K.Σ .= Σ end diff --git a/lib/ArgosCUDA.jl/src/reduction.jl b/lib/ArgosCUDA.jl/src/reduction.jl index 1a93b2a0..d4c6ae92 100644 --- a/lib/ArgosCUDA.jl/src/reduction.jl +++ b/lib/ArgosCUDA.jl/src/reduction.jl @@ -24,21 +24,21 @@ end function Argos.update!(K::Argos.HJDJ, A, D, Σ) m = size(A, 1) - ev = _scale_transpose_kernel!(CUDADevice())( + _scale_transpose_kernel!(CUDABackend())( K.Jt.nzVal, A.rowPtr, A.colVal, A.nzVal, D, K.transperm, ndrange=(m, 1), ) - wait(ev) + KA.synchronize(CUDABackend()) spgemm!('N', 'N', 1.0, K.Jt, A, 0.0, K.JtJ, 'O') K.Σ .= Σ end function Argos.update!(K::Argos.HJDJ, A, D) m = size(A, 1) - ev = _scale_transpose_kernel!(CUDADevice())( + _scale_transpose_kernel!(CUDABackend())( K.Jt.nzVal, A.rowPtr, A.colVal, A.nzVal, D, K.transperm, ndrange=(m, 1), ) - wait(ev) + KA.synchronize(CUDABackend()) spgemm!('N', 'N', 1.0, K.Jt, A, 0.0, K.JtJ, 'O') end @@ -47,7 +47,13 @@ function MadNLP.set_aug_diagonal!(kkt::Argos.BieglerKKTSystem{T, VI, VT, MT}, ip pr_diag_h = kkt.etc[:pr_diag_host]::Vector{T} # Broadcast is not working as MadNLP array are allocated on the CPU, # whereas pr_diag is allocated on the GPU - pr_diag_h .= ips.zl./(ips.x.-ips.xl) .+ ips.zu./(ips.xu.-ips.x) + x = MadNLP.full(ips.x) + xl = MadNLP.full(ips.xl) + xu = MadNLP.full(ips.xu) + zl = MadNLP.full(ips.zl) + zu = MadNLP.full(ips.zu) + + pr_diag_h .= zl ./ (x .- xl) .+ zu ./ (xu .- x) copyto!(kkt.pr_diag, pr_diag_h) fill!(kkt.du_diag, 0.0) end diff --git a/scripts/jump/scopf.jl b/scripts/jump/scopf.jl new file mode 100644 index 00000000..b48a91db --- /dev/null +++ b/scripts/jump/scopf.jl @@ -0,0 +1,239 @@ + +using ExaPF +using JuMP +using Ipopt +using LinearAlgebra +using LazyArtifacts +using DelimitedFiles + +using HSL_jll + +const PS = ExaPF.PowerSystem + +const DATA = joinpath(artifact"ExaData", "ExaData") +const SCENARIOS = joinpath(artifact"ExaData", "ExaData", "mp_demand") + +constraint_index(cons::Vector{NonlinearConstraintRef{ScalarShape}}) = getfield.(JuMP.index.(cons), :value) + +function PowerNetworkContingency(network::PS.PowerNetwork, remove_line) + data = Dict{String, Array}() + data["bus"] = copy(network.buses) + data["branch"] = copy(network.branches) + data["gen"] = copy(network.generators) + data["baseMVA"] = Float64[network.baseMVA] + data["cost"] = copy(network.costs) + return PS.PowerNetwork(data; remove_lines=Int[remove_line]) +end + +function build_scopf_model(polar, buffer, lines_id, solver) + nbus = ExaPF.get(polar, PS.NumberOfBuses()) + ngen = ExaPF.get(polar, PS.NumberOfGenerators()) + nlines = ExaPF.get(polar, PS.NumberOfLines()) + + pf = polar.network + baseMVA = pf.baseMVA + + pg_min, pg_max = PS.bounds(pf, PS.Generators(), PS.ActivePower()) + qg_min, qg_max = PS.bounds(pf, PS.Generators(), PS.ReactivePower()) + vm_min, vm_max = PS.bounds(pf, PS.Buses(), PS.VoltageMagnitude()) + + flow_min, flow_max = PS.bounds(pf, PS.Lines(), PS.ActivePower()) + flow_max = min.(1e5, flow_max) + + vm0 = buffer.vmag + va0 = buffer.vang + pg0 = buffer.pgen + + Pd = buffer.pload + Qd = buffer.qload + + cost_coefs = PS.get_costs_coefficients(pf) + + bus2gen = PS.get_bus_generators(pf.buses, pf.generators, pf.bus_to_indexes) + + # Power flow data + nscen = length(lines_id) + + #= + Build model + =# + + opfmodel = Model(solver) + + # VARIABLES + @variable(opfmodel, pg_min[i] <= Pg[i=1:ngen, k=1:nscen+1] <= pg_max[i], start=pg0[i]) + @variable(opfmodel, qg_min[i] <= Qg[i=1:ngen, k=1:nscen+1] <= qg_max[i]) + @variable(opfmodel, vm_min[i] <= Vm[i=1:nbus, k=1:nscen+1] <= vm_max[i], start=vm0[i]) + @variable(opfmodel, Va[i=1:nbus, k=1:nscen+1], start=va0[i]) + + # Preventive constraints + # u0 = uk + gref = findfirst(isequal(pf.ref[1]), pf.gen2bus) + for k in 1:nscen + for g in 1:ngen + (g == gref) && continue # power generation at ref is a slack + @constraint(opfmodel, Pg[g, 1] == Pg[g, k+1]) + end + for b in [pf.pv; pf.ref] + @constraint(opfmodel, Vm[b, 1] == Vm[b, k+1]) + end + end + + # Reference angle + for k in 1:nscen+1 + for b in pf.ref + @constraint(opfmodel, Va[b, k] == 0) + end + end + + # Power-flow constraints + for k in 1:nscen+1 + # TODO: this is inefficient and should be removed in a + # future version to avoid creating a new admittance matrix + # for each contingency. + network = if k >= 2 + PowerNetworkContingency(pf, lines_id[k-1]) + else + pf + end + Ybus = network.Ybus + rows = Ybus.rowval + yvals = Ybus.nzval + g_ij = real.(yvals) + b_ij = imag.(yvals) + + ## active + opfmodel.ext[:active_pf] = @NLconstraint( + opfmodel, [b=1:nbus], + Vm[b, k] * sum( + Vm[rows[c], k] * (g_ij[c] * cos(Va[b, k] - Va[rows[c], k]) + b_ij[c] * sin(Va[b, k] - Va[rows[c], k])) + for c in (Ybus.colptr[b]):(Ybus.colptr[b+1]-1) + ) == (sum(Pg[g, k] for g in get(bus2gen, b, Int[])) - Pd[b]) + ) + ## reactive + opfmodel.ext[:reactive_pf] = @NLconstraint( + opfmodel, [b=1:nbus], + Vm[b, k] * sum( + Vm[rows[c], k] * (g_ij[c] * sin(Va[b, k] - Va[rows[c], k]) - b_ij[c] * cos(Va[b, k] - Va[rows[c], k])) + for c in (Ybus.colptr[b]):(Ybus.colptr[b+1]-1)) == (sum(Qg[g, k] for g in get(bus2gen, b, Int[])) - Qd[b]) + ) + end + + # Line constraints + f = pf.lines.from_buses + t = pf.lines.to_buses + # Data + yff_re = real.(pf.lines.Yff) + yff_im = imag.(pf.lines.Yff) + yft_re = real.(pf.lines.Yft) + yft_im = imag.(pf.lines.Yft) + ytf_re = real.(pf.lines.Ytf) + ytf_im = imag.(pf.lines.Ytf) + ytt_re = real.(pf.lines.Ytt) + ytt_im = imag.(pf.lines.Ytt) + ## from lines + yff_abs = yff_re.^2 .+ yff_im.^2 + yft_abs = yft_re.^2 .+ yft_im.^2 + yre_fr = yff_re .* yft_re .+ yff_im .* yft_im + yim_fr = - yff_re .* yft_im .+ yff_im .* yft_re + + ## to lines + ytf_abs = ytf_re.^2 .+ ytf_im.^2 + ytt_abs = ytt_re.^2 .+ ytt_im.^2 + yre_to = ytf_re .* ytt_re .+ ytf_im .* ytt_im + yim_to = - ytf_re .* ytt_im .+ ytf_im .* ytt_re + + for k in 1:nscen+1 + for ℓ in 1:nlines + # Remove the lines except if base case k == 1 + if k >= 2 && ℓ == lines_id[k-1] + continue + end + @NLconstraint( + opfmodel, + Vm[f[ℓ], k]^2 * ( + yff_abs[ℓ] * Vm[f[ℓ], k]^2 + yft_abs[ℓ] * Vm[t[ℓ], k]^2 + + 2.0 * Vm[f[ℓ], k] * Vm[t[ℓ], k] * (yre_fr[ℓ] * cos(Va[f[ℓ], k]-Va[t[ℓ], k]) - yim_fr[ℓ] * sin(Va[f[ℓ], k]-Va[t[ℓ], k])) + ) <= flow_max[ℓ] + ) + @NLconstraint( + opfmodel, + Vm[t[ℓ], k]^2 * ( + ytf_abs[ℓ] * Vm[f[ℓ], k]^2 + ytt_abs[ℓ] * Vm[t[ℓ], k]^2 + + 2.0 * Vm[f[ℓ], k] * Vm[t[ℓ], k] * (yre_to[ℓ] * cos(Va[f[ℓ], k]-Va[t[ℓ], k]) - yim_to[ℓ] * sin(Va[f[ℓ], k]-Va[t[ℓ], k])) + ) <= flow_max[ℓ] + ) + end + end + + # Objective : costs in base case + @objective( + opfmodel, + Min, + sum( + cost_coefs[g, 4] * Pg[g, 1]^2 + cost_coefs[g, 3] * Pg[g, 1] + cost_coefs[g, 2] + for g in 1:ngen + ) + ) + + opfmodel.ext[:exapf] = pf + + return opfmodel +end + +# Need a quick fix in Ipopt.jl to work properly! +function attach_callback!(opfmodel::Model) + pf_active = constraint_index(opfmodel.ext[:active_pf]) + pf_reactive = constraint_index(opfmodel.ext[:reactive_pf]) + index_line_fr = constraint_index(opfmodel.ext[:line_fr]) + index_line_to = constraint_index(opfmodel.ext[:line_to]) + line_flow = [index_line_fr; index_line_to] + + opfmodel.ext[:feas_active_pf] = Float64[] + opfmodel.ext[:feas_reactive_pf] = Float64[] + opfmodel.ext[:feas_line_flow] = Float64[] + + function my_callback( + prob::IpoptProblem, + alg_mod::Cint, + iter_count::Cint, + obj_value::Float64, + inf_pr::Float64, + inf_du::Float64, + mu::Float64, + d_norm::Float64, + regularization_size::Float64, + alpha_du::Float64, + alpha_pr::Float64, + ls_trials::Cint, + ) + push!(opfmodel.ext[:feas_active_pf], norm(prob.g[pf_active], Inf)) + push!(opfmodel.ext[:feas_reactive_pf], norm(prob.g[pf_reactive], Inf)) + push!(opfmodel.ext[:feas_line_flow], max(maximum(prob.g[line_flow]), 0.0)) + return true + end + + MOI.set(opfmodel, Ipopt.CallbackFunction(), my_callback) + return opfmodel +end + +function scopf(case::String; max_contingencies=32) + datafile = joinpath(DATA, case) + polar = ExaPF.PolarForm(datafile) + stack = ExaPF.NetworkStack(polar) + instance = split(case, ".")[1] + lines_id = readdlm(joinpath(SCENARIOS, "$(instance).Ctgs"), ',', Int)[:] + + if length(lines_id) >= max_contingencies + @info "Capping the number of contingencies to $(max_contingencies)" + lines_id = lines_id[1:max_contingencies] + end + + model = build_scopf_model(polar, stack, lines_id, Ipopt.Optimizer) + JuMP.set_attribute(model, "tol", 1e-5) + JuMP.set_attribute(model, "hsllib", HSL_jll.libhsl_path) + JuMP.set_attribute(model, "linear_solver", "ma27") + JuMP.optimize!(model) + return model +end + diff --git a/scripts/scopf/Artifacts.toml b/scripts/scopf/Artifacts.toml new file mode 100644 index 00000000..1674f490 --- /dev/null +++ b/scripts/scopf/Artifacts.toml @@ -0,0 +1,6 @@ +[ExaData] +git-tree-sha1 = "d0421dd6bbdf0749990fb237f7017eb34b3f90c8" +lazy = true + [[ExaData.download]] + url = "https://web.cels.anl.gov/~mschanen/ExaData-affbca6.tar.gz" + sha256 = "96126184b0d5ec5c2924855766752f090a5a2e644491d9f5ff2a74c7334fa0aa" diff --git a/scripts/scopf/LocalPreferences.toml b/scripts/scopf/LocalPreferences.toml new file mode 100644 index 00000000..5da06c74 --- /dev/null +++ b/scripts/scopf/LocalPreferences.toml @@ -0,0 +1,2 @@ +[CUDA_Runtime_jll] +version = "11.8" diff --git a/scripts/scopf/Manifest.toml b/scripts/scopf/Manifest.toml new file mode 100644 index 00000000..955d0738 --- /dev/null +++ b/scripts/scopf/Manifest.toml @@ -0,0 +1,793 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.9.0" +manifest_format = "2.0" +project_hash = "c304673fa7bba6cab782a8c82da3e2ed8f75b020" + +[[deps.AMD]] +deps = ["Libdl", "LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "00163dc02b882ca5ec032400b919e5f5011dbd31" +uuid = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e" +version = "0.5.0" + +[[deps.ASL_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "6252039f98492252f9e47c312c8ffda0e3b9e78d" +uuid = "ae81ac8f-d209-56e5-92de-9978fef736f9" +version = "0.1.3+0" + +[[deps.AbstractFFTs]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "16b6dbc4cf7caee4e1e75c49485ec67b667098a0" +uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" +version = "1.3.1" + + [deps.AbstractFFTs.extensions] + AbstractFFTsChainRulesCoreExt = "ChainRulesCore" + + [deps.AbstractFFTs.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + +[[deps.Adapt]] +deps = ["LinearAlgebra", "Requires"] +git-tree-sha1 = "76289dc51920fdc6e0013c872ba9551d54961c24" +uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +version = "3.6.2" +weakdeps = ["StaticArrays"] + + [deps.Adapt.extensions] + AdaptStaticArraysExt = "StaticArrays" + +[[deps.ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +version = "1.1.1" + +[[deps.Argos]] +deps = ["ExaPF", "KernelAbstractions", "LinearAlgebra", "MadNLP", "MathOptInterface", "NLPModels", "Printf", "SparseArrays"] +path = "../.." +uuid = "ef244971-cf80-42b0-9762-2c2c832df5d5" +version = "0.3.2" + +[[deps.ArgosCUDA]] +deps = ["Argos", "CUDA", "CUSOLVERRF", "ExaPF", "KernelAbstractions", "LinearAlgebra", "MadNLP", "MadNLPGPU", "SparseArrays"] +path = "../../lib/ArgosCUDA.jl" +uuid = "8946db8d-321b-4174-84ad-48e2f9b69c56" +version = "0.1.0" + +[[deps.ArnoldiMethod]] +deps = ["LinearAlgebra", "Random", "StaticArrays"] +git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9" +uuid = "ec485272-7323-5ecc-a04f-4719b315124d" +version = "0.1.0" + +[[deps.ArrayInterface]] +deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"] +git-tree-sha1 = "d3f758863a47ceef2248d136657cb9c033603641" +uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +version = "7.4.8" + + [deps.ArrayInterface.extensions] + ArrayInterfaceBandedMatricesExt = "BandedMatrices" + ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices" + ArrayInterfaceCUDAExt = "CUDA" + ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore" + ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore" + ArrayInterfaceTrackerExt = "Tracker" + + [deps.ArrayInterface.weakdeps] + BandedMatrices = "aae01518-5342-5314-be14-df237901396f" + BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" + CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" + StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" + Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.Atomix]] +deps = ["UnsafeAtomics"] +git-tree-sha1 = "c06a868224ecba914baa6942988e2f2aade419be" +uuid = "a9b6321e-bd34-4604-b9c9-b65b8de01458" +version = "0.1.0" + +[[deps.BFloat16s]] +deps = ["LinearAlgebra", "Printf", "Random", "Test"] +git-tree-sha1 = "dbf84058d0a8cbbadee18d25cf606934b22d7c66" +uuid = "ab4f0b2a-ad5b-11e8-123f-65d77653426b" +version = "0.4.2" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[deps.BenchmarkTools]] +deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"] +git-tree-sha1 = "d9a9701b899b30332bbcb3e1679c41cce81fb0e8" +uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +version = "1.3.2" + +[[deps.BinaryProvider]] +deps = ["Libdl", "Logging", "SHA"] +git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058" +uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232" +version = "0.5.10" + +[[deps.Bzip2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2" +uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" +version = "1.0.8+0" + +[[deps.CEnum]] +git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90" +uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" +version = "0.4.2" + +[[deps.CUDA]] +deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "CompilerSupportLibraries_jll", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "Preferences", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "SpecialFunctions", "UnsafeAtomicsLLVM"] +git-tree-sha1 = "442d989978ed3ff4e174c928ee879dc09d1ef693" +uuid = "052768ef-5323-5732-b1bb-66c8b64840ba" +version = "4.3.2" + +[[deps.CUDA_Driver_jll]] +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] +git-tree-sha1 = "498f45593f6ddc0adff64a9310bb6710e851781b" +uuid = "4ee394cb-3365-5eb0-8335-949819d2adfc" +version = "0.5.0+1" + +[[deps.CUDA_Runtime_Discovery]] +deps = ["Libdl"] +git-tree-sha1 = "bcc4a23cbbd99c8535a5318455dcf0f2546ec536" +uuid = "1af6417a-86b4-443c-805f-a4643ffb695f" +version = "0.2.2" + +[[deps.CUDA_Runtime_jll]] +deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] +git-tree-sha1 = "5248d9c45712e51e27ba9b30eebec65658c6ce29" +uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2" +version = "0.6.0+0" + +[[deps.CUSOLVERRF]] +deps = ["CEnum", "CUDA", "KLU", "LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "79b418b0b42d34e0bb91f8ba3317cfdf35b8b28c" +uuid = "a8cc9031-bad2-4722-94f5-40deabb4245c" +version = "0.2.1" + +[[deps.CodecBzip2]] +deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"] +git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7" +uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" +version = "0.7.2" + +[[deps.CodecZlib]] +deps = ["TranscodingStreams", "Zlib_jll"] +git-tree-sha1 = "9c209fb7536406834aa938fb149964b985de6c83" +uuid = "944b1d66-785c-5afd-91f1-9de20f533193" +version = "0.7.1" + +[[deps.CommonSubexpressions]] +deps = ["MacroTools", "Test"] +git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" +uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" +version = "0.3.0" + +[[deps.Compat]] +deps = ["UUIDs"] +git-tree-sha1 = "7a60c856b9fa189eb34f5f8a6f6b5529b7942957" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.6.1" +weakdeps = ["Dates", "LinearAlgebra"] + + [deps.Compat.extensions] + CompatLinearAlgebraExt = "LinearAlgebra" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.0.2+0" + +[[deps.ConstructionBase]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "738fec4d684a9a6ee9598a8bfee305b26831f28c" +uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +version = "1.5.2" + + [deps.ConstructionBase.extensions] + ConstructionBaseIntervalSetsExt = "IntervalSets" + ConstructionBaseStaticArraysExt = "StaticArrays" + + [deps.ConstructionBase.weakdeps] + IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[[deps.DataStructures]] +deps = ["Compat", "InteractiveUtils", "OrderedCollections"] +git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0" +uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +version = "0.18.13" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[deps.DiffResults]] +deps = ["StaticArraysCore"] +git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621" +uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +version = "1.1.0" + +[[deps.DiffRules]] +deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] +git-tree-sha1 = "3d6d95e5f728741d74f121b9d0948c7372eee4b7" +uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" +version = "1.15.0" + +[[deps.Distributed]] +deps = ["Random", "Serialization", "Sockets"] +uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" + +[[deps.DocStringExtensions]] +deps = ["LibGit2"] +git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.9.3" + +[[deps.Downloads]] +deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +version = "1.6.0" + +[[deps.ExaPF]] +deps = ["CUDA", "ForwardDiff", "KernelAbstractions", "Krylov", "LazyArtifacts", "LightGraphs", "LinearAlgebra", "Metis", "Printf", "SparseArrays", "SparseDiffTools"] +git-tree-sha1 = "d44e4979913b644775f54c61f40da90ead99c0d5" +repo-rev = "fp/scopf" +repo-url = "https://github.com/exanauts/ExaPF.jl.git" +uuid = "0cf0e50c-a82e-488f-ac7e-41ffdff1b8aa" +version = "0.9.1" + +[[deps.ExprTools]] +git-tree-sha1 = "c1d06d129da9f55715c6c212866f5b1bddc5fa00" +uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" +version = "0.1.9" + +[[deps.FastClosures]] +git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef" +uuid = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" +version = "0.3.2" + +[[deps.FileWatching]] +uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" + +[[deps.FiniteDiff]] +deps = ["ArrayInterface", "LinearAlgebra", "Requires", "Setfield", "SparseArrays"] +git-tree-sha1 = "c6e4a1fbe73b31a3dea94b1da449503b8830c306" +uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" +version = "2.21.1" + + [deps.FiniteDiff.extensions] + FiniteDiffBandedMatricesExt = "BandedMatrices" + FiniteDiffBlockBandedMatricesExt = "BlockBandedMatrices" + FiniteDiffStaticArraysExt = "StaticArrays" + + [deps.FiniteDiff.weakdeps] + BandedMatrices = "aae01518-5342-5314-be14-df237901396f" + BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[[deps.ForwardDiff]] +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"] +git-tree-sha1 = "00e252f4d706b3d55a8863432e742bf5717b498d" +uuid = "f6369f11-7733-5829-9624-2563aa707210" +version = "0.10.35" +weakdeps = ["StaticArrays"] + + [deps.ForwardDiff.extensions] + ForwardDiffStaticArraysExt = "StaticArrays" + +[[deps.Future]] +deps = ["Random"] +uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" + +[[deps.GPUArrays]] +deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"] +git-tree-sha1 = "a3351bc577a6b49297248aadc23a4add1097c2ac" +uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" +version = "8.7.1" + +[[deps.GPUArraysCore]] +deps = ["Adapt"] +git-tree-sha1 = "2d6ca471a6c7b536127afccfa7564b5b39227fe0" +uuid = "46192b85-c4d5-4398-a991-12ede77f4527" +version = "0.1.5" + +[[deps.GPUCompiler]] +deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "TimerOutputs", "UUIDs"] +git-tree-sha1 = "cb090aea21c6ca78d59672a7e7d13bd56d09de64" +uuid = "61eb1bfa-7361-4325-ad38-22787b887f55" +version = "0.20.3" + +[[deps.Graphs]] +deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] +git-tree-sha1 = "1cf1d7dcb4bc32d7b4a5add4232db3750c27ecb4" +uuid = "86223c79-3864-5bf0-83f7-82e725a168b6" +version = "1.8.0" + +[[deps.Inflate]] +git-tree-sha1 = "5cd07aab533df5170988219191dfad0519391428" +uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" +version = "0.1.3" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[deps.Ipopt]] +deps = ["Ipopt_jll", "LinearAlgebra", "MathOptInterface", "OpenBLAS32_jll", "PrecompileTools"] +git-tree-sha1 = "d99b87f31175f7b4908f8d3719397792498b63de" +uuid = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +version = "1.4.0" + +[[deps.Ipopt_jll]] +deps = ["ASL_jll", "Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "MUMPS_seq_jll", "OpenBLAS32_jll", "libblastrampoline_jll"] +git-tree-sha1 = "0d3939fb672b84082f3ee930b51de03df935be31" +uuid = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7" +version = "300.1400.1200+0" + +[[deps.IrrationalConstants]] +git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" +uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" +version = "0.2.2" + +[[deps.JLLWrappers]] +deps = ["Preferences"] +git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.4.1" + +[[deps.JSON]] +deps = ["Dates", "Mmap", "Parsers", "Unicode"] +git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "0.21.4" + +[[deps.KLU]] +deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse_jll"] +git-tree-sha1 = "764164ed65c30738750965d55652db9c94c59bfe" +uuid = "ef3ab10e-7fda-4108-b977-705223b18434" +version = "0.4.0" + +[[deps.KernelAbstractions]] +deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"] +git-tree-sha1 = "47be64f040a7ece575c2b5f53ca6da7b548d69f4" +uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c" +version = "0.9.4" + +[[deps.Krylov]] +deps = ["LinearAlgebra", "Printf", "SparseArrays"] +git-tree-sha1 = "0356a64062656b0cbb43c504ad5de338251f4bda" +uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" +version = "0.9.1" + +[[deps.LDLFactorizations]] +deps = ["AMD", "LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "cbf4b646f82bfc58bb48bcca9dcce2eb88da4cd1" +uuid = "40e66cde-538c-5869-a4ad-c39174c6795b" +version = "0.10.0" + +[[deps.LLVM]] +deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"] +git-tree-sha1 = "26a31cdd9f1f4ea74f649a7bf249703c687a953d" +uuid = "929cbde3-209d-540e-8aea-75f648917ca0" +version = "5.1.0" + +[[deps.LLVMExtra_jll]] +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] +git-tree-sha1 = "09b7505cc0b1cee87e5d4a26eea61d2e1b0dcd35" +uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab" +version = "0.0.21+0" + +[[deps.LazyArtifacts]] +deps = ["Artifacts", "Pkg"] +uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" + +[[deps.LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.6.3" + +[[deps.LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" +version = "7.84.0+0" + +[[deps.LibGit2]] +deps = ["Base64", "NetworkOptions", "Printf", "SHA"] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[deps.LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" +version = "1.10.2+0" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[deps.LightGraphs]] +deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] +git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6" +uuid = "093fc24a-ae57-5d10-9952-331d41423f4d" +version = "1.3.5" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[deps.LinearOperators]] +deps = ["FastClosures", "LDLFactorizations", "LinearAlgebra", "Printf", "SparseArrays", "TimerOutputs"] +git-tree-sha1 = "a58ab1d18efa0bcf9f0868c6d387e4126dad3e72" +uuid = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125" +version = "2.5.2" + +[[deps.LogExpFunctions]] +deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "c3ce8e7420b3a6e071e0fe4745f5d4300e37b13f" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "0.3.24" + + [deps.LogExpFunctions.extensions] + LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" + LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" + LogExpFunctionsInverseFunctionsExt = "InverseFunctions" + + [deps.LogExpFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[deps.METIS_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "1fd0a97409e418b78c53fac671cf4622efdf0f21" +uuid = "d00139f3-1899-568f-a2f0-47f597d42d70" +version = "5.1.2+0" + +[[deps.MUMPS_seq_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "METIS_jll", "OpenBLAS32_jll", "Pkg", "libblastrampoline_jll"] +git-tree-sha1 = "f429d6bbe9ad015a2477077c9e89b978b8c26558" +uuid = "d7ed1dd3-d0ae-5e8e-bfb4-87a502085b8d" +version = "500.500.101+0" + +[[deps.MacroTools]] +deps = ["Markdown", "Random"] +git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.10" + +[[deps.MadNLP]] +deps = ["Libdl", "LinearAlgebra", "Logging", "MathOptInterface", "NLPModels", "Pkg", "Printf", "SolverCore", "SparseArrays", "SuiteSparse"] +git-tree-sha1 = "9e90600dd892820f4458ce5e031738a23a464074" +repo-rev = "ss/v1.9" +repo-url = "https://github.com/MadNLP/MadNLP.jl.git" +uuid = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6" +version = "0.6.0" + +[[deps.MadNLPGPU]] +deps = ["CUDA", "KernelAbstractions", "LinearAlgebra", "MadNLP"] +git-tree-sha1 = "b1e6d544a7c0c003152bb43d37ebc5b165ee31d5" +repo-rev = "fp/ka0.9" +repo-subdir = "lib/MadNLPGPU" +repo-url = "https://github.com/MadNLP/MadNLP.jl.git" +uuid = "d72a61cc-809d-412f-99be-fd81f4b8a598" +version = "0.5.0" + +[[deps.MadNLPHSL]] +deps = ["BinaryProvider", "Libdl", "MadNLP", "Pkg"] +git-tree-sha1 = "ca372dbb2521665b170444ce653358364ba4199f" +uuid = "7fb6135f-58fe-4112-84ca-653cf5be0c77" +version = "0.3.1" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[deps.MathOptInterface]] +deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test", "Unicode"] +git-tree-sha1 = "b5d162b61be1ab67de1311ae651b56729fd1b858" +uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +version = "1.17.0" + +[[deps.MbedTLS_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" +version = "2.28.2+0" + +[[deps.Metis]] +deps = ["CEnum", "LinearAlgebra", "METIS_jll", "SparseArrays"] +git-tree-sha1 = "66a4f74edb3ac5f28c74de60f9acc2a541fbbe28" +uuid = "2679e427-3c69-5b7f-982b-ece356f1e94b" +version = "1.4.0" + + [deps.Metis.extensions] + MetisGraphs = "Graphs" + MetisLightGraphs = "LightGraphs" + MetisSimpleWeightedGraphs = ["SimpleWeightedGraphs", "Graphs"] + + [deps.Metis.weakdeps] + Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" + LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" + SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622" + +[[deps.Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[deps.MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +version = "2022.10.11" + +[[deps.MutableArithmetics]] +deps = ["LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "964cb1a7069723727025ae295408747a0b36a854" +uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +version = "1.3.0" + +[[deps.NLPModels]] +deps = ["FastClosures", "LinearAlgebra", "LinearOperators", "Printf", "SparseArrays"] +git-tree-sha1 = "8205e052c2d67960a55dbe06cf31b577d19d76ad" +uuid = "a4795742-8479-5a88-8948-cc11e1c8c1a6" +version = "0.18.5" + +[[deps.NaNMath]] +deps = ["OpenLibm_jll"] +git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "1.0.2" + +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" + +[[deps.OpenBLAS32_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "2fb9ee2dc14d555a6df2a714b86b7125178344c2" +uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2" +version = "0.3.21+0" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.21+4" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.1+0" + +[[deps.OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.5+0" + +[[deps.OrderedCollections]] +git-tree-sha1 = "d321bf2de576bf25ec4d3e4360faca399afca282" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.6.0" + +[[deps.Parsers]] +deps = ["Dates", "PrecompileTools", "UUIDs"] +git-tree-sha1 = "a5aef8d4a6e8d81f171b2bd4be5265b01384c74c" +uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +version = "2.5.10" + +[[deps.Pkg]] +deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +version = "1.9.0" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "9673d39decc5feece56ef3940e5dafba15ba0f81" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.1.2" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "7eb1686b4f04b82f96ed7a4ea5890a4f0c7a09f1" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.4.0" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[deps.Profile]] +deps = ["Printf"] +uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" + +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[deps.Random]] +deps = ["SHA", "Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.Random123]] +deps = ["Random", "RandomNumbers"] +git-tree-sha1 = "552f30e847641591ba3f39fd1bed559b9deb0ef3" +uuid = "74087812-796a-5b5d-8853-05524746bad3" +version = "1.6.1" + +[[deps.RandomNumbers]] +deps = ["Random", "Requires"] +git-tree-sha1 = "043da614cc7e95c703498a491e2c21f58a2b8111" +uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143" +version = "1.5.3" + +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "1.2.2" + +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.0" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.Scratch]] +deps = ["Dates"] +git-tree-sha1 = "30449ee12237627992a99d5e30ae63e4d78cd24a" +uuid = "6c6a2e73-6563-6170-7368-637461726353" +version = "1.2.0" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.Setfield]] +deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"] +git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac" +uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" +version = "1.1.1" + +[[deps.SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + +[[deps.SimpleTraits]] +deps = ["InteractiveUtils", "MacroTools"] +git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231" +uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d" +version = "0.9.4" + +[[deps.Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[deps.SolverCore]] +deps = ["LinearAlgebra", "NLPModels", "Printf"] +git-tree-sha1 = "9fb0712d597d6598857ae50b7744df17b1137b38" +uuid = "ff4d7338-4cf1-434d-91df-b86cb86fb843" +version = "0.3.7" + +[[deps.SparseArrays]] +deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.SparseDiffTools]] +deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "Graphs", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays", "VertexSafeGraphs"] +git-tree-sha1 = "e19ac47477c9a8fcca06dab5e5471417d5d9d723" +uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804" +version = "1.31.0" + +[[deps.SpecialFunctions]] +deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "ef28127915f4229c971eb43f3fc075dd3fe91880" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "2.2.0" + + [deps.SpecialFunctions.extensions] + SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" + + [deps.SpecialFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"] +git-tree-sha1 = "8982b3607a212b070a5e46eea83eb62b4744ae12" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.5.25" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.0" + +[[deps.Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +version = "1.9.0" + +[[deps.SuiteSparse]] +deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] +uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" + +[[deps.SuiteSparse_jll]] +deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"] +uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" +version = "5.10.1+6" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" +version = "1.10.0" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[deps.TimerOutputs]] +deps = ["ExprTools", "Printf"] +git-tree-sha1 = "f548a9e9c490030e545f72074a41edfd0e5bcdd7" +uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +version = "0.5.23" + +[[deps.TranscodingStreams]] +deps = ["Random", "Test"] +git-tree-sha1 = "9a6ae7ed916312b41236fcef7e0af564ef934769" +uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" +version = "0.9.13" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[deps.UnsafeAtomics]] +git-tree-sha1 = "6331ac3440856ea1988316b46045303bef658278" +uuid = "013be700-e6cd-48c3-b4a1-df204f14c38f" +version = "0.2.1" + +[[deps.UnsafeAtomicsLLVM]] +deps = ["LLVM", "UnsafeAtomics"] +git-tree-sha1 = "ea37e6066bf194ab78f4e747f5245261f17a7175" +uuid = "d80eeb9a-aca5-4d75-85e5-170c8b632249" +version = "0.1.2" + +[[deps.VertexSafeGraphs]] +deps = ["Graphs"] +git-tree-sha1 = "8351f8d73d7e880bfc042a8b6922684ebeafb35c" +uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" +version = "0.2.0" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.13+0" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.7.0+0" + +[[deps.nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" +version = "1.48.0+0" + +[[deps.p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" +version = "17.4.0+0" diff --git a/scripts/scopf/Project.toml b/scripts/scopf/Project.toml new file mode 100644 index 00000000..6a72bd65 --- /dev/null +++ b/scripts/scopf/Project.toml @@ -0,0 +1,17 @@ +[deps] +Argos = "ef244971-cf80-42b0-9762-2c2c832df5d5" +ArgosCUDA = "8946db8d-321b-4174-84ad-48e2f9b69c56" +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +ExaPF = "0cf0e50c-a82e-488f-ac7e-41ffdff1b8aa" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" +MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6" +MadNLPGPU = "d72a61cc-809d-412f-99be-fd81f4b8a598" +MadNLPHSL = "7fb6135f-58fe-4112-84ca-653cf5be0c77" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" + +[extras] +CUDA_Runtime_jll = "76a88914-d11a-5bdc-97e0-2f5a05c973a2" diff --git a/scripts/scopf/benchmark.jl b/scripts/scopf/benchmark.jl new file mode 100644 index 00000000..22f5e6eb --- /dev/null +++ b/scripts/scopf/benchmark.jl @@ -0,0 +1,102 @@ +include("common.jl") +include("config.jl") +include("screening.jl") + +function benchmark(model, lines, kkt; use_gpu=false, ntrials=3, options...) + blk = build_scopf_model(model, lines; use_gpu=use_gpu) + + ## Warm-up + solver = build_madnlp(blk, kkt; max_iter=1) + MadNLP.solve!(solver) + + ## Benchmark + t_total, t_callbacks, t_linear_solver = (0.0, 0.0, 0.0) + n_it, obj = 0, 0.0 + for _ in 1:ntrials + Argos.reset!(Argos.backend(blk)) + # Solve + solver = build_madnlp(blk, kkt; options...) + MadNLP.solve!(solver) + # Save results + t_total += solver.cnt.total_time + t_callbacks += solver.cnt.eval_function_time + t_linear_solver += solver.cnt.linear_solver_time + n_it += solver.cnt.k + obj += solver.obj_val + # Clean memory + use_gpu && refresh_memory() + end + + return ( + iters = n_it / ntrials, + obj = obj / ntrials, + total = t_total / ntrials, + callbacks = t_callbacks / ntrials, + linear_solver = t_linear_solver / ntrials, + ) +end + +function main(cases, kkt, max_contingencies, ntrials, save_results; use_gpu=false, src_contingencies=:exadata, options...) + # Setup + dev = use_gpu ? :cuda : :cpu + form = isa(kkt, Argos.BieglerReduction) ? :biegler : :full + + nexp = length(cases) + results = zeros(nexp, 7) + + for (i, case) in enumerate(cases) + @info "Case: $case" + datafile = joinpath(DATA, case) + model = ExaPF.PolarForm(datafile) + nbus = PS.get(model, PS.NumberOfBuses()) + nlines = PS.get(model, PS.NumberOfLines()) + + if src_contingencies == :exadata + instance = split(case, ".")[1] + lines = readdlm(joinpath(SCENARIOS, "$(instance)_onehour_60.Ctgs"), ',', Int)[:] + lines = filter_islanding(model.network, lines) + elseif src_contingencies == :generated + all_lines = collect(1:nlines) + # Basic screening to remove islanding + lines = filter_islanding(model.network, all_lines) + else + error("Only :exadata and :generated are supported for contingencies") + end + + if max_contingencies >= 1 && length(lines) > max_contingencies + lines = lines[1:max_contingencies] + end + nscens = length(lines) + println("Num contingencies: ", nscens) + + try + r = benchmark(model, lines, kkt; ntrials=ntrials, use_gpu=use_gpu, options...) + results[i, :] .= (nbus, nscens, r.iters, r.obj, r.total, r.callbacks, r.linear_solver) + catch err + println("Error when solving problem $case: $(err)") + end + end + + if save_results + output_dir = joinpath(dirname(@__FILE__), RESULTS_DIR) + if !isdir(output_dir) + mkdir(output_dir) + end + output_file = joinpath(output_dir, "benchmark_opf_$(form)_$(dev).txt") + writedlm(output_file, results) + end + return results +end + +main( + cases, + kkt, + max_contingencies, + ntrials, + save_results; + print_level=print_level, + src_contingencies=src_contingencies, + linear_solver=linear_solver, + use_gpu=use_gpu, +) + diff --git a/scripts/scopf/benchmark_biegler.jl b/scripts/scopf/benchmark_biegler.jl new file mode 100644 index 00000000..70cff2f5 --- /dev/null +++ b/scripts/scopf/benchmark_biegler.jl @@ -0,0 +1,93 @@ + +include("common.jl") + +function _build_madnlp_gpu(blk::Argos.OPFModel) + madnlp_options = Dict{Symbol, Any}() + madnlp_options[:linear_solver] = LapackGPUSolver + madnlp_options[:lapack_algorithm] = MadNLP.CHOLESKY + madnlp_options[:dual_initialized] = true + madnlp_options[:max_wall_time] = 1800.0 + madnlp_options[:max_iter] = 1000 + madnlp_options[:print_level] = MadNLP.DEBUG + madnlp_options[:tol] = 1e-5 + opt_ipm, opt_linear, logger = MadNLP.load_options(; madnlp_options...) + KKT = Argos.BieglerKKTSystem{Float64, CuVector{Int}, CuVector{Float64}, CuMatrix{Float64}} + return MadNLP.MadNLPSolver{Float64, KKT}(blk, opt_ipm, opt_linear; logger=logger) +end + +function benchmark_biegler(model, lines; ntrials=3) + blk = build_model(model, lines; use_gpu=true) + + t_total, t_callbacks, t_linear_solver = (0.0, 0.0, 0.0) + n_it = 0 + obj = 0 + + solver = _build_madnlp_gpu(blk) + # Warm-up + MadNLP.solve!(solver; max_iter=1) + + for _ in 1:ntrials + solver = _build_madnlp_gpu(blk) + MadNLP.solve!(solver) + t_total += solver.cnt.total_time + t_callbacks += solver.cnt.eval_function_time + t_linear_solver += solver.cnt.linear_solver_time + n_it += solver.cnt.k + obj += solver.obj_val + end + return ( + iters = n_it / ntrials, + obj = obj / ntrials, + total = t_total / ntrials, + callbacks = t_callbacks / ntrials, + linear_solver = t_linear_solver / ntrials, + ) +end + +function main(cases, max_contingencies, ntrials, save_results) + @assert CUDA.has_cuda() + + nexp = length(cases) + results = zeros(nexp, 7) + for (i, case) in enumerate(cases) + @info "Case: $case" + datafile = joinpath(DATA, case) + model = ExaPF.PolarForm("../matpower/data/case_ACTIVSg2000.m") + nbus = PS.get(model, PS.NumberOfBuses()) + + instance = split(case, ".")[1] + lines = readdlm(joinpath(SCENARIOS, "$(instance).Ctgs"), ',', Int)[:] + lines = screen_contingencies(model.network, lines) + if max_contingencies >= 1 && length(lines) > max_contingencies + lines = lines[1:max_contingencies] + end + nscens = length(lines) + println("Num contingencies: ", nscens) + + try + r = benchmark_biegler(model, lines; ntrials=ntrials) + results[i, :] .= (nbus, nscens, r.iters, r.obj, r.total, r.callbacks, r.linear_solver) + catch + println("fail to solve problem $case.") + end + refresh_memory() + end + + if save_results + output_dir = joinpath(dirname(@__FILE__), RESULTS_DIR) + if !isdir(output_dir) + mkdir(output_dir) + end + time_now = now() + hash = "$(today())-$(hour(time_now))h$(minute(time_now))" + output_file = joinpath(output_dir, "benchmark_biegler_scopf_$(hash).txt") + writedlm(output_file, results) + end +end + +cases = [ + "case_ACTIVSg2000.m", +] + +main(cases, 1, 3, true) + diff --git a/scripts/scopf/benchmark_full_space.jl b/scripts/scopf/benchmark_full_space.jl new file mode 100644 index 00000000..1ae6fbec --- /dev/null +++ b/scripts/scopf/benchmark_full_space.jl @@ -0,0 +1,91 @@ + +include("common.jl") + +function _build_madnlp_cpu(blk::Argos.OPFModel) + return MadNLP.MadNLPSolver( + blk; + linear_solver=Ma27Solver, + dual_initialized=true, + max_wall_time=1800.0, + max_iter=1000, + print_level=MadNLP.DEBUG, + tol=1e-5, + ) +end + +function benchmark_extensive(model, lines; ntrials=3, gpu_ad=false) + blk = build_model(model, lines; use_gpu=gpu_ad) + + t_total, t_callbacks, t_linear_solver = (0.0, 0.0, 0.0) + n_it = 0 + obj = 0 + + solver = _build_madnlp_cpu(blk) + # Warm-up + MadNLP.solve!(solver; max_iter=1) + + for _ in 1:ntrials + solver = _build_madnlp_cpu(blk) + MadNLP.solve!(solver) + t_total += solver.cnt.total_time + t_callbacks += solver.cnt.eval_function_time + t_linear_solver += solver.cnt.linear_solver_time + n_it += solver.cnt.k + obj += solver.obj_val + end + return ( + iters = n_it / ntrials, + obj = obj / ntrials, + total = t_total / ntrials, + callbacks = t_callbacks / ntrials, + linear_solver = t_linear_solver / ntrials, + ) +end + +function main(cases, max_contingencies, ntrials, save_results) + gpu_ad = (CUDA.has_cuda()) ? true : false + + nexp = length(cases) + results = zeros(nexp, 7) + for (i, case) in enumerate(cases) + @info "Case: $case" + datafile = joinpath(DATA, case) + model = ExaPF.PolarForm(datafile) + nbus = PS.get(model, PS.NumberOfBuses()) + + instance = split(case, ".")[1] + lines = readdlm(joinpath(SCENARIOS, "$(instance).Ctgs"), ',', Int)[:] + lines = screen_contingencies(model.network, lines) + if max_contingencies >= 1 && length(lines) > max_contingencies + lines = lines[1:max_contingencies] + end + nscens = length(lines) + println("Num contingencies: ", nscens) + + try + r = benchmark_extensive(model, lines; ntrials=ntrials, gpu_ad=gpu_ad) + results[i, :] .= (nbus, nscens, r.iters, r.obj, r.total, r.callbacks, r.linear_solver) + catch + println("fail to solve problem $case.") + end + refresh_memory() + end + + if save_results + output_dir = joinpath(dirname(@__FILE__), RESULTS_DIR) + if !isdir(output_dir) + mkdir(output_dir) + end + time_now = now() + hash = "$(today())-$(hour(time_now))h$(minute(time_now))" + output_file = joinpath(output_dir, "benchmark_ma27_scopf_$(gpu_ad)_$(hash).txt") + writedlm(output_file, results) + end +end + +cases = [ + "case_ACTIVSg2000.m", +] + +main(cases, 1, 1, true) + diff --git a/scripts/scopf/common.jl b/scripts/scopf/common.jl new file mode 100644 index 00000000..e49048ed --- /dev/null +++ b/scripts/scopf/common.jl @@ -0,0 +1,119 @@ +using Dates +using DelimitedFiles +using Graphs +using LazyArtifacts +using LinearAlgebra +using Printf +using Random + +using NLPModels +using Argos +using ExaPF +using MadNLP + +# HSL +using MadNLPHSL + +# GPU +using CUDA +using KernelAbstractions +using ArgosCUDA +using MadNLPGPU + +const PS = ExaPF.PowerSystem + +const DATA = joinpath(artifact"ExaData", "ExaData") +const MATPOWER = "/home/fpacaud/dev/matpower/data" +const SCENARIOS = joinpath(artifact"ExaData", "ExaData", "mp_demand") +RESULTS_DIR = "results" + +if CUDA.has_cuda() + CUDA.allowscalar(false) +end + +function refresh_memory() + GC.gc(true) + CUDA.has_cuda() && CUDA.reclaim() + return +end + +function init_model!(blk) + x0 = NLPModels.get_x0(blk) + nnzj = NLPModels.get_nnzj(blk) + jac = zeros(nnzj) + NLPModels.jac_coord!(blk, x0, jac) + return +end + +function build_scopf_model(model, lines; use_gpu=false) + contingencies = [ExaPF.LineContingency(l) for l in lines] + nblks = length(contingencies) + 1 + if use_gpu + model_gpu = PolarForm(model, CUDABackend()) + nlp = Argos.StochEvaluator(model_gpu, nblks; contingencies=contingencies) + blk = Argos.OPFModel(Argos.bridge(nlp)) + else + nlp = Argos.StochEvaluator(model, nblks; contingencies=contingencies) + blk = Argos.OPFModel(nlp) + end + init_model!(blk) + return blk +end + +function build_madnlp( + blk::Argos.OPFModel, + ::Argos.FullSpace; + max_iter=max_iter, + dual_initialized=true, + tol=1e-5, + print_level=MadNLP.ERROR, + linear_solver=Ma27Solver, +) + return MadNLP.MadNLPSolver(blk; max_iter=max_iter, dual_initialized=dual_initialized, tol=tol, print_level=print_level, linear_solver=linear_solver) +end + +function build_madnlp( + blk::Argos.OPFModel, + ::Argos.BieglerReduction; + max_iter=max_iter, + dual_initialized=true, + tol=1e-5, + print_level=MadNLP.ERROR, + linear_solver=nothing, +) + madnlp_options = Dict{Symbol, Any}() + madnlp_options[:linear_solver] = LapackGPUSolver + madnlp_options[:lapack_algorithm] = MadNLP.CHOLESKY + madnlp_options[:dual_initialized] = dual_initialized + madnlp_options[:max_iter] = max_iter + madnlp_options[:print_level] = print_level + madnlp_options[:tol] = tol + opt_ipm, opt_linear, logger = MadNLP.load_options(; madnlp_options...) + KKT = Argos.BieglerKKTSystem{Float64, CuVector{Int}, CuVector{Float64}, CuMatrix{Float64}} + return MadNLP.MadNLPSolver{Float64, KKT}(blk, opt_ipm, opt_linear; logger=logger) +end + +# # Custom scaling +function MadNLP.scale_objective(nlp::Argos.OPFModel, grad::AbstractVector; max_gradient=1e-8) + return min(1, max_gradient / norm(grad, Inf)) +end + +function MadNLP.scale_constraints!( + nlp::Argos.OPFModel, + con_scale::AbstractVector, + jac::AbstractMatrix; + max_gradient=1e-8, +) + blk = Argos.backend(nlp) + ncons = length.(blk.constraints.exprs) + cnt = cumsum(ncons) + + # Powerflow + con_scale[1:cnt[1]] .= 1e-0 + # Power generation + con_scale[cnt[1]+1:cnt[2]] .= 1e0 + # Line flows + con_scale[cnt[2]+1:cnt[3]] .= 1e0 + return +end + diff --git a/scripts/scopf/config.jl b/scripts/scopf/config.jl new file mode 100644 index 00000000..19349ca8 --- /dev/null +++ b/scripts/scopf/config.jl @@ -0,0 +1,53 @@ +#= + Config file to run SCOPF benchmark. +=# + +# OPF instances +cases = [ + # "case9.m", + "case118.m", + # "case300.m", + "case1354pegase.m", + # "case2869pegase.m", + # "case9241pegase.m", + # "case1951rte.m", + # "case_ACTIVSg500.m", + # "case_ACTIVSg2000.m", + # # "case_ACTIVSg10k.m", +] + +# Source of contingencies (:exadata or :generated) +src_contingencies = :exadata + +# Maximum contingencies allowable. +max_contingencies = 10 + +# Contingencies +case_contingencies = [4:4:16, 4:4:16] + +# Number of trial runs to estimate running time. +ntrials = 3 + +# Save results on disk? +save_results = true + +# Should we use the GPU to evaluate the derivatives? +use_gpu = true + +# KKT system (FullSpace() or BieglerReduction()) +kkt = Argos.BieglerReduction() + +# KKT linear solver (Ma27Solver, Ma57Solver, LapackGPUSolver) +linear_solver = LapackGPUSolver + +# max iterations +max_iter = 100 + +# Verbose level +verbose = true +print_level = if verbose + MadNLP.DEBUG +else + MadNLP.ERROR +end + diff --git a/scripts/scopf/corrective_scopf.jl b/scripts/scopf/corrective_scopf.jl new file mode 100644 index 00000000..9b1b15c9 --- /dev/null +++ b/scripts/scopf/corrective_scopf.jl @@ -0,0 +1,119 @@ + +using Revise +using Random +using NLPModels +using ExaPF +using Argos +using MadNLP +using MadNLPHSL + +function generate_loads(model, nscen, magnitude) + nbus = get(model, ExaPF.PS.NumberOfBuses()) + stack = ExaPF.NetworkStack(model) + pload_det = stack.pload |> Array + qload_det = stack.qload |> Array + + # pload_det .*= 2.12 + # qload_det .*= 2.12 + + has_load = (pload_det .> 0) + + Random.seed!(1) + pload = magnitude .* (randn(nbus, nscen) .* has_load) .+ pload_det + qload = magnitude .* (randn(nbus, nscen) .* has_load) .+ qload_det + return pload, qload +end + +function instantiate!(blk::Argos.OPFModel) + x0 = NLPModels.get_x0(blk) + nnzj = NLPModels.get_nnzj(blk) + jac = zeros(nnzj) + NLPModels.jac_coord!(blk, x0, jac) + return +end + +function solve_corrective_opf(case, nscen) + model = ExaPF.PolarForm(case) + pl, ql = generate_loads(model, nscen, 0.05) + ev = Argos.CorrectiveEvaluator(model, pl, ql; line_constraints=true, epsilon=1e-4) + blk = Argos.OPFModel(ev) + instantiate!(blk) + solver = MadNLP.MadNLPSolver( + blk; + dual_initialized=true, + linear_solver=Ma27Solver, + max_iter=250, + print_level=MadNLP.DEBUG, + tol=1e-5, + ) + MadNLP.solve!(solver) + return ( + solver=solver, + model=ev, + ) +end + +function solve_corrective_tracking_opf(case, nscen, stack_ref) + model = ExaPF.PolarForm(case) + pl, ql = generate_loads(model, nscen, 0.05) + ev = Argos.CorrectiveEvaluator(model, pl, ql; line_constraints=true, tracking=true, stack_ref=stack_ref) + blk = Argos.OPFModel(ev) + instantiate!(blk) + solver = MadNLP.MadNLPSolver( + blk; + dual_initialized=true, + linear_solver=Ma27Solver, + max_iter=250, + print_level=MadNLP.DEBUG, + tol=1e-5, + ) + MadNLP.solve!(solver) + return ( + solver=solver, + model=ev, + ) +end + +function analyse_solution(sol) + nbus = get(Argos.model(sol.model), ExaPF.PS.NumberOfBuses()) + ngen = get(Argos.model(sol.model), ExaPF.PS.NumberOfGenerators()) + nx = ExaPF.number(Argos.model(sol.model), State()) + + y = sol.solver.y + active_cons = [1:nx*nscen; nx*nscen .+ findall(abs.(y[1+nx*nscen:end]) .> 1e-5)] + + # Jacobian + J = sol.model.jac.J[active_cons, :] + Gx = sol.model.jac.J[1:nx*nscen, 1:nx*nscen] + Gxu = J[1:nx*nscen, :] + + sJ = svd(Array(J)).S + sgx = svd(Array(Gx)).S + + println() + println("* Min singular value Gx: ", sgx[end]) + println("* Min singular value J: ", sJ[end]) +end + + +nscen = 1 +mode = :corrective + +case = "../matpower/data/case9.m" +sol = if mode == :corrective + solve_corrective_opf(case, nscen) +elseif mode == :tracking + # For this example, we took as a reference MATPOWER's base solution + model = ExaPF.PolarForm(case) + stack_ref = ExaPF.NetworkStack(model) + solve_corrective_tracking_opf(case, nscen, stack_ref) +end + +# Fetch solution +stack = sol.model.stack +ngen = get(Argos.model(sol.model), ExaPF.PS.NumberOfGenerators()) +# Adjustment +delta = stack.vuser[1:nscen] +# Reference setpoint for generators +pg_ref = stack.vuser[nscen+1:nscen+ngen] + diff --git a/scripts/scopf/ipopt.jl b/scripts/scopf/ipopt.jl new file mode 100644 index 00000000..fce8c34e --- /dev/null +++ b/scripts/scopf/ipopt.jl @@ -0,0 +1,27 @@ + +include("common.jl") + +using Ipopt +using MathOptInterface +const MOI = MathOptInterface + +# using HSL_jll + +datafile = "/home/fpacaud/dev/matpower/data/case118.m" +lines = readdlm(joinpath(SCENARIOS, "case118_onehour_60.Ctgs"), ',', Int)[:] + +model = ExaPF.PolarForm(datafile) +nbus = PS.get(model, PS.NumberOfBuses()) +nlines = PS.get(model, PS.NumberOfLines()) + +lines = filter_islanding(model.network, lines)[1:8] +blk = build_scopf_model(model, lines; use_gpu=true) + +optimizer = Ipopt.Optimizer() +MOI.set(optimizer, MOI.RawOptimizerAttribute("print_level"), 5) +MOI.set(optimizer, MOI.RawOptimizerAttribute("tol"), 1e-5) +# MOI.set(optimizer, MOI.RawOptimizerAttribute("hsllib"), HSL_jll.libhsl_path) +# MOI.set(optimizer, MOI.RawOptimizerAttribute("linear_solver"), "ma57") + +solution = Argos.optimize!(optimizer, blk.nlp) +MOI.empty!(optimizer) diff --git a/scripts/scopf/kkt.jl b/scripts/scopf/kkt.jl new file mode 100644 index 00000000..69443a7c --- /dev/null +++ b/scripts/scopf/kkt.jl @@ -0,0 +1,173 @@ +include("common.jl") +include("config.jl") +include("screening.jl") + +function analyse_memory(A::CUSPARSE.CuSparseMatrixCSR) + tot_mem = 0 + tot_mem += sizeof(A.nzVal) + tot_mem += sizeof(A.rowPtr) + tot_mem += sizeof(A.colVal) + return tot_mem +end + +function analyse_memory(F::ArgosCUDA.CUSOLVERRF.RFLU) + tot_mem = 0 + tot_mem += analyse_memory(F.M) + tot_mem += analyse_memory(F.P) + tot_mem += analyse_memory(F.Q) + tot_mem += sizeof(F.T) + tot_mem += sizeof(F.dsm.buffer) + tot_mem += sizeof(F.tsm.buffer) + tot_mem += sizeof(F.tsv.buffer) + return tot_mem +end + +function analyse_memory(kkt::Argos.BieglerKKTSystem) + tot_mem = 0 + # Buffers + tot_mem += sizeof(kkt._wj1) + tot_mem += sizeof(kkt._wx1) + tot_mem += sizeof(kkt._wx2) + tot_mem += sizeof(kkt._wxu1) + tot_mem += sizeof(kkt._wxu2) + tot_mem += sizeof(kkt._wxu3) + + # Reduced matrix + tot_mem += sizeof(kkt.aug_com) + + # Sparse transformation + tot_mem += sizeof(kkt.h_V) + tot_mem += sizeof(kkt.j_V) + tot_mem += analyse_memory(kkt.Gx) + tot_mem += analyse_memory(kkt.Gu) + tot_mem += analyse_memory(kkt.A) + tot_mem += sizeof(kkt.j_V) + + # IPM + tot_mem += sizeof(kkt.pr_diag) + tot_mem += sizeof(kkt.du_diag) + tot_mem += sizeof(kkt.con_scale) + tot_mem += sizeof(kkt.jacobian_scaling) + + # Condensed + tot_mem += sizeof(kkt.K.Σ) + tot_mem += sizeof(kkt.K.transperm) + tot_mem += analyse_memory(kkt.K.JtJ) + tot_mem += analyse_memory(kkt.K.Jt) + tot_mem += analyse_memory(kkt.K.W) + + # + tot_mem += analyse_memory(kkt.G_fac) + tot_mem += sizeof(kkt.reduction.z) + tot_mem += sizeof(kkt.reduction.ψ) + tot_mem += sizeof(kkt.reduction.tangents) + + return tot_mem +end + +function benchmark_kkt(model, lines, kkt; use_gpu=false, ntrials=3, options...) + use_gpu && refresh_memory() + blk = build_scopf_model(model, lines; use_gpu=use_gpu) + + ## Warm-up + solver = build_madnlp(blk, kkt; max_iter=1) + MadNLP.solve!(solver) + + ## Benchmark + t_build, t_factorization, t_backsolve = (0.0, 0.0, 0.0) + delta_err = 0.0 + n_it, obj = 0, 0.0 + for _ in 1:ntrials + t_build += CUDA.@elapsed begin + MadNLP.build_kkt!(solver.kkt) + end + t_factorization += CUDA.@elapsed begin + MadNLP.factorize!(solver.linear_solver) + end + t_backsolve += CUDA.@elapsed begin + MadNLP.solve_refine_wrapper!(solver, solver.d, solver.p) + end + + dsol = MadNLP.primal_dual(solver.d) + n = length(dsol) + psol = zeros(n) + + mul!(psol, solver.kkt, dsol) + + delta_err += norm(psol .- MadNLP.primal_dual(solver.p), Inf) + end + + memory_usage = if isa(kkt, Argos.BieglerReduction) + analyse_memory(solver.kkt) + else + 0 + end + println(memory_usage) + + return ( + build=t_build / ntrials, + factorization=t_factorization / ntrials, + backsolve=t_backsolve / ntrials, + accuracy=delta_err / ntrials, + memory=memory_usage, + ) +end + +function launch_kkt(cases, case_contingencies, kkt, ntrials, save_results; use_gpu=false, src_contingencies=:exadata, options...) + # Setup + dev = use_gpu ? :cuda : :cpu + form = isa(kkt, Argos.BieglerReduction) ? :biegler : :full + + nexp = sum(length.(case_contingencies)) + results = zeros(nexp, 7) + + i = 0 + for (case, contingencies) in zip(cases, case_contingencies) + datafile = joinpath(DATA, case) + model = ExaPF.PolarForm(datafile) + nbus = PS.get(model, PS.NumberOfBuses()) + nlines = PS.get(model, PS.NumberOfLines()) + if src_contingencies == :exadata + instance = split(case, ".")[1] + all_lines = readdlm(joinpath(SCENARIOS, "$(instance)_onehour_60.Ctgs"), ',', Int)[:] + all_lines = filter_islanding(model.network, all_lines) + elseif src_contingencies == :generated + all_lines = collect(1:nlines) + # Basic screening to remove islanding + all_lines = filter_islanding(model.network, all_lines) + else + error("Only :exadata and :generated are supported for contingencies") + end + + for ncont in contingencies + i += 1 + @info "Case: $case N: $ncont" + lines = all_lines[1:ncont] + r = benchmark_kkt(model, lines, kkt; ntrials=ntrials, use_gpu=use_gpu, options...) + results[i, :] .= (nbus, ncont, r.build, r.factorization, r.backsolve, r.accuracy, r.memory / 1024^2) + end + end + + if save_results + output_dir = joinpath(dirname(@__FILE__), RESULTS_DIR) + if !isdir(output_dir) + mkdir(output_dir) + end + output_file = joinpath(output_dir, "benchmark_kkt_$(form)_$(dev).txt") + writedlm(output_file, results) + end + return results +end + +launch_kkt( + cases, + case_contingencies, + kkt, + ntrials, + save_results; + print_level=print_level, + src_contingencies=src_contingencies, + linear_solver=linear_solver, + use_gpu=use_gpu, +) + diff --git a/scripts/scopf/preventive_scopf.jl b/scripts/scopf/preventive_scopf.jl new file mode 100644 index 00000000..1118676e --- /dev/null +++ b/scripts/scopf/preventive_scopf.jl @@ -0,0 +1,81 @@ + +using Revise + +include("common.jl") + +case = "1354pegase" +datafile = "../matpower/data/case$(case).m" +ctgs = readdlm(joinpath(SCENARIOS, "case$(case)_onehour_60.Ctgs"), ',', Int)[:] +# ctgs = readdlm(joinpath(SCENARIOS, "case$(case).Ctgs"), ',', Int)[:] + +function _build_madnlp_gpu(blk::Argos.OPFModel) + madnlp_options = Dict{Symbol, Any}() + madnlp_options[:linear_solver] = LapackGPUSolver + madnlp_options[:lapack_algorithm] = MadNLP.CHOLESKY + madnlp_options[:dual_initialized] = true + madnlp_options[:max_iter] = 250 + madnlp_options[:print_level] = MadNLP.DEBUG + madnlp_options[:tol] = 1e-5 + opt_ipm, opt_linear, logger = MadNLP.load_options(; madnlp_options...) + KKT = Argos.BieglerKKTSystem{Float64, CuVector{Int}, CuVector{Float64}, CuMatrix{Float64}} + return MadNLP.MadNLPSolver{Float64, KKT}(blk, opt_ipm, opt_linear; logger=logger) +end + +function test_1(datafile, line) + network = PS.PowerNetwork(datafile; remove_lines=[lines]) + network = PS.PowerNetwork(datafile) + println(line) + network.lines.Yff[line] = 0.0im + network.lines.Yft[line] = 0.0im + network.lines.Ytf[line] = 0.0im + network.lines.Ytt[line] = 0.0im + + model = ExaPF.PolarForm(network, ExaPF.CPU()) + stack = ExaPF.NetworkStack(model) + + conv = ExaPF.run_pf(model, stack; verbose=1) + return conv +end + +function test_2(model, lines) + nlp = build_model(model, lines; use_gpu=false) + ips = MadNLP.MadNLPSolver(nlp; linear_solver=Ma27Solver, tol=1e-5, print_level=MadNLP.DEBUG) + MadNLP.solve!(ips) + return ips +end + +function test_3(model, lines) + model = ExaPF.PolarForm(datafile, CPU()) + nlp = build_model(model, lines; use_gpu=true) + ips = _build_madnlp_gpu(nlp) + MadNLP.solve!(ips) + return ips +end + +function test_4(datafile, lines) + contingencies = [ExaPF.LineContingency(l) for l in lines] + nblk = length(contingencies) + 1 + blk = ExaPF.BlockPolarForm(datafile, nblk) + stack = ExaPF.NetworkStack(blk) + + basis = ExaPF.PolarBasis(blk) + pf = ExaPF.PowerFlowBalance(blk, contingencies) ∘ basis + + jac_pf = ExaPF.ArrowheadJacobian(blk, pf, State()) + ExaPF.jacobian!(jac_pf, stack) + # Solve power flow equations + solver = ExaPF.NewtonRaphson(; verbose=1) + conv = ExaPF.nlsolve!(solver, jac_pf, stack) + return conv +end + +model = ExaPF.PolarForm(datafile, CPU()) +nlines = PS.get(model, PS.NumberOfLines()) +@info "Screen contingencies" + + +# sctgs = g[1:4] +@info "Solve" +# conv = test_2(model, sctgs) +# conv = test_3(model, sctgs) +# ips = test_1(datafile, 4) diff --git a/scripts/scopf/screening.jl b/scripts/scopf/screening.jl new file mode 100644 index 00000000..1c42428a --- /dev/null +++ b/scripts/scopf/screening.jl @@ -0,0 +1,74 @@ + +const CONTINGENCIES_DIR = joinpath(dirname(@__FILE__), "contingencies") + +function filter_islanding(network, contingencies) + nbus, nlines = PS.get(network, PS.NumberOfBuses()), PS.get(network, PS.NumberOfLines()) + + graph = Graph() + add_vertices!(graph, nbus) + + f, t = network.lines.from_buses, network.lines.to_buses + for (i, j) in zip(f, t) + add_edge!(graph, i, j) + end + @assert is_connected(graph) + screened = Int[] + for (k, c) in enumerate(contingencies) + i, j = f[c], t[c] + rem_edge!(graph, i, j) + if is_connected(graph) + push!(screened, c) + end + add_edge!(graph, i, j) + end + return screened +end + +function filter_infeasible(network, contingencies) + is_feasible = Int[] + for line in contingencies + println("Screen contingency $(line) [Total: $(length(contingencies))]") + # Drop line from power network + post = PS.PowerNetwork(network) + post.lines.Yff[line] = 0.0im + post.lines.Yft[line] = 0.0im + post.lines.Ytf[line] = 0.0im + post.lines.Ytt[line] = 0.0im + + model = ExaPF.PolarForm(post, ExaPF.CPU()) + nlp = Argos.FullSpaceEvaluator(model) + + ips = MadNLP.MadNLPSolver( + Argos.OPFModel(nlp); + linear_solver=Ma27Solver, + max_iter=500, tol=1e-5, + print_level=MadNLP.ERROR, + ) + MadNLP.solve!(ips) + + if ips.status == MadNLP.SOLVE_SUCCEEDED + push!(is_feasible, line) + end + end + return is_feasible +end + +function screen_contingencies(model::PolarForm) + nlines = PS.get(model, PS.NumberOfLines()) + + all_lines = 1:nlines + cont1 = filter_islanding(model.network, all_lines) + cont2 = filter_infeasible(model.network, cont1) + ncont = length(cont2) + println("Keep $(ncont) contingencies out of $(nlines).") + return cont2 +end + +function generate_contingencies(case) + instance = split(case, ".")[1] + datafile = joinpath("/home/fpacaud/dev/matpower/data", case) + model = ExaPF.PolarForm(datafile, CPU()) + contingencies = screen_contingencies(model) + writedlm(joinpath(CONTINGENCIES_DIR, "$(instance).Ctgs"), contingencies) +end + diff --git a/scripts/scopf/tracking.jl b/scripts/scopf/tracking.jl new file mode 100644 index 00000000..4abf68df --- /dev/null +++ b/scripts/scopf/tracking.jl @@ -0,0 +1,47 @@ + + +using Printf +using MadNLP +using ExaPF +using Argos + +include("common.jl") +include("screening.jl") + +function tracking(model::ExaPF.PolarForm, lines, kkt, horizon; use_gpu=false) + blk = build_scopf_model(model, lines; use_gpu=use_gpu) + ev = Argos.backend(blk) + + # Warm-up + solver = build_madnlp(blk, kkt; print_level=MadNLP.DEBUG) + MadNLP.solve!(solver; max_iter=1) + + # Initial solve + solver = build_madnlp(blk, kkt; print_level=MadNLP.ERROR) + MadNLP.solve!(solver) + + cumulative_times = Float64[] + cumulative_iterations = Int[] + # Track solution (every 5mn) + for t in 1:horizon + # Update loads + ev.stack.pload .*= .99 + ev.stack.qload .*= .99 + x = Argos.initial(ev) + Argos.update!(ev, x) + copyto!(blk.meta.x0, x) + MadNLP.solve!(solver; mu_init=1e-6) + pg = sum(ev.stack.pgen) + push!(cumulative_times, solver.cnt.total_time) + push!(cumulative_iterations, solver.cnt.k) + @printf(" %4d %4d %12.5e %12.5e %10.4f %3d\n", t, solver.cnt.k, solver.obj_val, pg, solver.cnt.total_time, Int(solver.status)) + end + println("TOTAL TIME (s): ", solver.cnt.total_time) + return solver, cumulative_times, cumulative_iterations +end + +model = PolarForm("/home/fpacaud/dev/matpower/data/case1354pegase.m") +kkt = Argos.FullSpace() +lines = readdlm(joinpath(SCENARIOS, "case1354pegase_onehour_60.Ctgs"), ',', Int)[:] +lines = filter_islanding(model.network, lines)[1:16] +solver, cum_times, cum_it = tracking(model, lines, kkt, 13; use_gpu=true) diff --git a/src/Evaluators/Evaluators.jl b/src/Evaluators/Evaluators.jl index 70f63333..3861faf2 100644 --- a/src/Evaluators/Evaluators.jl +++ b/src/Evaluators/Evaluators.jl @@ -308,6 +308,7 @@ include("scalers.jl") include("reduced_evaluator.jl") include("full_space_evaluator.jl") include("stoch_evaluator.jl") +include("corrective_evaluator.jl") include("slack_evaluator.jl") include("feasibility_evaluator.jl") include("bridge_evaluator.jl") diff --git a/src/Evaluators/corrective_evaluator.jl b/src/Evaluators/corrective_evaluator.jl new file mode 100644 index 00000000..155f5852 --- /dev/null +++ b/src/Evaluators/corrective_evaluator.jl @@ -0,0 +1,280 @@ +mutable struct CorrectiveEvaluator{T, VI, VT, MT, JacCons, HessLag} <: AbstractNLPEvaluator + model::ExaPF.PolarFormRecourse{T, VI, VT, MT} + nscen::Int + nx::Int + nu::Int + blk_mapu::VI + blk_mapx::VI + mapz::VI + + # Expressions + basis::AutoDiff.AbstractExpression + costs::AutoDiff.AbstractExpression + constraints::AutoDiff.AbstractExpression + + _obj::VT + _cons::VT + _grad_control::VT + _multipliers::VT + + x_min::VT + x_max::VT + g_min::VT + g_max::VT + + # Cache + stack::ExaPF.NetworkStack + ∂stack::ExaPF.NetworkStack + + jac::JacCons + hess::HessLag + + map2tril::VI +end + +#= + Ordering: [x1, x2, ..., xN, u] +=# +function CorrectiveEvaluator( + model::PolarForm{T, VI, VT, MT}, nscen::Int; + contingencies=ExaPF.AbstractContingency[], + relax_eps=1e-6, line_constraints=true, epsilon=1e-4, + tracking=false, stack_ref=nothing, +) where {T, VI, VT, MT} + blk_model = ExaPF.PolarFormRecourse(model, nscen) + + nx = ExaPF.number(blk_model, State()) + nu = ExaPF.number(blk_model, Control()) + + # Block mappings + blk_mapxu = ExaPF.mapping(blk_model, ExaPF.AllVariables(), nscen) + blk_mapx = ExaPF.mapping(blk_model, State(), nscen) + blk_mapu = ExaPF.mapping(blk_model, Control(), nscen) + + mapz = [blk_mapx; blk_mapu[1:nu]] + + # Expressions + basis = ExaPF.PolarBasis(blk_model) + costs = if tracking + @assert isa(stack_ref, ExaPF.NetworkStack) + ExaPF.TrackingCost(blk_model, stack_ref) + else + ExaPF.QuadraticCost(blk_model) + end + + if length(contingencies) >= 1 + constraints_expr = Any[ + ExaPF.PowerFlowRecourse(blk_model, contingencies; epsilon=epsilon), + ExaPF.ReactivePowerBounds(blk_model, contingencies), + ExaPF.LineFlows(blk_model, contingencies), + ] + else + constraints_expr = Any[ + ExaPF.PowerFlowRecourse(blk_model; epsilon=epsilon), + ExaPF.ReactivePowerBounds(blk_model), + ExaPF.LineFlows(blk_model), + ] + end + + constraints = ExaPF.MultiExpressions(constraints_expr) + m = length(constraints) + + stack = ExaPF.NetworkStack(blk_model) + ∂stack = ExaPF.NetworkStack(blk_model) + + # Buffers + obj = VT(undef, nscen) + cons = VT(undef, m) + grad_control = VT(undef, nu * nscen) + y = VT(undef, m + nscen) + + s_min, s_max = ExaPF.bounds(blk_model, stack) + x_min, x_max = s_min[mapz], s_max[mapz] + + g_min, g_max = ExaPF.bounds(blk_model, constraints) + + # Remove bounds below a given threshold + g_max = min.(g_max, 1e5) + # Remove equalities + ggl = @view g_min[nx*nscen+1:end] + ggu = @view g_max[nx*nscen+1:end] + idx_eq = findall(ggl .== ggu) + if length(idx_eq) > 0 + println("[Argos] Elastic relaxation of $(length(idx_eq)) operational eq. constraints") + ggu[idx_eq] .+= 1e-6 + ggl[idx_eq] .-= 1e-6 + end + + jac = ExaPF.ArrowheadJacobian(blk_model, constraints ∘ basis, ExaPF.AllVariables()) + lagrangian_expr = [costs; constraints_expr] + lagrangian = ExaPF.MultiExpressions(lagrangian_expr) + hess = ExaPF.ArrowheadHessian(blk_model, lagrangian ∘ basis, ExaPF.AllVariables()) + + map2tril = tril_mapping(hess.H) + + return CorrectiveEvaluator( + blk_model, nscen, nx*nscen, nu, VI(blk_mapu), VI(blk_mapx), VI(mapz), + basis, costs, constraints, + obj, cons, grad_control, y, x_min, x_max, g_min, g_max, + stack, ∂stack, jac, hess, map2tril, + ) +end + +function CorrectiveEvaluator( + model::PolarForm, + ploads::Array{Float64, 2}, + qloads::Array{Float64, 2}; + options... +) + @assert size(ploads, 2) == size(qloads, 2) + nscen = size(ploads, 2) + evaluator = CorrectiveEvaluator(model, nscen; options...) + # Set loads as parameters in the model + ExaPF.set_params!(evaluator.stack, ploads, qloads) + ExaPF.set_params!(evaluator.∂stack, ploads, qloads) + ExaPF.set_params!(evaluator.jac, evaluator.stack) + ExaPF.set_params!(evaluator.hess, evaluator.stack) + return evaluator +end + +function CorrectiveEvaluator( + datafile::String, + pload::Array{Float64, 2}, + qload::Array{Float64, 2}; + device=ExaPF.CPU(), + options... +) + return CorrectiveEvaluator(ExaPF.PolarForm(datafile, device), pload, qload; options...) +end + +model(nlp::CorrectiveEvaluator) = nlp.model +backend(nlp::CorrectiveEvaluator) = nlp + +n_variables(nlp::CorrectiveEvaluator) = nlp.nx + nlp.nu +n_constraints(nlp::CorrectiveEvaluator) = length(nlp.g_min) + +constraints_type(::CorrectiveEvaluator) = :inequality +has_hessian(nlp::CorrectiveEvaluator) = true +has_hessian_lagrangian(nlp::CorrectiveEvaluator) = true + +# Initial position +function initial(nlp::CorrectiveEvaluator{T,VI,VT,MT}) where {T,VI,VT,MT} + x = VT(undef, n_variables(nlp)) + copyto!(x, nlp.stack, nlp.mapz) + return x +end + +# Bounds +bounds(nlp::CorrectiveEvaluator, ::Variables) = (nlp.x_min, nlp.x_max) +bounds(nlp::CorrectiveEvaluator, ::Constraints) = (nlp.g_min, nlp.g_max) + +## Callbacks +# Update buffer with new state and control +function update!(nlp::CorrectiveEvaluator, x) + copyto!(nlp.stack, nlp.mapz, x) + u = view(x, nlp.nx+1:nlp.nx+nlp.nu) + ExaPF.blockcopy!(nlp.stack, nlp.blk_mapu, u) + # Full forward pass + nlp.basis(nlp.stack.ψ, nlp.stack) + nlp.costs(nlp._obj, nlp.stack) + nlp.constraints(nlp._cons, nlp.stack) + return true +end + +objective(nlp::CorrectiveEvaluator, x) = 1.0 / nlp.nscen * sum(nlp._obj) +constraint!(nlp::CorrectiveEvaluator, c, x) = copyto!(c, nlp._cons) + +### +# First-order code +#### +function gradient!(nlp::CorrectiveEvaluator, g, x) + ExaPF.empty!(nlp.∂stack) + ExaPF.adjoint!(nlp.costs, nlp.∂stack, nlp.stack, 1.0) + ExaPF.adjoint!(nlp.basis, nlp.∂stack, nlp.stack, nlp.∂stack.ψ) + copyto!(g, nlp.∂stack, nlp.mapz) + copyto!(nlp._grad_control, nlp.∂stack, nlp.blk_mapu) + + gu = view(g, nlp.nx+1:nlp.nx+nlp.nu) + sum!(gu, reshape(nlp._grad_control, nlp.nu, nlp.nscen)) + g .*= 1.0 / nlp.nscen + return +end + +function jprod!(nlp::CorrectiveEvaluator, jv, x, v) + ExaPF.jacobian!(nlp.jac, nlp.stack) + mul!(jv, nlp.jac.J, v) + return +end + +function jtprod!(nlp::CorrectiveEvaluator, jv, x, v) + ExaPF.jacobian!(nlp.jac, nlp.stack) + mul!(jv, nlp.jac.J', v) + return +end + +function jacobian_structure(nlp::CorrectiveEvaluator) + i, j, _ = findnz(nlp.jac.J) + return i, j +end + +function jacobian_coo!(nlp::CorrectiveEvaluator, jacval::AbstractVector, x) + ExaPF.jacobian!(nlp.jac, nlp.stack) + copyto!(jacval, nonzeros(nlp.jac.J)) + return +end + +### +# Second-order code +#### +function hessian_lagrangian_prod!( + nlp::CorrectiveEvaluator, hessvec, x, y, σ, w, +) + m = n_constraints(nlp)::Int + nlp._multipliers[1:nlp.nscen] .= σ / nlp.nscen + nlp._multipliers[nlp.nscen+1:m+nlp.nscen] .= y + + ExaPF.hessian!(nlp.hess, nlp.stack, nlp._multipliers) + mul!(hessvec, nlp.hess.H, w) + return +end + +function hessian_lagrangian_coo!(nlp::CorrectiveEvaluator, hess, x, y, σ) + n = n_variables(nlp)::Int + m = n_constraints(nlp) + nlp._multipliers[1:nlp.nscen] .= σ / nlp.nscen * 0.01 + copyto!(nlp._multipliers, nlp.nscen + 1, y, 1, m) + ExaPF.hessian!(nlp.hess, nlp.stack, nlp._multipliers) + # Keep only lower-triangular part + transfer2tril!(hess, nlp.hess.H, nlp.map2tril) + return +end + +# Return lower-triangular matrix +function hessian_structure(nlp::CorrectiveEvaluator) + i, j, _ = findnz(nlp.hess.H) + rows = Int[ix for (ix, jx) in zip(i, j) if jx <= ix] + cols = Int[jx for (ix, jx) in zip(i, j) if jx <= ix] + return rows, cols +end + +function Base.show(io::IO, nlp::CorrectiveEvaluator) + n = n_variables(nlp) + m = n_constraints(nlp) + println(io, "A CorrectiveEvaluator object") + println(io, " * device: ", nlp.model.device) + println(io, " * #vars: ", n) + println(io, " * #cons: ", m) +end + +function reset!(nlp::CorrectiveEvaluator) + # Reset buffer + fill!(nlp._obj, 0) + fill!(nlp._multipliers, 0) + fill!(nlp._cons, 0) + empty!(nlp.stack) + empty!(nlp.∂stack) + ExaPF.init!(nlp.model, nlp.stack; update_loads=false) + fill!(nonzeros(nlp.hess.H), 0) + return +end + diff --git a/src/Evaluators/stoch_evaluator.jl b/src/Evaluators/stoch_evaluator.jl index 05f0248d..ea3c049d 100644 --- a/src/Evaluators/stoch_evaluator.jl +++ b/src/Evaluators/stoch_evaluator.jl @@ -36,12 +36,14 @@ end Ordering: [x1, x2, ..., xN, u] =# function StochEvaluator( - model::PolarForm{T, VI, VT, MT}, ploads::Array{T, 2}, qloads::Array{T, 2}; + model::PolarForm{T, VI, VT, MT}, + nscen::Int; + contingencies=ExaPF.AbstractContingency[], relax_eps=1e-6, ) where {T, VI, VT, MT} - @assert size(ploads, 2) == size(qloads, 2) - nscen = size(ploads, 2) - + if length(contingencies) >= 1 + @assert nscen == length(contingencies) + 1 + end blk_model = ExaPF.BlockPolarForm(model, nscen) nx = ExaPF.number(model, State()) @@ -57,18 +59,24 @@ function StochEvaluator( # Expressions basis = ExaPF.PolarBasis(blk_model) costs = ExaPF.CostFunction(blk_model) - constraints_expr = [ - ExaPF.PowerFlowBalance(blk_model), - ExaPF.PowerGenerationBounds(blk_model), - ExaPF.LineFlows(blk_model), - ] + if length(contingencies) >= 1 + constraints_expr = [ + ExaPF.PowerFlowBalance(blk_model, contingencies), + ExaPF.PowerGenerationBounds(blk_model, contingencies), + ExaPF.LineFlows(blk_model, contingencies), + ] + else + constraints_expr = [ + ExaPF.PowerFlowBalance(blk_model), + ExaPF.PowerGenerationBounds(blk_model), + ExaPF.LineFlows(blk_model), + ] + end constraints = ExaPF.MultiExpressions(constraints_expr) m = length(constraints) stack = ExaPF.NetworkStack(blk_model) ∂stack = ExaPF.NetworkStack(blk_model) - ExaPF.set_params!(stack, ploads, qloads) - ExaPF.set_params!(∂stack, ploads, qloads) # Buffers obj = VT(undef, nscen) @@ -87,13 +95,12 @@ function StochEvaluator( ggu = @view g_max[nx*nscen+1:end] idx_eq = findall(ggl .== ggu) if length(idx_eq) > 0 - println("[Argos] Elastic relaxation of operational eq. constraints $(idx_eq)") + println("[Argos] Elastic relaxation of $(length(idx_eq)) operational eq. constraints") ggu[idx_eq] .+= 1e-6 ggl[idx_eq] .-= 1e-6 end jac = ExaPF.ArrowheadJacobian(blk_model, constraints ∘ basis, ExaPF.AllVariables()) - ExaPF.set_params!(jac, stack) lagrangian_expr = [costs; constraints_expr] lagrangian = ExaPF.MultiExpressions(lagrangian_expr) hess = ExaPF.ArrowheadHessian(blk_model, lagrangian ∘ basis, ExaPF.AllVariables()) @@ -108,9 +115,31 @@ function StochEvaluator( stack, ∂stack, jac, hess, map2tril, ) end -function StochEvaluator(datafile::String, pload, qload; device=ExaPF.CPU(), options...) +function StochEvaluator( + datafile::String, + pload::Array{Float64, 2}, + qload::Array{Float64, 2}; + device=ExaPF.CPU(), + options... +) return StochEvaluator(ExaPF.PolarForm(datafile, device), pload, qload; options...) end +function StochEvaluator( + model::PolarForm{T, VI, VT, MT}, + ploads::Array{T, 2}, + qloads::Array{T, 2}; + options... +) where {T, VI, VT, MT} + @assert size(ploads, 2) == size(qloads, 2) + nscen = size(ploads, 2) + evaluator = StochEvaluator(model, nscen; options...) + # Set loads as parameters in the model + ExaPF.set_params!(evaluator.stack, ploads, qloads) + ExaPF.set_params!(evaluator.∂stack, ploads, qloads) + ExaPF.set_params!(evaluator.jac, evaluator.stack) + ExaPF.set_params!(evaluator.hess, evaluator.stack) + return evaluator +end model(nlp::StochEvaluator) = nlp.model backend(nlp::StochEvaluator) = nlp diff --git a/src/KKT/auglag_kkt.jl b/src/KKT/auglag_kkt.jl index 5f438024..03c7a655 100644 --- a/src/KKT/auglag_kkt.jl +++ b/src/KKT/auglag_kkt.jl @@ -35,7 +35,7 @@ Supports only bound-constrained optimization problem (so no Jacobian). [PMSSA2022] Pacaud, François, Daniel Adrian Maldonado, Sungho Shin, Michel Schanen, and Mihai Anitescu. "A feasible reduced space method for real-time optimal power flow." Electric Power Systems Research 212 (2022): 108268. """ -struct MixedAuglagKKTSystem{T, VT, MT} <: MadNLP.AbstractKKTSystem{T, VT, MT} +struct MixedAuglagKKTSystem{T, VT, MT} <: MadNLP.AbstractKKTSystem{T, VT, MT, MadNLP.ExactHessian{T, VT}} aug::AbstractNLPEvaluator # for Auglag information n::Int m::Int diff --git a/src/KKT/reduced_newton.jl b/src/KKT/reduced_newton.jl index 6195da93..bd17095b 100644 --- a/src/KKT/reduced_newton.jl +++ b/src/KKT/reduced_newton.jl @@ -40,7 +40,7 @@ in the reduction algorithm. [PSSMA2022] Pacaud, François, Sungho Shin, Michel Schanen, Daniel Adrian Maldonado, and Mihai Anitescu. "Condensed interior-point methods: porting reduced-space approaches on GPU hardware." arXiv preprint arXiv:2203.11875 (2022). """ -struct BieglerKKTSystem{T, VI, VT, MT, SMT} <: MadNLP.AbstractReducedKKTSystem{T, VT, MT} +struct BieglerKKTSystem{T, VI, VT, MT, SMT} <: MadNLP.AbstractReducedKKTSystem{T, VT, MT, MadNLP.ExactHessian{T, VT}} K::HJDJ{VI,VT,SMT} Wref::SMT W::SMT @@ -164,7 +164,7 @@ function BieglerKKTSystem{T, VI, VT, MT}(nlp::OPFModel, ind_cons=MadNLP.get_inde ind_W_fixed, ind_W_fixed_diag = get_fixed_nnz(W, ind_fixed, true) # Buffers - etc = Dict{Symbol, Any}(:reduction_time=>0.0) + etc = Dict{Symbol, Any}(:reduction_time=>0.0, :cond=>Float64[]) return BieglerKKTSystem{T, VI, VT, MT, SMT}( K, Wref, W, J, A, Gx, Gu, mapA, mapGx, mapGu, @@ -344,6 +344,8 @@ function MadNLP.build_kkt!(kkt::BieglerKKTSystem{T, VI, VT, MT}) where {T, VI, V # Regularize final reduced matrix to ensure it is full-rank fixed_diag!(kkt.aug_com, kkt.ind_fixed .- nx, 1.0) + + return end function MadNLP.solve_refine_wrapper!( @@ -386,8 +388,6 @@ function MadNLP.solve_refine_wrapper!( Σd = view(kkt.du_diag, nx+1:m) α = view(kkt.con_scale, nx+1:m) - Λ = 1.0 ./ (Σd .- α.^2 ./ Σₛ) - # RHS r₁₂ = view(b, 1:nx+nu) r₁ = view(b, 1:nx) # / state @@ -403,13 +403,20 @@ function MadNLP.solve_refine_wrapper!( dλ = view(x, ips.n+1:ips.n+nx) # / equality cons dy = view(x, ips.n+nx+1:ips.n+m) # / inequality cons + Λ = Σₛ ./ (Σd .* Σₛ .- α.^2) + # Reduction (1) --- Condensed - vj .= Λ .* (r₅ .+ α .* r₃ ./ Σₛ) # v = (α Σₛ⁻¹ α)⁻¹ * (r₅ + α Σₛ⁻¹ r₃) - mul!(jv, kkt.A', vj, -1.0, 0.0) # jᵥ = Aᵀ v + # vj .= Λ .* (r₅ .+ α .* r₃ ./ Σₛ) # v = (α Σₛ⁻¹ α)⁻¹ * (r₅ + α Σₛ⁻¹ r₃) + vj .= (Σₛ .* r₅ .+ r₃) # v = (Σₛ r₅ + α r₃) + mul!(jv, kkt.A', vj, 1.0, 0.0) # jᵥ = Aᵀ v jv .+= r₁₂ # r₁₂ - Aᵀv # Reduction (2) --- Biegler sx1 .= r₄ # r₄ ldiv!(Gxi, sx1) # Gₓ⁻¹ r₄ + + tt = similar(sx1) + mul!(tt, kkt.Gx, sx1) + sx2 .= tx # tx = jv[1:nx] kvx .= sx1 ; kvu .= 0.0 mul!(kh, K, kv) # [Kₓₓ Gₓ⁻¹ r₄ ; Kᵤₓ Gₓ⁻¹ r₄ ] @@ -431,18 +438,29 @@ function MadNLP.solve_refine_wrapper!( dλ .= tx # tₓ mul!(kh, K, dxu) # Kₓₓ dₓ + Kₓᵤ dᵤ axpy!(-1.0, khx, dλ) # tₓ - Kₓₓ dₓ + Kₓᵤ dᵤ + + # TODO: SEGFAULT ldiv!(Gxi', dλ) # dₗ = Gₓ⁻ᵀ(tₓ - Kₓₓ dₓ + Kₓᵤ dᵤ) + # (2) Extract Condensed mul!(vj, kkt.A, dxu) # Aₓ dₓ + Aᵤ dᵤ - dy .= Λ .* (r₅ .- vj .+ α .* r₃ ./ Σₛ) - ds .= (r₃ .+ α .* dy) ./ Σₛ + # dy .= Λ .* (r₅ .- vj .+ α .* r₃ ./ Σₛ) + # ds .= (r₃ .+ α .* dy) ./ Σₛ + ds .= (vj .- r₅) + dy .= Σₛ .* ds .- r₃ - x[ips.ind_fixed] .= 0 + x[ips.ind_fixed] .= 0.0 copyto!(x_h, x) return solve_status end function MadNLP.set_aug_RR!(kkt::BieglerKKTSystem, ips::MadNLP.MadNLPSolver, RR::MadNLP.RobustRestorer) - copyto!(kkt.pr_diag, ips.zl./(ips.x.-ips.xl) .+ ips.zu./(ips.xu.-ips.x) .+ RR.zeta.*RR.D_R.^2) + x = MadNLP.full(ips.x) + xl = MadNLP.full(ips.xl) + xu = MadNLP.full(ips.xu) + zl = MadNLP.full(ips.zl) + zu = MadNLP.full(ips.zu) + copyto!(kkt.pr_diag, zl./(x.-xl) .+ zu./(xu.-x) .+ RR.zeta.*RR.D_R.^2) copyto!(kkt.du_diag, .-RR.pp./RR.zp .- RR.nn./RR.zn) end + diff --git a/test/Algorithms/MadNLP_wrapper.jl b/test/Algorithms/MadNLP_wrapper.jl index 539e4d76..ab4738b5 100644 --- a/test/Algorithms/MadNLP_wrapper.jl +++ b/test/Algorithms/MadNLP_wrapper.jl @@ -70,10 +70,11 @@ end "case30.m", "case57.m", ] + tol = 1e-6 datafile = joinpath(INSTANCES_DIR, case) options = Dict{Symbol, Any}( :dual_initialized=>true, - :tol=>1e-6, + :tol=>tol, :print_level=>MadNLP.ERROR, ) @testset "Reduce-then-linearize" begin @@ -81,16 +82,16 @@ end ips = _madnlp_default(nlp; options...) @test ips.status == MadNLP.SOLVE_SUCCEEDED ipd = _madnlp_dense_kkt(nlp; options...) - _test_results_match(ips, ipd; atol=1e-8) + _test_results_match(ips, ipd; atol=tol) ipc = _madnlp_condensed_kkt(nlp; options...) - _test_results_match(ips, ipc; atol=1e-8) + _test_results_match(ips, ipc; atol=tol) end @testset "Linearize-then-reduce" begin flp = Argos.FullSpaceEvaluator(datafile) ips = _madnlp_default(flp; options...) @test ips.status == MadNLP.SOLVE_SUCCEEDED ipb = _madnlp_biegler_kkt(flp; options...) - _test_results_match(ips, ipb; atol=1e-8) + _test_results_match(ips, ipb; atol=tol) @test ipb.kkt.Wref === flp.hess.H end end diff --git a/test/Evaluators/TestEvaluators.jl b/test/Evaluators/TestEvaluators.jl index c42e5a83..21927475 100644 --- a/test/Evaluators/TestEvaluators.jl +++ b/test/Evaluators/TestEvaluators.jl @@ -97,6 +97,20 @@ function runtests(datafile, device, AT) test_evaluator_hessian_lagrangian(stoch, device, AT) test_evaluator_sparse_callbacks(stoch, device, AT) end + @testset "Argos.CorrectiveEvaluator Interface" begin + nblocks = 5 + case_name = split(split(datafile, '/')[end], '.')[1] + demands = joinpath(artifact"ExaData", "ExaData", "mp_demand") + pload = readdlm(joinpath(demands, "$(case_name)_oneweek_168.Pd"))[:, 1:nblocks] ./ 100 + qload = readdlm(joinpath(demands, "$(case_name)_oneweek_168.Qd"))[:, 1:nblocks] ./ 100 + + stoch = Argos.CorrectiveEvaluator(datafile, pload, qload) + test_evaluator_api(stoch, device, AT) + test_evaluator_callbacks(stoch, device, AT) + # TODO: currently broken + # test_evaluator_hessian_lagrangian(stoch, device, AT) + # test_evaluator_sparse_callbacks(stoch, device, AT) + end end end diff --git a/test/Evaluators/api.jl b/test/Evaluators/api.jl index 30642d8d..cec3146d 100644 --- a/test/Evaluators/api.jl +++ b/test/Evaluators/api.jl @@ -132,7 +132,7 @@ function test_evaluator_hessian(nlp, device, M; rtol=1e-6) u0 = u |> Array hess_fd = FiniteDiff.finite_difference_hessian(reduced_cost, u0) - @test H * w == hv + @test myisapprox(H * w, hv, rtol=rtol) @test myisapprox(H, hess_fd.data, rtol=rtol) end @@ -279,7 +279,7 @@ function test_evaluator_sparse_callbacks(nlp, device, M; rtol=1e-6) @test J_s ≈ J_d rtol=rtol end -function test_evaluator_bridged(bdg, device, M; rtol=1e-6) +function test_evaluator_bridged(bdg, device, M; rtol=1e-8) nlp = Argos.backend(bdg) (n, m) = Argos.n_variables(nlp), Argos.n_constraints(nlp) @@ -303,7 +303,8 @@ function test_evaluator_bridged(bdg, device, M; rtol=1e-6) @test isa(u2, M) @test M(u1) == u2 - @test M(g1) == g2 @test M(c1) == c2 + # Gradient is less accurate as it involves the solution of a linear system + @test myisapprox(g1, g2; rtol=rtol) end diff --git a/test/runtests.jl b/test/runtests.jl index b7ae1721..2e87bd3a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,9 +19,8 @@ const CASES = ["case9.m", "case30.m"] ARCHS = Any[(CPU(), Array, SparseMatrixCSC)] if has_cuda_gpu() - using CUDAKernels using ArgosCUDA - CUDA_ARCH = (CUDADevice(), CuArray, nothing) + CUDA_ARCH = (CUDABackend(), CuArray, nothing) push!(ARCHS, CUDA_ARCH) end