Add prototype functions for cascade algorithm

src/Data.jl

@@ -215,6 +215,25 @@ end
 
 #*------------------------------------------------------------------------------
 
+#=
+Returns an array of tuples containing branch indices of branches with a minimum loading of min_loading and their actual loading (considering the power flow type pf_type). Choosing min_loading = 1. yields overloaded branches.
+=#
+function get_branches(
+        network_data::Dict{String,<:Any}; 
+        min_loading = 0.0, # minimum loading of branches
+        pf_type = "MVA-loading" # or "MW-loading", "Mvar-loading"
+    )
+
+    branches = [
+        (i, br[pf_type]) for (i, br) in network_data["branch"]
+        if br["br_status"] == 1 && br[pf_type] >= min_loading
+    ]
+
+    return branches
+end
+
+#*------------------------------------------------------------------------------
+
 #=
 Deactivates branches that have been overloaded due to a power flow redistribution (flow/capacity > 1.). The branches are identified by their string indices in branch_ids. The status of these branches is set to 0 in network_data. Parallel branches are not automatically deactivated.
 =#

src/ITCPG.jl

@@ -25,7 +25,7 @@ export disable_branches!, destroy_tl!
 
 include("PowerFlow.jl")
 export update_pf_data!, calc_init_op, calc_ac_pf!, calc_branchloads!
-export calc_total_impact!
+export calc_total_impact!, restore_p_balance!
 
 include("PlotUtils.jl")
 

src/PowerFlow.jl

@@ -53,21 +53,29 @@ end
 Calculates the AC power flow solution for the given NDD and updates it in the NDD.
 CAUTION: Overwrites any previous AC solution contained in the NDD!
 =#
-function calc_ac_pf!(network_data::Dict{String,<:Any}, method=:JuMP)
+function calc_ac_pf!(
+        network_data::Dict{String,<:Any}, 
+        method = :JuMP; 
+        print_level = 0
+    )
+
     if method == :JuMP # calculate AC-PF solution using a JuMP model
         pf_result = run_ac_pf(
             network_data, 
-            optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
+            optimizer_with_attributes(
+                Ipopt.Optimizer, 
+                "print_level" => print_level
+            )
         )
     elseif method == :NLsolve # calculate AC-PF solution using NLsolve.jl
         pf_result = Dict{String,Any}()
-        pf_result["solution"] = compute_ac_pf(network_data)
+        pf_result["solution"] = compute_ac_pf(network_data, show_trace=true)
     else
         throw(ArgumentError("Unknown method $method."))
     end
 
     update_pf_data!(network_data, pf_result, model=:ac) # update PF in NDD
-    return network_data
+    return nothing
 end
 
 #=
@@ -133,3 +141,234 @@ end
 
 #*------------------------------------------------------------------------------
 
+#=
+Calculates the total impact of a hurricane given the time series of damaged transmission lines. Each time a transmission line is destroyed, the AC power flow equations are solved and failure cascades may unfold. The output of the solver is saved in a .txt file.
+=#
+function calc_total_impact!(
+        network_data::Dict{String,<:Any},
+        tl_damage::Array{Tuple{String,Int64},1}, # time series of tl damage
+        savepath::String;
+        method = :JuMP, # solver to use
+        print_level = 5 # Ipopt output details
+    )
+
+    output = @capture_out _calc_total_impact!(
+        network_data, tl_damage, method=method, print_level=print_level
+    )
+    ### Write terminal output of solver into file
+    open(savepath, "w") do io
+        write(io, output)
+    end
+
+    return nothing
+end
+
+function _calc_total_impact!(
+        network_data::Dict{String,<:Any},
+        tl_damage::Array{Tuple{String,Int64},1}; # time series of tl damage
+        method = :JuMP, # solver to use
+        print_level = 0 # Ipopt output
+    )
+
+    ### Get unique frames in which transmission lines were damaged
+    unique_frames = sort(unique([tl_damage[i][2] for i in 1:length(tl_damage)]))
+
+    ### Go through frames in which transmission lines were destroyed
+    for f in unique_frames
+        destroyed_tl = [dmg[1] for dmg in tl_damage if dmg[2] == f]
+        destroy_tl!(network_data, destroyed_tl)
+        ### Resolve AC-PF and see whether a cascade unfolds
+        _calc_cascade!(network_data, method, print_level)
+    end
+    
+    return nothing
+end
+
+#=
+Calculates whether a cascading failure unfolds in a damaged network. If any network contained in network_data misses a slack bus, the largest generator is selected in each connected component. Additionally, the active power balance is restored for each connected component in every step, if it is violated.
+=#
+function _calc_cascade!(
+        network_data::Dict{String,<:Any},
+        method = :JuMP,
+        print_level = 0
+    )
+
+    ### Update network topology (remove dangling buses etc.)
+    simplify_network!(network_data)
+
+    ### Define new slacks, if necessary
+    correct_reference_buses!(network_data)
+    
+    ### Restore active power balance, if violated
+    restore_p_balance!(network_data)
+    
+    # connected_components = calc_connected_components(network_data)
+    # for cc in connected_components
+    #     restore_active_power_balance!(network_data, cc)
+    # end
+
+    ### Calculate new AC-PF solution
+    set_ac_pf_start_values!(network_data) # last state as initial condition
+    calc_ac_pf!(network_data, method, print_level=print_level)
+    
+    ### Check for overloaded branches
+    overloaded_br = get_branches(network_data, min_loading=1.)
+
+    if isempty(overloaded_br) == false
+        ### Print overloaded branches as output
+        println("-------------------------------------------------------------")
+        println("Overloaded branches:", overloaded_br)
+        println("-------------------------------------------------------------")
+        
+        ### Deactivate overloaded branches
+        disable_branches!(network_data, [br[1] for br in overloaded_br])
+        ### Repeat cascade algorithm
+        _calc_cascade!(network_data, method, print_level)
+    end
+
+    return nothing
+end
+
+function restore_p_balance!(network_data::Dict{String,<:Any})
+    ### Lookup which generators and loads are connected to which buses
+    gen_lookup = bus_gen_lookup(network_data["gen"], network_data["bus"])
+    load_lookup = bus_load_lookup(network_data["load"], network_data["bus"])
+
+    ### Calculate connected components
+    connected_components = calc_connected_components(network_data)
+    for cc in connected_components
+        _restore_p_balance!(network_data, gen_lookup, load_lookup, cc)
+    end
+
+    return nothing
+end
+
+function _restore_p_balance!(
+        network_data::Dict{String,<:Any},
+        gen_lookup::Dict{Int64,Array{Any,1}},
+        load_lookup::Dict{Int64,Array{Any,1}},
+        cc::Set{Int64} # indices of buses in connected component
+    )
+
+    ### Find generators in connected component 
+    c_gens = [
+        gens for (i, gens) in gen_lookup 
+        if i ∈ cc && !isempty(gens)
+    ]
+    ### Find loads in connected component
+    c_loads = [
+        loads for (i, loads) in load_lookup 
+        if i ∈ cc && !isempty(loads)
+    ]
+    
+    ### Calculate active power mismatches
+    c_Δp, c_Δp_lim = _calc_p_mismatch(
+        network_data, c_gens, c_loads, cc
+    )
+    
+    println(c_Δp)
+    println(c_Δp_lim)
+    println(length(c_gens))
+    # if c_Δp < 0 # overproduction
+    #     if isapprox(c_Δp_lim, 0, atol=1e-10) # "sufficient" balance
+    #         ### Set all dispatches to their minimum
+    #         for gen in c_gens
+    #             str_id = string(gen["index"])
+    #             network_data["gen"][str_id]["pg"] = gen["pmin"]
+    #         end
+    #     elseif c_Δp_lim > 0 # adjustment possible
+    #         ### Reduce all dispatches uniformly
+            
+    # end
+
+
+    # get total active power load of loads connected to buses in cc
+    # get total active power capacity of generators connected to buses in cc
+    # check whether generators can match the total active power load
+    # -> if yes: Increase generation uniformly (check if a generator hits capacity)
+    # -> if not: Reduce active power loads uniformly until total load is equal to total active power capacity of generators
+    # update generators (and loads) in network data dictionary
+
+end
+
+function _calc_p_mismatch(
+        network_data::Dict{String,<:Any},
+        c_gens::Array{Array{Any,1},1}, 
+        c_loads::Array{Array{Any,1},1},
+        cc::Set{Int64}
+    )
+
+    ### Calculate active power demand in connected component
+    c_pd = 0.
+    for load_list in c_loads
+        for load in load_list
+            if load["status"] != 0 # load is active
+                c_pd += load["pd"]
+            end
+        end
+    end
+    
+    ### Calculate total active power loss in branches in connected component
+    plosses = [
+        br["pt"] + br["pf"] for br in values(network_data["branch"]) 
+        if br["f_bus"] ∈ cc && br["t_bus"] in cc && br["br_status"] == 1
+    ]
+    c_ploss = sum(plosses)
+
+    ### Calculate active power generation in connected component
+    c_pg, c_pmin, c_pmax = 0., 0., 0.
+    for gen_list in c_gens
+        for gen in gen_list
+            if gen["gen_status"] != 0 # generator is active
+                c_pg += gen["pg"]
+                c_pmin += gen["pmin"]
+                c_pmax += gen["pmax"]
+            end
+        end
+    end
+
+    ### Calculate relevant mismatches
+    c_Δp = c_pd + c_ploss - c_pg
+    if isapprox(c_Δp, 0, atol=1e-10) # "sufficient" power balance
+        c_Δp = 0.0
+        c_Δp_lim = NaN
+    elseif c_Δp < 0 # overproduction
+        c_Δp_lim = c_pd + c_ploss - c_pmin
+    elseif c_Δp > 0 # underproduction
+        c_Δp_lim = c_pd + c_ploss - c_pmax
+    end
+    
+    # println("max. capacity: ", c_pmax)
+    # println("generation: ", c_pg)
+    # println("demand: ", c_pd)
+    # println("generation - demand: ", c_pg - c_pd)
+    # println("losses: ", c_ploss)
+    # println("mismatch: ", c_pd + c_ploss - c_pg)
+    # println("mismatch lim: ", c_pd + c_ploss - c_pmax)
+    # println(c_Δp)
+
+    return c_Δp, c_Δp_lim
+end
+
+#=
+
+Ideas for treating grid partitioning during cascading failures
+---
+
+-> Start by deactivating transmission line(s) that are the first destroyed by the hurricane
+
+-> Run simplify_network! to update the grid topology 
+- save any deactivated load buses, since they contribute to the overall power outage!
+
+-> Calculate connected components via calc_connected_components and check whether there are more than one connected component
+-> If there is only one connected component: Jump to power flow calculation
+-> If there are more than one connected component: Use correct_reference_buses! to define new slack buses for all components (uses "pmax" to find generator with largest capacity). Check whether active power balance is violated in any component and if so correct it! Jump to power flow calculation
+
+-> Calculate AC power flow using run_ac_pf (JuMP model), since it tries to solve all connected components in parallel (NLsolve does not do that).
+
+Question: Still use previous voltage magnitudes and angles as starting values after a graph partition happened??
+
+MONITOR VOLTAGES DURING A CASCADE!? (reactive power imbalance may lead to voltage instability)
+CALCULATE TOTAL FLOW CHANGES IN BRANCHES DURING CASCADES
+
+=#
\ No newline at end of file