GitHub

Analysis at Scale

To follow along, you can download this tutorial as a Julia script (.jl) or Jupyter notebook (.ipynb).

A dense PTDF or LODF is an O(N²) array of Float64, which on a real interconnection is tens of gigabytes, and a screen reads each row exactly once. This tutorial makes use of virtual matrices that compute rows on demand to answer the following question:

Across every single-line outage, which contingencies drive a surviving line closest to — or past — its limit?

Note

We use a small network (73 buses) for demonstration purposes.

using PowerNetworkMatrices
import PowerSystems
import PowerSystemCaseBuilder
using LinearAlgebra: dot
using DataFrames

Load an example test system with PowerSystemCaseBuilder.build_system:

sys = PowerSystemCaseBuilder.build_system(
    PowerSystemCaseBuilder.PSISystems,
    "RTS_GMLC_DA_sys",
);
[ Info: Loaded time series from storage file existing=/home/runner/.julia/packages/PowerSystemCaseBuilder/zuCHR/data/serialized_system/e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855/RTS_GMLC_DA_sys_time_series_storage.h5 new=/tmp/jl_JtYkD4 compression=InfrastructureSystems.CompressionSettings(false, InfrastructureSystems.CompressionTypesModule.CompressionTypes.DEFLATE = 1, 3, true)
[ Info: Unit System changed to InfrastructureSystems.UnitSystemModule.UnitSystem.SYSTEM_BASE = 0

Step 1 — Set up the study: base-case flows and limits

A post-contingency flow requires the flow each line carries at present and the limit it must stay under.

Base-case flows. Under the DC approximation a line's flow is its PTDF row dotted with the vector of net bus injections, that is, generation minus load. VirtualPTDF supplies those rows one at a time instead of materializing the whole matrix.

vptdf = VirtualPTDF(sys)
DC_vPTDF_Matrix{PowerNetworkMatrices.KLUWrapper.KLULinSolveCache{Float64, Int64}}
    Dimension 1, [(110, 111), (315, 316), (113, 215), (111, 114), (106, 110), (305, 310), (204, 209), (208, 209), (103, 124), (217, 222)  …  (318, 321), (320, 323), (119, 120), (218, 221), (115, 121), (120, 123), (220, 223), (118, 121), (315, 321), (319, 320)]
    Dimension 2, [101, 102, 103, 104, 105, 106, 107, 108, 109, 110  …  316, 317, 318, 319, 320, 321, 322, 323, 324, 325]
And data with size (73, 108)

The injection vector has to line up with the matrix's bus dimension, so get_bus_lookup is what maps a bus number to its position in that vector:

bus_lookup = get_bus_lookup(vptdf)
Dict{Int64, Int64} with 73 entries:
  111 => 11
  107 => 7
  103 => 3
  115 => 15
  119 => 19
  205 => 29
  306 => 54
  302 => 50
  209 => 33
  124 => 24
  213 => 37
  222 => 46
  319 => 67
  310 => 58
  311 => 59
  201 => 25
  318 => 66
  214 => 38
  123 => 23
  ⋮   => ⋮

We build the injection vector from the system, accumulating per bus because a bus can host several generators and loads, ordered to match the matrix's bus axis.

injection = zeros(Float64, length(bus_lookup))
for gen in PowerSystems.get_components(
    d -> !isa(
        d,
        Union{PowerSystems.ElectricLoad, PowerSystems.SynchronousCondenser},
    ),
    PowerSystems.StaticInjection,
    sys,
)
    PowerSystems.get_available(gen) || continue
    bus_number = PowerSystems.get_number(PowerSystems.get_bus(gen))
    injection[bus_lookup[bus_number]] += PowerSystems.get_active_power(gen)
end
for load in PowerSystems.get_components(
    d -> !isa(d, PowerSystems.FixedAdmittance),
    PowerSystems.ElectricLoad,
    sys,
)
    PowerSystems.get_available(load) || continue
    bus_number = PowerSystems.get_number(PowerSystems.get_bus(load))
    injection[bus_lookup[bus_number]] -= PowerSystems.get_active_power(load)
end
injection
73-element Vector{Float64}:
  0.5999999999999999
  0.7100000000000002
 -1.8
 -0.74
 -0.71
 -1.36
  2.3
 -1.71
 -1.75
 -1.95
  ⋮
  0.0
  0.21999999999999975
 -1.81
 -1.28
  3.55
  3.0999999999999996
  7.1
  0.0
  0.0

The system defaults to its per-unit system base, so these injections and the ratings read below are already on the same 100-MVA base and are directly comparable. The injections do not need to sum to zero because the reference bus balances the difference.

The base flow on every line follows. get_arc_axis enumerates the arcs to sweep, and each vptdf[arc, :] computes that line's PTDF row on first access and caches it, so one pass touches each row once.

arcs = get_arc_axis(vptdf)
base_flow = Dict(arc => dot(vptdf[arc, :], injection) for arc in arcs)
Dict{Tuple{Int64, Int64}, Float64} with 108 entries:
  (110, 111) => -1.26186
  (315, 316) => 0.328842
  (113, 215) => -1.69168
  (111, 114) => -0.916528
  (106, 110) => -0.898257
  (305, 310) => -0.412374
  (204, 209) => -0.444902
  (215, 221) => -3.76631
  (119, 120) => -1.1865
  (208, 209) => -1.04594
  (103, 124) => -1.98655
  (217, 222) => -1.33286
  (209, 211) => -1.24719
  (109, 112) => -1.21264
  (107, 203) => 0.530554
  (116, 117) => -2.61428
  (311, 313) => -1.97578
  (115, 121) => -3.8223
  (211, 214) => -0.706016
  ⋮          => ⋮
Note

Building every base flow this way touches the whole PTDF, one row at a time. When only the base flows are needed, the sparse ABA/BA DC solve (ABA_Matrix, BA_Matrix) obtains them in a single factorization; the row-at-a-time route here is what the N-1 screen below requires anyway.

Line limits. Each branch's rating is read and keyed by arc. Parallel branches share an arc, so their ratings are summed into the combined corridor limit.

line_rating = Dict{Tuple{Int, Int}, Float64}()
for branch in PowerSystems.get_components(PowerSystems.ACTransmission, sys)
    arc = PowerSystems.get_arc(branch)
    key = (
        PowerSystems.get_number(PowerSystems.get_from(arc)),
        PowerSystems.get_number(PowerSystems.get_to(arc)),
    )
    line_rating[key] = get(line_rating, key, 0.0) + PowerSystems.get_rating(branch)
end

Step 2 — Screen every contingency on a VirtualLODF

The LODF gives the redistribution: when outaged trips, line monitored picks up LODF[monitored, outaged] of the outaged line's pre-trip flow. The post-contingency flow on monitored is therefore base_flow[monitored] + LODF[monitored, outaged] · base_flow[outaged], and its loading is that quantity over its rating.

VirtualLODF is constructed and indexed exactly like the dense LODF, with the same constructor and the same [monitored, outaged] indexing, but never forms the whole matrix. max_cache_size caps the row cache in MiB.

vlodf = VirtualLODF(sys; max_cache_size = 100)
VirtualLODF{Tuple{Vector{Tuple{Int64, Int64}}, Vector{Tuple{Int64, Int64}}}, Tuple{Dict{Tuple{Int64, Int64}, Int64}, Dict{Tuple{Int64, Int64}, Int64}}, PowerNetworkMatrices.KLUWrapper.KLULinSolveCache{Float64, Int64}}
    Dimension 1, [(110, 111), (315, 316), (113, 215), (111, 114), (106, 110), (305, 310), (204, 209), (208, 209), (103, 124), (217, 222)  …  (318, 321), (320, 323), (119, 120), (218, 221), (115, 121), (120, 123), (220, 223), (118, 121), (315, 321), (319, 320)]
    Dimension 2, [(110, 111), (315, 316), (113, 215), (111, 114), (106, 110), (305, 310), (204, 209), (208, 209), (103, 124), (217, 222)  …  (318, 321), (320, 323), (119, 120), (218, 221), (115, 121), (120, 123), (220, 223), (118, 121), (315, 321), (319, 320)]
And data with size (108, 108)

A row, meaning one monitored line's factors against every outage, is the unit the cache stores, so the sweep proceeds row by row: compute each monitored line's row once, then score it against every outage. Both dimensions of a LODF are arcs, so get_arc_axis lists the lines to sweep and get_arc_lookup turns an outaged arc into its position within a row. For each outage we keep the single worst-loaded survivor.

lines = get_arc_axis(vlodf)
outage_col = get_arc_lookup(vlodf)

worst = Dict{Tuple{Int, Int}, @NamedTuple{line::Tuple{Int, Int}, loading::Float64}}()
for monitored in lines
    factors = vlodf[monitored, :]
    rating = line_rating[monitored]
    f_m = base_flow[monitored]
    for outaged in lines
        outaged == monitored && continue
        post = f_m + factors[outage_col[outaged]] * base_flow[outaged]
        loading = abs(post) / rating
        if !haskey(worst, outaged) || loading > worst[outaged].loading
            worst[outaged] = (line = monitored, loading = loading)
        end
    end
end

Ranking the outages by the worst loading they cause gives the screen's result.

screen = sort(
    DataFrame(;
        outaged_line = [o for o in keys(worst)],
        worst_monitored = [w.line for w in values(worst)],
        loading = [round(w.loading; digits = 3) for w in values(worst)],
    ),
    :loading; rev = true,
)
98 rows omitted
outaged_lineworst_monitoredloading
Tuple{Int64, Int64}Tuple{Int64, Int64}Float64
(107, 203)(107, 108)1.314
(107, 108)(107, 203)1.314
(315, 321)(316, 317)1.3
(115, 121)(116, 117)1.196
(113, 215)(107, 108)1.194
(103, 124)(107, 108)1.116
(115, 124)(107, 108)1.116
(114, 116)(107, 108)1.099
(112, 123)(107, 108)1.096
(110, 112)(107, 108)1.094
⋮⋮⋮

The day-ahead schedule is not N-1 secure: most single-line outages here push some surviving line over its rating.

count(>(1.0), screen.loading)
86

The worst contingency drives a line to roughly 1.3× its limit. The top pair is reciprocal — tripping (107, 108) overloads (107, 203) and vice versa — because they form a tightly coupled corridor in which each inherits essentially the entire flow of the other. The numbers are identical to those of a dense LODF; the difference is that no dense matrix was ever stored.

The sweep visited every row, so all of them are now cached. For a virtual matrix, get_lodf_data returns exactly the rows the cache is holding, keyed by row index:

length(get_lodf_data(vlodf))
108

On RTS-GMLC that is approximately 108 short rows, a negligible memory footprint, which is why nothing was evicted. At realistic scale the situation is the opposite: the rows do not all fit, and max_cache_size is a hard ceiling. Once the cache is full the least-recently-used row is dropped. A full screen still completes, trading a bounded memory footprint for the recomputation of an evicted row should the screen return to it. That trade is what allows an N-1 screen to run on a grid whose dense LODF would not fit.

Step 3 — Reuse rows across operating points

This screen is not run once. It reruns at every operating point, but the LODF is a property of the topology rather than the dispatch: the factors do not change from hour to hour, only the base flows they multiply. A study that reruns the screen therefore wants its rows to remain resident rather than be recomputed on each pass.

In practice a defined set of facilities is monitored every cycle, here the inter-area tie corridors. Declaring them up front as the VirtualLODF constructor's persistent_arcs holds those rows in the cache and makes them exempt from eviction, so no amount of churn from the rest of a screen can force them to be re-solved.

tie_lines = [(107, 203), (113, 215), (123, 217)]
vlodf_watch = VirtualLODF(sys; persistent_arcs = tie_lines, max_cache_size = 100)
VirtualLODF{Tuple{Vector{Tuple{Int64, Int64}}, Vector{Tuple{Int64, Int64}}}, Tuple{Dict{Tuple{Int64, Int64}, Int64}, Dict{Tuple{Int64, Int64}, Int64}}, PowerNetworkMatrices.KLUWrapper.KLULinSolveCache{Float64, Int64}}
    Dimension 1, [(110, 111), (315, 316), (113, 215), (111, 114), (106, 110), (305, 310), (204, 209), (208, 209), (103, 124), (217, 222)  …  (318, 321), (320, 323), (119, 120), (218, 221), (115, 121), (120, 123), (220, 223), (118, 121), (315, 321), (319, 320)]
    Dimension 2, [(110, 111), (315, 316), (113, 215), (111, 114), (106, 110), (305, 310), (204, 209), (208, 209), (103, 124), (217, 222)  …  (318, 321), (320, 323), (119, 120), (218, 221), (115, 121), (120, 123), (220, 223), (118, 121), (315, 321), (319, 320)]
And data with size (108, 108)

Across repeated operating points only the base flows are recomputed, which is cheap, and the tie-line factors are read straight from the pinned rows. VirtualPTDF takes the same keyword. Pinned rows still count against the budget, so the constructor errors if the pinned set alone would exceed max_cache_size.

Step 4 — Reclaim the memory

For VirtualPTDF and VirtualLODF, the row cache is emptied in place to free it. Pinned rows are removed too, with a warning.

empty!(vlodf.cache)

Where to go next