Add option to save time series of supply in calc_secondary_dmg!

src/Data.jl

@@ -403,7 +403,7 @@ end
 """
     deactivate_overloads!(network_data::Dict{String,<:Any})
 """
-function deactivate_overloads!(network_data::Dict{String,<:Any})
+function deactivate_overloads!(network_data::Dict{String,<:Any}, f=0)
     ### Get overloaded active branches and deactivate them
     overloaded_br = get_branches(network_data, min_loading=1.)
     if isempty(overloaded_br) == false
@@ -420,6 +420,13 @@ function deactivate_overloads!(network_data::Dict{String,<:Any})
 
     simplify_network!(network_data)
 
+    ### Add currently supplied load to history
+    if haskey(network_data, "history")
+        network_data["history"][f][length(network_data["history"][f])+1] = (
+            get_total_load(network_data)[1], "overload"
+        )
+    end
+
     return overloaded_br
 end
 
@@ -427,8 +434,7 @@ 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-ID's of destroyed transmission lines
+        network_data::Dict{String,<:Any}, tl_ids::Array{String,1}, f = 0
     )
 
     ### Go through transmission lines to destroy
@@ -449,6 +455,13 @@ function destroy_tl!(
 
     simplify_network!(network_data) # remove dangling buses etc.
 
+    ### Add currently supplied load to history
+    if haskey(network_data, "history")
+        network_data["history"][f][length(network_data["history"][f])+1] = (
+            get_total_load(network_data)[1], "wind"
+        )
+    end
+
     return nothing
 end
 

src/PowerFlow.jl

@@ -249,6 +249,7 @@ There are several additional keyword arguments that control the way calculations
 function calc_secondary_dmg!(
         network_data::Dict{String,<:Any},
         primary_dmg::Array{Tuple{String,Int64},1}; # time series of wind damage
+        save_supply_history = true, # whether to save time series of supply
         logpath = "", # optional path to .log file for saving logs of interest
         kwargs... # other keyword arguments (see documentation)
     )
@@ -268,22 +269,29 @@ function calc_secondary_dmg!(
         unique([primary_dmg[i][2] for i in 1:length(primary_dmg)])
     )
 
+    ### Save the time series (history) of supplied load
+    if save_supply_history == true
+        network_data["history"] = Dict{Int64,Dict{Int64,Tuple{Float64,String}}}(
+            f => Dict{Int64,Tuple{Float64,String}}() for f in dmg_frames
+        )
+    end
+
     pf_iter = 0 # increased by 1 every time the power flow is recalculated
     ### Go through wind frames and calculate cascading failures
     for f in dmg_frames
         ### Deactivate transmission lines destroyed by the wind in frame f 
-        destroyed_tl = [dmg[1] for dmg in primary_dmg if dmg[2] == f]
+        destroyed_tl = String[dmg[1] for dmg in primary_dmg if dmg[2] == f]
         info(LOGGER, 
             "Wind frame $f: Destroying overhead transmission lines $destroyed_tl"
         ) # log destroyed transmission lines
-        destroy_tl!(network_data, destroyed_tl) # update network topology
+        destroy_tl!(network_data, destroyed_tl, f) # update network topology
 
         ### Calculate new state of the system
         pf_iter = calc_new_state!(network_data, pf_iter, f; kwargs...)
 
         ### Deactivate any overloaded branches and calculate cascade
-        deactivated_br = deactivate_overloads!(network_data)
         cascade_iter = 0 # count every time overloaded branches fail
+        deactivated_br = deactivate_overloads!(network_data, f)
         while isempty(deactivated_br) == false
             cascade_iter += 1
             info(LOGGER, 
@@ -292,7 +300,7 @@ function calc_secondary_dmg!(
             ### Find new state of the system
             pf_iter = calc_new_state!(network_data, pf_iter, f; kwargs...)
             ### Deactivate new overloaded branches
-            deactivated_br = deactivate_overloads!(network_data)
+            deactivated_br = deactivate_overloads!(network_data, f)
         end
     end
 
@@ -365,7 +373,7 @@ function calc_new_state!(
             save(save_pf_iter[2], "network_data", network_data)
         end
         ### Restore active power balance and calculate DC-PF
-        restore_p_balance!(network_data)
+        restore_p_balance!(network_data, f)
         correct_reference_buses!(network_data)
         info(LOGGER,
             "PF-iteration $pf_iter: Calculating a new DC-PF solution."
@@ -766,13 +774,17 @@ end
 
 Restores active power balance (APB) in the network described by `network_data`.
 """
-function restore_p_balance!(network_data::Dict{String,<:Any}, apb_tol=1e-10)
+function restore_p_balance!(
+        network_data::Dict{String,<:Any}, f = 0; 
+        apb_tol = 1e-10
+    )
+
     ### Calculate connected components
     connected_components = calc_connected_components(network_data)
 
     ### Restore active power balance for every connected component
     for (cc_i, cc) in enumerate(connected_components)
-        _restore_p_balance!(network_data, cc, cc_i, apb_tol)
+        _restore_p_balance!(network_data, cc, cc_i, f, apb_tol)
     end
 
     ### Update network topology (remove dangling buses etc.)
@@ -791,6 +803,7 @@ function _restore_p_balance!(
         network_data::Dict{String,<:Any}, 
         cc::Set{Int64}, # set of bus indices in connected component
         cc_i::Int64, # index/number of connected component
+        f = 0, # wind frame for saving currently supplied load
         apb_tol = 1e-10 # mismatch tolerated for active power balance
     )
     
@@ -843,6 +856,12 @@ function _restore_p_balance!(
                 "Grid #$cc_i: Overproduction is unavoidable => #$cc_i fails!"
             )
             _deactivate_cc!(network_data, cc, active_cc_gens, active_cc_loads)
+            if haskey(network_data, "history")
+                supply_f = network_data["history"][f]
+                supply_f[length(supply_f)+1] = (
+                    get_total_load(network_data)[1], "overproduction"
+                )
+            end
         end
     elseif cc_Δp < -apb_tol # underproduction
         info(LOGGER, "Grid #$cc_i: Detected underproduction: $cc_Δp")
@@ -876,6 +895,12 @@ function _restore_p_balance!(
             info(LOGGER, 
                "Grid #$cc_i: Global active power mismatch after load shedding: $cc_Δp"
             )
+            if haskey(network_data, "history")
+                supply_f = network_data["history"][f]
+                supply_f[length(supply_f)+1] = (
+                    get_total_load(network_data)[1], "load_shedding"
+                )
+            end
         else # sufficient dispatch increase possible
             info(LOGGER, 
                 "Grid #$cc_i: Global active power mismatch requires an increase of generation!"