ECE 467: Power Systems Analysis, Operations, and Markets

Claudio Cañizares

Estimated study time: 1 hr 20 min

Table of contents

Sources and References

Primary textbook — A. Gómez-Expósito, A. J. Conejo, and C. A. Cañizares (Eds.), Electric Energy Systems: Analysis and Operation, 2nd ed., CRC Press, 2018. Supplementary texts — J. J. Grainger and W. D. Stevenson, Power Systems Analysis, McGraw-Hill, 1994; J. D. Glover, T. J. Overbye, and M. S. Sarma, Power Systems Analysis and Design, 6th ed., Cengage Learning, 2017; A. J. Wood, B. F. Wollenberg, and G. B. Sheblé, Power Generation, Operation, and Control, 3rd ed., Wiley, 2013.


Chapter 1: Power Systems Overview and Market Structure

1.1 Introduction

Modern electric power systems are among the most complex engineered systems in existence. They span geographical areas of thousands of square kilometres and must continuously match generation to load in real time. Understanding their structure, operation, and economic context is the foundation for every topic in this course.

1.2 Generation, Transmission, and Distribution

A power system consists of three principal subsystems: generation, which converts primary energy sources into electrical energy; transmission, which transports bulk power at high voltage over long distances with minimal losses; and distribution, which delivers power to end consumers at utilisation voltages.

Generation covers a wide portfolio of technologies. Synchronous generators driven by steam turbines (coal, natural gas, nuclear) or hydro turbines account for the majority of installed capacity in most grids. Inverter-based resources — wind turbines and photovoltaic plants interfaced through power electronics — now constitute a rapidly growing share. Each generation technology exhibits distinct cost structures, ramp rates, minimum stable output levels, and environmental characteristics that collectively shape system operations and markets.

Transmission networks operate at voltages ranging from 115 kV to 765 kV (extra-high voltage) or beyond. High voltage reduces the current required to transfer a given power, thereby reducing resistive losses proportional to \(I^2 R\). Transmission grids are typically meshed, providing multiple power-flow paths and improving reliability through redundancy. Key components include overhead lines, underground cables, transformers, shunt reactors and capacitors, and flexible AC transmission system (FACTS) devices.

Distribution networks operate at medium voltage (4–35 kV) and low voltage (120/240 V or 230/400 V) and are predominantly radial in configuration. Distribution systems connect millions of customers and increasingly host distributed energy resources (DER), including rooftop solar, battery storage, and electric vehicles, creating bidirectional power flows not present in traditional designs.

1.3 History of Electric Power Industry Restructuring

For most of the 20th century, electric utilities were vertically integrated monopolies regulated by government agencies. A single utility owned generation, transmission, and distribution assets and provided service within an exclusive franchise territory. Regulators approved rates based on cost of service, allowing recovery of prudent capital costs plus a regulated rate of return.

Beginning in the 1970s and accelerating through the 1990s, several forces pushed toward restructuring:

  • The Public Utility Regulatory Policies Act (PURPA, US, 1978) required utilities to purchase power from qualifying independent generators, introducing competition in generation.
  • The Energy Policy Act of 1992 (US) mandated open-access transmission, allowing third parties to use the transmission grid at cost-based rates.
  • FERC Orders 888 and 889 (1996) required non-discriminatory open access to transmission, separating transmission services from generation.
  • California restructured its market in 1998; the 2000–2001 energy crisis exposed significant market design flaws.
  • European liberalisation proceeded through EU Directives (1996, 2003, 2009), unbundling transmission from generation and retail.
Restructuring separates the competitive elements (generation, retail supply) from the natural monopoly elements (transmission, distribution wires). The latter remain regulated because the economics of replication do not support competing infrastructure.

1.4 Structure of Restructured Electricity Markets

In restructured environments, several institutional entities coordinate system operation and market clearing:

Independent System Operator (ISO) / Regional Transmission Organization (RTO): A neutral entity that operates the transmission grid, administers energy and ancillary service markets, and ensures reliable operation. Examples include CAISO (California), PJM (Mid-Atlantic), IESO (Ontario), and AESO (Alberta).

Market participants: Generators, load-serving entities, retailers, and financial traders submit bids and offers into market auctions.

Market timeline: Most markets operate a day-ahead (DA) market that clears 24 hourly intervals the day before real time, followed by a real-time (RT) or balancing market that adjusts for deviations. Ancillary service markets (regulation, spinning reserve, non-spinning reserve) clear simultaneously with energy.

Locational Marginal Pricing (LMP): In nodal markets, the clearing price at each node reflects the marginal cost of serving an additional increment of load at that location, including energy, congestion, and losses components:

\[ \text{LMP}_k = \lambda + \mu_k^{\text{cong}} + \mu_k^{\text{loss}} \]

where \(\lambda\) is the system energy price, \(\mu_k^{\text{cong}}\) captures transmission congestion, and \(\mu_k^{\text{loss}}\) captures marginal losses at node \(k\).


Chapter 2: Power Systems Modeling

2.1 Introduction to System Modeling

Analytical study of a power system requires mathematical models of every component. This chapter develops models for the per-unit normalisation framework, transmission lines, transformers, synchronous machines, induction machines, and loads.

2.2 Per-Unit System

Motivation and Base Quantities

Transformers change voltage and current levels throughout a power system. Comparing impedances across different voltage levels is cumbersome in physical units. The per-unit (pu) system normalises all quantities by base values, collapsing transformer turns ratios and facilitating uniform analysis.

For a single-phase system, choose a base apparent power \(S_{\text{base}}\) (VA) and a base voltage \(V_{\text{base}}\) (V) at a chosen point in the network. The remaining base quantities follow: \[ I_{\text{base}} = \frac{S_{\text{base}}}{V_{\text{base}}}, \qquad Z_{\text{base}} = \frac{V_{\text{base}}^2}{S_{\text{base}}} \]

The per-unit value of any quantity \(X\) is \(x_{\text{pu}} = X / X_{\text{base}}\).

For three-phase systems with line voltage \(V_{\text{base,LL}}\) and three-phase apparent power \(S_{\text{base,3\phi}}\):

\[ Z_{\text{base}} = \frac{(V_{\text{base,LL}})^2}{S_{\text{base,3\phi}}} \]

Base Selection Across Transformer Banks

The system base \(S_{\text{base}}\) is constant throughout the network. When a transformer connects two zones, the voltage base changes according to the turns ratio. If the transformer has turns ratio \(a = N_1/N_2\), then:

\[ V_{\text{base,2}} = \frac{V_{\text{base,1}}}{a} \]

and the impedance base transforms as:

\[ Z_{\text{base,2}} = \frac{V_{\text{base,2}}^2}{S_{\text{base}}} \]

Impedance Conversion Between Bases

Equipment nameplates give impedances in pu on the equipment’s own base (\(S_{\text{rated}}\), \(V_{\text{rated}}\)). To convert to system base:

\[ z_{\text{pu,new}} = z_{\text{pu,old}} \times \frac{S_{\text{base,new}}}{S_{\text{base,old}}} \times \left(\frac{V_{\text{base,old}}}{V_{\text{base,new}}}\right)^2 \]
Example — Per-unit impedance conversion: A transformer rated 100 MVA, 138/13.8 kV has a leakage reactance of 0.10 pu on its own base. The system base is 200 MVA, with voltage bases of 138 kV and 13.8 kV (consistent with the transformer's turns ratio). Find the transformer reactance on system base.

Since the voltage bases are already consistent:

\[ x_{\text{pu,sys}} = 0.10 \times \frac{200\,\text{MVA}}{100\,\text{MVA}} \times \left(\frac{138\,\text{kV}}{138\,\text{kV}}\right)^2 = 0.10 \times 2 \times 1 = 0.20\,\text{pu} \]

2.3 Transmission Line Models

Distributed Parameter Model

A transmission line of length \(\ell\) is characterised by distributed series resistance \(r\) (Ω/km), series inductance \(l\) (H/km), shunt conductance \(g\) (S/km), and shunt capacitance \(c\) (F/km). Define \(\gamma = \sqrt{(r + j\omega l)(g + j\omega c)}\) as the propagation constant and \(Z_c = \sqrt{(r + j\omega l)/(g + j\omega c)}\) as the characteristic impedance.

The exact two-port ABCD model is:

\[ \begin{pmatrix} V_S \\ I_S \end{pmatrix} = \begin{pmatrix} \cosh(\gamma\ell) & Z_c \sinh(\gamma\ell) \\ \sinh(\gamma\ell)/Z_c & \cosh(\gamma\ell) \end{pmatrix} \begin{pmatrix} V_R \\ I_R \end{pmatrix} \]

Nominal-\(\pi\) Equivalent (Medium Lines)

For lines up to roughly 300 km, the nominal-\(\pi\) model adequately represents the line. The series impedance is \(Z = (r + j\omega l)\ell\) and total shunt admittance is \(Y = j\omega c \ell\). Half the shunt admittance is placed at each terminal:

\[ Y_{\text{shunt}} = Y/2 \text{ at each end} \]

This model is used throughout load-flow analysis in this course.

Short-Line Model

For lines under 80 km, shunt capacitance is negligible and the model reduces to a simple series impedance \(Z = R + jX\).

2.4 Transformer Models

The ideal transformer model introduces a turns ratio \(a\) that scales voltage and current. For an off-nominal turns ratio transformer (e.g., a tap-changing transformer), the circuit is represented as an ideal transformer of ratio \(a : 1\) in series with the leakage impedance \(Z_t\). In the admittance form used for the Y-bus:

\[ Y_{ii}^{\text{tap}} = \frac{Y_t}{a^2}, \quad Y_{ij}^{\text{tap}} = -\frac{Y_t}{a}, \quad Y_{ji}^{\text{tap}} = -\frac{Y_t}{a}, \quad Y_{jj}^{\text{tap}} = Y_t \]

where \(Y_t = 1/Z_t\) is the transformer series admittance.

Phase-shifting transformers introduce a complex turns ratio \(a e^{j\phi}\) that shifts the voltage phase angle and are used to control active power flow on transmission lines.

2.5 Synchronous Machine Models

Steady-State Equivalent Circuit

Under balanced steady-state conditions, a synchronous generator is modelled by an internal EMF \(\bar{E}\) behind a synchronous reactance \(X_s\) (round-rotor approximation):

\[ \bar{E} = \bar{V}_t + jX_s \bar{I}_a \]

where \(\bar{V}_t\) is the terminal voltage and \(\bar{I}_a\) is the armature current. The internal EMF magnitude is controlled by field excitation and determines the reactive power output; the power angle \(\delta\) between \(\bar{E}\) and \(\bar{V}_t\) determines active power output.

Salient-pole machines have different direct-axis and quadrature-axis reactances \(X_d \neq X_q\), requiring decomposition of currents into \(d\) and \(q\) components for accurate modelling.

Dynamic Models

For stability studies, the swing equation (Chapter 10) governs rotor dynamics. The machine presents sub-transient reactance \(X_d''\) immediately after a disturbance, transitioning to transient reactance \(X_d'\), and finally to steady-state \(X_d\), with corresponding time constants \(T_{d0}''\) and \(T_{d0}'\).

2.6 Induction Machine Models

An induction motor is modelled by the standard equivalent circuit with slip \(s = (\omega_s - \omega_r)/\omega_s\). The rotor branch impedance is \(R_r/s + jX_r\). At slip \(s = 0\) the machine runs at synchronous speed (no torque); at \(s = 1\) the rotor is stalled. For load-flow studies, induction motors are commonly represented as constant-power or constant-impedance loads, though detailed dynamic models track the rotor flux states for transient stability simulation.

2.7 Load Models

Load models relate active and reactive power consumption to bus voltage and frequency:

Constant power (P): \(P = P_0\), \(Q = Q_0\) — independent of voltage.

Constant current (I): \(P = P_0 \lvert V \rvert / V_0\), \(Q = Q_0 \lvert V \rvert / V_0\).

Constant impedance (Z): \(P = P_0 (\lvert V \rvert / V_0)^2\), \(Q = Q_0 (\lvert V \rvert / V_0)^2\).

Composite (ZIP): \(P = P_0 [a_Z (\lvert V \rvert/V_0)^2 + a_I (\lvert V \rvert/V_0) + a_P]\), with \(a_Z + a_I + a_P = 1\).

Voltage-dependent load models are important for voltage stability analysis. Frequency-dependent load models include a term \(K_f \Delta f\) and are relevant for frequency stability and control studies.


Chapter 3: Power Flow Analysis

3.1 Introduction

Power flow (load flow) analysis determines the steady-state bus voltages, line currents, and power flows for a given operating condition. It is the fundamental computational tool for planning, operations, and contingency analysis.

3.2 Power Flow Equations

Nodal Admittance Matrix (Y-bus)

For an \(n\)-bus network, the Y-bus is defined by:

\[ I_{\text{bus}} = Y_{\text{bus}} \cdot V_{\text{bus}} \]

where \(\mathbf{Y}_{\text{bus}} \in \mathbb{C}^{n \times n}\). The diagonal element \(Y_{ii}\) equals the sum of all admittances connected to bus \(i\); the off-diagonal element \(Y_{ij} = Y_{ji}\) equals the negative of the admittance of the branch directly connecting buses \(i\) and \(j\) (zero if no direct branch).

Injection Equations

The complex power injection at bus \(i\) is:

\[ S_i = P_i + jQ_i = V_i I_i^* = V_i \left(\sum_{k=1}^{n} Y_{ik} V_k\right)^* \]

Writing \(V_i = V_i e^{j\theta_i}\) and \(Y_{ik} = G_{ik} + jB_{ik}\):

\[ P_i = V_i \sum_{k=1}^{n} V_k (G_{ik} \cos\theta_{ik} + B_{ik} \sin\theta_{ik}) \]\[ Q_i = V_i \sum_{k=1}^{n} V_k (G_{ik} \sin\theta_{ik} - B_{ik} \cos\theta_{ik}) \]

where \(\theta_{ik} = \theta_i - \theta_k\).

3.3 Bus Classification

Each bus in the power flow problem has four variables: \(P_i, Q_i, V_i, \theta_i\). Two are specified (known), two are unknowns to be solved.

Bus TypeSpecifiedUnknown
Slack (reference)\(\lvert V \rvert\), \(\theta = 0\)\(P\), \(Q\)
PV (generator)\(P\), \(\lvert V \rvert\)\(Q\), \(\theta\)
PQ (load)\(P\), \(Q\)\(\lvert V \rvert\), \(\theta\)

The slack bus absorbs the mismatch between total generation and total load plus losses and provides the phase reference.

3.4 Newton-Raphson Power Flow

Mismatch Equations

Define the power mismatch at each non-slack bus:

\[ \Delta P_i = P_i^{\text{sch}} - P_i^{\text{calc}}(\mathbf{V}, \boldsymbol{\theta}) \]\[ \Delta Q_i = Q_i^{\text{sch}} - Q_i^{\text{calc}}(\mathbf{V}, \boldsymbol{\theta}) \]

For PQ buses both \(\Delta P_i\) and \(\Delta Q_i\) appear; for PV buses only \(\Delta P_i\) (since \(V_i\) is specified, \(Q_i\) is free and adjusted to hold voltage).

Jacobian Matrix Structure

The Newton-Raphson update at iteration \(\nu\) is:

\[ \begin{pmatrix} \Delta \mathbf{P} \\ \Delta \mathbf{Q} \end{pmatrix} = \begin{pmatrix} J_1 & J_2 \\ J_3 & J_4 \end{pmatrix} \begin{pmatrix} \Delta \boldsymbol{\theta} \\ \Delta \mathbf{V}/\mathbf{V} \end{pmatrix} \]

where \(\Delta \mathbf{V}/\mathbf{V}\) denotes the vector of fractional voltage magnitude corrections. The Jacobian submatrices contain:

Off-diagonal elements (\(i \neq k\)):

\[ J_{1,ik} = \frac{\partial P_i}{\partial \theta_k} = V_i V_k (G_{ik} \sin\theta_{ik} - B_{ik}\cos\theta_{ik}) \]\[ J_{2,ik} = \frac{\partial P_i}{\partial V_k} V_k = V_i V_k (G_{ik}\cos\theta_{ik} + B_{ik}\sin\theta_{ik}) \]\[ J_{3,ik} = \frac{\partial Q_i}{\partial \theta_k} = V_i V_k (-G_{ik}\cos\theta_{ik} - B_{ik}\sin\theta_{ik}) \]\[ J_{4,ik} = \frac{\partial Q_i}{\partial V_k}V_k = V_i V_k (G_{ik}\sin\theta_{ik} - B_{ik}\cos\theta_{ik}) \]

Diagonal elements (\(i = k\)):

\[ J_{1,ii} = -Q_i - B_{ii}V_i^2, \quad J_{2,ii} = P_i + G_{ii}V_i^2 \]\[ J_{3,ii} = P_i - G_{ii}V_i^2, \quad J_{4,ii} = Q_i - B_{ii}V_i^2 \]

Convergence Criterion

Iteration continues until the maximum mismatch satisfies a tolerance:

\[ \max_i (\lvert \Delta P_i \rvert, \lvert \Delta Q_i \rvert) \leq \epsilon \]

Typical tolerance: \(\epsilon = 10^{-5}\) pu. NR converges quadratically near the solution, requiring typically 3–7 iterations regardless of system size.

Example — Newton-Raphson iteration step: Consider a 3-bus system (bus 1 slack, bus 2 PV, bus 3 PQ). At the flat-start (\(V_i = 1, \theta_i = 0\)):

Calculate \(P_3^{\text{calc}}\) using the power flow equations, form the Jacobian (dimension \(3 \times 3\) for 1 PV + 1 PQ bus), solve for \([\Delta\theta_2, \Delta\theta_3, \Delta V_3/V_3]^T\), update \(\theta_2 \mathrel{+}= \Delta\theta_2\), etc., and re-evaluate mismatches. Repeat to convergence.

3.5 Gauss-Seidel Power Flow

The Gauss-Seidel (GS) method iterates bus voltages one at a time. For a PQ bus:

\[ V_i^{(\nu+1)} = \frac{1}{Y_{ii}}\left[\frac{P_i - jQ_i}{(V_i^{(\nu)})^*} - \sum_{k \neq i} Y_{ik} V_k^{(\nu)}\right] \]

GS is simple to implement and has low memory requirements but converges linearly (slowly) and may require hundreds of iterations for large systems. Acceleration factors \(\alpha \approx 1.6\) can improve convergence. NR is preferred in practice.

3.6 Fast Decoupled Power Flow

Decoupling Approximation

In high-voltage transmission networks, active power \(P\) is primarily determined by voltage angles \(\boldsymbol{\theta}\), while reactive power \(Q\) is primarily determined by voltage magnitudes \(\mathbf{V}\). This physical observation motivates the fast decoupled power flow (FDPF), which sets \(J_2 = J_3 = 0\) and simplifies the Jacobian submatrices.

Under the assumptions that \(G_{ik} \ll B_{ik}\), \(\cos\theta_{ik} \approx 1\), \(\sin\theta_{ik} \approx 0\), \(B_{ii} \gg G_{ii}\), and \(Q_i \ll B_{ii}V_i^2\):

\[ J_1 \approx -\text{diag}(V_i) \cdot B' \cdot \text{diag}(V_i) \]\[ J_4 \approx -\text{diag}(V_i) \cdot B'' \cdot \text{diag}(V_i) \]

Decoupled B’ and B’’ Matrices

The FDPF equations become:

\[ \frac{\Delta \mathbf{P}}{\mathbf{V}} = -B' \Delta \boldsymbol{\theta} \]\[ \frac{\Delta \mathbf{Q}}{\mathbf{V}} = -B'' \Delta \mathbf{V} \]

B’ matrix derivation: \(B'\) is formed from the network susceptance matrix, omitting shunt elements (capacitors, reactors) and the effects of off-nominal transformer taps. Element \(B'_{ij} = -1/x_{ij}\) for a purely reactive line of reactance \(x_{ij}\) between buses \(i\) and \(j\); diagonal \(B'_{ii} = \sum_k 1/x_{ik}\).

B’’ matrix derivation: \(B''\) is formed similarly but includes shunt susceptances and the reactive portions of the transformer model. PV bus rows and columns are eliminated from \(B''\) because voltage magnitudes at PV buses are fixed.

Both \(B'\) and \(B''\) are constant real matrices, so their LU factorizations are computed once and reused each iteration, giving very fast iteration cycles. Convergence is linear but acceptable for operations applications.

3.8 DC Power Flow Model

Derivation

The DC power flow (DCPF) is a linearised lossless approximation obtained from the AC flow equations under:

  1. \(V_i \approx 1.0\,\text{pu}\) for all buses
  2. \(\theta_{ik}\) is small: \(\sin\theta_{ik} \approx \theta_{ik}\), \(\cos\theta_{ik} \approx 1\)
  3. Lines are lossless: \(G_{ik} = 0\), so \(Y_{ik} = -jB_{ik}\)

Then \(P_i = \sum_k B_{ik}(\theta_i - \theta_k)\), or in matrix form:

\[ \mathbf{P} = \mathbf{B}_{\text{DC}} \boldsymbol{\theta} \]

where \(\mathbf{B}_{\text{DC}}\) is the DC susceptance matrix. Deleting the slack row/column gives a reduced system solvable for bus angles.

Line flow on branch \(i\text{-}k\) is:

\[ P_{ik} = \frac{\theta_i - \theta_k}{x_{ik}} \]

Power Transfer Distribution Factors (PTDF)

Define the PTDF as the sensitivity of line flow to power injections. For a unit injection at bus \(m\) and withdrawal at the slack bus:

\[ \text{PTDF}_{l,m} = \frac{\partial P_l}{\partial P_m} \]

Derivation: From \(\boldsymbol{\theta} = \mathbf{B}_{\text{DC}}^{-1} \mathbf{P}\) and \(P_l = (1/x_l)(\mathbf{e}_i - \mathbf{e}_k)^T \boldsymbol{\theta}\):

\[ \text{PTDF}_{l,m} = \frac{1}{x_l}(\mathbf{e}_i - \mathbf{e}_k)^T \mathbf{B}_{\text{DC}}^{-1} \mathbf{e}_m \]

PTDFs are constant (independent of operating point) and are used in contingency analysis and electricity market congestion management. The full PTDF matrix \(\mathbf{H} \in \mathbb{R}^{L \times n-1}\) (lines × non-slack buses) allows efficient computation of line flow changes:

\[ \Delta \mathbf{P}_{\text{lines}} = \mathbf{H} \Delta \mathbf{P}_{\text{inj}} \]

3.9 Y-Bus and Z-Bus Matrix Analysis

The Z-bus (impedance matrix) is the inverse of the Y-bus: \(\mathbf{Z}_{\text{bus}} = \mathbf{Y}_{\text{bus}}^{-1}\). The driving-point impedance \(Z_{ii}\) (diagonal element) is the Thevenin impedance seen from bus \(i\). The transfer impedance \(Z_{ij}\) relates voltage at bus \(i\) to current injected at bus \(j\). The Z-bus is the key tool in fault analysis.


Chapter 5 & Appendix B: Power System Economic Operations

5.2 Economic Load Dispatch

Problem Formulation

Economic load dispatch (ELD) determines the optimal allocation of total load \(P_D\) among \(N_g\) online generators to minimise total operating cost while satisfying power balance and generator limits.

Objective:

\[ \min \sum_{i=1}^{N_g} C_i(P_i) \]

Equality constraint (power balance):

\[ \sum_{i=1}^{N_g} P_i = P_D + P_L \]

where \(P_L\) is total transmission loss.

Inequality constraints (generator limits):

\[ P_i^{\min} \leq P_i \leq P_i^{\max}, \quad i = 1, \ldots, N_g \]

Cost curves \(C_i(P_i)\) are typically quadratic approximations to the heat rate curve of a thermal unit times the fuel cost:

\[ C_i(P_i) = a_i + b_i P_i + c_i P_i^2 \quad (\$/\text{hr}) \]

The incremental (marginal) cost is:

\[ \frac{dC_i}{dP_i} = b_i + 2c_i P_i \]

Equal Incremental Cost Criterion (Lossless Case)

Using Lagrangian optimisation with multiplier \(\lambda\) for the equality constraint:

\[ \mathcal{L} = \sum_{i=1}^{N_g} C_i(P_i) + \lambda\left(P_D - \sum_{i=1}^{N_g} P_i\right) \]

KKT stationarity conditions give:

\[ \frac{\partial \mathcal{L}}{\partial P_i} = \frac{dC_i}{dP_i} - \lambda = 0 \quad \Rightarrow \quad \frac{dC_i}{dP_i} = \lambda \quad \forall i \]
Equal Incremental Cost (Merit Order) Principle: At the economic dispatch optimum (ignoring limits and losses), all online generators operate at the same incremental cost \(\lambda\), the system lambda or locational marginal price (in the lossless case). This \(\lambda\) is determined by the power balance constraint.

For unconstrained generators with quadratic costs:

\[ P_i = \frac{\lambda - b_i}{2c_i} \]

Lambda Iteration Algorithm

Since \(\sum P_i(\lambda) = P_D\) must be satisfied, iterate \(\lambda\):

  1. Initialise \(\lambda\) (e.g., interpolate from marginal cost curves)
  2. Compute \(P_i = (\lambda - b_i)/(2c_i)\) for each generator
  3. Enforce limits: if \(P_i < P_i^{\min}\), set \(P_i = P_i^{\min}\); if \(P_i > P_i^{\max}\), set \(P_i = P_i^{\max}\)
  4. Compute \(P_{\text{total}} = \sum P_i\)
  5. If \(P_{\text{total}} < P_D\), increase \(\lambda\); if \(P_{\text{total}} > P_D\), decrease \(\lambda\)
  6. Repeat until \(\lvert P_{\text{total}} - P_D \rvert \leq \epsilon\)

When a generator hits a limit, it is fixed at that limit and its incremental cost no longer participates in the \(\lambda\) search. The optimality condition then requires that the binding generator’s marginal cost equal \(\lambda\) only if the limit were not active; otherwise \(\lambda\) lies above or below its actual marginal cost.

Kuhn-Tucker Conditions

The full KKT conditions for the ELD problem with limits are:

\[ \frac{dC_i}{dP_i} - \lambda + \mu_i^{\max} - \mu_i^{\min} = 0 \]\[ \mu_i^{\max} \geq 0, \quad \mu_i^{\max}(P_i - P_i^{\max}) = 0 \]\[ \mu_i^{\min} \geq 0, \quad \mu_i^{\min}(P_i^{\min} - P_i) = 0 \]

where \(\mu_i^{\max}\) and \(\mu_i^{\min}\) are non-negative multipliers for the upper and lower limits respectively. For an unconstrained generator: \(\mu_i^{\max} = \mu_i^{\min} = 0\) and \(dC_i/dP_i = \lambda\). For a generator at its maximum: \(\mu_i^{\max} = \lambda - dC_i/dP_i > 0\).

Inclusion of Transmission Losses: B-Coefficient Method

When losses are non-negligible, the power balance becomes \(\sum P_i = P_D + P_L\). Using the Kron reduction, total losses are approximated as:

\[ P_L = \sum_{i=1}^{N_g} \sum_{j=1}^{N_g} P_i B_{ij} P_j + \sum_{i=1}^{N_g} B_{0i} P_i + B_{00} \]

where \(B_{ij}\) are the loss B-coefficients derived from the network admittance and assumed operating point. The modified optimality condition becomes the coordination equation:

\[ \frac{dC_i}{dP_i} \cdot \frac{1}{1 - \partial P_L/\partial P_i} = \lambda \]

The penalty factor \(\beta_i = 1/(1 - \partial P_L/\partial P_i) \geq 1\) penalises generators electrically distant from load centres. The lambda iteration is repeated with penalty factors updated at each step.

5.3 Unit Commitment Problem

Problem Formulation

ELD assumes which units are online is given. The unit commitment (UC) problem determines the optimal subset of generators to start up and shut down over a planning horizon (typically 24–168 hours) to minimise total cost including start-up costs, subject to minimum up/down time constraints and spinning reserve requirements.

Decision variables: \(u_{i,t} \in \{0, 1\}\) (commitment status) and \(P_{i,t}\) (output) for each unit \(i\) and time period \(t\).

Objective (simplified):

\[ \min \sum_{t=1}^{T} \sum_{i=1}^{N_g} \left[C_i(P_{i,t}) u_{i,t} + SC_i \cdot z_{i,t}\right] \]

where \(SC_i\) is the start-up cost of unit \(i\) and \(z_{i,t} = \max(u_{i,t} - u_{i,t-1}, 0)\) is the start-up binary.

Key constraints:

  • Power balance: \(\sum_i P_{i,t} u_{i,t} = P_{D,t}\) for each \(t\)
  • Generator limits: \(P_i^{\min} u_{i,t} \leq P_{i,t} \leq P_i^{\max} u_{i,t}\)
  • Minimum up time: once committed, unit must remain on for at least \(MUT_i\) hours
  • Minimum down time: once de-committed, unit must remain off for at least \(MDT_i\) hours
  • Spinning reserve: \(\sum_i P_i^{\max} u_{i,t} \geq P_{D,t} + R_t\)

B.2 Unit Commitment Solution Methods

Priority List Method

The simplest heuristic orders generators by their average production cost (at some characteristic loading level). Generators are committed in merit order until total available capacity exceeds the load plus reserve requirement. This method ignores start-up costs and temporal constraints but is computationally fast.

Dynamic Programming

Dynamic programming (DP) solves UC exactly (within modelling approximations) by decomposing the horizon into stages. Let \(F(t, X)\) be the minimum cost from period 1 to \(t\) with unit commitment state \(X\) (a binary vector of length \(N_g\)) at period \(t\).

Recursion:

\[ F(t, X_t) = \min_{X_{t-1}} \left[F(t-1, X_{t-1}) + \text{Transition}(X_{t-1}, X_t) + \text{ELD}(P_{D,t}, X_t)\right] \]

where \(\text{Transition}(X_{t-1}, X_t)\) sums start-up costs for units transitioning from 0 to 1, and \(\text{ELD}(P_{D,t}, X_t)\) is the optimal dispatch cost for the active set \(X_t\) delivering \(P_{D,t}\).

The number of states grows as \(2^{N_g}\), so DP is tractable only for small systems. For large systems, Lagrangian relaxation, mixed-integer programming (MIP), or heuristic approaches are used.

Modern ISOs solve UC as a large-scale mixed-integer linear program (MILP), exploiting the linear programming relaxation with branch-and-cut algorithms. Commercial solvers (CPLEX, Gurobi) handle systems with hundreds of units and 48–168 time periods in minutes.

Chapter 6: Optimal Power Flow and Contingency Analysis

6.3 Contingency Analysis

N-1 Security Criterion

Power systems are designed and operated to remain secure under any single-element outage (generator or transmission line). This “N-1” criterion is checked by simulating each credible contingency and verifying that no thermal limits are exceeded and bus voltages remain within acceptable bounds.

For a system with \(N\) lines, the brute-force approach requires \(N\) full power flow solutions — computationally intensive for large systems. DC power flow with PTDFs offers a rapid screening method.

Contingency Analysis Using PTDFs

When line \(c\) is lost, the change in line flow on any remaining line \(l\) is estimated by the line outage distribution factor (LODF):

\[ \text{LODF}_{l,c} = \frac{\Delta P_l}{P_c^0} \]

where \(P_c^0\) is the pre-contingency flow on the outaged line. LODFs are derived from PTDFs:

\[ \text{LODF}_{l,c} = \frac{\text{PTDF}_{l,m} - \text{PTDF}_{l,n}}{1 - (\text{PTDF}_{c,m} - \text{PTDF}_{c,n})} \]

where line \(c\) connects buses \(m\) and \(n\). Post-contingency flow on line \(l\):

\[ P_l^{\text{post}} = P_l^0 + \text{LODF}_{l,c} \cdot P_c^0 \]

This linear approximation is very fast (no power flow solution needed) and enables real-time security assessment.

6.4 Optimal Power Flow

Formulation

The Optimal Power Flow (OPF) extends ELD by incorporating the full AC power flow equations as equality constraints, plus transmission thermal limits and voltage limits as inequality constraints.

General OPF:

\[ \min_{P_g, Q_g, V, \theta} \quad \sum_i C_i(P_{g,i}) \]

subject to:

  • AC power flow equations: \(P_i(\mathbf{V}, \boldsymbol{\theta}) = P_{g,i} - P_{d,i}\) and \(Q_i(\mathbf{V}, \boldsymbol{\theta}) = Q_{g,i} - Q_{d,i}\)
  • Voltage limits: \(V_i^{\min} \leq V_i \leq V_i^{\max}\)
  • Thermal limits: \(S_{l} \leq S_l^{\max}\)
  • Generator limits: \(P_{g,i}^{\min} \leq P_{g,i} \leq P_{g,i}^{\max}\), \(Q_{g,i}^{\min} \leq Q_{g,i} \leq Q_{g,i}^{\max}\)

Solution Methods (Appendix B.3–B.5)

Interior Point Methods (B.3): Transform the problem by introducing slack variables for inequality constraints and solving the KKT system with a barrier function. These methods have polynomial complexity and are the dominant approach in commercial OPF solvers.

Sequential Quadratic Programming (SQP, B.4): Solve a sequence of quadratic subproblems (QP) approximating the nonlinear OPF, updating the operating point each iteration. SQP has good local convergence properties.

Lagrangian Relaxation (B.5): Dualise complicating constraints (e.g., coupling constraints in multi-period problems) and solve the resulting Lagrangian subproblems, updating dual variables via subgradient or augmented Lagrangian methods.


Chapter 5.4–5.5: Electricity Markets

5.4 Fundamentals of Electricity Markets

Unique Characteristics of Electricity

Electricity markets differ from most commodity markets in fundamental ways:

  • Non-storability: Electricity cannot be economically stored at large scale; generation must equal load at every instant.
  • Network externalities: Power flows obey Kirchhoff’s laws, not contracts. An injection anywhere affects all line flows simultaneously.
  • Inelastic demand: Short-run price elasticity of demand is very low, creating the potential for extreme price spikes during scarcity.
  • Real-time balance requirement: System frequency deviates if generation and load are not balanced, requiring continuous control.

Supply and Demand in Electricity Markets

The supply curve in a spot electricity market is constructed by stacking generator offers in order of increasing marginal cost — the supply (or “merit order”) stack. The intersection of supply and demand determines the market clearing price (MCP) and quantity.

Short-run marginal cost (SRMC) for a thermal generator: fuel cost per MWh (heat rate × fuel price) plus variable operating cost. The system MCP equals the SRMC of the marginal (most expensive dispatched) unit.

Long-run equilibrium: In competitive markets, the MCP must average above the levelised cost of new entry (CONE) to incentivise investment. Capacity markets or energy price caps affect this equilibrium.

Market Power

Market power arises when one or more generators can profitably withhold capacity to raise prices above competitive levels. Mitigation measures include offer caps, market power screens, and structural remedies (divestiture).

5.5 Auction Mechanisms

Single-Period Energy Auction

In a single settlement period, generators submit price-quantity offers \((c_i, P_i^{\text{max}})\). The ISO ranks offers from lowest to highest (merit order) and dispatches in order until load \(P_D\) is satisfied. The clearing price \(\lambda\) is the offer price of the last (marginal) unit dispatched.

Uniform-price auction: All dispatched generators receive \(\lambda\). This incentivises truthful bidding if generators are price-takers (competitive). With market power, generators may shade offers above SRMC.

Pay-as-bid auction: Each dispatched generator receives its own offer price. This reduces the clearing price but distorts bidding incentives (generators withhold true cost information).

Day-Ahead Market

The DA market clears a 24-period auction the day before real time. Generators submit hourly energy offers, start-up costs, minimum generation levels, and ramp rate constraints. The ISO solves a UC + economic dispatch optimisation (a MILP) to minimise total cost. Cleared generators receive energy payments at the hourly LMP plus start-up uplift if applicable.

Locational Marginal Pricing in Auctions

When the network is congested (a transmission limit binds), the optimal dual variable of the line constraint contributes to LMPs differently at each bus:

\[ \text{LMP}_k = \lambda - \sum_l \mu_l \cdot \text{PTDF}_{l,k} \]

where \(\mu_l \geq 0\) is the shadow price (congestion rent) of line \(l\), and \(\text{PTDF}_{l,k}\) is the sensitivity of line \(l\) to injection at bus \(k\) relative to the slack. This disaggregation ensures that generators and loads facing their true marginal cost, including network costs, have efficient dispatch incentives.

Ancillary Service Markets

Alongside energy, ISOs procure ancillary services co-optimised with energy:

  • Spinning reserve: Online capacity that can respond within 10 minutes; compensates for unexpected generation loss.
  • Non-spinning (supplemental) reserve: Offline capacity startable within 10–30 minutes.
  • Regulation (AGC): Fast-responding resources (seconds to minutes) that follow automatic generation control signals to balance real-time deviations.

Co-optimisation ensures that a unit providing reserve is compensated for the opportunity cost of not generating energy.


Chapter 7–8: Fault Analysis

7.2 Symmetrical Components

Motivation

Unbalanced faults (single-line-to-ground, line-to-line, double-line-to-ground) produce unbalanced three-phase voltages and currents even in an otherwise balanced system. Symmetrical components decompose unbalanced phasors into three balanced sets, enabling separate analysis of each sequence network.

Transformation Matrix

For a set of three-phase phasors \(\mathbf{F}_{abc} = [F_a, F_b, F_c]^T\), the symmetrical components \(\mathbf{F}_{012} = [F_0, F_1, F_2]^T\) (zero, positive, negative sequence) are related by:

\[ \mathbf{F}_{abc} = \mathbf{A} \mathbf{F}_{012} \]

where

\[ \mathbf{A} = \begin{pmatrix} 1 & 1 & 1 \\ 1 & a^2 & a \\ 1 & a & a^2 \end{pmatrix}, \qquad a = e^{j2\pi/3} = -\frac{1}{2} + j\frac{\sqrt{3}}{2} \]

The inverse transformation:

\[ \mathbf{F}_{012} = \mathbf{A}^{-1} \mathbf{F}_{abc} = \frac{1}{3}\begin{pmatrix} 1 & 1 & 1 \\ 1 & a & a^2 \\ 1 & a^2 & a \end{pmatrix} \mathbf{F}_{abc} \]

Properties of \(a\): \(a^3 = 1\), \(1 + a + a^2 = 0\), \(\lvert a \rvert = 1\).

Physical interpretation:

  • Positive sequence (\(F_1\)): balanced set with normal \(a\text{-}b\text{-}c\) phase sequence.
  • Negative sequence (\(F_2\)): balanced set with reverse phase sequence (\(a\text{-}c\text{-}b\)).
  • Zero sequence (\(F_0\)): all three phasors in phase with each other.

8.2 Sequence Network Models

For a balanced three-phase network, positive, negative, and zero sequence networks are decoupled. Each network is a single-phase equivalent:

Positive sequence network: Contains all the machines and lines at their positive-sequence impedances. This is the normal (prefault) network.

Negative sequence network: Machines have negative-sequence impedance \(Z_2 \approx Z_d'' = X_d''\) (subtransient). Transmission lines and transformers have the same impedance as in the positive sequence (they are balanced devices). Rotating machines have \(Z_2 \neq Z_1\).

Zero sequence network: Network topology changes. Lines carry zero-sequence current only if there is a ground return path. Transformer connections determine zero-sequence current flow:

  • Grounded wye / grounded wye: zero sequence passes through both windings.
  • Grounded wye / delta: zero sequence circulates in the delta, not transmitted.
  • Ungrounded wye: zero sequence cannot flow.
The Thevenin equivalent of each sequence network seen from the faulted bus \(f\) is characterised by:
  • Positive sequence: open-circuit voltage \(V_f\) (prefault voltage), Thevenin impedance \(Z_{1,ff}\)
  • Negative sequence: zero source voltage, Thevenin impedance \(Z_{2,ff}\)
  • Zero sequence: zero source voltage, Thevenin impedance \(Z_{0,ff}\)

8.4 Balanced (Three-Phase) Faults

A three-phase (3Φ) fault at bus \(f\) involves all three phases simultaneously to ground with fault impedance \(Z_f\). Only the positive sequence network is involved:

\[ I_1 = \frac{V_f}{Z_{1,ff} + Z_f} \]

The fault current in phase \(a\) is:

\[ I_{fa} = I_1 = \frac{V_f}{Z_{1,ff} + Z_f} \]

(since \(I_2 = I_0 = 0\) for a balanced fault). The fault current magnitude is \(\lvert I_{fa} \rvert\) and sets the interrupting duty for circuit breakers.

Post-fault bus voltages are found from:

\[ V_k^{\text{post}} = V_k^{\text{pre}} - Z_{k,f} I_{fa} \]

where \(Z_{k,f}\) is the element of the Z-bus relating voltage at bus \(k\) to current at bus \(f\).

8.5 Unbalanced Faults

Single Line-to-Ground (LG) Fault

For a phase-\(a\) to ground fault with fault impedance \(Z_f\):

Boundary conditions: \(I_b = I_c = 0\), \(V_a = Z_f I_a\).

These give \(I_0 = I_1 = I_2 = I_a/3\) and the connection of sequence networks in series:

\[ I_1 = \frac{V_f}{Z_{1,ff} + Z_{2,ff} + Z_{0,ff} + 3Z_f} \]\[ I_a = 3 I_1 = \frac{3 V_f}{Z_{1,ff} + Z_{2,ff} + Z_{0,ff} + 3Z_f} \]

Line-to-Line (LL) Fault

For a phase-\(b\) to phase-\(c\) fault with \(Z_f\):

Boundary conditions: \(I_a = 0\), \(I_b = -I_c\), \(V_b - V_c = Z_f I_b\).

These give \(I_0 = 0\) and the positive and negative sequence networks in parallel:

\[ I_1 = \frac{V_f}{Z_{1,ff} + Z_{2,ff} + Z_f}, \quad I_2 = -I_1 \]\[ I_b = (a^2 - a) I_1 = -j\sqrt{3} I_1 \]

Double Line-to-Ground (LLG) Fault

For phases \(b\) and \(c\) to ground with \(Z_f\) in each phase and \(Z_g\) in common ground:

Boundary conditions: \(I_a = 0\), \(V_b = Z_f I_b + Z_g I_g\), \(V_c = Z_f I_c + Z_g I_g\).

Sequence networks: positive sequence drives the current; negative and zero sequences are connected in parallel across the positive sequence:

\[ I_1 = \frac{V_f}{Z_{1,ff} + \dfrac{Z_{2,ff}(Z_{0,ff} + 3Z_g + 3Z_f)}{Z_{2,ff} + Z_{0,ff} + 3Z_g + 3Z_f}} \]
The magnitudes of fault currents generally satisfy: LLG > LL > LG (for solidly grounded systems with \(Z_{0,ff} < Z_{1,ff}\)). However, if \(Z_{0,ff} > Z_{1,ff}\), LG faults may be smaller than LL faults. The practitioner must check for each specific system.

Chapter 9–10: Power System Stability and Control

10.1 Stability Definitions

Power system stability is the ability of a power system, for a given initial operating condition, to regain a state of operating equilibrium after being subjected to a physical disturbance, with most system variables bounded and the system remaining intact.

The IEEE/CIGRE classification (2004, updated 2021) categorises stability by the physical phenomenon:

  • Rotor angle stability: Ability to maintain synchronism under torque imbalances.
    • Transient (large-disturbance) angle stability: response to severe contingencies (faults, large load steps).
    • Small-signal angle stability: damping of oscillatory modes under small disturbances.
  • Voltage stability: Ability to maintain steady voltages at all buses following disturbances; failures lead to voltage collapse.
  • Frequency stability: Ability to maintain frequency in the nominal range following a large imbalance between generation and load.

10.2 Angle Stability — Swing Equation and Equal Area Criterion

Swing Equation

The rotor dynamics of a synchronous generator are governed by Newton’s second law applied to rotational motion:

\[ \frac{2H}{\omega_s} \frac{d^2\delta}{dt^2} = P_m - P_e(\delta) - D\frac{d\delta}{dt} \]

where:

  • \(H\) = inertia constant (seconds, stored kinetic energy at rated speed / rated MVA)
  • \(\omega_s = 2\pi f_0\) = synchronous angular frequency (rad/s)
  • \(\delta\) = rotor angle (rad), measuring electrical angle between machine’s EMF and an infinite bus reference
  • \(P_m\) = mechanical power input (pu)
  • \(P_e(\delta)\) = electrical power output (pu)
  • \(D\) = damping coefficient

The electrical power from a machine connected to an infinite bus through total reactance \(X\) is:

\[ P_e(\delta) = \frac{E' V_\infty}{X} \sin\delta = P_{\max} \sin\delta \]

At the pre-disturbance equilibrium: \(P_m = P_e(\delta_0) = P_{\max}\sin\delta_0\).

Equal Area Criterion (EAC)

For a lossless single-machine infinite-bus (SMIB) system without damping, the swing equation has the form of a conservative system. Stability can be assessed graphically without solving the differential equation.

Accelerating area (A1): Area between the \(P_m\) line and the pre-fault \(P_e\) curve from \(\delta_0\) to \(\delta_c\) (fault clearing angle). This represents kinetic energy gained by the rotor during the fault.

Decelerating area (A2): Area between the post-fault \(P_e\) curve and the \(P_m\) line from \(\delta_c\) to the maximum allowable angle \(\delta_{\max}\) (the supplement of the post-fault unstable equilibrium angle). This represents maximum kinetic energy that can be absorbed.

Equal Area Criterion: The system is transiently stable if and only if a decelerating area \(A_2 \geq A_1\) can be found, i.e., the rotor can decelerate and return to a stable equilibrium before reaching the unstable equilibrium point. Mathematically: \[ \int_{\delta_0}^{\delta_c} (P_m - P_e^{\text{fault}}(\delta))\,d\delta \leq \int_{\delta_c}^{\delta_{\max}} (P_e^{\text{post}}(\delta) - P_m)\,d\delta \]

Critical Clearing Time (CCT)

The CCT is the maximum fault duration for which the system remains transiently stable. It is found by equating \(A_1 = A_2\) and solving for the critical clearing angle \(\delta_{\text{cc}}\):

\[ \cos\delta_{\text{cc}} = \frac{P_m(\delta_{\max} - \delta_0) + P_{\max}^{\text{post}}\cos\delta_{\max} - P_{\max}^{\text{fault}}\cos\delta_0}{P_{\max}^{\text{post}} - P_{\max}^{\text{fault}}} \]

where \(\delta_{\max} = \pi - \delta_0^{\text{post}}\) and \(\delta_0^{\text{post}} = \arcsin(P_m / P_{\max}^{\text{post}})\).

The CCT in seconds is obtained from the swing equation by numerically (or for zero-fault-power, analytically) integrating from \(\delta_0\) to \(\delta_{\text{cc}}\). For the special case \(P_e^{\text{fault}} = 0\):

\[ t_{\text{cc}} = \sqrt{\frac{4H(\delta_{\text{cc}} - \delta_0)}{\omega_s P_m}} \]
Example — CCT calculation: A generator with \(H = 5\,\text{s}\), \(P_m = 0.8\,\text{pu}\), \(P_{\max}^{\text{pre}} = 2.0\,\text{pu}\), \(P_{\max}^{\text{fault}} = 0\), \(P_{\max}^{\text{post}} = 1.5\,\text{pu}\), \(f_0 = 60\,\text{Hz}\).

Pre-fault equilibrium: \(\delta_0 = \arcsin(0.8/2.0) = 23.58°\).

Post-fault equilibrium: \(\delta_0^{\text{post}} = \arcsin(0.8/1.5) = 32.23°\).

\(\delta_{\max} = 180° - 32.23° = 147.77°\).

\[ \cos\delta_{\text{cc}} = \frac{0.8(147.77° - 23.58°)\cdot\frac{\pi}{180} + 1.5\cos(147.77°) - 0}{1.5 - 0} \]

Numerically: \(\cos\delta_{\text{cc}} = [0.8 \times 2.165 - 1.5 \times 0.844] / 1.5 = [1.732 - 1.266]/1.5 = 0.311\), so \(\delta_{\text{cc}} = 71.9°\).

\[ t_{\text{cc}} = \sqrt{\frac{4 \times 5 \times (71.9° - 23.58°) \times \pi/180}{2\pi \times 60 \times 0.8}} = \sqrt{\frac{4 \times 5 \times 0.843}{301.6}} = \sqrt{0.0559} \approx 0.236\,\text{s} \]

Small-Signal Stability

Linearise the swing equation around the equilibrium \(\delta_0\):

\[ \frac{2H}{\omega_s}\Delta\ddot{\delta} + D\Delta\dot{\delta} + K_s \Delta\delta = 0 \]

where \(K_s = \partial P_e/\partial\delta\big|_{\delta_0} = P_{\max}\cos\delta_0\) is the synchronising torque coefficient. The characteristic roots are:

\[ s_{1,2} = -\frac{D\omega_s}{4H} \pm \sqrt{\left(\frac{D\omega_s}{4H}\right)^2 - \frac{K_s\omega_s}{2H}} \]

Oscillatory stability requires \(K_s > 0\) (met when \(\delta_0 < 90°\)) and \(D > 0\). Power system stabilisers (PSS) add damping torque to suppress inter-area oscillations.

10.3 Voltage Stability

PV and QV Curves

Voltage stability is analysed via the PV curve (nose curve) and QV curve, which characterise how bus voltages respond to increasing load.

PV curve: Plot of bus voltage magnitude \(\lvert V \rvert\) versus total active load power \(P\) (with constant power factor). As \(P\) increases, \(\lvert V \rvert\) decreases. The “nose point” (maximum loadability point or saddle-node bifurcation) represents the voltage collapse limit \(P_{\max}\). Beyond this point, no operating solution exists.

For a simple 2-bus system with source voltage \(E\), line reactance \(X\), and load with constant power factor \(\cos\phi\):

\[ P_{\max} = \frac{E^2}{2X(1 + \sin\phi)} \]

QV curve: Plot of bus reactive power injection \(Q\) versus voltage \(V\) at the bus, holding active power constant. The minimum of the QV curve indicates the critical reactive power deficit:

  • If the minimum \(Q\) is below zero, the bus can support the load without reactive compensation.
  • If the minimum \(Q > 0\), reactive compensation is needed to maintain stability.
The voltage stability margin is the distance in MW (or in reactive power, Mvar) between the current operating point and the nose point of the PV (or QV) curve. It quantifies how much additional load can be served before voltage collapse.

Continuation Power Flow

The continuation power flow (CPF) traces the PV curve through the nose point using predictor-corrector steps parameterised by a load multiplier \(\lambda\), avoiding the ill-conditioning of the standard power flow Jacobian near the bifurcation.

Reactive Power Compensation

Voltage support is provided by:

  • Shunt capacitors: Fixed or switched, provide reactive power proportional to \(V^2\).
  • Static VAr compensators (SVC): Thyristor-controlled, provide variable reactive support.
  • STATCOM: VSC-based, provides reactive support independent of voltage.
  • Synchronous condensers: Rotating machines, provide inertia and continuous reactive control.
  • Under-load tap-changing (ULTC) transformers: Adjust voltage by changing turns ratio under load.

9.4 Voltage Control

Voltage control is hierarchical:

  • Primary voltage control: Automatic voltage regulators (AVR) on generators maintain terminal voltage by adjusting field current. Response time: milliseconds to seconds.
  • Secondary voltage control: Coordinates multiple generators in a region to maintain pilot bus voltages. Response time: seconds to minutes.
  • Tertiary voltage control (OPF-based): Optimal scheduling of reactive resources. Response time: minutes.

The droop characteristic of a generator’s AVR relates terminal voltage deviation to reactive power output:

\[ Q = Q_0 - \frac{1}{R_Q}(V - V_{\text{ref}}) \]

where \(R_Q\) is the voltage regulation droop.

10.4 Frequency Stability

Frequency stability concerns the system’s ability to maintain nominal frequency following a large active power imbalance. A large generator trip creates a sudden deficit \(\Delta P\) that the inertia of all synchronous machines initially absorbs:

\[ \frac{2H_{\text{sys}}}{\omega_s} \frac{df}{dt} = P_m - P_e = -\Delta P \]

where \(H_{\text{sys}} = \sum_i H_i S_i / S_{\text{sys}}\) is the system inertia. The initial rate of change of frequency (RoCoF):

\[ \text{RoCoF} = \frac{df}{dt}\bigg|_{t=0} = -\frac{\Delta P \cdot f_0}{2 H_{\text{sys}} S_{\text{sys}}} \]

A high RoCoF — common in low-inertia systems with high renewable penetration — can trigger protection relays before governor responses activate.

The frequency nadir (minimum frequency) is reached when governor responses arrest the frequency decline and matches the loss:

\[ \Delta P = D \Delta f_{\text{nadir}} + \sum_i \Delta P_{m,i} \]

where \(D\) is the system load damping coefficient and \(\Delta P_{m,i}\) are governor responses.

9.3 Frequency Control

Primary Frequency Control (Governor Droop)

Speed governors on generators respond to frequency deviation by adjusting mechanical power. The droop characteristic is:

\[ \frac{\Delta P_m}{P_{\text{rated}}} = -\frac{1}{R} \frac{\Delta f}{f_0} \]

where \(R\) is the droop setting (typically 4–5%). A 5% droop means a 5% change in frequency causes a 100% change in generator output.

Static frequency characteristic (SFC): At steady state after a disturbance, the frequency settles at a new value determined by the aggregate droop of all contributing generators and the frequency sensitivity of loads:

\[ \Delta f_{\text{ss}} = -\frac{\Delta P_L}{D_L + \sum_i (1/R_i)} \]

where \(D_L\) is the load-frequency sensitivity coefficient and \(\sum_i(1/R_i)\) is the aggregate generator stiffness.

Example — Droop and frequency response: A 1000 MW system has two generators (G1: 600 MW rated, \(R_1 = 0.05\); G2: 400 MW rated, \(R_2 = 0.05\)) and load damping \(D_L = 1\,\text{pu MW/Hz}\) on a 1000 MVA base, \(f_0 = 60\,\text{Hz}\). A 100 MW (0.1 pu) load step occurs.

Generator stiffness: \(1/R_1 = 0.6/0.05 = 12\,\text{pu MW/Hz}\) on own base \(= 12 \times 0.6 = 7.2\,\text{pu MW/Hz}\) on system base; similarly for G2: \(8\,\text{pu MW/Hz}\).

\[ \Delta f_{\text{ss}} = -\frac{0.1}{1 + 7.2 + 8} \approx -\frac{0.1}{16.2} \approx -0.00617\,\text{pu} = -0.37\,\text{Hz} \]

The frequency settles at approximately 59.63 Hz.

Secondary Frequency Control (Automatic Generation Control — AGC)

Primary droop control arrests frequency deviation but leaves a steady-state error \(\Delta f_{\text{ss}} \neq 0\). Automatic generation control (AGC) eliminates this error by adjusting the reference set points of governors based on the Area Control Error (ACE):

\[ \text{ACE} = \Delta P_{\text{tie}} + B \Delta f \]

where \(\Delta P_{\text{tie}}\) is the deviation of actual tie-line power from scheduled value and \(B\) is the frequency bias coefficient. The AGC integrator drives ACE → 0, restoring both frequency and inter-area power exchange to scheduled values.

Tertiary Frequency Control (Economic Dispatch)

The optimal re-dispatch of generation every 5–15 minutes restores reserve margins depleted by primary and secondary control, preparing the system for the next potential disturbance.

In power systems with high penetration of inverter-based resources (wind, solar), total system inertia decreases because power-electronic converters do not inherently provide inertial response. Synthetic (virtual) inertia emulation — using the controls of wind turbines or battery storage to mimic inertial response — is an active area of research and deployment to maintain adequate frequency stability margins.

Summary of Key Formulas

Per-Unit System

\[ z_{\text{new}} = z_{\text{old}} \cdot \frac{S_{\text{base,new}}}{S_{\text{base,old}}} \cdot \left(\frac{V_{\text{base,old}}}{V_{\text{base,new}}}\right)^2 \]

Power Flow Equations

\[ P_i = V_i \sum_k V_k (G_{ik}\cos\theta_{ik} + B_{ik}\sin\theta_{ik}) \]\[ Q_i = V_i \sum_k V_k (G_{ik}\sin\theta_{ik} - B_{ik}\cos\theta_{ik}) \]

DC Power Flow and PTDF

\[ \mathbf{P}_{\text{inj}} = \mathbf{B}_{\text{DC}} \boldsymbol{\theta}, \qquad P_{ik} = \frac{\theta_i - \theta_k}{x_{ik}} \]\[ \text{PTDF}_{l,m} = \frac{1}{x_l}(\mathbf{e}_i - \mathbf{e}_k)^T \mathbf{B}_{\text{DC}}^{-1} \mathbf{e}_m \]

Economic Dispatch

\[ \frac{dC_i}{dP_i} = \lambda \quad \forall i \text{ (unconstrained, lossless)}, \qquad P_i = \frac{\lambda - b_i}{2c_i} \]

Fault Currents

Fault TypeSequence currentsTotal fault current
3-phase\(I_1 = V_f/(Z_1 + Z_f)\)\(I_a = I_1\)
LG\(I_0 = I_1 = I_2 = V_f/(Z_1+Z_2+Z_0+3Z_f)\)\(I_a = 3I_1\)
LL\(I_1 = V_f/(Z_1+Z_2+Z_f)\), \(I_2 = -I_1\)\(\lvert I_b \rvert = \sqrt{3}\lvert I_1\rvert\)

Swing Equation and CCT

\[ \frac{2H}{\omega_s}\ddot{\delta} = P_m - P_{\max}\sin\delta, \qquad t_{\text{cc}} = \sqrt{\frac{4H(\delta_{\text{cc}} - \delta_0)}{\omega_s P_m}} \text{ (zero-fault power)} \]

Frequency Response

\[ \text{RoCoF} = -\frac{\Delta P \cdot f_0}{2H_{\text{sys}}}, \qquad \Delta f_{\text{ss}} = -\frac{\Delta P_L}{D_L + \sum_i (1/R_i)} \]
Back to top