Notebooks · lib/power/study

Violations, distribution factors and the transfer limit in a generator-interconnection study

A transmission owner asked to connect a new generator at a point of interconnection (POI) runs a system impact study: the network is solved with and without the project, in the base case and under every single contingency, and the project is charged with the limits it breaks that were not broken before. study does the steady-state part:

the project is a generator at the POI bus of mw megawatts. With a reactive range (qmin < qmax) it holds the POI at its pre-project voltage as a PV bus; with qmin = qmax it injects that fixed Q;
its power is absorbed by the reference bus (sink 1) or by every other generator in proportion to its output (sink 2), the usual stand-ins for the generation the project displaces;
n1_ac (power.contingency) runs the base case and the contingencies before and after; a violation (contingency, kind, element) present after and not before is new; one present in both whose loading rose by more than worse_pct percentage points is worsened;
the distribution factors are the DC model's flow change on every branch per MW of the project's output, injected at the POI and taken out at the sink;
the transfer limit is the largest output with no new violation, found by bisection on the full AC N-1 between 0 and mw_max to within tol_mw.

The contingency numbering is the same before and after: the project's generator is appended after the case's own generators, so every branch and generator keeps its row.

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 (Study), 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 33–242

from power.case import BUS_I, BUS_TYPE, Case, GEN_STATUS, PG, PQ, PV, VM, read_case
from power.contingency import (
    Limits, N1, V_CONT, V_ELEM, V_KIND, V_PCT, all_n1, dc_factors, default_limits,
    n1_ac,
)
from power.pf import internal, runpf, NR

#: where the project's power goes
SINK_SLACK: Int = 1
SINK_GENS: Int = 2
#: MATPOWER's generator matrix width and the columns written here
GEN_W: Int = 21
G_QMAX: Int = 3
G_QMIN: Int = 4
G_VG: Int = 5
G_MBASE: Int = 6
G_PMAX: Int = 8
G_QG: Int = 2

record Study {
    pre: N1
    post: N1
    #: per violation row of post: 1 new, 2 worsened, 0 neither
    status: [Int]
    df: Arr
    limit_mw: Num
    limit_found: Bool
    #: the bisection: rows (MW, new violations)
    steps: Arr
}

#| The case with the project added: a generator of `mw` MW at bus row
#| `poi` (PV at `vset` within [qmin, qmax], or a fixed Q = qmin when
#| qmin == qmax), its power taken from the sink.
def with_injection(c: Case, poi: Int, mw: Num, vset: Num, qmin: Num, qmax: Num, sink: Int) -> Case {
    ng: Int = int(shape(c.gen)[0])
    nc: Int = int(shape(c.gen)[1])
    gen = zeros(ng + 1, max(nc, GEN_W))
    for g to ng {
        for j to nc {
            gen[g, j] = c.gen[g, j]
        }
    }
    bus = copy(c.bus)
    fixed_q = qmin == qmax
    gen[ng, 0] = bus[poi, BUS_I]
    gen[ng, PG] = mw
    gen[ng, G_QG] = qmin if fixed_q else 0.0
    gen[ng, G_QMAX] = qmax
    gen[ng, G_QMIN] = qmin
    gen[ng, G_VG] = vset
    gen[ng, G_MBASE] = c.base_mva
    gen[ng, GEN_STATUS] = 1.0
    gen[ng, G_PMAX] = mw
    if not fixed_q and int(bus[poi, BUS_TYPE]) == PQ {
        bus[poi, BUS_TYPE] = 1.0 * PV
    }
    if fixed_q and int(bus[poi, BUS_TYPE]) == PQ {
        #: a PQ bus keeps its type; the generator's fixed Q is part of its injection
        gen[ng, G_VG] = bus[poi, VM]
    }
    if sink == SINK_GENS {
        total = 0.0
        for g to ng {
            if gen[g, GEN_STATUS] > 0.0 and gen[g, PG] > 0.0 {
                total += gen[g, PG]
            }
        }
        for g to ng {
            if gen[g, GEN_STATUS] > 0.0 and gen[g, PG] > 0.0 {
                gen[g, PG] -= mw * gen[g, PG] / total
            }
        }
    }
    return Case { base_mva: c.base_mva, bus: bus, gen: gen, branch: c.branch }
}

#| The DC flow change on every branch per MW of the project (injected at
#| bus row poi, taken out at the sink).
def distribution_factors(c: Case, poi: Int, sink: Int) -> Arr {
    d = dc_factors(c)
    nl: Int = int(shape(d.ptdf)[0])
    df = zeros(nl)
    for l to nl {
        df[l] = d.ptdf[l, poi]
    }
    if sink == SINK_GENS {
        net = internal(c)
        total = 0.0
        for g to len(net.gbus) {
            if net.gen[g, GEN_STATUS] > 0.0 and net.gen[g, PG] > 0.0 {
                total += net.gen[g, PG]
            }
        }
        for g to len(net.gbus) {
            if net.gen[g, GEN_STATUS] > 0.0 and net.gen[g, PG] > 0.0 {
                w = net.gen[g, PG] / total
                for l to nl {
                    df[l] -= w * d.ptdf[l, net.gbus[g]]
                }
            }
        }
    }
    return df
}

#| For each violation row of `post`: 1 when (contingency, kind, element)
#| is not among `pre`'s, 2 when it is but its percentage rose by more
#| than worse_pct, else 0.
def classify(pre: Arr, post: Arr, worse_pct: Num) -> [Int] {
    known = map()
    npre: Int = int(shape(pre)[0])
    for i to npre {
        known[_key(pre, i)] = pre[i, V_PCT]
    }
    n: Int = int(shape(post)[0])
    out = int_zeros(n)
    for i to n {
        k = _key(post, i)
        if not has(known, k) {
            out[i] = 1
        }
        else if post[i, V_PCT] > known[k] + worse_pct {
            out[i] = 2
        }
    }
    return out
}

def _key(v: Arr, i: Int) -> Str {
    return f"{int(v[i, V_CONT])}:{int(v[i, V_KIND])}:{int(v[i, V_ELEM])}"
}

def _count_new(pre: Arr, post: Arr) -> Int {
    s = classify(pre, post, inf)
    n: Int = 0
    for k to len(s) {
        if s[k] == 1 {
            n += 1
        }
    }
    return n
}

#| The study of a project of `mw` MW at bus row `poi`, and its transfer
#| limit up to `mw_max`.
def study(c: Case, poi: Int, mw: Num, qmin: Num, qmax: Num, sink: Int, kinds: [Int], elems: [Int], lim: Limits,
          mw_max: Num, tol_mw: Num, worse_pct: Num) -> Study {
    base = runpf(c, NR, false, 1e-8, 10)
    vset = base.bus[poi, VM]
    pre = n1_ac(c, kinds, elems, lim, 1e-8, 10)
    post = n1_ac(with_injection(c, poi, mw, vset, qmin, qmax, sink), kinds, elems, lim, 1e-8, 10)
    status = classify(pre.viol, post.viol, worse_pct)
    df = distribution_factors(c, poi, sink)

    #: bisection on the output: lo has no new violation, hi has one
    steps = []
    hi_r = n1_ac(with_injection(c, poi, mw_max, vset, qmin, qmax, sink), kinds, elems, lim, 1e-8, 10)
    n_hi: Int = _count_new(pre.viol, hi_r.viol)
    steps = append(steps, [mw_max, 1.0 * n_hi])
    found = n_hi > 0
    lo = 0.0
    hi = mw_max
    if found {
        while hi - lo > tol_mw {
            mid = 0.5 * (lo + hi)
            r = n1_ac(with_injection(c, poi, mid, vset, qmin, qmax, sink), kinds, elems, lim, 1e-8, 10)
            n: Int = _count_new(pre.viol, r.viol)
            steps = append(steps, [mid, 1.0 * n])
            if n > 0 {
                hi = mid
            }
            else {
                lo = mid
            }
        }
    }
    limit = lo if found else mw_max
    return Study { pre: pre, post: post, status: status, df: df, limit_mw: limit, limit_found: found,
                   steps: steps }
}

test "case30, 60 MW at bus 30 sunk at the reference: the DC factors sum to one over the cut into bus 30" {
    c = read_case("lib/power/cases/case30.m")
    df = distribution_factors(c, 29, SINK_SLACK)
    #: bus 30 is fed by branches 27-30 and 29-30 (rows 37 and 38): together they carry the whole MW back
    assert_near(df[37] + df[38], -1.0, 1e-12)
}

test "case30, a project at bus 30: no new violation at the limit, one within tol_mw above it" {
    c = read_case("lib/power/cases/case30.m")
    kinds, elems = all_n1(c)
    s = study(c, 29, 30.0, -20.0, 20.0, SINK_SLACK, kinds, elems, default_limits(), 200.0, 8.0, 1.0)
    assert_true(s.limit_found)
    #: the bisection's last bracket: the smallest output tried that broke a limit is within tol_mw of the limit
    hi = 200.0
    for k to int(shape(s.steps)[0]) {
        if s.steps[k, 1] > 0.0 and s.steps[k, 0] > s.limit_mw {
            hi = min(hi, s.steps[k, 0])
        }
    }
    assert_true(hi - s.limit_mw <= 8.0)
    base = runpf(c, NR, false, 1e-8, 10)
    at = n1_ac(with_injection(c, 29, s.limit_mw, base.bus[29, VM], -20.0, 20.0, SINK_SLACK), kinds, elems,
               default_limits(), 1e-8, 10)
    assert_true(_count_new(s.pre.viol, at.viol) == 0)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

from power.case import BUS_I, BUS_TYPE, Case, GEN_STATUS, PG, PQ, PV, VM, read_case
from power.contingency import (
    Limits, N1, V_CONT, V_ELEM, V_KIND, V_PCT, all_n1, dc_factors, default_limits,
    n1_ac,
)
from power.pf import internal, runpf, NR

#: where the project's power goes
SINK_SLACK: Int = 1
SINK_GENS: Int = 2
#: MATPOWER's generator matrix width and the columns written here
GEN_W: Int = 21
G_QMAX: Int = 3
G_QMIN: Int = 4
G_VG: Int = 5
G_MBASE: Int = 6
G_PMAX: Int = 8
G_QG: Int = 2

record Study {
    pre: N1
    post: N1
    #: per violation row of post: 1 new, 2 worsened, 0 neither
    status: [Int]
    df: Arr
    limit_mw: Num
    limit_found: Bool
    #: the bisection: rows (MW, new violations)
    steps: Arr
}

#| The case with the project added: a generator of `mw` MW at bus row
#| `poi` (PV at `vset` within [qmin, qmax], or a fixed Q = qmin when
#| qmin == qmax), its power taken from the sink.
def with_injection(c: Case, poi: Int, mw: Num, vset: Num, qmin: Num, qmax: Num, sink: Int) -> Case {
    ng: Int = int(shape(c.gen)[0])
    nc: Int = int(shape(c.gen)[1])
    gen = zeros(ng + 1, max(nc, GEN_W))
    for g to ng {
        for j to nc {
            gen[g, j] = c.gen[g, j]
        }
    }
    bus = copy(c.bus)
    fixed_q = qmin == qmax
    gen[ng, 0] = bus[poi, BUS_I]
    gen[ng, PG] = mw
    gen[ng, G_QG] = qmin if fixed_q else 0.0
    gen[ng, G_QMAX] = qmax
    gen[ng, G_QMIN] = qmin
    gen[ng, G_VG] = vset
    gen[ng, G_MBASE] = c.base_mva
    gen[ng, GEN_STATUS] = 1.0
    gen[ng, G_PMAX] = mw
    if not fixed_q and int(bus[poi, BUS_TYPE]) == PQ {
        bus[poi, BUS_TYPE] = 1.0 * PV
    }
    if fixed_q and int(bus[poi, BUS_TYPE]) == PQ {
        #: a PQ bus keeps its type; the generator's fixed Q is part of its injection
        gen[ng, G_VG] = bus[poi, VM]
    }
    if sink == SINK_GENS {
        total = 0.0
        for g to ng {
            if gen[g, GEN_STATUS] > 0.0 and gen[g, PG] > 0.0 {
                total += gen[g, PG]
            }
        }
        for g to ng {
            if gen[g, GEN_STATUS] > 0.0 and gen[g, PG] > 0.0 {
                gen[g, PG] -= mw * gen[g, PG] / total
            }
        }
    }
    return Case { base_mva: c.base_mva, bus: bus, gen: gen, branch: c.branch }
}

#| The DC flow change on every branch per MW of the project (injected at
#| bus row poi, taken out at the sink).
def distribution_factors(c: Case, poi: Int, sink: Int) -> Arr {
    d = dc_factors(c)
    nl: Int = int(shape(d.ptdf)[0])
    df = zeros(nl)
    for l to nl {
        df[l] = d.ptdf[l, poi]
    }
    if sink == SINK_GENS {
        net = internal(c)
        total = 0.0
        for g to len(net.gbus) {
            if net.gen[g, GEN_STATUS] > 0.0 and net.gen[g, PG] > 0.0 {
                total += net.gen[g, PG]
            }
        }
        for g to len(net.gbus) {
            if net.gen[g, GEN_STATUS] > 0.0 and net.gen[g, PG] > 0.0 {
                w = net.gen[g, PG] / total
                for l to nl {
                    df[l] -= w * d.ptdf[l, net.gbus[g]]
                }
            }
        }
    }
    return df
}

#| For each violation row of `post`: 1 when (contingency, kind, element)
#| is not among `pre`'s, 2 when it is but its percentage rose by more
#| than worse_pct, else 0.
def classify(pre: Arr, post: Arr, worse_pct: Num) -> [Int] {
    known = map()
    npre: Int = int(shape(pre)[0])
    for i to npre {
        known[_key(pre, i)] = pre[i, V_PCT]
    }
    n: Int = int(shape(post)[0])
    out = int_zeros(n)
    for i to n {
        k = _key(post, i)
        if not has(known, k) {
            out[i] = 1
        }
        else if post[i, V_PCT] > known[k] + worse_pct {
            out[i] = 2
        }
    }
    return out
}

def _key(v: Arr, i: Int) -> Str {
    return f"{int(v[i, V_CONT])}:{int(v[i, V_KIND])}:{int(v[i, V_ELEM])}"
}

def _count_new(pre: Arr, post: Arr) -> Int {
    s = classify(pre, post, inf)
    n: Int = 0
    for k to len(s) {
        if s[k] == 1 {
            n += 1
        }
    }
    return n
}

#| The study of a project of `mw` MW at bus row `poi`, and its transfer
#| limit up to `mw_max`.
def study(c: Case, poi: Int, mw: Num, qmin: Num, qmax: Num, sink: Int, kinds: [Int], elems: [Int], lim: Limits,
          mw_max: Num, tol_mw: Num, worse_pct: Num) -> Study {
    base = runpf(c, NR, false, 1e-8, 10)
    vset = base.bus[poi, VM]
    pre = n1_ac(c, kinds, elems, lim, 1e-8, 10)
    post = n1_ac(with_injection(c, poi, mw, vset, qmin, qmax, sink), kinds, elems, lim, 1e-8, 10)
    status = classify(pre.viol, post.viol, worse_pct)
    df = distribution_factors(c, poi, sink)

#: bisection on the output: lo has no new violation, hi has one
    steps = []
    hi_r = n1_ac(with_injection(c, poi, mw_max, vset, qmin, qmax, sink), kinds, elems, lim, 1e-8, 10)
    n_hi: Int = _count_new(pre.viol, hi_r.viol)
    steps = append(steps, [mw_max, 1.0 * n_hi])
    found = n_hi > 0
    lo = 0.0
    hi = mw_max
    if found {
        while hi - lo > tol_mw {
            mid = 0.5 * (lo + hi)
            r = n1_ac(with_injection(c, poi, mid, vset, qmin, qmax, sink), kinds, elems, lim, 1e-8, 10)
            n: Int = _count_new(pre.viol, r.viol)
            steps = append(steps, [mid, 1.0 * n])
            if n > 0 {
                hi = mid
            }
            else {
                lo = mid
            }
        }
    }
    limit = lo if found else mw_max
    return Study { pre: pre, post: post, status: status, df: df, limit_mw: limit, limit_found: found,
                   steps: steps }
}

test "case30, 60 MW at bus 30 sunk at the reference: the DC factors sum to one over the cut into bus 30" {
    c = read_case("lib/power/cases/case30.m")
    df = distribution_factors(c, 29, SINK_SLACK)
    #: bus 30 is fed by branches 27-30 and 29-30 (rows 37 and 38): together they carry the whole MW back
    assert_near(df[37] + df[38], -1.0, 1e-12)
}

test "case30, a project at bus 30: no new violation at the limit, one within tol_mw above it" {
    c = read_case("lib/power/cases/case30.m")
    kinds, elems = all_n1(c)
    s = study(c, 29, 30.0, -20.0, 20.0, SINK_SLACK, kinds, elems, default_limits(), 200.0, 8.0, 1.0)
    assert_true(s.limit_found)
    #: the bisection's last bracket: the smallest output tried that broke a limit is within tol_mw of the limit
    hi = 200.0
    for k to int(shape(s.steps)[0]) {
        if s.steps[k, 1] > 0.0 and s.steps[k, 0] > s.limit_mw {
            hi = min(hi, s.steps[k, 0])
        }
    }
    assert_true(hi - s.limit_mw <= 8.0)
    base = runpf(c, NR, false, 1e-8, 10)
    at = n1_ac(with_injection(c, 29, s.limit_mw, base.bus[29, VM], -20.0, 20.0, SINK_SLACK), kinds, elems,
               default_limits(), 1e-8, 10)
    assert_true(_count_new(s.pre.viol, at.viol) == 0)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors