Rework the cascading failure algorithm

src/Analysis.jl

+#* Functions for analysing the state of a power grid given its data dictionary
+#*------------------------------------------------------------------------------
+
+
+
+
+
+
+
+#*------------------------------------------------------------------------------
+#* Functions that check AC-PF constraints
+
+#=
+Checks and logs the constraints considered in the JuMP model that PowerModels.jl builts for the AC power flow problem. The constraints are:
+1. Kirchoff's current law: AC power balance at each bus
+2. Slack bus constraints: Voltage angle and magnitude fixed
+3. PV-bus constraints: Voltage magnitudes fixed at setpoints
+4. Generator constraints (generators are connected to PV-buses): Active power dispatches fixed
+=#
+function check_ac_constraints(
+        new_ndd::Dict{String,<:Any}, # network data with new solution
+        old_ndd::Dict{String,<:Any} # network data with old solution
+    )
+    
+    ### Check nodal power balances (Kirchhoff's current law)
+    check_nodal_ac_power_balance(new_ndd)
+
+    ### Check constraints 2. - 4. for each connected component, respectively
+    connected_components = calc_connected_components(new_ndd)
+    for (cc_i, cc) in enumerate(connected_components)
+        check_ac_constraints(new_ndd, old_ndd, cc, cc_i)
+    end
+    
+    return nothing
+end
+
+#=
+Calculates and logs the maximum nodal active and reactive power mismatch. Both the magnitude of the mismatch and the corresponding bus are logged.
+=#
+function check_nodal_ac_power_balance(network_data::Dict{String,<:Any})
+    ### Calculate nodal AC power mismatches
+    nodal_balance = calc_power_balance(network_data)
+
+    ### Determine maximum active and reactive power mismatch
+    max_p_delta = 0.
+    max_p_delta_bus = ""
+    max_q_delta = 0.
+    max_q_delta_bus = ""
+    for (key, bus) in nodal_balance["bus"]
+        if abs(bus["p_delta"]) > max_p_delta
+            max_p_delta = abs(bus["p_delta"])
+            max_p_delta_bus = key
+        end
+        if abs(bus["q_delta"]) > max_q_delta
+            max_q_delta = abs(bus["q_delta"])
+            max_q_delta_bus = key
+        end
+    end
+
+    ### Log results
+    info(LOGGER, "Maximum nodal active power mismatch: $max_p_delta at bus $max_p_delta_bus")
+    info(LOGGER, "Maximum nodal reactive power mismatch: $max_q_delta at bus $max_q_delta_bus")
+
+    return nothing
+end
+
+function check_area_nodal_power_balance(
+        network_data::Dict{String,<:Any}, 
+        area = 7
+    )
+    
+    active_area_buses = String[
+        "$(b["index"])" for b in values(network_data["bus"]) if 
+        b["area"] == area && b["bus_type"] != 4
+    ]
+    nodal_balance = calc_power_balance(network_data)["bus"]
+
+    area_nodal_balance = Dict(
+        i => values for (i, values) in nodal_balance if i ∈ active_area_buses
+    )
+
+    max_p_delta = 0.
+    max_p_delta_bus = ""
+    max_q_delta = 0.
+    max_q_delta_bus = ""
+
+    for (key, bus) in area_nodal_balance
+        if abs(bus["p_delta"]) > max_p_delta
+            max_p_delta = abs(bus["p_delta"])
+            max_p_delta_bus = key
+        end
+        if abs(bus["q_delta"]) > max_q_delta
+            max_q_delta = abs(bus["q_delta"])
+            max_q_delta_bus = key
+        end
+    end
+
+    info(LOGGER, "Area $area: Maximum nodal active power mismatch: $max_p_delta at bus $max_p_delta_bus")
+    info(LOGGER, "Area $area: Maximum nodal reactive power mismatch: $max_q_delta at bus $max_q_delta_bus")
+
+    return nothing
+end
+
+function check_ac_constraints(
+        new_ndd::Dict{String,<:Any}, # network data with new solution
+        old_ndd::Dict{String,<:Any}, # network data with old solution
+        cc::Set{Int64}, # set of bus indices in connected component
+        cc_i::Int64 # index/number of connected component
+    )
+
+    ### Check slack constraints
+    cc_slack_new = [
+        b for b in values(new_ndd["bus"]) if b["bus_type"]==3 && b["index"]∈cc
+    ] # slack bus in new network data
+    cc_slack_old = [
+        b for b in values(old_ndd["bus"]) if b["bus_type"]==3 && b["index"]∈cc
+    ] # slack bus in old network data
+    if length(cc_slack_new) != 1
+        warn(LOGGER, "Grid #$cc_i has $(length(cc_slack_new)) slack buses!")
+    elseif cc_slack_new[1]["index"] != cc_slack_old[1]["index"]
+        warn(LOGGER, 
+            "Slack bus $(cc_slack_new[1]["index"]) in the new state of $cc_i does not agree with slack bus $(cc_slack_old[1]["index"]) in the old state!"
+        )
+    end
+    cc_slack_gen = [
+        gen["index"] for gen in values(new_ndd["gen"]) if 
+        gen["gen_bus"] == cc_slack_new[1]["index"]
+    ] # index of generator connected to the slack bus
+    info(LOGGER, 
+        "Grid #$cc_i: Slack bus $(cc_slack_new[1]["index"]) (with generator $cc_slack_gen) voltage angle constraint: va_new = $(cc_slack_new[1]["va"]), va_old = $(cc_slack_old[1]["va"])"
+    )
+    info(LOGGER, 
+        "Grid #$cc_i: Slack bus $(cc_slack_new[1]["index"]) (with generator $cc_slack_gen) voltage magnitude constraint: vm_new = $(cc_slack_new[1]["vm"]), vm_old = $(cc_slack_old[1]["vm"])"
+    )
+    info(LOGGER, 
+        "Grid #$cc_i: Slack bus $(cc_slack_new[1]["index"]) (with generator $cc_slack_gen) active power dispatch: pg_new = $(new_ndd["gen"]["$(cc_slack_gen[1])"]["pg"]), pg_old = $(old_ndd["gen"]["$(cc_slack_gen[1])"]["pg"])"
+    )
+
+    ### Check PV-bus constraints
+    cc_PV_new = [
+        b for b in values(new_ndd["bus"]) if b["bus_type"] == 2 && b["index"]∈cc
+    ]
+    max_vm_diff = 0.0 # maximum voltage magnitude mismatch
+    max_vm_diff_bus = 0 # index of corresponding PV-bus
+    min_vm_diff = 1e10 # minimum voltage magnitude mismatch
+    min_vm_diff_bus = 0 # index of corresponding PV-bus
+    # vm_mismatches = zeros(length(cc_PV_new))
+    for PV_bus in cc_PV_new
+        vm_diff = PV_bus["vm"] - old_ndd["bus"]["$(PV_bus["index"])"]["vm"]
+        if abs(vm_diff) >= abs(max_vm_diff)
+            max_vm_diff = vm_diff
+            max_vm_diff_bus = PV_bus["index"]
+        end
+        if abs(vm_diff) < abs(min_vm_diff)
+            min_vm_diff = vm_diff
+            min_vm_diff_bus = PV_bus["index"]
+        end
+    end
+    info(LOGGER, 
+        "Grid #$cc_i: Maximum voltage magnitude mismatch: vm_new - vm_old = $max_vm_diff at PV-bus $max_vm_diff_bus"
+    )
+    info(LOGGER, 
+        "Grid #$cc_i: Minimum voltage magnitude mismatch: vm_new - vm_old = $min_vm_diff at PV-bus $min_vm_diff_bus"
+    )
+
+    ### Check generator active power constraints
+    active_cc_gens = [
+        gen for gen in values(new_ndd["gen"]) if gen["gen_bus"] ∈ cc &&
+        gen["gen_status"] != 0 && gen["index"] ∉ cc_slack_gen
+    ] # active generators that are connected to buses in cc
+    max_pg_diff = 0.0 # maximum mismatch in active power dispatch
+    max_pg_diff_gen = 0 # index of corresponding generator
+    min_pg_diff = 1e10 # minimum mismatch in active power dispatch
+    min_pg_diff_gen = 0 # index of corresponding generator
+    for gen in active_cc_gens
+        pg_diff = gen["pg"] - old_ndd["gen"]["$(gen["index"])"]["pg"]
+        if abs(pg_diff) >= abs(max_pg_diff)
+            max_pg_diff = pg_diff
+            max_pg_diff_gen = gen["index"]
+        end
+        if abs(pg_diff) < abs(min_pg_diff)
+            min_pg_diff = pg_diff
+            min_pg_diff_gen = gen["index"]
+        end
+    end
+    info(LOGGER, 
+        "Grid #$cc_i: Maximum mismatch in active power dispatch: pg_new - pg_old = $max_pg_diff at generator $max_pg_diff_gen"
+    )
+    info(LOGGER, 
+        "Grid #$cc_i: Minimum mismatch in active power dispatch: pg_new - pg_old = $min_pg_diff at generator $min_pg_diff_gen"
+    )
+
+    ### Check maximum voltage angle difference along branches
+    cc_branches = [
+        (br["f_bus"], br["t_bus"]) for br in values(new_ndd["branch"]) if 
+        br["br_status"] == 1 && br["f_bus"] ∈ cc && br["t_bus"] ∈ cc
+    ]
+    max_va_diff = 0.0 # maximum voltage angle difference
+    max_va_diff_br = (0, 0) # from- and to-bus of corresponding branch
+    min_va_diff = 1e10 # minimum voltage angle difference
+    min_va_diff_br = (0, 0) # from- and to-bus of corresponding branch
+    for br in cc_branches
+        va_diff = new_ndd["bus"]["$(br[2])"]["va"] - new_ndd["bus"]["$(br[1])"]["va"]
+        if abs(va_diff) >= abs(max_va_diff)
+            max_va_diff = va_diff
+            max_va_diff_br = br
+        end
+        if abs(va_diff) < abs(min_va_diff)
+            min_va_diff = va_diff
+            min_va_diff_br = br
+        end
+    end
+    info(LOGGER, 
+        "Grid #$cc_i: Maximum branch voltage angle difference: va_to - va_fr = $(rad2deg(max_va_diff))° in branch $max_va_diff_br"
+    )
+    info(LOGGER, 
+        "Grid #$cc_i: Minimum branch voltage angle difference: va_to - va_fr = $(rad2deg(min_va_diff))° in branch $min_va_diff_br"
+    )
+    
+    return nothing
+end
+
+function check_area_power_balance(network_data::Dict{String,<:Any}, area=7)
+    active_area_buses = Int64[
+        b["index"] for b in values(network_data["bus"]) if 
+        b["area"] == area && b["bus_type"] != 4
+    ]
+    active_area_loads = Int64[
+        l["index"] for l in values(network_data["load"]) if
+        l["load_bus"] ∈ active_area_buses && l["status"] == 1
+    ]
+    active_area_gens = Int64[
+        g["index"] for g in values(network_data["gen"]) if
+        g["gen_bus"] ∈ active_area_buses && g["gen_status"] == 1
+    ]
+
+    ### Active power balance
+    area_pd = sum(
+        Float64[network_data["load"]["$l"]["pd"] for l in active_area_loads]
+    )
+    area_pg = sum(
+        Float64[network_data["gen"]["$g"]["pg"] for g in active_area_gens]
+    )
+    area_ploss = sum(
+        Float64[br["pt"] + br["pf"] for br in values(network_data["branch"]) if 
+        br["f_bus"] ∈ active_area_buses && br["t_bus"] in active_area_buses
+        && br["br_status"] == 1]
+    )
+
+    area_Δp = area_pg - area_pd - area_ploss
+    info(LOGGER, "Active power balance in area $area: $area_Δp")
+
+    ### Reactive power balance
+    area_qd = sum(
+        Float64[network_data["load"]["$l"]["qd"] for l in active_area_loads]
+    )
+    area_qg = sum(
+        Float64[network_data["gen"]["$g"]["qg"] for g in active_area_gens]
+    )
+    area_qloss = sum(
+        Float64[br["qt"] + br["qf"] for br in values(network_data["branch"]) if 
+        br["f_bus"] ∈ active_area_buses && br["t_bus"] in active_area_buses
+        && br["br_status"] == 1]
+    )
+    area_qshunt = 0.
+    for shunt in values(network_data["shunt"])
+        if shunt["status"] == 1 && shunt["shunt_bus"] ∈ active_area_buses
+            vm = network_data["bus"]["$(shunt["shunt_bus"])"]["vm"]
+            area_qshunt += shunt["bs"] * vm^2
+        end
+    end
+
+    area_Δq = area_qg - area_qd - area_qloss + area_qshunt
+    info(LOGGER, "Reactive power balance in area $area: $area_Δq")
+
+    return nothing
+end
\ No newline at end of file

src/Data.jl

@@ -9,24 +9,37 @@ busfile - .csv file that contains geographic locations of all buses
 branchfile - optional .csv file that contains transmission line lengths
 nddfile - .jld2 file in which the final network data dictionary will be saved
 =#
+
+"""
+    build_network_data(file::String; kwargs...)
+
+Reads a .RAW or .m `file` containing the data of a power grid and calculates the power flow according to the chosen model (see below). Returns a network data dictionary conform with the [PowerModels Network Data Format](https://lanl-ansi.github.io/PowerModels.jl/stable/network-data/) that also contains the calculated power flow solution.
+"""
 function build_network_data(
-        model = :ac; 
-        pgfile::String,  busfile::String, nddfile::String, branchfile = ""
+        file::String;
+        locfile = "", lenfile = "", outfile = "",
+        pf_model = :ac,
+        kwargs...
     )
 
-    ### Read grid data and calculate initial operation point
-    network_data = calc_init_op(pgfile, model=model)
+    ### Read grid data and calculate operation point
+    network_data = calc_pf(file, pf_model=pf_model; kwargs...)
     
-    ### Add additional data to dictionary
-    add_locs!(network_data, busfile) # add geographic bus locations
-    if isempty(branchfile) == false # check whether to add line lengths
-        add_tl_lengths!(network_data, branchfile) 
+    ### Add additional data to network_data, if given
+    if isempty(locfile) == false
+        add_locs!(network_data, locfile) # add geographic bus locations
+    end
+    if isempty(lenfile) == false
+        add_tl_lengths!(network_data, lenfile)
     end
+
     add_tl_voltages!(network_data) # add transmission line voltages
-    
-    save(nddfile, "network_data", network_data) # save network data
 
-    return nothing
+    if isempty(outfile) == false
+        save(outfile, "network_data", network_data) # save network data
+    end
+
+    return network_data
 end
 
 #*------------------------------------------------------------------------------
@@ -155,7 +168,6 @@ function get_bustypes(network_data::Dict{String,<:Any})
     filter!(i -> i ∉ slack, load)
     filter!(i -> i ∉ slack, gen)
     filter!(i -> i ∉ slack, load_and_gen)
-
     
     ### Empty buses with neither load nor generation
     empty = unique(
@@ -242,25 +254,97 @@ function get_MW_loads(network_data::Dict{String,<:Any})
     return load_dict
 end
 
+function get_total_load(network_data::Dict{String,<:Any})
+    total_MW_load = 0.
+    total_Mvar_load = 0.
+    total_MVA_load = 0.
+
+    for load in [l for l in values(network_data["load"]) if l["status"] == 1]
+        total_MW_load += load["pd"]
+        total_Mvar_load += load["qd"]
+        total_MVA_load += sqrt(load["pd"]^2 + load["qd"]^2)
+    end
+
+    return total_MW_load, total_Mvar_load, total_MVA_load
+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.
-=#
+# Returns an array of tuples that each contain information about specific branches of interest. The tuples contain the respective branch index, the loading of the branch and a boolean value set to `1`, if the branch is a transformer, or `0`, if the branch is a transmission line. 
+"""
+    get_branches(network_data::Dict{String,<:Any}; kwargs...)
+
+Returns the indices of branches in `network_data` that are loaded (flow/capacity) by at least `min_loading` (see **keyword arguments**).
+
+# Keyword Arguments
+
+`pf_type` can be set to `"MVA-loading"` (default), `"MW-loading"` or `"Mvar-loading"` and specifies whether the loading due to apparent (MVA), active (MW) or reactive power (Mvar) should be considered. If `network_data["pf_model"] => :dc` meaning that it was last updated using a DC power flow solution, `pf_type` is automatically set to `"MW-loading"`.
+
+`min_loading` is a `Float64` (default `0.0`) that specifies the minimum loading of branches to be returned (according to `pf_type`). Only branches with a loading equal or above `min_loading` will be returned. Hence, `min_loading = 1.0` can be used to obtain overloaded branches.
+
+`br_status` can be set to `1` (default) or `0` and specifies whether active (`1`) or inactive branches (`0`) should be returned.
+"""
 function get_branches(
         network_data::Dict{String,<:Any}; 
+        pf_type = "MVA-loading", # or "MW-loading", "Mvar-loading"
         min_loading = 0.0, # minimum loading of branches
-        pf_type = "MVA-loading" # or "MW-loading", "Mvar-loading"
+        br_status = 1 # active or inactive branches
     )
 
-    branches = [
-        (i, br[pf_type]) for (i, br) in network_data["branch"]
-        if br["br_status"] == 1 && br[pf_type] >= min_loading
-    ]
+    ### Check whether network_data contains a DC-PF solution
+    if network_data["pf_model"] == :dc
+        pf_type = "MW-loading" # use branch loadings due to active power
+    end
+
+    # branches = [
+    #     (i, br[pf_type], br["transformer"]) for (i, br) in network_data["branch"] if br["br_status"] == 1 && br[pf_type] >= min_loading
+    # ]
+    if br_status == 1
+        branches = [
+            br["index"] for br in values(network_data["branch"]) if 
+            br["br_status"] == br_status && br[pf_type] >= min_loading
+        ]
+    elseif br_status == 0
+        branches = [
+            br["index"] for br in values(network_data["branch"]) if 
+            br["br_status"] == br_status
+        ]
+    else
+        throw(ArgumentError("Unknown branch status $(br_status)!"))
+    end
 
     return branches
 end
 
+"""
+    get_cc_data(network_data::Dict{String,<:Any}, cc::Set{Int64})
+"""
+function get_cc_data(
+        network_data::Dict{String,<:Any}, cc::Set{Int64}; name = "cc"
+    )
+
+    cc_network_data = Dict{String,Any}(
+        "baseMVA" => network_data["baseMVA"],
+        "per_unit" => network_data["per_unit"],
+        "source_type" => network_data["source_type"],
+        "name" => name,
+        "source_version" => network_data["source_version"],
+        "bus" => Dict{String,Any}(),
+        "gen" => Dict{String,Any}(),
+        "load" => Dict{String,Any}(),
+        "shunt" => Dict{String,Any}(),
+        "branch" => Dict{String,Any}(),
+        "dcline" => Dict{String,Any}(),
+        "storage" => Dict{String,Any}(),
+        "switch" => Dict{String,Any}()
+    )
+    
+    for b in cc
+        cc_network_data["bus"]["$b"] = network_data["bus"]["$b"]
+        
+    end
+end
+
 #*------------------------------------------------------------------------------
 
 #=
@@ -268,7 +352,7 @@ Deactivates branches that have been overloaded due to a power flow redistributio
 =#
 function disable_branches!(
         network_data::Dict{String,<:Any}, 
-        branch_ids::Array
+        branch_ids::Array{Int64,1}
     )
     ### Go through branches to deactivate
     for id in branch_ids
@@ -280,12 +364,35 @@ function disable_branches!(
     return nothing
 end
 
+"""
+    deactivate_overloads!(network_data::Dict{String,<:Any})
+"""
+function deactivate_overloads!(network_data::Dict{String,<:Any})
+    ### Get overloaded active branches and deactivate them
+    overloaded_br = get_branches(network_data, min_loading=1.)
+    if isempty(overloaded_br) == false
+        info(LOGGER, 
+            "Deactivating the following $(length(overloaded_br)) overloaded branches: $overloaded_br"
+        )
+    else
+        info(LOGGER, "No branches are overloaded.")
+    end
+    for br in overloaded_br
+        network_data["branch"]["$br"]["br_status"] = 0
+        network_data["branch"]["$br"]["failure_reason"] = "overload"
+    end
+
+    simplify_network!(network_data)
+
+    return overloaded_br
+end
+
 #=
 Function that deactivates (overhead) transmission lines identified by their indices (string) in tl_ids that have been destroyed by a hurricane. For any given line index all parallel lines also become deactivated. The status of these lines is set to 0 in network_data. Afterwards PowerModels function simplify_network! is run in order to remove dangling buses etc. that might occur due to the failures.
 =#
 function destroy_tl!(
         network_data::Dict{String,<:Any}, 
-        tl_ids::Array # string ids of destroyed transmission lines
+        tl_ids::Array # string-ID's of destroyed transmission lines
     )
 
     ### Go through transmission lines to destroy
@@ -300,6 +407,7 @@ function destroy_tl!(
         )
         for idx in ids
             network_data["branch"][idx]["br_status"] = 0
+            network_data["branch"][idx]["failure_reason"] = "wind"
         end
     end
 
@@ -330,5 +438,127 @@ function _deactivate_cc!(
     ### Deactivate rest of connected component (branches, dangling buses etc.)
     simplify_network!(network_data)
     
+    return nothing
+end
+
+#*------------------------------------------------------------------------------
+
+"""
+    check_voltage_bounds(network_data::Dict{String,Any}; kwargs...)
+
+Checks whether active buses in `network_data` have voltage magnitudes outside of their tolerance interval (usually ranging from 0.9 p.u. to 1.1 p.u.) and returns a dictionary with the bus indices as keys and the respective active loads connected to the bus as values (PQ-generators are excluded). ITCPG's logger prints a `warn` message everytime a voltage violation is detected.
+
+# Keyword Arguments
+
+`limit` can be set to `"lower"`, `"upper"` or `"both"` (default) and specifies which voltage bound should be considered. If e.g. `limit = "lower"` is chosen, only buses with under-voltages are returned.
+
+`filter_active_loads` can be set to `true` or `false` (default) and specifies whether only buses with active loads connected to them should be returned (useful for under-voltage load shedding). 
+
+See also: [`run_outer_ac_pf!`](@ref), [`run_uvls!`](@ref)
+"""
+function check_voltage_bounds(
+        network_data::Dict{String,<:Any}; 
+        limit = "both", # considered voltage bounds
+        filter_active_loads = false # whether to filter buses with active loads
+    )
+
+    vv_buses = Dict{String,Array{String,1}}()
+    ### Go through buses and check their voltage magnitudes
+    for (key, bus) in network_data["bus"]
+        vm = bus["vm"]
+        if limit == "both" || limit == "lower"
+            if vm < bus["vmin"] && bus["bus_type"] != 4
+                warn(LOGGER, 
+                    "Bus $key (type $(network_data["bus"][key]["bus_type"])) has a voltage magnitude of |V|=$vm that is below its lower bound of $(bus["vmin"])" 
+                )
+                ### Get string indices of active loads connected to the bus
+                connected_loads = [
+                    l_id for (l_id, l) in network_data["load"] if 
+                    l["load_bus"] == parse(Int64, key) && l["status"] == 1 && 
+                    haskey(l, "is_gen") == false # exclude PQ-generators
+                ]
+                vv_buses[key] = connected_loads
+            end
+        end
+        if limit == "both" || limit == "upper"
+            if vm > bus["vmax"] && bus["bus_type"] != 4
+                warn(LOGGER, 
+                    "Bus $key (type $(network_data["bus"][key]["bus_type"])) has a voltage magnitude of |V|=$vm that is above its upper bound of $(bus["vmax"])" 
+                )
+                ### Get string indices of active loads connected to the bus
+                connected_loads = [
+                    l_id for (l_id, l) in network_data["load"] if 
+                    l["load_bus"] == parse(Int64, key) && l["status"] == 1 && 
+                    haskey(l, "is_gen") == false # exclude PQ-generators
+                ]
+                vv_buses[key] = connected_loads
+            end
+        end
+    end
+
+    ### Check whether to exclude buses without active loads
+    if filter_active_loads == true
+        filter!(x -> isempty(last(x)) == false, vv_buses)
+    end
+
+    return vv_buses
+end
+
+function check_gen_bounds(network_data::Dict{String,<:Any})
+    for (i, gen) in network_data["gen"]
+        if gen["gen_status"] == 1
+            if gen["pg"] > gen["pmax"] || gen["pg"] < gen["pmin"]
+                warn(LOGGER, 
+                    "Active power dispatch of generator $i lies outside of interval [$(gen["pmin"]), $(gen["pmax"])]: $(gen["pg"])"
+                )
+            end
+            if gen["qg"] > gen["qmax"] || gen["qg"] < gen["qmin"]
+                warn(LOGGER,
+                    "Reactive power dispatch of generator $i lies outside of interval [$(gen["qmin"]), $(gen["qmax"])]: $(gen["qg"])"
+                )
+            end
+        end
+    end
+
+    return nothing
+end
+
+function check_gen_bounds(
+        network_data::Dict{String,<:Any},
+        cc::Set{Int64}, # set of bus indices in connected component
+        cc_i::Int64; # index/number of connected components
+        power = "active" # "active": active power, "both": reactive power too
+    )
+
+    cc_slack = [
+        b["index"] for b in values(network_data["bus"]) if 
+        b["bus_type"] == 3 && b["index"] ∈ cc
+    ]
+    
+    if length(cc_slack) != 1
+        warn(LOGGER, "(Sub-)Grid #$cc_i has $(length(cc_slack)) slack buses!")
+    end
+
+    for (i, gen) in network_data["gen"]
+        if gen["gen_bus"] ∈ cc && gen["gen_status"] == 1
+            if gen["gen_bus"] ∈ cc_slack
+                info(LOGGER, 
+                    "Active power dispatch of slack bus $i in #$cc_i with interval [$(gen["pmin"]), $(gen["pmax"])]: $(gen["pg"])"
+                )
+            else
+                if gen["pg"] > gen["pmax"] || gen["pg"] < gen["pmin"]
+                    warn(LOGGER, "Active power dispatch of generator $i in #$cc_i lies outside of interval [$(gen["pmin"]), $(gen["pmax"])]: $(gen["pg"])")
+                end
+                if power == "both"
+                    if gen["qg"] > gen["qmax"] || gen["qg"] < gen["qmin"]
+                        warn(LOGGER,
+                            "Reactive power dispatch of generator $i lies outside of interval [$(gen["qmin"]), $(gen["qmax"])]: $(gen["qg"])"
+                        )
+                    end
+                end
+            end
+        end
+    end
+
     return nothing
 end
\ No newline at end of file

src/ITCPG.jl

@@ -45,20 +45,28 @@ export config_logger
 include("Data.jl")
 export build_network_data
 export add_locs!, add_tl_lengths!, add_tl_voltages!
-export get_bustypes, get_underground_tl, get_MW_loads, get_branches
-export disable_branches!, destroy_tl!
+export get_bustypes, get_underground_tl, get_MW_loads, get_total_load
+export get_branches, deactivate_overloads!, destroy_tl!
+export check_voltage_bounds, check_gen_bounds
+
+include("Analysis.jl")
+export check_ac_constraints, check_nodal_power_balance 
+export check_area_nodal_power_balance, check_area_power_balance
 
 include("PowerFlow.jl")
-export update_pf_data!, calc_init_op, calc_ac_pf!, calc_branchloads!
-export calc_total_impact!, calc_cascade!, restore_p_balance!, calc_Δp
+export calc_ac_pf!, calc_dc_pf!, calc_secondary_dmg!
+export update_pf_data!, calc_pf, calc_branchloads!
+export calc_new_state!, run_outer_ac_pf!, run_inner_ac_pf!
+export run_pvpq_switching!, run_uvls!
+export restore_p_balance!, calc_Δp
 
 include("PlotUtils.jl")
 
 include("PlotPG.jl")
-export plot_pg_map, plot_pg_overhead_tl_segments, plot_pg
+export plot_pg_map, plot_pg_overhead_tl_segments, plot_pg, plot_pg_areas
 
 include("Impact.jl")
-export sim_impact
+export sim_primary_dmg
 export calc_overhead_tl_windloads, get_windfield, calc_overhead_tl_segments
 
 include("PlotImpact.jl")

src/Impact.jl

 #* Functions for modelling the impact of hurricanes on power grids.
 #*------------------------------------------------------------------------------
 
-function sim_impact(seg_data::Dict{String,<:Any}, γ::Float64)
+function sim_primary_dmg(seg_data::Dict{String,<:Any}, γ::Float64)
     tl_failures = Array{Tuple{String,Int64},1}() # array for line failures
     
     ### Go through overhead transmission lines and simulate until failure