Power Flow Analysis and Voltage Stability Assessment of a Transmission Network

A dependency-light solver, cross-checked two ways

This report documents a custom, first-principles Newton-Raphson AC power flow solver, implemented in Python without commercial or third-party power system libraries, used to evaluate the steady-state performance, contingency security, and voltage stability of the WSCC 9-bus benchmark transmission network.

The primary validation method is a direct cross-check against pandapower: a dedicated script builds the identical network in pandapower and runs its own AC power flow solver for comparison. As a secondary, independently coded check, a Gauss-Seidel power flow implementation was run against the same network model, agreeing with the Newton-Raphson solution to within 0%–2% voltage magnitude across all buses.

Four analyses were carried out on the validated model:

• N-1 contingency screening across all six transmission lines and three generator step-up (GSU) transformers: finding that losing any GSU transformer isolates its generator entirely, while line outages redistribute flow onto parallel corridors by roughly 62%–253%.
• Meshed vs. radial reconfiguration: confirming that forcing radial operation increases system losses from a base case 4.64 MW up to 13.21 MW.
• Progressive loading scan: showing voltage drops of 5–10% at heavily loaded buses once demand reaches roughly 1.55×–1.9× nominal, with the solver failing to converge beyond about 2.2× nominal load.
• Transformer thermal loading: showing the most heavily loaded GSU unit exceeding its nameplate rating at roughly 1.55× nominal demand, ahead of the network-wide voltage limit.

A supplementary three-phase symmetrical fault study, performed with the bus impedance (Z_bus) method, characterised fault current magnitudes at representative buses for protection coordination context. Newton-Raphson convergence was consistently fast (4–8 iterations), owing to the small, well-conditioned nature of the benchmark. Overall, the study shows the custom solver reproduces standard power system behaviour with high numerical fidelity and provides an auditable, dependency-light platform for contingency and voltage-stability screening.

Objectives & scope

Modern transmission networks must remain secure not only under normal conditions but also following the unplanned loss of individual elements, the N-1 security criterion. Verifying this requires repeatedly solving the AC power flow equations under a range of network topologies and loading conditions. As demand grows, or as generation and transmission assets are lost, networks can also approach voltage instability, where bus voltages degrade sharply for comparatively small further increases in load.

This project implements the full analytical workflow needed to study these phenomena: a from-scratch Newton-Raphson AC power flow solver, an independent cross-validation solver, an automated N-1 contingency screening routine, a meshed-versus-radial reconfiguration comparison, a progressive loading voltage stability scan, a transformer thermal loading check, and a three-phase fault current calculation.

Objectives

• Implement a Newton-Raphson AC power flow solver from first principles (Y_bus construction, polar-form mismatch equations, analytical Jacobian) for slack, PV, and PQ buses.
• Validate solver accuracy primarily against pandapower, with an independently coded Gauss-Seidel solver as a secondary cross-check.
• Screen the network for N-1 contingency security, quantifying post-outage flow redistribution and identifying resulting radial (demeshed) configurations.
• Characterise voltage stability behaviour by scanning demand upward from nominal, tracking voltage decline at the most heavily loaded buses and GSU transformer thermal loading.
• Extend the steady-state solver to a basic fault response capability (three-phase symmetrical fault analysis).

Scope & test network

All analyses are performed on the WSCC 9-bus system (the Anderson & Fouad benchmark), consisting of three generating units (one slack, two PV), three load buses, three junction buses, three GSU transformers, and six 230 kV lines arranged in a single meshed ring. It was chosen because its topology and data are publicly documented and reproducible, keeping the study fully self-contained and auditable.

Figure 1: WSCC 9-bus test network diagram. G1/G2/G3 are the slack and PV generator buses; 5, 6, 8 are load buses; 4, 7, 9 are junction buses.

Network data, formulation & validation strategy

Test network data

The WSCC 9-bus system comprises 9 buses on a 100 MVA system base. Bus 1 (16.5 kV) is the slack bus; Buses 2 (18.0 kV) and 3 (13.8 kV) are PV buses; Buses 4–9 (230 kV) are junction and load buses. Three GSU transformers connect Buses 1, 2, and 3 to the 230 kV network at Buses 4, 7, and 9 respectively. Six 230 kV lines form a single meshed ring: 4–5, 4–6, 5–7, 6–9, 7–8, and 8–9.

Newton-Raphson AC power-flow formulation

For each bus i, injected active and reactive power are expressed in polar coordinates as functions of all bus voltage magnitudes V, angles θ, and the network admittance matrix Y_bus = G + jB. The slack bus fixes V and θ as the angle reference; PV buses specify P and V (Q floats); PQ buses specify both P and Q (V floats). The mismatch vector ΔF is driven to zero by Newton's method:

The Jacobian J is assembled analytically in its standard four-block form (∂P/∂θ, ∂P/∂V, ∂Q/∂θ, ∂Q/∂V) and re-evaluated every iteration. Y_bus is built from branch R, X, and shunt charging B using the standard π-model stamping procedure, with GSU transformers modelled as series reactances. The solver iterates until the largest absolute mismatch falls below a tolerance (1×10⁻⁶–1×10⁻⁹ pu depending on the study) or a maximum iteration count is reached.

Two initial guess conventions are used and compared: a "generator-informed" flat start (V = 1.0 pu, θ = 0 at PQ buses, PV/slack initialised at scheduled setpoints), and a fully flat start (V = 1.0 pu, θ = 0 everywhere, including generators).

Validation strategy

The primary validation method is a direct comparison against pandapower. A dedicated script, src/pandapower_validation.py, builds the identical WSCC 9-bus network in pandapower's data model and calls pandapower.runpp to obtain an independent solution for direct comparison against the custom solver's bus voltages and angles.

As a secondary cross-check, a Gauss-Seidel solver was independently implemented from a completely different iterative formulation (fixed-point bus voltage updates rather than Newton linearisation) and run against the identical Y_bus/injection model. Because Gauss-Seidel solves the same nonlinear equations by an unrelated numerical route, close agreement between it and the Newton-Raphson solution is meaningful evidence that both the network model and the Newton-Raphson implementation are correct.

The meshed-vs-radial comparison opens each of the six ring lines in turn and compares total real-power losses and bus voltage spread against the closed-ring base case. The loading scan scales all three load buses together from 1.00× to 2.20× nominal in 0.05× steps, holding generation dispatch fixed, and re-solves at each step, continuing until the solver fails to converge (used as a practical, non-continuation-method indicator of maximum loadability). At every step, GSU transformer apparent power is also expressed as a percentage of nameplate rating. The fault analysis applies a bolted three-phase fault at bus k using the bus impedance matrix Z_bus = Y_bus⁻¹ and the pre-fault voltage profile.

Base case, contingencies, stability & faults

4.1Base case power flow

With generator buses initialised at their scheduled voltage setpoints (the standard flat start), the Newton-Raphson solver converged in 4 iterations to a maximum mismatch tolerance of 1×10⁻⁶ pu. With a fully flat start, convergence took 5 iterations, both reflecting the well-known quadratic convergence property of Newton-Raphson on a small, well-conditioned network.

Figure 2: Base case bus voltage profile. All buses remain within the ±5% (0.95–1.05 pu) operating band.

Figure 3: Newton-Raphson convergence history for both initial guess conventions (log-scale mismatch vs. iteration), illustrating quadratic convergence.

4.2Validation

The independently coded Gauss-Seidel solver converged in 43 iterations (tolerance 1×10⁻⁸ pu). Across all nine buses, the maximum deviation between the Newton-Raphson, pandapower, and Gauss-Seidel voltage magnitude solutions was 3.00×10⁻⁶%, several orders of magnitude inside the 2% threshold, confirming the Y_bus model and Newton-Raphson implementation are correct.

4.3N-1 contingency screening

All six line contingencies converged. The mildest, loss of line 4–6, redistributes the worst-affected surviving branch's flow by 61.5%; the most severe, loss of line 7–8, redistributes flow by 253.3% and pushes line 5–7 to 155.0 MVA against its 150 MVA rating (103% of rating), the only thermal violation found in the base-case N-1 screen. Loss of line 4–5 also drives Bus 5 to 0.839 pu, a voltage violation of the ±5% band.

All three GSU transformer contingencies resulted in immediate islanding: each generator has only one radial connection into the 230 kV mesh, so losing its step-up transformer disconnects the unit entirely rather than producing a converged but stressed power flow case.

Figure 4: Maximum post-outage branch flow redistribution for each N-1 line contingency (%). All three GSU transformer contingencies caused islanding rather than redistribution.

4.4Meshed (ring) vs. radial configuration

The closed-ring base case sustains total transmission losses of 4.64 MW. Forcing radial operation by opening any one of the six ring lines increases losses in every case, from a modest +0.71 MW (opening 8–9) up to +8.57 MW (opening 5–7), a loss penalty of up to roughly 185% relative to meshed operation. Radial configurations also show a wider bus voltage spread, confirming that ring/meshed operation improves both loss performance and voltage profile relative to radial feeding, at the cost of the additional protection coordination a meshed topology requires.

Figure 5: Total real-power transmission losses (MW): meshed operation vs. each single-line-open radial reconfiguration.

4.5Voltage stability & increased demand loading scan

As system demand is scaled uniformly upward, voltage at the three load buses declines monotonically. Bus 5, the electrically weakest, furthest from the nearest generator along the highest-impedance corridor, first crosses a 5% drop at 1.55× nominal load and a 10% drop at 1.9× nominal load. Buses 6 and 8 follow the same trend but decline more gently, reflecting their stronger electrical connection to the nearest generator.

The power flow continued to converge up to 2.2× nominal demand; beyond that point (tested up to 2.30×) the solver failed to converge for every initial guess attempted, consistent with the network approaching its maximum loadability (voltage collapse) boundary in this region.

Figure 6: Voltage decline at the three load buses vs. increasing system demand (× nominal). Dashed line marks the network's practical convergence limit near 2.2×.

4.6Transformer thermal loading

At the base case, GSU transformer T1 (Gen 1, 250 MVA) carries only 30.6% of rating, while T2 (200 MVA) and T3 (100 MVA) are already comparatively heavily loaded at 81.6% and 86.3% respectively. As demand grows, T1's loading rises fastest, crossing 100% of its rating at approximately 1.55× nominal demand, noticeably before the network-wide voltage stability limits become binding. T3 also approaches its rating (100.1% at 1.90×) later in the scan, while T2 remains below 95% throughout the tested range. This identifies T1 thermal loading as the first equipment-level constraint the system encounters under sustained load growth, ahead of the voltage stability boundary itself.

Figure 7: GSU transformer loading (% of nameplate MVA rating) vs. increasing system demand. T1 is the first unit to exceed its thermal rating.

4.7Three-phase symmetrical fault analysis

A bolted three-phase fault was applied in turn at Bus 1 and at each of Buses 4, 5, 6, 7, 8, and 9. At every faulted bus the local voltage collapses to zero, as expected, while remaining buses retain a partial voltage depending on their electrical distance from the fault. The fault at Bus 1 produces the highest per-unit fault current (1.670 pu) due to the generator's low internal (transformer) impedance, while faults on the 230 kV network produce broadly similar per-unit currents (1.46–1.51 pu), reflecting the relatively uniform impedance of the meshed ring as seen from any one of its buses.

Figure 8: Three-phase bolted fault current magnitude (kA) by faulted bus location.

What the results mean operationally

5.1Convergence behaviour

The project scope originally targeted stable convergence within 10–12 Newton-Raphson iterations. In practice, the WSCC 9-bus benchmark converged in 4 iterations from a generator-informed start and 5 from a fully flat start, rising only to 6–8 iterations even under combined heavy loading and contingency stress, before the solver stopped converging altogether beyond the network's loadability limit. This isn't a solver shortfall, quadratic convergence on a compact, strongly meshed 9-bus system is expected to be fast, but the 10–12 iteration figure is more characteristic of larger transmission networks, more weakly conditioned systems, or industrial-scale EMS runs solved without a warm start.

5.2Contingency security implications

The N-1 screening results split cleanly into two risk categories. Line contingencies stress the network by redistributing flow within the surviving mesh, severely in some cases, but the network continues to operate, converging to a valid (if more stressed) solution every time. Transformer contingencies are categorically different: because each generator has only a single radial GSU transformer into the mesh, losing it removes the generator entirely rather than stressing the network. Line contingencies are network security events requiring flow and voltage limit checks; GSU transformer contingencies are generation adequacy events requiring reserve/redispatch planning, and the islanding detection automatically flags which category each contingency falls into.

5.3Voltage stability & equipment margins

Sections 4.5 and 4.6 together show that demand-growth headroom is bounded by two different mechanisms that bind at different points: transformer T1 reaches its thermal rating at roughly 1.55× nominal demand, well before bus voltage limits become binding, and well before the network approaches voltage collapse near 2.2×. In this network, equipment thermal limits, not voltage stability, are the first constraint encountered under sustained, uniform load growth. Voltage stability analysis alone would understate the actual growth headroom if performed without also checking equipment thermal ratings.

Meshed operation as a loss & voltage stability resource

The meshed-vs-radial comparison reinforces the contingency results: every radial configuration increases both losses and voltage spread relative to the closed-ring base case, by amounts that vary widely with which line is opened (+0.71 MW to +8.57 MW). The same physical event, a line outage, has a materially different impact depending on which corridor is affected, which is precisely why blanket N-1 rules of thumb are insufficient and case-by-case contingency screening remains necessary even on a network this small.

Where the results should be read carefully

Where this could go next

• Add PV bus reactive power limit enforcement (automatic PV → PQ switching at Q_max/Q_min), a standard feature of production-grade power flow solvers.
• Extend to a larger IEEE benchmark (30, 57, or 118 bus) to test iteration-count behaviour on a more weakly conditioned system and obtain results on a more realistically meshed topology.
• Implement a continuation power flow (predictor-corrector) method to trace the complete P-V and Q-V nose curves for a rigorous maximum loadability margin.
• Extend fault analysis to negative- and zero-sequence networks to support single-line-to-ground and line-to-line fault calculations.
• Source or design actual thermal ratings for branches and transformers if this workflow is applied to a real utility network.
• Add time-domain transient stability simulation (swing equation integration following a fault).

A complete, auditable power flow toolchain

This project delivered a complete power flow analysis toolchain: a Newton-Raphson AC solver, a primary pandapower validation script, a secondary independently coded Gauss-Seidel cross-check, N-1 contingency screening, meshed-versus-radial reconfiguration comparison, a voltage stability loading scan, transformer thermal loading checks, and a three-phase fault current calculation, applied to the WSCC 9-bus transmission benchmark.

The solver converged reliably across all tested operating conditions up to the network's maximum loadability boundary, with the pandapower check showing agreement to within 3.00×10⁻⁶%. The analyses produced several concrete, physically consistent engineering findings: N-1 line contingencies redistribute flow onto surviving corridors by tens to hundreds of percent depending on severity, while the network's three GSU transformers have no contingency redundancy and instead island their generators when lost; meshed operation reduces both losses and voltage spread relative to any radial configuration; heavily loaded buses see voltage drops in the operationally significant 5–10% range once demand grows to roughly 1.55×–1.9× nominal; and, on this network, a transformer thermal limit is reached before the voltage stability boundary itself becomes binding.

Taken together, the results demonstrate a sound, auditable, and extensible power system analysis capability built without dependence on external power flow libraries for its core solver, with a third-party pandapower cross-check confirming its functionality.

Software structure

The analysis toolchain is implemented in Python 3 (numpy, scipy, matplotlib):

• data/network_9bus.py: WSCC 9-bus network dataset (buses, generators, loads, lines, transformers).
• src/ybus.py: Admittance-matrix (Y_bus) builder from branch data, with optional branch outage exclusion for contingency studies.
• src/newton_raphson.py: Newton-Raphson AC power flow solver (polar form, analytical Jacobian) and branch power flow post-processing.
• src/pandapower_validation.py: builds the identical network in pandapower and runs pandapower's own AC power flow solver for comparison.
• src/gauss_seidel.py: independently coded Gauss-Seidel AC power flow solver, the project's secondary cross-check.
• src/network_setup.py: assembles solver-ready Y_bus/injection arrays from the raw network dataset, with load scaling and outage index options.
• src/contingency.py: automated N-1 contingency screening across all branches, including islanding detection and thermal/voltage violation checks.
• src/fault_analysis.py: three-phase symmetrical fault analysis via the bus impedance (Z_bus) method.
• src/main.py: orchestrates all analyses, saves results to results/results.json, and generates figures.

All numerical results and figures in this report were generated directly by this codebase; no numbers were entered by hand. The complete source code and results.json data file accompany this report.

Console output: results log generated by src/main.py on a full run of the toolchain.