Notebooks · lib/power/cpf

The continuation power flow, a PV curve up to the nose and past it

Push more and more power along one direction (a project's output at the POI, taken out at a sink) and the power flow's solution moves: the voltages sag, and at some transfer the solution stops existing. Near that point, the nose of the PV curve, Newton's method on the plain power flow fails, because its Jacobian becomes singular there: the curve turns back and two solutions (the normal one and a low-voltage one) merge.

The continuation power flow (V. Ajjarapu and C. Christy, "The continuation power flow: a tool for steady state voltage stability analysis", IEEE Trans. Power Systems 7 (1992) 416-423; the pseudo-arclength form of MATPOWER's runcpf) traces the curve through the nose by making the transfer lambda an unknown and adding one equation:

F(x, lambda) = S(V) - S0 - lambda d = 0       (the power flow, x = Va, Vm)
z' ([x; lambda] - [x_k; lambda_k]) = sigma    (the step along the curve)

where z is the unit tangent at the last point, found from [J, -d; z_prev'] z = [0; 1]. Each step predicts x_k + sigma z and corrects by Newton on the augmented system, whose matrix stays nonsingular at the nose. The step grows after an easy correction and halves after a failed one. When a step lands past the nose (lambda falls after rising), the trace returns to the point before the last one and approaches again with a quarter of the step, until the step is below nose_tol; the nose is then known to within about the square of that step. The trace stops once lambda has fallen back below stop_frac of its largest value, or after max_steps.

Here d is +1 per unit at the POI, made up by the reference bus (sink 1) or by the generators in proportion to their output (sink 2), so lambda is the project's output in per unit of baseMVA. Generator reactive limits are not enforced along the curve.

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 (Cpf), 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 38–336

from power.case import Case, GEN_STATUS, PG, VM, read_case
from power.pf import NR, bus_types, initial_v, internal, jacobian, make_sbus, mismatch, runpf
from power.splu import lu_solve, mindeg_order, splu
from power.ybus import Ybus, full_network, make_ybus
from power.study import SINK_GENS, SINK_SLACK

PIVOT: Num = 0.001
#: the step grows by this factor after a correction of at most EASY iterations
GROW: Num = 1.5
EASY: Int = 3
#: the smallest step tried before the trace gives up
SIGMA_MIN: Num = 1e-6
#: past the nose, the trace steps back and shrinks the step by this factor until it is below nose_tol
NOSE_SHRINK: Num = 4.0

record Cpf {
    #: per point: lambda (per unit), Vm at the POI, the smallest Vm, corrector iterations
    lam: Arr
    v_poi: Arr
    v_min: Arr
    #: the voltages at every point (points x buses)
    vm: Arr
    va: Arr
    lam_max: Num
    nose: Int
    #: the smallest singular value of the power-flow Jacobian at each point
    smin: Arr
}

#| The transfer direction: +1 at bus row poi, made up by the sink.
def direction(net_gen: Arr, gbus: [Int], nb: Int, poi: Int, sink: Int) -> Arr {
    d = zeros(nb)
    d[poi] = 1.0
    if sink == SINK_GENS {
        total = 0.0
        for g to len(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(gbus) {
            if net_gen[g, GEN_STATUS] > 0.0 and net_gen[g, PG] > 0.0 {
                d[gbus[g]] -= net_gen[g, PG] / total
            }
        }
    }
    return d
}

#| The augmented matrix [[J, col]; [row']] in compressed sparse columns.
def augment(n: Int, cp: [Int], ri: [Int], vals: Arr, col: Arr, row: Arr) {
    m: Int = n + 1
    acp = int_zeros(m + 1)
    ari: [Int] = []
    av = []
    for j to n {
        acp[j] = len(ari)
        for p from cp[j] to cp[j + 1] {
            ari = append(ari, ri[p])
            av = append(av, vals[p])
        }
        ari = append(ari, n)
        av = append(av, row[j])
    }
    acp[n] = len(ari)
    for i to n {
        if col[i] != 0.0 {
            ari = append(ari, i)
            av = append(av, col[i])
        }
    }
    ari = append(ari, n)
    av = append(av, row[n])
    acp[m] = len(ari)
    return acp, ari, av
}

def _dense(n: Int, cp: [Int], ri: [Int], vals: Arr) -> Arr {
    a = zeros(n, n)
    for j to n {
        for p from cp[j] to cp[j + 1] {
            a[ri[p], j] += vals[p]
        }
    }
    return a
}

#| The PV curve for a transfer into bus row poi, from the case's power
#| flow solution, with the first step sigma0 (per unit of the arc length).
def cpf(c: Case, poi: Int, sink: Int, sigma0: Num, sigma_max: Num, stop_frac: Num, max_steps: Int,
        tol: Num, nose_tol: Num, with_smin: Bool) -> Cpf {
    base = runpf(c, NR, false, tol, 10)
    if not base.success {
        stop("cpf: the base case does not converge")
    }
    net = internal(Case { base_mva: c.base_mva, bus: base.bus, gen: base.gen, branch: c.branch })
    nb = net.nb
    y = make_ybus(nb, net.fbus, net.tbus, net.bus, net.branch, net.base, full_network())
    ref, pv, pq = bus_types(net)
    npv: Int = len(pv)
    npq: Int = len(pq)
    n: Int = npv + 2 * npq
    s0 = make_sbus(net)
    dbus = direction(net.gen, net.gbus, nb, poi, sink)
    #: dF/dlambda = -d on the P rows; the reference's share is not an equation
    dsel = zeros(n)
    for k to npv {
        dsel[k] = -dbus[pv[k]]
    }
    for k to npq {
        dsel[npv + k] = -dbus[pq[k]]
    }
    v = initial_v(net)
    x = zeros(n + 1)
    for k to npv {
        x[k] = atan2(im(v[pv[k]]), re(v[pv[k]]))
    }
    for k to npq {
        x[npv + k] = atan2(im(v[pq[k]]), re(v[pq[k]]))
        x[npv + npq + k] = abs(v[pq[k]])
    }
    vm0 = zeros(nb)
    va0 = zeros(nb)
    for b to nb {
        vm0[b] = abs(v[b])
        va0[b] = atan2(im(v[b]), re(v[b]))
    }

    lams = [0.0]
    vps = [vm0[poi]]
    vmins = [amin(vm0)]
    vms = []
    vms = append(vms, vm0)
    vas = []
    vas = append(vas, va0)
    smins = [_smin(y, v, pv, pq, with_smin)]
    z = zeros(n + 1)
    z[n] = 1.0
    #: every accepted point's unknowns and tangent, to step back to when the nose is passed
    xs = []
    xs = append(xs, x)
    zs = []
    zs = append(zs, z)
    sigma = sigma0
    lam_max = 0.0
    nose: Int = 0
    q: [Int] = []
    steps: Int = 0
    while steps < max_steps {
        steps += 1
        #: the tangent: [J, -d; z_prev'] z = [0; 1]
        cp, ri, vals = jacobian(y, v, pv, pq)
        acp, ari, av = augment(n, cp, ri, vals, dsel, z)
        if len(q) == 0 {
            q = mindeg_order(n + 1, acp, ari)
        }
        rhs = zeros(n + 1)
        rhs[n] = 1.0
        tz = lu_solve(splu(n + 1, acp, ari, av, q, PIVOT), rhs)
        tz = tz / sqrt(sum(tz * tz))
        #: predict, then correct on the augmented system
        ok = false
        its: Int = 0
        xn = copy(x)
        vn = copy(v)
        while not ok and sigma > SIGMA_MIN {
            xn = x + sigma * tz
            vn = _volts(v, xn, pv, pq, npv, npq)
            ok, its, xn, vn = _correct(y, s0, dbus, dsel, xn, vn, x, tz, sigma, pv, pq, n, q, tol)
            if not ok {
                sigma = sigma / 2.0
            }
        }
        if not ok {
            break
        }
        x = xn
        v = vn
        z = tz
        lam = x[n]
        vmk = zeros(nb)
        vak = zeros(nb)
        for b to nb {
            vmk[b] = abs(v[b])
            vak[b] = atan2(im(v[b]), re(v[b]))
        }
        k: Int = len(lams)
        if k >= 2 and lam < lams[k - 1] and lams[k - 1] >= lams[k - 2] and sigma > nose_tol {
            #: the nose lies between points k - 2 and k: go back to k - 2 and come up on it with a quarter step
            back: Int = k - 2
            lams = lams[:back + 1]
            vps = vps[:back + 1]
            vmins = vmins[:back + 1]
            vms = vms[:back + 1]
            vas = vas[:back + 1]
            smins = smins[:back + 1]
            xs = xs[:back + 1]
            zs = zs[:back + 1]
            x = xs[back, :]
            z = zs[back, :]
            for b to nb {
                v[b] = vms[back, b] * (cos(vas[back, b]) + sin(vas[back, b]) * 1.0i)
            }
            sigma = sigma / NOSE_SHRINK
            lam_max = amax(lams)
            continue
        }
        lams = append(lams, lam)
        vps = append(vps, vmk[poi])
        vmins = append(vmins, amin(vmk))
        vms = append(vms, vmk)
        vas = append(vas, vak)
        smins = append(smins, _smin(y, v, pv, pq, with_smin))
        xs = append(xs, x)
        zs = append(zs, z)
        if lam > lam_max {
            lam_max = lam
            nose = len(lams) - 1
        }
        if lam < stop_frac * lam_max {
            break
        }
        if its <= EASY {
            sigma = min(sigma * GROW, sigma_max)
        }
    }
    return Cpf { lam: lams, v_poi: vps, v_min: vmins, vm: vms, va: vas, lam_max: lam_max, nose: nose,
                 smin: smins }
}

def _volts(v: CArr, x: Arr, pv: [Int], pq: [Int], npv: Int, npq: Int) -> CArr {
    out = copy(v)
    for k to npv {
        b = pv[k]
        out[b] = abs(v[b]) * (cos(x[k]) + sin(x[k]) * 1.0i)
    }
    for k to npq {
        b = pq[k]
        out[b] = x[npv + npq + k] * (cos(x[npv + k]) + sin(x[npv + k]) * 1.0i)
    }
    return out
}

#| Newton on [F(x, lambda); z'(x - x_k) - sigma] from the prediction.
def _correct(y: Ybus, s0: CArr, dbus: Arr, dsel: Arr, x0: Arr, v0: CArr, xk: Arr, z: Arr, sigma: Num,
             pv: [Int], pq: [Int], n: Int, q: [Int], tol: Num) {
    x = copy(x0)
    v = copy(v0)
    npv: Int = len(pv)
    npq: Int = len(pq)
    its: Int = 0
    max_it: Int = 10
    while its < max_it {
        sb = czeros(len(s0))
        for b to len(s0) {
            sb[b] = s0[b] + x[n] * dbus[b]
        }
        f = mismatch(y, v, sb, pv, pq)
        g = zeros(n + 1)
        for k to n {
            g[k] = f[k]
        }
        g[n] = sum(z * (x - xk)) - sigma
        if amax(abs(g)) < tol {
            return true, its, x, v
        }
        its += 1
        cp, ri, vals = jacobian(y, v, pv, pq)
        acp, ari, av = augment(n, cp, ri, vals, dsel, z)
        dx = lu_solve(splu(n + 1, acp, ari, av, q, PIVOT), g)
        x = x - dx
        v = _volts(v, x, pv, pq, npv, npq)
    }
    return false, its, x, v
}

def _smin(y: Ybus, v: CArr, pv: [Int], pq: [Int], on: Bool) -> Num {
    if not on {
        return 0.0
    }
    n: Int = len(pv) + 2 * len(pq)
    cp, ri, vals = jacobian(y, v, pv, pq)
    _, s, _ = svd(_dense(n, cp, ri, vals))
    return amin(s)
}

test "case9: the PV curve for 1 pu more at bus 5 has a nose, where the Jacobian is near singular" {
    c = read_case("lib/power/cases/case9.m")
    r = cpf(c, 4, SINK_SLACK, 0.5, 1.0, 0.7, 200, 1e-10, 1e-4, true)
    assert_true(r.lam_max > 1.0)
    assert_true(r.nose > 0 and r.nose < len(r.lam) - 1)
    #: the curve passes the nose: the last point is below it, at a lower voltage
    assert_true(r.lam[len(r.lam) - 1] < r.lam_max)
    assert_true(r.v_poi[len(r.lam) - 1] < r.v_poi[r.nose])
    assert_true(r.smin[r.nose] < 0.05 * r.smin[0])
}

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

from power.case import Case, GEN_STATUS, PG, VM, read_case
from power.pf import NR, bus_types, initial_v, internal, jacobian, make_sbus, mismatch, runpf
from power.splu import lu_solve, mindeg_order, splu
from power.ybus import Ybus, full_network, make_ybus
from power.study import SINK_GENS, SINK_SLACK

PIVOT: Num = 0.001
#: the step grows by this factor after a correction of at most EASY iterations
GROW: Num = 1.5
EASY: Int = 3
#: the smallest step tried before the trace gives up
SIGMA_MIN: Num = 1e-6
#: past the nose, the trace steps back and shrinks the step by this factor until it is below nose_tol
NOSE_SHRINK: Num = 4.0

record Cpf {
    #: per point: lambda (per unit), Vm at the POI, the smallest Vm, corrector iterations
    lam: Arr
    v_poi: Arr
    v_min: Arr
    #: the voltages at every point (points x buses)
    vm: Arr
    va: Arr
    lam_max: Num
    nose: Int
    #: the smallest singular value of the power-flow Jacobian at each point
    smin: Arr
}

#| The transfer direction: +1 at bus row poi, made up by the sink.
def direction(net_gen: Arr, gbus: [Int], nb: Int, poi: Int, sink: Int) -> Arr {
    d = zeros(nb)
    d[poi] = 1.0
    if sink == SINK_GENS {
        total = 0.0
        for g to len(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(gbus) {
            if net_gen[g, GEN_STATUS] > 0.0 and net_gen[g, PG] > 0.0 {
                d[gbus[g]] -= net_gen[g, PG] / total
            }
        }
    }
    return d
}

#| The augmented matrix [[J, col]; [row']] in compressed sparse columns.
def augment(n: Int, cp: [Int], ri: [Int], vals: Arr, col: Arr, row: Arr) {
    m: Int = n + 1
    acp = int_zeros(m + 1)
    ari: [Int] = []
    av = []
    for j to n {
        acp[j] = len(ari)
        for p from cp[j] to cp[j + 1] {
            ari = append(ari, ri[p])
            av = append(av, vals[p])
        }
        ari = append(ari, n)
        av = append(av, row[j])
    }
    acp[n] = len(ari)
    for i to n {
        if col[i] != 0.0 {
            ari = append(ari, i)
            av = append(av, col[i])
        }
    }
    ari = append(ari, n)
    av = append(av, row[n])
    acp[m] = len(ari)
    return acp, ari, av
}

def _dense(n: Int, cp: [Int], ri: [Int], vals: Arr) -> Arr {
    a = zeros(n, n)
    for j to n {
        for p from cp[j] to cp[j + 1] {
            a[ri[p], j] += vals[p]
        }
    }
    return a
}

#| The PV curve for a transfer into bus row poi, from the case's power
#| flow solution, with the first step sigma0 (per unit of the arc length).
def cpf(c: Case, poi: Int, sink: Int, sigma0: Num, sigma_max: Num, stop_frac: Num, max_steps: Int,
        tol: Num, nose_tol: Num, with_smin: Bool) -> Cpf {
    base = runpf(c, NR, false, tol, 10)
    if not base.success {
        stop("cpf: the base case does not converge")
    }
    net = internal(Case { base_mva: c.base_mva, bus: base.bus, gen: base.gen, branch: c.branch })
    nb = net.nb
    y = make_ybus(nb, net.fbus, net.tbus, net.bus, net.branch, net.base, full_network())
    ref, pv, pq = bus_types(net)
    npv: Int = len(pv)
    npq: Int = len(pq)
    n: Int = npv + 2 * npq
    s0 = make_sbus(net)
    dbus = direction(net.gen, net.gbus, nb, poi, sink)
    #: dF/dlambda = -d on the P rows; the reference's share is not an equation
    dsel = zeros(n)
    for k to npv {
        dsel[k] = -dbus[pv[k]]
    }
    for k to npq {
        dsel[npv + k] = -dbus[pq[k]]
    }
    v = initial_v(net)
    x = zeros(n + 1)
    for k to npv {
        x[k] = atan2(im(v[pv[k]]), re(v[pv[k]]))
    }
    for k to npq {
        x[npv + k] = atan2(im(v[pq[k]]), re(v[pq[k]]))
        x[npv + npq + k] = abs(v[pq[k]])
    }
    vm0 = zeros(nb)
    va0 = zeros(nb)
    for b to nb {
        vm0[b] = abs(v[b])
        va0[b] = atan2(im(v[b]), re(v[b]))
    }

lams = [0.0]
    vps = [vm0[poi]]
    vmins = [amin(vm0)]
    vms = []
    vms = append(vms, vm0)
    vas = []
    vas = append(vas, va0)
    smins = [_smin(y, v, pv, pq, with_smin)]
    z = zeros(n + 1)
    z[n] = 1.0
    #: every accepted point's unknowns and tangent, to step back to when the nose is passed
    xs = []
    xs = append(xs, x)
    zs = []
    zs = append(zs, z)
    sigma = sigma0
    lam_max = 0.0
    nose: Int = 0
    q: [Int] = []
    steps: Int = 0
    while steps < max_steps {
        steps += 1
        #: the tangent: [J, -d; z_prev'] z = [0; 1]
        cp, ri, vals = jacobian(y, v, pv, pq)
        acp, ari, av = augment(n, cp, ri, vals, dsel, z)
        if len(q) == 0 {
            q = mindeg_order(n + 1, acp, ari)
        }
        rhs = zeros(n + 1)
        rhs[n] = 1.0
        tz = lu_solve(splu(n + 1, acp, ari, av, q, PIVOT), rhs)
        tz = tz / sqrt(sum(tz * tz))
        #: predict, then correct on the augmented system
        ok = false
        its: Int = 0
        xn = copy(x)
        vn = copy(v)
        while not ok and sigma > SIGMA_MIN {
            xn = x + sigma * tz
            vn = _volts(v, xn, pv, pq, npv, npq)
            ok, its, xn, vn = _correct(y, s0, dbus, dsel, xn, vn, x, tz, sigma, pv, pq, n, q, tol)
            if not ok {
                sigma = sigma / 2.0
            }
        }
        if not ok {
            break
        }
        x = xn
        v = vn
        z = tz
        lam = x[n]
        vmk = zeros(nb)
        vak = zeros(nb)
        for b to nb {
            vmk[b] = abs(v[b])
            vak[b] = atan2(im(v[b]), re(v[b]))
        }
        k: Int = len(lams)
        if k >= 2 and lam < lams[k - 1] and lams[k - 1] >= lams[k - 2] and sigma > nose_tol {
            #: the nose lies between points k - 2 and k: go back to k - 2 and come up on it with a quarter step
            back: Int = k - 2
            lams = lams[:back + 1]
            vps = vps[:back + 1]
            vmins = vmins[:back + 1]
            vms = vms[:back + 1]
            vas = vas[:back + 1]
            smins = smins[:back + 1]
            xs = xs[:back + 1]
            zs = zs[:back + 1]
            x = xs[back, :]
            z = zs[back, :]
            for b to nb {
                v[b] = vms[back, b] * (cos(vas[back, b]) + sin(vas[back, b]) * 1.0i)
            }
            sigma = sigma / NOSE_SHRINK
            lam_max = amax(lams)
            continue
        }
        lams = append(lams, lam)
        vps = append(vps, vmk[poi])
        vmins = append(vmins, amin(vmk))
        vms = append(vms, vmk)
        vas = append(vas, vak)
        smins = append(smins, _smin(y, v, pv, pq, with_smin))
        xs = append(xs, x)
        zs = append(zs, z)
        if lam > lam_max {
            lam_max = lam
            nose = len(lams) - 1
        }
        if lam < stop_frac * lam_max {
            break
        }
        if its <= EASY {
            sigma = min(sigma * GROW, sigma_max)
        }
    }
    return Cpf { lam: lams, v_poi: vps, v_min: vmins, vm: vms, va: vas, lam_max: lam_max, nose: nose,
                 smin: smins }
}

def _volts(v: CArr, x: Arr, pv: [Int], pq: [Int], npv: Int, npq: Int) -> CArr {
    out = copy(v)
    for k to npv {
        b = pv[k]
        out[b] = abs(v[b]) * (cos(x[k]) + sin(x[k]) * 1.0i)
    }
    for k to npq {
        b = pq[k]
        out[b] = x[npv + npq + k] * (cos(x[npv + k]) + sin(x[npv + k]) * 1.0i)
    }
    return out
}

#| Newton on [F(x, lambda); z'(x - x_k) - sigma] from the prediction.
def _correct(y: Ybus, s0: CArr, dbus: Arr, dsel: Arr, x0: Arr, v0: CArr, xk: Arr, z: Arr, sigma: Num,
             pv: [Int], pq: [Int], n: Int, q: [Int], tol: Num) {
    x = copy(x0)
    v = copy(v0)
    npv: Int = len(pv)
    npq: Int = len(pq)
    its: Int = 0
    max_it: Int = 10
    while its < max_it {
        sb = czeros(len(s0))
        for b to len(s0) {
            sb[b] = s0[b] + x[n] * dbus[b]
        }
        f = mismatch(y, v, sb, pv, pq)
        g = zeros(n + 1)
        for k to n {
            g[k] = f[k]
        }
        g[n] = sum(z * (x - xk)) - sigma
        if amax(abs(g)) < tol {
            return true, its, x, v
        }
        its += 1
        cp, ri, vals = jacobian(y, v, pv, pq)
        acp, ari, av = augment(n, cp, ri, vals, dsel, z)
        dx = lu_solve(splu(n + 1, acp, ari, av, q, PIVOT), g)
        x = x - dx
        v = _volts(v, x, pv, pq, npv, npq)
    }
    return false, its, x, v
}

def _smin(y: Ybus, v: CArr, pv: [Int], pq: [Int], on: Bool) -> Num {
    if not on {
        return 0.0
    }
    n: Int = len(pv) + 2 * len(pq)
    cp, ri, vals = jacobian(y, v, pv, pq)
    _, s, _ = svd(_dense(n, cp, ri, vals))
    return amin(s)
}

test "case9: the PV curve for 1 pu more at bus 5 has a nose, where the Jacobian is near singular" {
    c = read_case("lib/power/cases/case9.m")
    r = cpf(c, 4, SINK_SLACK, 0.5, 1.0, 0.7, 200, 1e-10, 1e-4, true)
    assert_true(r.lam_max > 1.0)
    assert_true(r.nose > 0 and r.nose < len(r.lam) - 1)
    #: the curve passes the nose: the last point is below it, at a lower voltage
    assert_true(r.lam[len(r.lam) - 1] < r.lam_max)
    assert_true(r.v_poi[len(r.lam) - 1] < r.v_poi[r.nose])
    assert_true(r.smin[r.nose] < 0.05 * r.smin[0])
}

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