Introduction

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

This tutorial uses a PTDF to examine current power flow, a LODF to examine how a tripped line's power redirects, and network reductions to make the analysis faster to repeat. Together this will answer:

If a key transmission line trips, which other lines are most at risk of overloading?

using PowerNetworkMatrices
import PowerSystemCaseBuilder
using Logging

Note

PowerSystemCaseBuilder.jl only supplies the ready-made example systems used throughout this documentation. To build a System from your own data, see the PowerSystems.jl documentation.

Step 1 — Load the network

Network matrices are built from a PowerSystems.System. Load an example test system with PowerSystemCaseBuilder.build_system:

sys = with_logger(NullLogger()) do
    PowerSystemCaseBuilder.build_system(
        PowerSystemCaseBuilder.PSSEParsingTestSystems,
        "psse_14_network_reduction_test_system",
    )
end;

Step 2 — Look at the topology

Every matrix below is indexed by an arc tuple(from_bus, to_bus), so the network's topology is exactly the set of arcs the matrices carry. The IncidenceMatrix is that topology on its own — signed bus-to-arc connectivity and nothing else — and get_arc_axis reads out its arcs.

incidence = IncidenceMatrix(sys)
get_arc_axis(incidence)

23-element Vector{Tuple{Int64, Int64}}:
 (110, 111)
 (107, 108)
 (109, 110)
 (109, 114)
 (106, 112)
 (104, 109)
 (118, 104)
 (103, 104)
 (104, 106)
 (102, 103)
 ⋮
 (106, 111)
 (101, 117)
 (102, 104)
 (109, 1001)
 (112, 1002)
 (104, 1001)
 (114, 1002)
 (110, 1002)
 (107, 1001)

Collapsing those arcs per bus turns the same information into a map of the network. get_arc_axis supplies the pairs and get_bus_axis lists the buses in the order the matrices use them.

neighbors = Dict(bus => Int[] for bus in get_bus_axis(incidence))
for (from, to) in get_arc_axis(incidence)
    push!(neighbors[from], to)
    push!(neighbors[to], from)
end

for bus in get_bus_axis(incidence)
    println(bus, " → ", join(sort(neighbors[bus]), ", "))
end

101 → 115, 117
102 → 103, 104, 115
103 → 102, 104, 116
104 → 102, 103, 106, 109, 118, 1001
106 → 104, 111, 112
107 → 108, 1001
108 → 107
109 → 104, 110, 114, 1001
110 → 109, 111, 1002
111 → 106, 110
112 → 106, 1002
114 → 109, 1002
115 → 101, 102
116 → 103
117 → 101, 118
118 → 104, 117
1001 → 104, 107, 109
1002 → 110, 112, 114

Two details of that listing are worth naming before going further. Buses 1001 and 1002 are not in the source data: they are the star points PowerSystems.jl creates for the case's two Transformer3W components, and each appears here as three arcs meeting at a point. Buses 105 and 113 are missing for the opposite reason — each sits behind a zero-impedance branch and has been merged into its neighbor (105 into 104, 113 into 112), which get_bus_reduction_map reports for the get_network_reduction_data the matrix carries:

get_bus_reduction_map(get_network_reduction_data(incidence))

Dict{Int64, Set{Int64}} with 18 entries:
  110  => Set()
  108  => Set()
  111  => Set()
  104  => Set([105])
  114  => Set()
  118  => Set()
  103  => Set()
  109  => Set()
  107  => Set()
  116  => Set()
  106  => Set()
  101  => Set()
  115  => Set()
  117  => Set()
  1001 => Set()
  1002 => Set()
  102  => Set()
  112  => Set([113])

The rest of the tutorial focuses on the line between buses 103 and 104. Reading it off the map above, 103 also reaches 102 and 116, while 104 reaches 102, 106, 109, 118, and the star bus 1001 — so a trip of (103, 104) has several places to push its flow.

Step 3 — Determine important lines with the PTDF

A PTDF (Power Transfer Distribution Factor) matrix is indexed by an arc tuple(from_bus, to_bus) and a bus number, answering what proportion of power injected at a given bus flows through a given branch (when withdrawn at the reference bus).

ptdf = PTDF(sys)

DC_PTDF_Matrix
    Dimension 1, [(110, 111), (107, 108), (109, 110), (109, 114), (106, 112), (104, 109), (118, 104), (103, 104), (104, 106), (102, 103)  …  (103, 116), (106, 111), (101, 117), (102, 104), (109, 1001), (112, 1002), (104, 1001), (114, 1002), (110, 1002), (107, 1001)]
    Dimension 2, [101, 102, 103, 104, 106, 107, 108, 109, 110, 111, 112, 114, 115, 116, 117, 118, 1001, 1002]
Note!! The data shown below corresponds to the transposed matrix.
And data with size (18, 23):
  0.0           0.0   0.0          …   0.0           0.0          0.0
  3.44558e-18   0.0  -1.61695e-20     -3.62198e-17   8.32819e-17  0.0
  1.4042e-17    0.0  -6.58967e-20     -1.47609e-16   3.39404e-16  0.0
  1.81628e-17   0.0  -8.5235e-20      -1.90926e-16   4.39007e-16  0.0
 -0.129584      0.0  -0.0057158       -0.87568       0.123868     0.0
  0.000423477   0.0   1.86791e-5   …   0.0028617    -0.000404798  1.0
  0.000423477  -1.0   1.86791e-5       0.0028617    -0.000404798  1.0
  0.000846954   0.0   3.73582e-5       0.0057234    -0.000809596  0.0
  0.00549369    0.0  -0.00840114      -0.977836      0.986105     0.0
 -0.135703      0.0  -0.0057348       -0.876403      0.129968     0.0
 -0.0011903     0.0  -0.00604335   …  -0.978101     -0.00485305   0.0
  0.00152037    0.0  -0.00370953       0.0141406    -0.0052299    0.0
  1.72279e-18   0.0  -8.08477e-21     -1.81099e-17   4.1641e-17   0.0
  1.4042e-17    0.0  -6.58967e-20     -1.47609e-16   3.39404e-16  0.0
  6.05427e-18   0.0  -2.84117e-20     -6.36422e-17   1.46336e-16  0.0
  1.21085e-17   0.0  -5.68234e-20  …  -1.27284e-16   2.92671e-16  0.0
  0.000423477   0.0   1.86791e-5       0.0028617    -0.000404798  0.0
  0.00194125    0.0  -0.00605134      -0.980599     -0.00799259   0.0

One column of the PTDF corresponds to one bus's influence on every branch. See bus 103's influence on arc(103, 104):

ptdf[(103, 104), 103]

0.7110286823280983

Compare a bus that barely touches this line:

ptdf[(103, 104), 102]

0.011214160544562514

Since power is withdrawn at reference buses, power injected at a reference bus does not flow. get_ref_bus returns those buses, one per electrical island:

ref_buses = get_ref_bus(ptdf)
@assert ptdf[(103, 104), only(ref_buses)] == 0.0
@assert ptdf[(102, 103), only(ref_buses)] == 0.0

Step 4 — Determine post-contingency power flow with the LODF

The LODF (Line Outage Distribution Factor) is indexed by two arc tuples, the monitored arc and the outaged arc, answering how flow is redistributed in an outage.

lodf = LODF(sys)

LODF{Tuple{Vector{Tuple{Int64, Int64}}, Vector{Tuple{Int64, Int64}}}, Tuple{Dict{Tuple{Int64, Int64}, Int64}, Dict{Tuple{Int64, Int64}, Int64}}, Matrix{Float64}}
    Dimension 1, [(110, 111), (107, 108), (109, 110), (109, 114), (106, 112), (104, 109), (118, 104), (103, 104), (104, 106), (102, 103)  …  (103, 116), (106, 111), (101, 117), (102, 104), (109, 1001), (112, 1002), (104, 1001), (114, 1002), (110, 1002), (107, 1001)]
    Dimension 2, [(110, 111), (107, 108), (109, 110), (109, 114), (106, 112), (104, 109), (118, 104), (103, 104), (104, 106), (102, 103)  …  (103, 116), (106, 111), (101, 117), (102, 104), (109, 1001), (112, 1002), (104, 1001), (114, 1002), (110, 1002), (107, 1001)]
Note!! The data shown below corresponds to the transposed matrix.
And data with size (23, 23):
 -1.0           1.34534e-21  -0.00310471   …   0.996895      3.31395e-18
  0.0          -1.0           0.0              0.0           0.0
 -0.00468628    1.48812e-21  -1.0             -0.995314      4.66416e-19
 -0.0800045    -7.79065e-20   0.445148         0.525152      2.61641e-18
 -0.555477      4.28637e-20   0.0014171        0.556894     -1.02022e-17
 -0.000863794   8.09319e-20  -3.8101e-5    …   0.000825693  -3.05189e-17
 -7.47149e-17   1.08007e-18   5.68252e-18      6.31369e-16   6.31369e-16
 -1.45567e-17   7.10212e-20   9.60995e-19     -9.84008e-17  -1.88436e-16
  0.147021     -3.71709e-20   0.00648494      -0.140536      2.09116e-17
 -1.5966e-17   -1.24928e-19  -1.2086e-19      -3.05516e-16   8.18015e-17
  ⋮                                        ⋱                
  1.0           1.38677e-20   0.00310471      -0.996895      3.28787e-18
 -7.64574e-17   3.53249e-19   4.24918e-18  …  -1.53584e-15   2.99533e-16
 -2.42526e-17   2.68139e-19   1.66934e-18     -7.10918e-16   3.75207e-17
  0.0335352    -1.08797e-20   0.0014792       -0.032056     -4.85635e-17
 -0.555477      4.48054e-20   0.0014171        0.556894      1.05173e-18
 -0.0335352     1.08553e-17  -0.0014792        0.032056     -3.62322e-15
 -0.0800045    -6.87383e-20   0.445148     …   0.525152      7.8895e-17
  0.60188       5.97418e-21  -0.39812         -1.0          -6.0992e-19
 -8.60878e-20   7.78701e-20  -3.79723e-21      8.22906e-20  -1.0

We can filter for monitored arcs that absorb flow past a given threshold. First, the monitored arcs, taken from the LODF's arc axis with get_arc_axis:

outaged = (103, 104)
# LODF for monitored arc `arc` when outage arc `outage` is outaged.
monitored = [(arc, lodf[arc, outaged]) for arc in get_arc_axis(lodf) if arc != outaged]

22-element Vector{Tuple{Tuple{Int64, Int64}, Float64}}:
 ((110, 111), -1.455672426015247e-17)
 ((107, 108), 7.102122596786162e-20)
 ((109, 110), 9.609945062769983e-19)
 ((109, 114), 1.4886516198960635e-16)
 ((106, 112), 8.11626105250763e-17)
 ((104, 109), -6.4737388673502484e-18)
 ((118, 104), 0.6489912669795742)
 ((104, 106), -1.2672211256780222e-17)
 ((102, 103), -1.0000000000000002)
 ((117, 118), 0.6489912669795743)
 ⋮
 ((106, 111), 8.908670980981309e-17)
 ((101, 117), 0.6489912669795743)
 ((102, 104), 0.3510087330204253)
 ((109, 1001), -1.8843559903117955e-16)
 ((112, 1002), 3.095443138419165e-16)
 ((104, 1001), -5.963807436568199e-16)
 ((114, 1002), 3.095443138419165e-16)
 ((110, 1002), -9.840083078372377e-17)
 ((107, 1001), -1.8843559903117955e-16)

Then the filter:

filter!(pair -> abs(pair[2]) > 0.05, monitored)
sort!(monitored; by = pair -> -abs(pair[2]))
monitored

7-element Vector{Tuple{Tuple{Int64, Int64}, Float64}}:
 ((102, 103), -1.0000000000000002)
 ((115, 102), -0.648991266979575)
 ((101, 115), -0.648991266979575)
 ((117, 118), 0.6489912669795743)
 ((101, 117), 0.6489912669795743)
 ((118, 104), 0.6489912669795742)
 ((102, 104), 0.3510087330204253)

Branch (102, 103) inherits the entire flow of the outaged line, a factor of -1.0, because it is the series partner on the far side of bus 103. The two parallel paths 101–115–102 and 101–117–118–104 each absorb roughly 65%, and branch (102, 104) takes the remaining 35%. Every other branch is unaffected.

A negative factor means the redistributed flow runs against the monitored branch's (from, to) orientation. Individual factors can also be read directly.

# How much does the arc(102,103) get affected by an outage on arc(103, 104)?
lodf[(102, 103), (103, 104)]

-1.0000000000000002

Step 5 — Reduce the network to repeat the analysis cheaply

Reliability studies are rerun across every credible outage and every operating point, so the matrices should be as small as possible. Networks typically contain buses that do not affect a study of this kind: dead-end (radial) buses and pass-through (degree-two) buses. PowerNetworkMatrices.jl can reduce them away, shrinking the matrices while leaving the surviving sensitivities unchanged.

Reductions are supplied to any constructor — PTDF included — through the network_reductions keyword, as a vector of NetworkReduction specifications. Here we combine RadialReduction, which drops dangling buses, with DegreeTwoReduction, which fuses pass-through chains.

reductions = NetworkReduction[RadialReduction(), DegreeTwoReduction()]
ptdf_reduced = PTDF(sys; network_reductions = reductions)
size(ptdf_reduced), size(ptdf)

((12, 17), (18, 23))

The reduced matrix has fewer branch rows and fewer bus columns, but the sensitivity from Step 3 is unchanged.

ptdf_reduced[(103, 104), 103], ptdf[(103, 104), 103]

(0.7110286823280985, 0.7110286823280983)

The same holds for the LODF screen.

lodf_reduced = LODF(sys; network_reductions = reductions)
lodf_reduced[(102, 103), (103, 104)], lodf[(102, 103), (103, 104)]

(-0.9999999999999992, -1.0000000000000002)

A reduction therefore gives a smaller and faster pair of matrices that answer identically for every element that survives it.

Warning

When screening a reduced network against ratings, a degree-two merge fuses several branches into one equivalent arc whose limit is an aggregate. Use get_single_element_contingency_rating and the related aggregated-rating accessors rather than a single branch's PowerSystems.get_rating.

Note that every lookup above used an arc tuple rather than a branch name. Arc tuples are an unambiguous identifier and they survive reductions, whereas named branches may be merged away, so they are the recommended identifier to use in code that runs both before and after a reduction.

Where to go next

Analysis at Scale — the second tutorial: screening many contingencies on a large network with the memory-light virtual matrices and cache control.
Matrix overview & indexing — the reference for every matrix type, its axes, indexing, and accessors.
The DC Power Flow Approximation — the theory behind these sensitivities.
How to Define and Apply Contingencies — post-contingency factors beyond the single-outage LODF.