Support for N-x contingencies

Artifacts.toml

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

src/Evaluators/ProxALEvalutor.jl

@@ -400,13 +400,22 @@ function optimize!(
         iteration()
 
         # Check convergence
-        if max(
+        iteration_limit = 0
+        for blk in opfBlockData.blkModel
+            solution = get_solution(blk)
+            if solution.status ∉ MOI_OPTIMAL_STATUSES
+                iteration_limit = 1
+            end
+        end
+        giteration_limit = comm_sum(iteration_limit, comm)
+
+        if (max(
             runinfo.maxviol_t[end],
             runinfo.maxviol_c[end],
             runinfo.maxviol_t_actual[end],
             runinfo.maxviol_c_actual[end],
             runinfo.maxviol_d[end]
-        ) <= algparams.tol
+        ) <= algparams.tol) && (giteration_limit == 0)
         # minviol = max(
         #     runinfo.maxviol_t[end],
         #     runinfo.maxviol_c[end],

src/OPF/opfdata.jl

@@ -240,7 +240,7 @@ function opf_loaddata(raw::RawData;
                       load_scale::Float64=1.0,
                       ramp_scale::Float64=0.0,
                       corr_scale::Float64=0.1,
-                      lineOff=Line(),
+                      lineOff=nothing,
                       genOff::Vector{Int}=Int[])
     #
     # load buses
@@ -278,15 +278,19 @@ function opf_loaddata(raw::RawData;
     #
     branch_arr = raw.branch_arr
     num_lines = size(branch_arr, 1)
-    lines_on = findall((branch_arr[:,11] .> 0) .& ((branch_arr[:,1].!=lineOff.from) .| (branch_arr[:,2].!=lineOff.to)) )
-    num_on   = length(lines_on)
-
-    if lineOff.from> 0 && lineOff.to > 0
-        # println("opf_loaddata: was asked to remove line from,to=", lineOff.from, ",", lineOff.to)
-    end
-    if length(findall(branch_arr[:,11].==0))>0
-        # println("opf_loaddata: ", num_lines-length(findall(branch_arr[:,11].>0)), " lines are off and will be discarded (out of ", num_lines, ")")
+    lines_on = findall((branch_arr[:,11] .> 0)) 
+    if !isa(lineOff, Nothing)
+        for line in lineOff
+            intersect!(lines_on, findall((branch_arr[:,11] .> 0) .& ((branch_arr[:,1].!=line.from) .| (branch_arr[:,2].!=line.to)) ))
+            if line.from> 0 && line.to > 0
+                # println("opf_loaddata: was asked to remove line from,to=", lineOff.from, ",", lineOff.to)
+            end
+            if length(findall(branch_arr[:,11].==0))>0
+                # println("opf_loaddata: ", num_lines-length(findall(branch_arr[:,11].>0)), " lines are off and will be discarded (out of ", num_lines, ")")
+            end
+        end
     end
+    num_on   = length(lines_on)
 
     lines = Array{Line}(undef, num_on)
     ncols_lines = size(branch_arr, 2)

src/OPF/opfmodel.jl

@@ -112,7 +112,7 @@ function opf_model_add_block_constraints(opfmodel::JuMP.Model, opfdata::OPFData,
         @views for t=1:T, k=1:K
             opfdata_c = (k == 1) ? opfdata :
                 opf_loaddata(rawdata;
-                    lineOff = opfdata.lines[rawdata.ctgs_arr[k - 1]],
+                    lineOff = opfdata.lines[rawdata.ctgs_arr[k - 1],:],
                     time_horizon_start = modelinfo.time_horizon_start + t - 1,
                     time_horizon_end = modelinfo.time_horizon_start + t - 1,
                     load_scale = modelinfo.load_scale,
@@ -130,7 +130,7 @@ function opf_model_add_block_constraints(opfmodel::JuMP.Model, opfdata::OPFData,
         @views for t=1:T, k=1:K
             opfdata_c = (k == 1) ? opfdata :
                 opf_loaddata(rawdata;
-                    lineOff = opfdata.lines[rawdata.ctgs_arr[k - 1]],
+                    lineOff = opfdata.lines[rawdata.ctgs_arr[k - 1],:],
                     time_horizon_start = modelinfo.time_horizon_start + t - 1,
                     time_horizon_end = modelinfo.time_horizon_start + t - 1,
                     load_scale = modelinfo.load_scale,

src/OPF/opfsolution.jl

@@ -303,10 +303,10 @@ function write(runinfo::ProxALProblem, nlp::AbstractNLPEvaluator, filename::Stri
     modelinfo = nlp.modelinfo
     x = runinfo.x
     λ = runinfo.λ
-    if isa(comm, MPI.Comm)
+    if isa(comm, MPI.Comm) && comm_ranks(comm) != 1
         info = MPI.Info()
         ff = h5open(filename, "w", comm, info)
-    elseif isa(comm, Nothing)
+    elseif isa(comm, Nothing) || comm_ranks(comm) == 1 
         ff = h5open(filename, "w")
     else
         error("Wrong type of comm")

src/backends.jl

@@ -112,7 +112,13 @@ function init!(block::EmptyBlockModel, algparams::AlgParams) end
 function set_objective!(block::EmptyBlockModel, algparams::AlgParams,
                         primal::AbstractPrimalSolution, dual::AbstractDualSolution)
 end
-function get_solution(block::EmptyBlockModel) end
+function get_solution(block::EmptyBlockModel)
+    # Empty block is always optimal
+    solution = (
+        status=MOI.OPTIMAL,
+    )
+    return solution
+end
 
 ### Implementation of JuMPBlockBackend
 struct JuMPBlockBackend <: AbstractBlockModel
@@ -210,7 +216,7 @@ function init!(block::JuMPBlockBackend, algparams::AlgParams)
         end
         opfdata_c = (j == 1) ? opfdata :
             opf_loaddata(block.rawdata;
-            lineOff = opfdata.lines[block.rawdata.ctgs_arr[j - 1]],
+            lineOff = opfdata.lines[block.rawdata.ctgs_arr[j - 1],:],
             time_horizon_start = modelinfo.time_horizon_start + t - 1,
             time_horizon_end = modelinfo.time_horizon_start + t - 1,
             load_scale = modelinfo.load_scale,
@@ -768,9 +774,9 @@ function get_solution(block::AdmmBlockBackend)
     k = block.k
     t = block.t
 
-    if status ∉ MOI_OPTIMAL_STATUSES
-        @warn("Block $blk[$k, $t] subproblem not solved to optimality. status: $status")
-    end
+    # if status ∉ MOI_OPTIMAL_STATUSES
+    #     @warn("Block $blk[$k, $t] subproblem not solved to optimality. status: $status")
+    # end
 
     s = ExaAdmmBackend.slack_values(model) |> Array
     solution = (

src/blocks.jl

@@ -75,13 +75,13 @@ function OPFBlocks(
         t, k;
         decompCtgs=false,
     )
-        lineOff = Line()
+        lineOff = nothing
         if length(rawdata.ctgs_arr) < k - 1
             error("Not enough contingencies in .ctg file while trying to read contingency $(k-1).")
         end
         if decompCtgs
             if k > 1
-                lineOff = opfdata.lines[rawdata.ctgs_arr[k - 1]]
+                lineOff = opfdata.lines[rawdata.ctgs_arr[k - 1,:]]
             end
         end
         data = opf_loaddata(

src/communication.jl

@@ -168,11 +168,11 @@ end
 Collective to reduce and return the sum over scalar `data`.
 
 """
-function comm_sum(data::Float64, comm::MPI.Comm)
+function comm_sum(data::Number, comm::MPI.Comm)
     return MPI.Allreduce(data, MPI.SUM, comm)
 end
 
-function comm_sum(data::Float64, comm::Nothing)
+function comm_sum(data::Number, comm::Nothing)
     return data
 end
 

test/n-2.jl

@@ -0,0 +1,111 @@
+using Test
+using ProxAL
+using DelimitedFiles, Printf
+using LinearAlgebra, JuMP
+using CatViews
+using CUDA
+using AMDGPU
+using MPI
+using LazyArtifacts
+
+use_MPI = !isempty(ARGS) && (parse(Bool, ARGS[1]) == true)
+use_MPI && MPI.Init()
+const DATA_DIR = joinpath(artifact"ExaData", "ExaData")
+# N-2 case9 files
+case = "case9_2"
+T = 2
+ramp_scale = 0.5
+load_scale = 1.0
+rtol = 1e-4
+
+# Load case
+case_file = joinpath(DATA_DIR, "mp_demand", "$(case).m")
+load_file = joinpath(DATA_DIR, "mp_demand", "$(case)")
+
+# Model/formulation settings
+modelinfo = ModelInfo()
+modelinfo.num_time_periods = T
+modelinfo.load_scale = load_scale
+modelinfo.ramp_scale = ramp_scale
+modelinfo.allow_obj_gencost = true
+modelinfo.allow_constr_infeas = false
+modelinfo.weight_freq_ctrl = 0.1
+modelinfo.time_link_constr_type = :penalty
+modelinfo.obj_scale = 1e-4
+
+# Algorithm settings
+algparams = AlgParams()
+algparams.verbose = 0
+
+solver_list = ["ExaAdmmCPU"]
+
+@testset "Test ProxAL on $(case)" begin
+    modelinfo.case_name = case
+
+    using Ipopt
+    backend = JuMPBackend()
+    algparams.optimizer =
+        optimizer_with_attributes(Ipopt.Optimizer,
+            "print_level" => Int64(algparams.verbose > 0)*5)
+    K = 1
+    modelinfo.num_ctgs = K
+    # ctgs_link = :preventive_penalty
+    OPTIMAL_OBJVALUE = []
+    OPTIMAL_PG = []
+    OPTIMAL_WT = []
+
+    algparams.iterlim = 10000
+    algparams.decompCtgs = true
+    modelinfo.ctgs_link_constr_type = :corrective_penalty
+    algparams.mode = :coldstart
+    algparams.tol = 1e-3
+    algparams.verbose = 0
+    nlp = ProxALEvaluator(case_file, load_file, modelinfo, algparams, backend, use_MPI ? MPI.COMM_WORLD : nothing)
+    runinfo = ProxAL.optimize!(nlp)
+    OPTIMAL_OBJVALUE = runinfo.objvalue[end]
+    OPTIMAL_PG = runinfo.x.Pg[:,1,:][:]
+    OPTIMAL_WT = runinfo.x.ωt[:]
+    @test runinfo.maxviol_c[end] <= algparams.tol
+    @test runinfo.maxviol_t[end] <= algparams.tol
+    @test runinfo.maxviol_c_actual[end] <= algparams.tol
+    @test runinfo.maxviol_t_actual[end] <= algparams.tol
+    @test runinfo.maxviol_d[end] <= algparams.tol
+    @test runinfo.iter <= algparams.iterlim
+    algparams.decompCtgs = false
+    algparams.verbose = 0
+    nlp = ProxALEvaluator(case_file, load_file, modelinfo, algparams, backend, use_MPI ? MPI.COMM_WORLD : nothing)
+    runinfo = ProxAL.optimize!(nlp)
+    @test_broken isapprox(runinfo.objvalue[end], OPTIMAL_OBJVALUE, rtol = 1e-2)
+    @test_broken isapprox(runinfo.x.Pg[:,1,:][:], OPTIMAL_PG, rtol = 1e-1)
+    @test isapprox(runinfo.maxviol_c[end], 0.0)
+    @test isapprox(runinfo.maxviol_c_actual[end], 0.0)
+    @test runinfo.maxviol_c_actual[end] <= algparams.tol
+    @test runinfo.maxviol_t_actual[end] <= algparams.tol
+    @test runinfo.maxviol_d[end] <= algparams.tol
+
+    backend = AdmmBackend()
+    algparams.tron_outer_iterlim=20
+    algparams.tron_inner_iterlim=800
+    algparams.tron_outer_eps=1e-6
+    algparams.decompCtgs = true
+    algparams.mode = :coldstart
+    algparams.iterlim = 10000
+    algparams.verbose = 0
+    algparams.init_opf = false
+    algparams.θ_c = 1.0
+    algparams.ρ_c = 1.0
+    # Reduce tolerance for N-2
+    algparams.tol = 1e-3
+    nlp = ProxALEvaluator(case_file, load_file, modelinfo, algparams, backend, use_MPI ? MPI.COMM_WORLD : nothing)
+    runinfo = ProxAL.optimize!(nlp)
+    @test isapprox(runinfo.objvalue[end], OPTIMAL_OBJVALUE, rtol = 1e-2)
+    @test isapprox(runinfo.x.Pg[:,1,:][:], OPTIMAL_PG, rtol = 1e-1)
+    @test runinfo.maxviol_c[end] <= algparams.tol
+    @test runinfo.maxviol_t[end] <= algparams.tol
+    @test runinfo.maxviol_c_actual[end] <= algparams.tol
+    @test runinfo.maxviol_t_actual[end] <= algparams.tol
+    @test runinfo.maxviol_d[end] <= algparams.tol
+    @test runinfo.iter <= algparams.iterlim
+end
+
+use_MPI && MPI.Finalize()

test/runtests.jl

@@ -15,6 +15,10 @@ testdir = @__DIR__
     @testset "Sequential tests" begin
         include("convergence.jl")
     end
+    # Testing N-2
+    @testset "Sequential tests" begin
+        include("n-2.jl")
+    end
 
     # Testing using 2 processes