Notebooks · lib/power/contingency

N-1 contingency analysis by AC re-solves and DC distribution factors

A transmission network is planned to survive the loss of any one element (NERC TPL-001, category P1): with a line, a transformer or a generator out of service, every branch must stay within its emergency rating and every bus voltage within its emergency band. n1_ac checks that the slow way and the exact way: for each outage it takes the element out, solves the AC power flow again by Newton's method starting from the base-case voltages (a warm start, two or three iterations each), and lists every limit the solution breaks.

Each branch's loading is the larger of the apparent powers entering its two ends, |S| in MVA, held to rate A in the base case and rate B after an outage (MATPOWER's columns RATE_A and RATE_B; a rating of 0 is no limit). Voltages are held to [vmin, vmax] before and [vmin_post, vmax_post] after.

An outage that splits the network leaves the buses no longer joined to the reference bus without supply: they are taken out (type 4), with their loads and generators, and counted as islanded. A lost generator's power is made up by the reference bus's generator, as in a power flow; when the lost generator was the reference's, the first remaining PV bus takes over.

dc_factors gives the linear sensitivities of the DC model (MATPOWER's makePTDF and makeLODF):

PTDF[l, j]: the change in the flow of branch l per unit of power injected at bus j and withdrawn at the reference bus;
LODF[l, k]: the share of branch k's pre-outage flow that moves onto branch l when k is taken out, H[l, k] / (1 - H[k, k]) with H[l, k] = PTDF[l, f_k] - PTDF[l, t_k] (and -1 on the diagonal).

With them, dc_screen estimates every post-outage flow in one matrix product and flags the ones above a fraction of rate B, the screening that decides which outages earn a full AC solve. A branch whose outage islands part of the network has 1 - H[k, k] = 0; its LODF column is left at zero and the branch is reported as a bridge.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

This program declares a record type (Limits), and records run only in compiled C for now, so the Run buttons here are off. Open it in the workbench and use its C tab to compile it.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 41–573

from power.case import (
    BR_STATUS, BR_X, BUS_TYPE, Case, GEN_STATUS, GS, NONE, PD, PF, PG, PT, QD, QF, QT, SHIFT, TAP, VA, VG, VM,
    read_case,
)
from power.pf import DC, Net, bus_types, internal, make_sbus, newton_pf, pf_solution, runpf
from power.splu import lu_solve, mindeg_order, splu
from power.ybus import full_network, make_ybus

#: what a contingency takes out
BRANCH_OUT: Int = 1
GEN_OUT: Int = 2
#: the kinds of violation
THERMAL: Int = 1
LOW_V: Int = 2
HIGH_V: Int = 3
DIVERGED: Int = 4
#: the columns of a violation row
V_CONT: Int = 0
V_KIND: Int = 1
V_ELEM: Int = 2
V_VALUE: Int = 3
V_LIMIT: Int = 4
V_PCT: Int = 5
VIOL_W: Int = 6
#: the row of the base case in the result matrices, and its contingency number in violation rows
BASE_ROW: Int = 0
BASE_CASE: Num = -1.0
#: the columns of a screening row: contingency, branch, DC flow, rating
SCREEN_W: Int = 4
#: MATPOWER's rating columns
RATE_A: Int = 5
RATE_B: Int = 6
PIVOT: Num = 0.001
PERCENT: Num = 100.0
#: 1 - H[k, k] below this: the branch is a bridge
BRIDGE_TOL: Num = 1e-10

record Limits {
    vmin: Num
    vmax: Num
    vmin_post: Num
    vmax_post: Num
    rate_pre: Int
    rate_post: Int
}

record N1 {
    #: per row (row 0 the base case, then one per contingency)
    converged: [Int]
    iterations: [Int]
    islanded: [Int]
    lost_mw: Arr
    #: |S| at the more loaded end of each branch, MVA (rows x branches)
    flows: Arr
    #: bus voltage magnitudes, per unit; 0 at a bus left without supply
    vm: Arr
    #: one row per violation: contingency (-1: base), kind, element, value, limit, percent of the limit
    viol: Arr
}

from power.case import (
    BR_STATUS, BR_X, BUS_TYPE, Case, GEN_STATUS, GS, NONE, PD, PF, PG, PT, QD, QF, QT, SHIFT, TAP, VA, VG, VM,
    read_case,
)
from power.pf import DC, Net, bus_types, internal, make_sbus, newton_pf, pf_solution, runpf
from power.splu import lu_solve, mindeg_order, splu
from power.ybus import full_network, make_ybus

#: what a contingency takes out
BRANCH_OUT: Int = 1
GEN_OUT: Int = 2
#: the kinds of violation
THERMAL: Int = 1
LOW_V: Int = 2
HIGH_V: Int = 3
DIVERGED: Int = 4
#: the columns of a violation row
V_CONT: Int = 0
V_KIND: Int = 1
V_ELEM: Int = 2
V_VALUE: Int = 3
V_LIMIT: Int = 4
V_PCT: Int = 5
VIOL_W: Int = 6
#: the row of the base case in the result matrices, and its contingency number in violation rows
BASE_ROW: Int = 0
BASE_CASE: Num = -1.0
#: the columns of a screening row: contingency, branch, DC flow, rating
SCREEN_W: Int = 4
#: MATPOWER's rating columns
RATE_A: Int = 5
RATE_B: Int = 6
PIVOT: Num = 0.001
PERCENT: Num = 100.0
#: 1 - H[k, k] below this: the branch is a bridge
BRIDGE_TOL: Num = 1e-10

record Limits {
    vmin: Num
    vmax: Num
    vmin_post: Num
    vmax_post: Num
    rate_pre: Int
    rate_post: Int
}

record N1 {
    #: per row (row 0 the base case, then one per contingency)
    converged: [Int]
    iterations: [Int]
    islanded: [Int]
    lost_mw: Arr
    #: |S| at the more loaded end of each branch, MVA (rows x branches)
    flows: Arr
    #: bus voltage magnitudes, per unit; 0 at a bus left without supply
    vm: Arr
    #: one row per violation: contingency (-1: base), kind, element, value, limit, percent of the limit
    viol: Arr
}