12.4 Power Flow Analysis & System Stability (Steady-State, Transient & Swing Equation)

12.4 Power Flow Analysis & System Stability (Steady-State, Transient & Swing Equation)

Executive Overview: Power system stability encompasses the ability of interconnected synchronous machines to remain in synchronism following normal operational adjustments and severe transient disturbances. On the PE Power exam, stability analysis spans three critical areas: formulating power flow equations and classifying network buses (Slack, PV, PQ); calculating steady-state power-angle limits ($P = \frac{V_1 V_2}{X}\sin\delta$) and synchronizing stiffness; and analyzing rotor electromechanical transients via the Swing Equation and Equal Area Criterion (EAC) to calculate Critical Clearing Angle ($\delta_{cr}$) and Critical Clearing Time ($t_{cr}$).

---

1. Power Flow Problem Formulation & Bus Classification

In an $N$-bus power system, each bus $i$ is characterized by four electrical state variables: active power ($P_i$), reactive power ($Q_i$), voltage magnitude ($|V_i|$), and voltage phase angle ($\delta_i$).

Nodal Power Balance at Bus i:

         (Generation: P_Gi + j*Q_Gi) 
                    |
                    v
       =============o=============  Bus i (|V_i| /_ delta_i)
                    |       \
                    v        ------> To interconnected buses k via Y_ik
           (Load: P_Li + j*Q_Li)

Bus Classification Architecture

At each bus, two variables are specified as known inputs, while the remaining two variables are calculated via iterative numerical solution:

Nodal Admittance & Newton-Raphson Formulation

The complex power injected into bus $i$ is governed by the bus admittance matrix $\mathbf{Y}_{bus}$:

$$P_i - j Q_i = V_i^* \sum_{k=1}^N Y_{ik} V_k = |V_i| \sum_{k=1}^N |V_k| |Y_{ik}| \angle(\theta_{ik} - \delta_i + \delta_k)$$

Separating into real and reactive power equations:

$$P_i = |V_i| \sum_{k=1}^N |V_k| |Y_{ik}| \cos(\theta_{ik} - \delta_i + \delta_k), \qquad Q_i = -|V_i| \sum_{k=1}^N |V_k| |Y_{ik}| \sin(\theta_{ik} - \delta_i + \delta_k)$$

The non-linear power mismatch equations are solved iteratively using the Newton-Raphson Jacobian matrix $\mathbf{J}$:

$$\begin{bmatrix} \Delta \mathbf{P} \\ \Delta \mathbf{Q} \end{bmatrix} = \begin{bmatrix} \mathbf{J}_{11} & \mathbf{J}_{12} \\ \mathbf{J}_{21} & \mathbf{J}_{22} \end{bmatrix} \begin{bmatrix} \Delta \boldsymbol{\delta} \\ \Delta |\mathbf{V}| / |\mathbf{V}| \end{bmatrix} = \begin{bmatrix} \frac{\partial \mathbf{P}}{\partial \boldsymbol{\delta}} & \frac{\partial \mathbf{P}}{\partial |\mathbf{V}|} \\ \frac{\partial \mathbf{Q}}{\partial \boldsymbol{\delta}} & \frac{\partial \mathbf{Q}}{\partial |\mathbf{V}|} \end{bmatrix} \begin{bmatrix} \Delta \boldsymbol{\delta} \\ \Delta |\mathbf{V}| \end{bmatrix}$$

• Fast Decoupled Power Flow (FDLF): Leverages the high $X/R$ ratio of transmission lines ($P$ strongly couples to $\delta$; $Q$ strongly couples to $|V|$), setting $\mathbf{J}_{12} \approx 0$ and $\mathbf{J}_{21} \approx 0$ to solve two independent, constant sub-matrices ($[\mathbf{B}']$ and $[\mathbf{B}'']$).

---

2. Steady-State Power-Angle Relationship & Synchronizing Power

For a synchronous generator delivering power across total transfer reactance $X_{total} = X_d' + X_{tr} + X_{line}$ to an infinite bus ($V_\infty = |V_\infty|\angle 0^\circ$):

$$P_e = \frac{|E'||V_\infty|}{X_{total}} \sin\delta = P_{max} \sin\delta$$ $$Q_e = \frac{|E'|^2 - |E'||V_\infty|\cos\delta}{X_{total}}$$

Steady-State Power-Angle Curve P(delta):

 Active Power P |
                |                 P_max (Peak at delta = 90 deg)
          P_max |                    .---.
                |                  /       \
            P_m |----------------o           o (Static Limit delta = 180 - delta_0)
                |               /|           |\
                |              / |           | \
                |             /  |           |  \
              0 +------------o---+-----------+---o-------------> Rotor Angle delta
                            0   delta_0     90   180

Steady-State Stability Margin (SSM)

$$SSM = \frac{P_{max} - P_0}{P_{max}} \times 100\% = \left( 1 - \sin\delta_0 \right) \times 100\%$$

Synchronizing Power Coefficient ($P_{syn}$)

The synchronizing power coefficient represents the electrical "spring stiffness" holding the machine in synchronism. It is the derivative of active power with respect to power angle:

$$P_{syn} = \frac{dP_e}{d\delta} = \frac{|E'||V_\infty|}{X_{total}} \cos\delta = P_{max} \cos\delta \quad [\text{MW/rad or pu/rad}]$$

• When $0^\circ \le \delta < 90^\circ$: $P_{syn} > 0$ (System is stable; an incremental increase in $\delta$ increases $P_e$, counteracting acceleration).
• When $\delta = 90^\circ$: $P_{syn} = 0$ (Steady-State Stability Limit).
• When $90^\circ < \delta \le 180^\circ$: $P_{syn} < 0$ (System is unstable; increasing $\delta$ reduces $P_e$, causing runaway rotor acceleration).

---

3. Rotor Dynamics & The Swing Equation

During a system fault, mechanical input power $P_m$ from the prime mover remains temporarily constant (due to governor inertia), while electrical power output $P_e$ drops abruptly. The net accelerating power $P_a = P_m - P_e$ accelerates the rotor mass.

Inertia Constant ($H$) and Angular Momentum ($M$)

The inertia constant $H$ normalizes stored kinetic energy at synchronous speed to the machine MVA base:

$$H = \frac{\text{Stored Kinetic Energy at Synchronous Speed } (E_k)}{\text{Machine MVA Rating } (S_{base})} = \frac{\frac{1}{2} J \omega_{sm}^2}{S_{base}} \quad [\text{MW}\cdot\text{s/MVA} \text{ or seconds}]$$

Angular momentum $M$ is expressed as:

$$M = \frac{2 H S_{base}}{\omega_s} = \frac{H}{\pi f_0} \quad [\text{MJ}\cdot\text{s/elec rad}], \qquad M = \frac{H}{180 f_0} \quad [\text{MJ}\cdot\text{s/elec degree}]$$

The Per-Unit Swing Equation

$$\frac{2 H}{\omega_s} \frac{d^2 \delta}{dt^2} = P_m - P_e = P_a \quad [\text{per unit}]$$

In terms of electrical degrees:

$$\frac{H}{180 f_0} \frac{d^2 \delta}{dt^2} = P_m - P_e = P_a \quad [\text{per unit}]$$

Initial rotor acceleration immediately following a fault (at $t = 0^+$ where $P_e = P_{e,fault}$):

$$\alpha = \left. \frac{d^2 \delta}{dt^2} \right|_{t=0^+} = \frac{180 f_0}{H} (P_m - P_{e,fault}) \quad [^\circ/\text{s}^2] = \frac{\pi f_0}{H} (P_m - P_{e,fault}) \quad [\text{rad/s}^2]$$

---

4. Transient Stability & The Equal Area Criterion (EAC)

For a Single-Machine Infinite-Bus (SMIB) system, the Equal Area Criterion evaluates transient stability graphically by equating kinetic energy stored during fault acceleration to potential energy absorbed during post-fault deceleration.

Equal Area Criterion Power-Angle Curves:

 Power P |
         |                  Curve 1: Pre-Fault (P_max1)
   P_max1|                    .---.
         |                  /   |   \
   P_max3|----------Curve 3: Post-Fault (P_max3)----
         |          /     /  A_2|     \      \
     P_m |---------o=====o======|======o------\--- P_m
         |        /| A_1 |      |      |\      \
   P_max2|----Curve 2: During-Fault (P_max2)-----
         |      /  |     |      |      |  \      \
       0 +-----o---+-----+------+------+---+------o-------> delta
              0  delta_0 delta_cr    delta_max   180

The Three Operating Curves

1. Pre-Fault Curve (1): $P_{e1} = P_{max1} \sin\delta$, initial operating angle $\delta_0 = \arcsin(P_m / P_{max1})$.
2. During-Fault Curve (2): $P_{e2} = P_{max2} \sin\delta$ ($P_{max2} \ll P_{max1}$; if bolted 3-phase fault at bus, $P_{max2} = 0$).
3. Post-Fault Curve (3): $P_{e3} = P_{max3} \sin\delta$ (Fault cleared by isolating one parallel line; $P_{max3} < P_{max1}$).

Energy Balance Formulation

1. Accelerating Area ($A_1$): Kinetic energy gained from $\delta_0$ to clearing angle $\delta_{cr}$: $$A_1 = \int_{\delta_0}^{\delta_{cr}} (P_m - P_{e2}) \, d\delta = P_m (\delta_{cr} - \delta_0) + P_{max2}(\cos\delta_{cr} - \cos\delta_0)$$
2. Decelerating Area ($A_2$): Maximum kinetic energy returned to grid from $\delta_{cr}$ to $\delta_{max}$: $$A_2 = \int_{\delta_{cr}}^{\delta_{max}} (P_{e3} - P_m) \, d\delta = P_{max3}(\cos\delta_{cr} - \cos\delta_{max}) - P_m (\delta_{max} - \delta_{cr})$$ where the maximum allowable swing angle is: $\delta_{max} = \pi - \arcsin\left( \frac{P_m}{P_{max3}} \right) \text{ [radians]}$.

Critical Clearing Angle ($\delta_{cr}$)

Setting accelerating area equal to maximum decelerating area ($A_1 = A_2$) yields the exact equation for Critical Clearing Angle:

$$\cos\delta_{cr} = \frac{P_m (\delta_{max} - \delta_0) + P_{max3} \cos\delta_{max} - P_{max2} \cos\delta_0}{P_{max3} - P_{max2}}$$

Critical Clearing Time ($t_{cr}$) for Bolted Fault ($P_{max2} = 0$)

When electrical power is completely zero during the fault ($P_{e2} = 0$), acceleration is constant: $\frac{d^2 \delta}{dt^2} = \frac{\omega_s P_m}{2 H}$. Integrating twice from $\delta(0) = \delta_0$:

$$\delta(t) = \delta_0 + \frac{\omega_s P_m}{4 H} t^2 \implies t_{cr} = \sqrt{\frac{4 H (\delta_{cr} - \delta_0)}{\omega_s P_m}} = \sqrt{\frac{2 H (\delta_{cr} - \delta_0)}{\pi f_0 P_m}} \quad [\text{seconds}]$$

(where $\delta_{cr}$ and $\delta_0$ are in radians).

---

5. Comprehensive Worked Calculation: Swing Equation & Critical Clearing

Problem Statement

A $60\text{ Hz}$, $500\text{ MVA}$ synchronous generator with inertia constant $H = 3.5\text{ MJ/MVA}$ is connected to an infinite bus through two identical parallel transmission lines. The system is operating in steady-state delivering rated active power $P_m = 1.0\text{ pu}$ at $\mathbf{V}_\infty = 1.0\angle 0^\circ\text{ pu}$ with generator internal transient voltage $|E'| = 1.25\text{ pu}$.

The power-angle relationships for the three system states are:

• Pre-Fault (both lines in service): $P_{e1} = 2.00 \sin\delta\;\text{pu}$
• During-Fault (3-phase bolted fault at line terminal): $P_{e2} = 0.50 \sin\delta\;\text{pu}$
• Post-Fault (fault cleared by tripping one faulted line): $P_{e3} = 1.50 \sin\delta\;\text{pu}$

Calculate:

1. The initial steady-state operating power angle $\delta_0$ (in degrees and radians).
2. The initial rotor acceleration $\alpha$ at the instant of fault inception.
3. The maximum rotor angle $\delta_{max}$ for stability (in degrees and radians).
4. The Critical Clearing Angle $\delta_{cr}$ (in degrees).
5. The accelerating area $A_1$ and verify $A_1 = A_2$.

============================== STEP-BY-STEP SOLUTION ==============================

Step 1: Compute Initial Operating Angle delta_0
  P_m = 1.0 pu = P_max1 * sin(delta_0) = 2.00 * sin(delta_0)
  sin(delta_0) = 1.0 / 2.00 = 0.5000
  delta_0 = arcsin(0.5000) = 30.00 deg
  delta_0_rad = 30.00 * (pi / 180) = 0.52360 rad

Step 2: Compute Initial Rotor Acceleration at t = 0+
  At fault inception (delta = 30.0 deg):
  P_e2(0+) = 0.50 * sin(30.0 deg) = 0.50 * 0.50 = 0.250 pu
  Net accelerating power P_a = P_m - P_e2 = 1.00 - 0.250 = 0.750 pu
  
  Using Swing Equation acceleration:
  alpha = (180 * f_0 / H) * P_a = (180 * 60 / 3.5) * 0.750
        = (10,800 / 3.5) * 0.750 = 3,085.71 * 0.750 = 2,314.29 deg/s^2
  
  In electrical radians:
  alpha_rad = (pi * 60 / 3.5) * 0.750 = 53.856 * 0.750 = 40.392 rad/s^2

Step 3: Compute Maximum Allowable Swing Angle delta_max
  Post-fault peak power P_max3 = 1.50 pu
  delta_post_nominal = arcsin(P_m / P_max3) = arcsin(1.00 / 1.50) = arcsin(0.6667)
                     = 41.810 deg = 0.72973 rad
  
  delta_max = 180 deg - delta_post_nominal = 180 deg - 41.810 deg = 138.190 deg
  delta_max_rad = 138.190 * (pi / 180) = 2.41187 rad

Step 4: Compute Critical Clearing Angle delta_cr
  Formula:
  cos(delta_cr) = [ P_m * (delta_max_rad - delta_0_rad) + P_max3 * cos(delta_max) 
                    - P_max2 * cos(delta_0) ] / (P_max3 - P_max2)
  
  Evaluate each term in numerator:
  Term 1: P_m * (delta_max_rad - delta_0_rad) = 1.00 * (2.41187 - 0.52360) = 1.88827
  Term 2: P_max3 * cos(delta_max) = 1.50 * cos(138.190 deg) = 1.50 * (-0.74536) = -1.11803
  Term 3: -P_max2 * cos(delta_0) = -0.50 * cos(30.0 deg) = -0.50 * (0.86603) = -0.43301
  
  Numerator = 1.88827 - 1.11803 - 0.43301 = 0.33723
  Denominator = P_max3 - P_max2 = 1.50 - 0.50 = 1.00
  
  cos(delta_cr) = 0.33723 / 1.00 = 0.33723
  delta_cr = arccos(0.33723) = 70.293 deg = 1.22685 rad

Step 5: Verify Equal Area Criterion Energy Balance (A1 = A2)
  Accelerating Area A1:
  A1 = P_m * (delta_cr_rad - delta_0_rad) + P_max2 * [cos(delta_cr) - cos(delta_0)]
     = 1.00 * (1.22685 - 0.52360) + 0.50 * [0.33723 - 0.86603]
     = 0.70325 + 0.50 * (-0.52880) = 0.70325 - 0.26440 = 0.43885 pu-rad
  
  Decelerating Area A2:
  A2 = P_max3 * [cos(delta_cr) - cos(delta_max)] - P_m * (delta_max_rad - delta_cr_rad)
     = 1.50 * [0.33723 - (-0.74536)] - 1.00 * (2.41187 - 1.22685)
     = 1.50 * (1.08259) - 1.18502 = 1.62388 - 1.18502 = 0.43886 pu-rad
  
  A1 == A2 = 0.43885 pu-rad (Exact match confirms mathematical precision!)
===================================================================================

---

6. Common Exam Traps & Strategic Pitfalls

• The Radians vs. Degrees Mismatch in EAC: Forgetting that the linear angle subtraction in the EAC equation ($P_m(\delta_{max} - \delta_0)$) must be in radians, while arguments to trigonometric functions ($\sin, \cos$) are in degrees. Multiplying $P_m$ by $(138.19^\circ - 30.0^\circ) = 108.19$ instead of $1.888\text{ rad}$ yields impossible values for $\cos\delta_{cr} > 1$.
• The Steady-State Stability Limit Misconception: Believing the steady-state stability limit is $\delta = 180^\circ$. The theoretical steady-state limit occurs at $\delta = 90^\circ$ ($P_{syn} = 0$). Beyond $90^\circ$, synchronizing torque is negative, leading to spontaneous pull-out.
• Generator Inertia Base Conversion Blindspot: Failing to adjust $H$ when changing system MVA base. $H_{new} = H_{old} \times \left( \frac{S_{base,old}}{S_{base,new}} \right)$.
• Incorrect Post-Fault Peak Angle Calculation: Setting $\delta_{max} = 180^\circ - \delta_0$ instead of $\delta_{max} = 180^\circ - \arcsin(P_m / P_{max3})$. The maximum angle is governed by the post-fault power curve, not the pre-fault curve.