CadCraft™ v0.8.1
Use it now
Menu

Resources Numerical accuracy Calculation policy

Load flow and short-circuit currents

Bus voltages, branch flows and losses of an AC network, and fault currents by symmetrical components.

How the precision is kept

Per-unit arithmetic in 64-bit. The load flow is a Newton iteration stopped when the largest power mismatch is below 1e-9 per unit. Fault currents are direct complex arithmetic on inverted bus matrices; the peak factor is the IEC 60909 formula as the file gives it.

Open item. Clamp the power factor to [0, 1] before acos in voltageDrop, and return the exact drop its comment names.

Calculations under this policy

12 calculations in 3 source files. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.

power.js Shared shell, every app 64-bit 7 calculations Reviewed

Tests bin/test-engines.mjs, tools/network-test/test.mjs

Bus admittance matrix

64-bit

The admittances that link every pair of buses.

y = 1 ÷ (r + jx); Yff = (y + jb ÷ 2) ÷ |a|²; Yft = −y ÷ conj(a); Ytt = y + jb ÷ 2

The series admittance is r ÷ (r² + x²) and −x ÷ (r² + x²), computed once per branch. Tap phase is converted from degrees with π ÷ 180. Entries are summed into dense arrays in branch order.

  • power.js
  • lines 18 to 39
/** Bus admittance matrix as real parts G, B (n×n) for the given system; also the per-branch four-element admittances. */
export function buildYbus(sys) {
  const n = sys.buses.length, idx = new Map(sys.buses.map((b, i) => [b.id, i]));
  const G = Array.from({ length: n }, () => new Array(n).fill(0)), B = Array.from({ length: n }, () => new Array(n).fill(0));
  const br = [];
  sys.buses.forEach((b, i) => { G[i][i] += b.Gs || 0; B[i][i] += b.Bs || 0; });
  for (const l of sys.branches) {
    const f = idx.get(l.from), t = idx.get(l.to);
    if (f == null || t == null) throw new Error(`branch ${l.id}: unknown bus`);
    const den = l.r * l.r + l.x * l.x; if (!(den > 0)) throw new Error(`branch ${l.id}: zero impedance`);
    const y = [l.r / den, -l.x / den], bc = (l.b || 0) / 2, tap = l.tap || 1, ph = (l.shift || 0) * DEG;
    const a = [tap * Math.cos(ph), tap * Math.sin(ph)], a2 = tap * tap;
    const yff = [y[0] / a2, (y[1] + bc) / a2], ytt = [y[0], y[1] + bc];
    const yA = [(y[0] * a[0] - y[1] * a[1]) / a2, (y[0] * a[1] + y[1] * a[0]) / a2];   // y·a/|a|² = y / conj(a)
    const yB = [(y[0] * a[0] + y[1] * a[1]) / a2, (y[1] * a[0] - y[0] * a[1]) / a2];   // y·conj(a)/|a|² = y / a
    const Yft = [-yA[0], -yA[1]], Ytf = [-yB[0], -yB[1]];
    G[f][f] += yff[0]; B[f][f] += yff[1]; G[t][t] += ytt[0]; B[t][t] += ytt[1];
    G[f][t] += Yft[0]; B[f][t] += Yft[1]; G[t][f] += Ytf[0]; B[t][f] += Ytf[1];
    br.push({ id: l.id, f, t, Yff: yff, Yft, Ytf, Ytt: ytt, rating: l.rating });
  }
  return { G, B, branches: br, idx };
}

Power injected at each bus

64-bit

Real and reactive power from the voltages.

Pᵢ = Σ Vᵢ Vⱼ (G cos δ + B sin δ); Qᵢ = Σ Vᵢ Vⱼ (G sin δ − B cos δ)

A plain sum over the buses in 64-bit. The mismatch is the difference between this sum and the scheduled power, so near the solution it is a small difference of larger terms, with rounding near 1e-16 times the size of those terms.

  • power.js
  • lines 41 to 49
function inject(G, B, V, th) {
  const n = V.length, P = new Array(n).fill(0), Q = new Array(n).fill(0);
  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) {
    if (!G[i][j] && !B[i][j]) continue;
    const d = th[i] - th[j], c = Math.cos(d), s = Math.sin(d);
    P[i] += V[i] * V[j] * (G[i][j] * c + B[i][j] * s); Q[i] += V[i] * V[j] * (G[i][j] * s - B[i][j] * c);
  }
  return { P, Q };
}

Load flow iteration

64-bit

Corrects angles and voltages until the power balances.

J Δx = (P_sched − P, Q_sched − Q); stop when max |mismatch| < 1e-9 pu

The Jacobian is the analytic one, rebuilt each iteration and solved by dense elimination. At most 30 iterations per round. A voltage that leaves 0.05 to 5 pu ends the run as a collapse.

  • power.js
  • lines 64 to 92
for (let it = 0; it <= maxIter; it++) {
  const { P, Q } = inject(G, B, V, th);
  const f = [...ang.map((i) => Psp[i] - P[i]), ...mag.map((i) => Qsp[i] - Q[i])];
  mismatch = f.length ? Math.max(...f.map(Math.abs)) : 0;
  if (mismatch < tol) { ok = true; iterations += it; break; }
  if (it === maxIter) { iterations += it; break; }
  const m = f.length, J = Array.from({ length: m }, () => new Array(m).fill(0));
  const row = (i, isQ) => (isQ ? ang.length + mag.indexOf(i) : ang.indexOf(i));
  for (const i of ang) for (const j of ang) {
    const d = th[i] - th[j], c = Math.cos(d), s = Math.sin(d);
    J[row(i, false)][ang.indexOf(j)] = i === j ? -Q[i] - B[i][i] * V[i] * V[i] : V[i] * V[j] * (G[i][j] * s - B[i][j] * c);
  }
  for (const i of ang) for (const j of mag) {
    const d = th[i] - th[j], c = Math.cos(d), s = Math.sin(d);
    J[row(i, false)][ang.length + mag.indexOf(j)] = i === j ? P[i] / V[i] + G[i][i] * V[i] : V[i] * (G[i][j] * c + B[i][j] * s);
  }
  for (const i of mag) for (const j of ang) {
    const d = th[i] - th[j], c = Math.cos(d), s = Math.sin(d);
    J[row(i, true)][ang.indexOf(j)] = i === j ? P[i] - G[i][i] * V[i] * V[i] : -V[i] * V[j] * (G[i][j] * c + B[i][j] * s);
  }
  for (const i of mag) for (const j of mag) {
    const d = th[i] - th[j], c = Math.cos(d), s = Math.sin(d);
    J[row(i, true)][ang.length + mag.indexOf(j)] = i === j ? Q[i] / V[i] - B[i][i] * V[i] : V[i] * (G[i][j] * s - B[i][j] * c);
  }
  const dx = solveLinear(J, f);
  if (!dx) { problems.push(problem('loadflow.singular', 'error', 'The load-flow Jacobian is singular: the network is islanded or has a zero-impedance loop.', [{ kind: 'system', id: 'all' }])); return { ok: false, problems, buses: {}, branches: {}, iterations, mismatch }; }
  ang.forEach((i, k) => { th[i] += dx[k]; }); mag.forEach((i, k) => { V[i] += dx[ang.length + k]; });
  if (V.some((v) => !(v > 0.05 && v < 5))) { problems.push(problem('loadflow.collapse', 'error', 'Voltage collapse: the load is beyond what the network can deliver (no real solution).', [{ kind: 'system', id: 'all' }])); return { ok: false, problems, buses: {}, branches: {}, iterations: iterations + it, mismatch }; }
}

Reactive limits of generator buses

64-bit

Switches a voltage-controlled bus to a fixed reactive power when it passes its limit.

switch when Q_g > Qmax + 1e-9 or Q_g < Qmin − 1e-9

The comparison has an absolute allowance of 1e-9 pu. A bus that was switched is not switched back in a later round.

  • power.js
  • lines 94 to 102
// Q limits on PV buses
const { Q } = inject(G, B, V, th); let changed = false;
sys.buses.forEach((b, i) => {
  if (type[i] !== 'pv' || opts.noQLimits) return;
  const Qg = Q[i] + (b.Qd || 0);
  if (b.Qmax != null && Qg > b.Qmax + 1e-9) { type[i] = 'pq'; Qsp[i] = b.Qmax - (b.Qd || 0); changed = true; switched.push({ bus: b.id, limit: 'Qmax' }); }
  else if (b.Qmin != null && Qg < b.Qmin - 1e-9) { type[i] = 'pq'; Qsp[i] = b.Qmin - (b.Qd || 0); changed = true; switched.push({ bus: b.id, limit: 'Qmin' }); }
});
if (!changed) break;

Branch flows and losses

64-bit

Power at both ends of each branch and the loss in it.

S_f = V_f conj(I_f); loss = P_f + P_t

The loss of a branch is the sum of the powers entering at its two ends, which are nearly equal and opposite. For a low-loss branch that sum loses digits in the ratio of the flow to the loss.

  • power.js
  • lines 109 to 115
for (const k of branches) {
  const cm = (a, b) => [a[0] * b[0] - a[1] * b[1], a[0] * b[1] + a[1] * b[0]];
  const If = [cm(k.Yff, Vc[k.f])[0] + cm(k.Yft, Vc[k.t])[0], cm(k.Yff, Vc[k.f])[1] + cm(k.Yft, Vc[k.t])[1]];
  const It = [cm(k.Ytf, Vc[k.f])[0] + cm(k.Ytt, Vc[k.t])[0], cm(k.Ytf, Vc[k.f])[1] + cm(k.Ytt, Vc[k.t])[1]];
  const Sf = cm(Vc[k.f], [If[0], -If[1]]), St = cm(Vc[k.t], [It[0], -It[1]]);
  outL[k.id] = { Pf: Sf[0], Qf: Sf[1], Pt: St[0], Qt: St[1], Sf: cabs(Sf) * base, St: cabs(St) * base, loss: Sf[0] + St[0], If: cabs(If) };
  lossP += Sf[0] + St[0]; lossQ += Sf[1] + St[1];

Per-unit base

64-bit

Ohms to per unit on a voltage and power base.

Z_base = kV² ÷ MVA; z_pu = z ÷ Z_base

One multiplication and two divisions in 64-bit, with no rounding step and no stored constant.

  • power.js
  • lines 125 to 126
export const zBase = (kV, MVA) => (kV * kV) / MVA;
export const ohmToPu = (z, kV, MVA) => z / zBase(kV, MVA);

Feeder voltage drop

64-bit

The three-phase voltage drop of a feeder.

ΔV ≈ √3 I (R cos φ + X sin φ), φ = acos(pf)

The power factor goes to acos without a clamp, so a value above 1 by rounding gives NaN. Only the usual approximation is returned, although the comment above the function also names the exact drop. The percent is taken against kV × 1000.

  • power.js
  • lines 127 to 130
export function voltageDrop({ R, X, I, pf, kV }) {
  const phi = Math.acos(pf), approx = Math.sqrt(3) * I * (R * pf + X * Math.sin(phi));
  return { approx, percent: kV ? 100 * approx / (kV * 1000) : null };
}
shortcircuit.js Shared shell, every app 64-bit 3 calculations Reviewed

Tests bin/test-engines.mjs, tools/network-test/test.mjs

Sequence impedance matrices

64-bit

The impedance seen from every bus in the positive, negative and zero sequence.

Z = Y⁻¹ for each sequence

Each admittance is 1 ÷ z by complex division. The three matrices are inverted in full by Gauss–Jordan elimination; the rounding error of an explicit inverse grows with the condition number of Y.

  • shortcircuit.js
  • lines 22 to 41
/** The three Zbus matrices. A zero-sequence matrix of a system with no ground path is flagged (`noGround0`). */
export function buildZbus(sys) {
  const n = sys.buses.length, idx = new Map(sys.buses.map((b, i) => [b.id, i]));
  const mk = () => Array.from({ length: n }, () => Array.from({ length: n }, () => [0, 0]));
  const Y = [mk(), mk(), mk()], hasGround = [false, false, false];
  const add = (s, i, j, y) => { Y[s][i][i] = cadd(Y[s][i][i], y); if (j != null) { Y[s][j][j] = cadd(Y[s][j][j], y); Y[s][i][j] = csub(Y[s][i][j], y); Y[s][j][i] = csub(Y[s][j][i], y); } };
  const get = (z, dflt) => (z ? z : dflt);
  for (const l of sys.branches) {
    const f = idx.get(l.from), t = idx.get(l.to); if (f == null || t == null) throw new Error(`branch ${l.id}: unknown bus`);
    add(0, f, t, inv(l.z1)); add(1, f, t, inv(get(l.z2, l.z1))); if (l.z0) add(2, f, t, inv(l.z0));
  }
  for (const s of sys.sources || []) {
    const i = idx.get(s.bus); if (i == null) throw new Error(`source: unknown bus ${s.bus}`);
    add(0, i, null, inv(s.z1)); add(1, i, null, inv(get(s.z2, s.z1))); hasGround[0] = hasGround[1] = true;
    if (s.z0) { add(2, i, null, inv(s.z0)); hasGround[2] = true; }
  }
  sys.buses.forEach((b, i) => { if (b.g0) { add(2, i, null, inv(b.g0)); hasGround[2] = true; } });
  const Z = Y.map((m, s) => (s === 2 && !hasGround[2] ? null : cinvert(m)));
  return { Z1: Z[0], Z2: Z[1], Z0: Z[2], noGround0: !hasGround[2], idx, n };
}

Sequence to phase values

64-bit

Phase currents and voltages from the sequence components.

I_a = I₀ + I₁ + I₂; I_b = I₀ + a² I₁ + a I₂; I_c = I₀ + a I₁ + a² I₂, a = −½ + j√3 ÷ 2

a is held as −0.5 and Math.sqrt(3) ÷ 2. For a balanced three-phase fault the three magnitudes agree only to rounding.

  • shortcircuit.js
  • line 43
const toPhase = (i0, i1, i2) => [cadd(cadd(i0, i1), i2), cadd(cadd(i0, cmul(A2, i1)), cmul(A, i2)), cadd(cadd(i0, cmul(A, i1)), cmul(A2, i2))];

Fault current

64-bit

Initial and peak fault current, fault power, and the voltages and branch currents during the fault.

3-phase: I₁ = c ÷ (Z₁ + Z_f); line to ground: I₁ = c ÷ (Z₁ + Z₂ + Z₀ + 3Z_f); i_p = κ √2 I_k, κ = 1.02 + 0.98 e^(−3R/X)

Direct complex arithmetic on the diagonal entries of the impedance matrices. κ uses the coefficients the file attributes to IEC 60909-0; R/X is taken from the positive-sequence impedance at the fault. Base current is baseMVA × 10⁶ ÷ (√3 × kV × 10³); without a bus kV the currents stay in per unit.

  • shortcircuit.js
  • lines 45 to 83
/**
 * One fault. type: '3ph' | 'lg' | 'll' | 'llg'; at: bus id; zf: [r, x] fault impedance (default 0); c: pre-fault voltage factor (default 1).
 * → { ok, I (A per phase [a,b,c] magnitudes) , pu: { Ia, Ib, Ic, I0, I1, I2 } complex, Ik (A, the largest phase), ip (A), MVA, xr, buses: { id: { Va, Vb, Vc (pu mags) } },
 *     branches: { id: { Ia, Ib, Ic (A mags) } }, problems }
 */
export function fault(sys, zb, { type = '3ph', at, zf = [0, 0], c = 1 }) {
  const k = zb.idx.get(at), problems = [];
  if (k == null) return { ok: false, problems: [problem('fault.bus', 'error', `Fault bus ${at} is not in the system.`, [{ kind: 'bus', id: String(at) }])] };
  const Z1 = zb.Z1[k][k], Z2 = zb.Z2[k][k], Z0 = zb.Z0 ? zb.Z0[k][k] : null, V = [c, 0];
  if ((type === 'lg' || type === 'llg') && (!Z0 || cabs(Z0) > 1e6)) return { ok: false, problems: [problem('fault.noground', 'error', `Bus ${at} has no zero-sequence path to ground, so a fault to ground there has no current in this model.`, [{ kind: 'bus', id: String(at) }])] };
  const Zf3 = cscale(zf, 3);
  let i0 = ZERO, i1, i2 = ZERO;
  if (type === '3ph') i1 = cdiv(V, cadd(Z1, zf));
  else if (type === 'lg') { i1 = cdiv(V, cadd(cadd(cadd(Z1, Z2), Z0), Zf3)); i2 = i1; i0 = i1; }
  else if (type === 'll') { i1 = cdiv(V, cadd(cadd(Z1, Z2), zf)); i2 = cscale(i1, -1); }
  else if (type === 'llg') {
    const Z0e = cadd(Z0, Zf3), par = cdiv(cmul(Z2, Z0e), cadd(Z2, Z0e));
    i1 = cdiv(V, cadd(Z1, par)); i2 = cscale(cmul(i1, cdiv(Z0e, cadd(Z2, Z0e))), -1); i0 = cscale(cmul(i1, cdiv(Z2, cadd(Z2, Z0e))), -1);
  } else return { ok: false, problems: [problem('fault.type', 'error', `Unknown fault type ${type}.`, [])] };
  const [Ia, Ib, Ic] = toPhase(i0, i1, i2), Ibase = (sys.baseMVA * 1e6) / (Math.sqrt(3) * (sys.buses[k].kV ? sys.buses[k].kV * 1e3 : NaN));
  const mags = [cabs(Ia), cabs(Ib), cabs(Ic)], pu = Math.max(...mags);
  const r = Z1[0], x = Z1[1], xr = r > 1e-12 ? x / r : Infinity, kappa = Number.isFinite(xr) ? 1.02 + 0.98 * Math.exp(-3 / xr) : 1.02;
  const buses = {}, branches = {};
  sys.buses.forEach((b, j) => {
    const v1 = csub(V, cmul(zb.Z1[j][k], i1)), v2 = cscale(cmul(zb.Z2[j][k], i2), -1), v0 = zb.Z0 ? cscale(cmul(zb.Z0[j][k], i0), -1) : ZERO;
    const [va, vb, vc] = toPhase(v0, v1, v2); buses[b.id] = { Va: cabs(va), Vb: cabs(vb), Vc: cabs(vc), seq: [v0, v1, v2] };
  });
  const vs = (j) => buses[sys.buses[j].id].seq;
  for (const l of sys.branches) {
    const f = zb.idx.get(l.from), t = zb.idx.get(l.to), a = vs(f), b = vs(t);
    const s1 = cdiv(csub(a[1], b[1]), l.z1), s2 = cdiv(csub(a[2], b[2]), l.z2 || l.z1), s0 = l.z0 ? cdiv(csub(a[0], b[0]), l.z0) : ZERO;
    const [ia, ib, ic] = toPhase(s0, s1, s2), sc = Number.isFinite(Ibase) ? Ibase : 1;
    branches[l.id] = { Ia: cabs(ia) * sc, Ib: cabs(ib) * sc, Ic: cabs(ic) * sc };
  }
  const unit = Number.isFinite(Ibase) ? Ibase : null;
  const out = { ok: true, type, at, pu: { Ia, Ib, Ic, I0: i0, I1: i1, I2: i2 }, Ik_pu: pu, Ik: unit ? pu * unit : null, ip: unit ? kappa * Math.SQRT2 * pu * unit : null, kappa, xr, MVA: sys.buses[k].kV ? Math.sqrt(3) * sys.buses[k].kV * pu * unit / 1e6 : null, buses, branches, problems };
  out.angle = carg(Ia) * 180 / Math.PI;
  return out;
}
numeric.js Shared shell, every app 64-bit 2 calculations Reviewed

Tests bin/test-engines.mjs, tools/network-test/test.mjs

Complex division

64-bit

a ÷ b for the impedances and currents of the fault study.

a ÷ b = a conj(b) ÷ (b_re² + b_im²)

The textbook form. The squared magnitude can overflow or underflow for values near the ends of the 64-bit range.

  • numeric.js
  • line 27
export const cdiv = (a, b) => { const d = b[0] * b[0] + b[1] * b[1]; return [(a[0] * b[0] + a[1] * b[1]) / d, (a[1] * b[0] - a[0] * b[1]) / d]; };

Complex matrix inverse

64-bit

Inverts a sequence admittance matrix.

Gauss–Jordan elimination with pivoting by magnitude

Singular when the pivot is at most 1e-14 of the largest entry, so the test scales with the matrix.

  • numeric.js
  • lines 34 to 47
/** Invert a complex n×n matrix (rows of [re, im] pairs) by Gauss–Jordan; returns null when singular. */
export function cinvert(A) {
  const n = A.length, M = A.map((r, i) => [...r.map((v) => [v[0], v[1]]), ...Array.from({ length: n }, (_, j) => (i === j ? [1, 0] : [0, 0]))]);
  let scale = 0; for (const r of A) for (const v of r) scale = Math.max(scale, cabs(v));
  const tiny = Math.max(scale, 1e-300) * 1e-14;
  for (let c = 0; c < n; c++) {
    let p = c; for (let r = c + 1; r < n; r++) if (cabs(M[r][c]) > cabs(M[p][c])) p = r;
    if (cabs(M[p][c]) <= tiny) return null;
    [M[c], M[p]] = [M[p], M[c]];
    const d = M[c][c]; for (let k = 0; k < 2 * n; k++) M[c][k] = cdiv(M[c][k], d);
    for (let r = 0; r < n; r++) if (r !== c) { const f = M[r][c]; if (f[0] || f[1]) for (let k = 0; k < 2 * n; k++) M[r][k] = csub(M[r][k], cmul(f, M[c][k])); }
  }
  return M.map((r) => r.slice(n));
}

In detail

Method
Newton–Raphson in polar form (the file cites Tinney and Hart 1967); π branch model with off-nominal tap; Zbus by Gauss–Jordan inversion of the sequence admittance matrices; sequence connections for each fault type; κ = 1.02 + 0.98 e^(−3R/X).
Number types
64-bit numbers; complex values as [re, im] pairs.
Tolerances
Mismatch < 1e-9 pu, 30 iterations, 12 rounds of reactive-limit switching with a 1e-9 allowance on the limits. A voltage outside 0.05 to 5 pu stops the run. Inversion: a pivot at most 1e-14 of the largest entry is singular. No zero-sequence path when |Z0| > 1e6 pu. X/R is infinite when R ≤ 1e-12.
Checked by
The engines test runs a small system and compares the three-phase fault current with 10 pu within 1e-5.
Limits
The fault study ignores load currents, line charging and taps, as the file states. voltageDrop returns only the approximate formula.

Approach: 64-bit, stated limit

The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.

Every policy that takes this approach