AC Branch Pi-Model + Transformer Handling

Implement the exact branch power flow equations in acopf-math-model.md using MATPOWER branch data:

[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, TAP, SHIFT, BR_STATUS, ANGMIN, ANGMAX]

Quick start

• Use scripts/branch_flows.py to compute per-unit branch flows.
• Treat the results as power leaving the “from” bus and power leaving the “to” bus (i.e., compute both directions explicitly).

Example:

import json
import numpy as np

from scripts.branch_flows import compute_branch_flows_pu, build_bus_id_to_idx

data = json.load(open("/root/network.json"))
baseMVA = float(data["baseMVA"])
buses = np.array(data["bus"], dtype=float)
branches = np.array(data["branch"], dtype=float)

bus_id_to_idx = build_bus_id_to_idx(buses)

Vm = buses[:, 7]  # initial guess VM
Va = np.deg2rad(buses[:, 8])  # initial guess VA

br = branches[0]
P_ij, Q_ij, P_ji, Q_ji = compute_branch_flows_pu(Vm, Va, br, bus_id_to_idx)

S_ij_MVA = (P_ij**2 + Q_ij**2) ** 0.5 * baseMVA
S_ji_MVA = (P_ji**2 + Q_ji**2) ** 0.5 * baseMVA
print(S_ij_MVA, S_ji_MVA)

Model details (match the task formulation)

Per-unit conventions

• Work in per-unit internally.
• Convert with baseMVA:

- \(P{pu} = P{MW} / baseMVA\) - \(Q{pu} = Q{MVAr} / baseMVA\) - \(|S|{MVA} = |S|{pu} \cdot baseMVA\)

Transformer handling (MATPOWER TAP + SHIFT)

• Use \(T_{ij} = tap \cdot e^{j \cdot shift}\).
• Implementation shortcut (real tap + phase shift):

- If abs(TAP) < 1e-12, treat tap = 1.0 (no transformer). - Convert SHIFT from degrees to radians. - Use the angle shift by modifying the angle difference: - \(\delta{ij} = \thetai - \thetaj - shift\) - \(\delta{ji} = \thetaj - \thetai + shift\)

Series admittance

Given BRR = r, BRX = x:

• If r == 0 and x == 0, set g = 0, b = 0 (avoid divide-by-zero).
• Else:

- \(y = 1/(r + jx) = g + jb\) - \(g = r/(r^2 + x^2)\) - \(b = -x/(r^2 + x^2)\)

Line charging susceptance

• BRB is the total line charging susceptance \(bc\) (per unit).
• Each end gets \(b_c/2\) in the standard pi model.

Power flow equations (use these exactly)

Let:

• \(Vi = |Vi| e^{j\thetai}\), \(Vj = |Vj| e^{j\thetaj}\)
• tap is real, shift is radians
• invt = 1/tap, invt2 = inv_t^2

Then the real/reactive power flow from i→j is:

• \(P{ij} = g |Vi|^2 inv\t2 - |Vi||Vj| inv\t (g\cos\delta{ij} + b\sin\delta{ij})\)
• \(Q{ij} = -(b + bc/2)|Vi|^2 inv\t2 - |Vi||Vj| inv\t (g\sin\delta{ij} - b\cos\delta_{ij})\)

And from j→i is:

• \(P{ji} = g |Vj|^2 - |Vi||Vj| inv\t (g\cos\delta{ji} + b\sin\delta_{ji})\)
• \(Q{ji} = -(b + bc/2)|Vj|^2 - |Vi||Vj| inv\t (g\sin\delta{ji} - b\cos\delta{ji})\)

Compute apparent power:

• \(|S{ij}| = \sqrt{P{ij}^2 + Q_{ij}^2}\)
• \(|S{ji}| = \sqrt{P{ji}^2 + Q_{ji}^2}\)

Common uses

Enforce MVA limits (rateA)

• RATE_A is an MVA limit (may be 0 meaning “no limit”).
• Enforce in both directions:

- \(|S{ij}| \le RATEA\) - \(|S{ji}| \le RATEA\)

Compute branch loading %

For reporting “most loaded branches”:

• loadingpct = 100 * max(|Sij|, |Sji|) / RATEA if RATE_A > 0, else 0.

Aggregate bus injections for nodal balance

To build the branch flow sum for each bus \(i\):

• Add \(P{ij}, Q{ij}\) to bus i
• Add \(P{ji}, Q{ji}\) to bus j

This yields arrays Pout[i], Qout[i] such that the nodal balance can be written as:

• \(P^g - P^d - G^s|V|^2 = P_{out}\)
• \(Q^g - Q^d + B^s|V|^2 = Q_{out}\)

Sanity checks (fast debug)

• With SHIFT=0 and TAP=1, if \(Vi = Vj\) and \(\thetai=\thetaj\), then \(P{ij}\approx 0\) and \(P{ji}\approx 0\) (lossless only if r=0).
• For a pure transformer (r=x=0) you should not get meaningful flows; treat as g=b=0 (no series element).