After completing this lesson, the student should be able to:

• derive an aggregate turbine-governor-frequency model for one control area;
• form the area-control error for an interconnected system;
• construct a matched, linearly parameterized adaptive error model;
• derive a Lyapunov update law for frequency and voltage channels;
• explain the role of projection, leakage, normalization, saturation, and rate limits;
• distinguish conceptual adaptive regulation from protection, dispatch, and transient-stability studies.

The lesson does not replace nonlinear differential-algebraic power-system models, electromagnetic-transient models, relay studies, or operator procedures. It isolates the adaptive-control reasoning that can be embedded inside a carefully validated secondary-control or excitation-control layer.

Frequency is mainly coupled to the active-power balance, whereas voltage is mainly coupled to reactive-power balance and excitation. This separation is approximate, especially in weak grids, converter-dominated systems, and stressed transmission networks.

flowchart TD
  A["Load, renewable, and topology changes"] --> B["Primary frequency response: \ndroop and fast reserves"]
  B --> C["Secondary frequency control: \nACE restoration"]
  C --> D["Economic dispatch and \nsupervisory optimization"]
  A --> E["Excitation and converter \nvoltage control"]
  E --> F["Secondary voltage and \nreactive-power coordination"]
  F --> D
  C --> G["Adaptive frequency \nparameters or gains"]
  F --> H["Adaptive voltage \nparameters or gains"]
  G --> I["Limits, projection, \nmonitoring, fallback"]
  H --> I
        

Adaptation should normally modify a bounded set of controller coefficients or uncertainty estimates. It should not bypass protection, reserve limits, turbine constraints, field-current limits, converter current limits, or supervisory authorization.

Let \( \Delta f \) be the per-unit frequency deviation, \( \Delta P_m \) the mechanical-power increment, \( \Delta P_v \) the governor-valve command, and \( \Delta P_L \) the active-load disturbance. A standard reduced model is

\[ M\dot{\Delta f}=\Delta P_m-\Delta P_L-D\Delta f-\Delta P_{\mathrm{tie}}, \]

\[ T_t\dot{\Delta P_m}=-\Delta P_m+\Delta P_v,\qquad T_g\dot{\Delta P_v}=-\Delta P_v-\frac{1}{R}\Delta f+u_f . \]

Here \( M \) is the aggregate inertia coefficient, \( D \) is load-frequency damping, \( R \) is speed droop, and \( T_t,T_g \) are turbine and governor time constants. The secondary command \( u_f \) changes the governor reference so that primary droop error is removed.

In an isolated area with constant load and only primary droop, the quasi-steady balance gives

\[ 0=-\Delta P_L-\left(D+\frac{1}{R}\right)\Delta f_{\mathrm{ss}}, \qquad \Delta f_{\mathrm{ss}} =-\frac{\Delta P_L}{D+1/R}. \]

Thus primary control arrests frequency but generally does not restore it exactly. Secondary integral or adaptive action is required for restoration.

For two coherent areas, the incremental tie-line power is commonly represented by

\[ \dot{\Delta P}_{12}=2\pi T_{12}\left(\Delta f_1-\Delta f_2\right), \]

where \( T_{12} \) is the synchronizing coefficient. With the export-positive convention,

\[ \mathrm{ACE}_1=B_1\Delta f_1+\Delta P_{12},\qquad \mathrm{ACE}_2=B_2\Delta f_2-\Delta P_{12}. \]

The frequency-bias coefficients \( B_i \) allocate responsibility between frequency restoration and scheduled interchange. If the secondary loop drives both ACE signals to zero and \( B_1+B_2\neq 0 \), then

\[ \left(B_1+B_2\right)\Delta f_{\mathrm{ss}}=0 \quad\Longrightarrow\quad \Delta f_{\mathrm{ss}}=0,\qquad \Delta P_{12,\mathrm{ss}}=0. \]

For teaching, the adaptive error variable may be chosen as \( e_f=\mathrm{ACE} \), a filtered ACE, or a frequency tracking error. Real implementations must preserve sign conventions, metering quality, communication timing, and balancing-authority rules.

After stable inner turbine-governor dynamics have been accounted for, a reduced secondary-loop model can be written as

\[ \dot e_f=\boldsymbol{\theta}_f^{*\mathsf T}\boldsymbol{\phi}_f +b_f u_f+d_f,\qquad b_f>0, \]

with regressor \( \boldsymbol{\phi}_f=[e_f,\xi_f,\Delta P_L]^{\mathsf T} \) and \( \dot{\xi}_f=e_f \). Unknown coefficients can represent aggregate damping, effective reserve response, operating-point dependence, or a local disturbance map. Assuming the dominant uncertainty is matched, use

\[ u_f=\frac{1}{b_f} \left(-\hat{\boldsymbol{\theta}}_f^{\mathsf T}\boldsymbol{\phi}_f -k_f e_f\right),\qquad k_f>0. \]

Define \( \tilde{\boldsymbol{\theta}}_f= \boldsymbol{\theta}_f^*-\hat{\boldsymbol{\theta}}_f \). With \( d_f=0 \), the closed-loop error equation is

\[ \dot e_f=-k_f e_f+ \tilde{\boldsymbol{\theta}}_f^{\mathsf T}\boldsymbol{\phi}_f. \]

Choose the update law and Lyapunov function

\[ \dot{\hat{\boldsymbol{\theta}}}_f =\boldsymbol{\Gamma}_f\boldsymbol{\phi}_f e_f,\qquad V_f=\frac{1}{2}e_f^2+ \frac{1}{2}\tilde{\boldsymbol{\theta}}_f^{\mathsf T} \boldsymbol{\Gamma}_f^{-1} \tilde{\boldsymbol{\theta}}_f , \]

where \( \boldsymbol{\Gamma}_f=\boldsymbol{\Gamma}_f^{\mathsf T} \succ 0 \). Since \( \dot{\tilde{\boldsymbol{\theta}}}_f =-\dot{\hat{\boldsymbol{\theta}}}_f \),

\[ \begin{aligned} \dot V_f &=e_f\dot e_f+ \tilde{\boldsymbol{\theta}}_f^{\mathsf T} \boldsymbol{\Gamma}_f^{-1} \dot{\tilde{\boldsymbol{\theta}}}_f\\ &=-k_f e_f^2+ e_f\tilde{\boldsymbol{\theta}}_f^{\mathsf T}\boldsymbol{\phi}_f -\tilde{\boldsymbol{\theta}}_f^{\mathsf T} \boldsymbol{\phi}_f e_f\\ &=-k_f e_f^2. \end{aligned} \]

Therefore all Lyapunov states are bounded and \( e_f\in L_2\cap L_\infty \). Under bounded \( \dot e_f \), Barbalat's lemma gives \( e_f(t)\to 0 \). Parameter convergence requires persistent excitation and is not implied by frequency restoration alone.

Around an operating point, a reduced voltage channel may be represented by

\[ T_v\dot{\Delta V} =-a_v\Delta V+b_v\Delta E_f-c_q\Delta Q_L, \]

where \( \Delta E_f \) is the exciter or converter voltage command and \( \Delta Q_L \) is a reactive-load disturbance. Define \( e_v=\Delta V-\Delta V_m \) and a linearly parameterized model

\[ \dot e_v= \boldsymbol{\theta}_v^{*\mathsf T}\boldsymbol{\phi}_v+ b_vu_v-\dot{\Delta V}_m+d_v, \quad \boldsymbol{\phi}_v=[e_v,\xi_v,\Delta Q_L]^{\mathsf T}, \quad \dot\xi_v=e_v. \]

The certainty-equivalent voltage command is

\[ u_v=\frac{1}{b_v} \left( -\hat{\boldsymbol{\theta}}_v^{\mathsf T}\boldsymbol{\phi}_v -k_v e_v+\dot{\Delta V}_m \right), \]

\[ \dot{\hat{\boldsymbol{\theta}}}_v= \boldsymbol{\Gamma}_v\boldsymbol{\phi}_v e_v. \]

The same cancellation argument gives \( \dot V_v=-k_v e_v^2 \) in the ideal matched case. In a physical exciter or inverter, \( u_v \) must be limited by field voltage, field current, stator current, converter current, reactive-power capability, and anti-windup logic.

flowchart TD
  A["Measure frequency, voltage, tie-line, P, Q"] --> B["Filter and validate signals"]
  B --> C["Build regressors and errors"]
  C --> D["Compute nominal plus adaptive command"]
  D --> E["Apply saturation and rate limits"]
  E --> F["Plant or grid-facing actuator"]
  C --> G["Normalized update law"]
  G --> H["Projection or leakage"]
  H --> I["Parameter estimates"]
  I --> D
  E --> J["Constraint and anomaly monitor"]
  J --> K["Freeze adaptation or revert to certified fallback"]
        

A normalized projected law has the conceptual form

\[ \dot{\hat{\boldsymbol{\theta}}} =\operatorname{Proj}_{\Omega} \left( \hat{\boldsymbol{\theta}}, \frac{\boldsymbol{\Gamma}\boldsymbol{\phi}e} {1+\nu\boldsymbol{\phi}^{\mathsf T}\boldsymbol{\phi}} \right), \qquad \nu>0. \]

Leakage can instead be introduced by

\[ \dot{\hat{\boldsymbol{\theta}}} =\boldsymbol{\Gamma}\boldsymbol{\phi}e -\sigma\boldsymbol{\Gamma} \left(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0\right), \qquad \sigma>0. \]

Projection enforces hard parameter bounds; leakage discourages drift but generally changes asymptotic convergence into uniform ultimate boundedness. Dead zones can stop adaptation when the regulation error is comparable with sensor noise or unavoidable ripple.

If the reduced frequency and voltage channels are weakly coupled, a composite Lyapunov function can be proposed:

\[ V=V_f+\rho V_v,\qquad \rho>0. \]

Let coupling terms satisfy \( |g_f(e_f,e_v)|\leq c_{ff}|e_f|+c_{fv}|e_v| \) and \( |g_v(e_f,e_v)|\leq c_{vf}|e_f|+c_{vv}|e_v| \). Young's inequality gives cross-term bounds such as

\[ |e_f e_v|\leq \frac{\epsilon}{2}e_f^2+ \frac{1}{2\epsilon}e_v^2,\qquad \epsilon>0. \]

The nominal gains \( k_f,k_v \) can dominate sufficiently small coupling bounds. This is a local reduced-model argument, not a proof of transient stability for a full network differential-algebraic model. Strong active-reactive coupling requires a MIMO adaptive design rather than two independently tuned scalar channels.

Chapter27_Lesson3.py

"""
Chapter27_Lesson3.py
Conceptual adaptive frequency and voltage regulation using matched,
linearly parameterized reduced-order error models.

Dependencies:
    numpy
    matplotlib

Run:
    python Chapter27_Lesson3.py
"""

from __future__ import annotations

import csv
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt


def project(theta: np.ndarray, lower: np.ndarray, upper: np.ndarray) -> np.ndarray:
    """Componentwise projection used to keep adaptive estimates bounded."""
    return np.minimum(np.maximum(theta, lower), upper)


def simulate(
    final_time: float = 80.0,
    dt: float = 0.01,
) -> dict[str, np.ndarray]:
    steps = int(final_time / dt) + 1
    time = np.linspace(0.0, final_time, steps)

    # Reduced matched models:
    # e_dot = theta_true^T phi + b*u + d_unmatched
    theta_f_true = np.array([-0.55, -0.20, 1.00])
    theta_v_true = np.array([-0.75, -0.28, 0.85])
    b_f, b_v = 1.20, 1.05

    # Stabilizing feedback and adaptation gains.
    k_f, k_v = 1.80, 2.10
    gamma_f = np.diag([2.0, 0.8, 3.0])
    gamma_v = np.diag([1.6, 0.7, 2.4])

    lower = np.array([-3.0, -3.0, -3.0])
    upper = np.array([3.0, 3.0, 3.0])

    e_f = np.zeros(steps)
    xi_f = np.zeros(steps)
    e_v = np.zeros(steps)
    xi_v = np.zeros(steps)
    u_f = np.zeros(steps)
    u_v = np.zeros(steps)
    load_p = np.zeros(steps)
    load_q = np.zeros(steps)
    theta_f_hat = np.zeros((steps, 3))
    theta_v_hat = np.zeros((steps, 3))

    # Nonzero initial regulation errors.
    e_f[0] = 0.025
    e_v[0] = -0.030
    theta_f_hat[0] = np.array([-0.10, -0.05, 0.20])
    theta_v_hat[0] = np.array([-0.15, -0.05, 0.10])

    for k in range(steps - 1):
        t = time[k]

        # Active- and reactive-power disturbances.
        load_p[k] = 0.0 if t < 8.0 else (0.035 if t < 42.0 else 0.055)
        load_q[k] = 0.0 if t < 15.0 else (-0.025 if t < 50.0 else 0.040)

        phi_f = np.array([e_f[k], xi_f[k], load_p[k]])
        phi_v = np.array([e_v[k], xi_v[k], load_q[k]])

        # Certainty-equivalent adaptive controls.
        u_f[k] = (-theta_f_hat[k] @ phi_f - k_f * e_f[k]) / b_f
        u_v[k] = (-theta_v_hat[k] @ phi_v - k_v * e_v[k]) / b_v

        # Small unmatched periodic terms illustrate why exact cancellation is not promised.
        d_f = 0.0015 * np.sin(0.35 * t)
        d_v = 0.0010 * np.cos(0.27 * t)

        e_f_dot = theta_f_true @ phi_f + b_f * u_f[k] + d_f
        e_v_dot = theta_v_true @ phi_v + b_v * u_v[k] + d_v
        xi_f_dot = e_f[k]
        xi_v_dot = e_v[k]

        theta_f_dot = gamma_f @ (phi_f * e_f[k])
        theta_v_dot = gamma_v @ (phi_v * e_v[k])

        e_f[k + 1] = e_f[k] + dt * e_f_dot
        e_v[k + 1] = e_v[k] + dt * e_v_dot
        xi_f[k + 1] = xi_f[k] + dt * xi_f_dot
        xi_v[k + 1] = xi_v[k] + dt * xi_v_dot
        theta_f_hat[k + 1] = project(
            theta_f_hat[k] + dt * theta_f_dot, lower, upper
        )
        theta_v_hat[k + 1] = project(
            theta_v_hat[k] + dt * theta_v_dot, lower, upper
        )

    load_p[-1] = load_p[-2]
    load_q[-1] = load_q[-2]
    u_f[-1] = u_f[-2]
    u_v[-1] = u_v[-2]

    return {
        "time": time,
        "frequency_error": e_f,
        "voltage_error": e_v,
        "frequency_control": u_f,
        "voltage_control": u_v,
        "active_load": load_p,
        "reactive_load": load_q,
        "theta_f_hat": theta_f_hat,
        "theta_v_hat": theta_v_hat,
    }


def save_csv(result: dict[str, np.ndarray], path: Path) -> None:
    with path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(
            [
                "time_s",
                "frequency_error_pu",
                "voltage_error_pu",
                "frequency_control_pu",
                "voltage_control_pu",
                "active_load_step_pu",
                "reactive_load_step_pu",
            ]
        )
        for row in zip(
            result["time"],
            result["frequency_error"],
            result["voltage_error"],
            result["frequency_control"],
            result["voltage_control"],
            result["active_load"],
            result["reactive_load"],
        ):
            writer.writerow(row)


def plot_result(result: dict[str, np.ndarray], output: Path) -> None:
    time = result["time"]
    fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True)
    axes[0].plot(time, result["frequency_error"], label="frequency error")
    axes[0].plot(time, result["active_load"], "--", label="active-load disturbance")
    axes[0].set_ylabel("Per-unit")
    axes[0].grid(True)
    axes[0].legend()

    axes[1].plot(time, result["voltage_error"], label="voltage error")
    axes[1].plot(time, result["reactive_load"], "--", label="reactive-load disturbance")
    axes[1].set_xlabel("Time (s)")
    axes[1].set_ylabel("Per-unit")
    axes[1].grid(True)
    axes[1].legend()

    fig.suptitle("Conceptual Adaptive Frequency and Voltage Regulation")
    fig.tight_layout()
    fig.savefig(output, dpi=180)
    plt.close(fig)


def main() -> None:
    result = simulate()
    save_csv(result, Path("Chapter27_Lesson3_python_results.csv"))
    plot_result(result, Path("Chapter27_Lesson3_python_results.png"))

    print(f"Final frequency error: {result['frequency_error'][-1]: .6e}")
    print(f"Final voltage error:   {result['voltage_error'][-1]: .6e}")
    print("Created CSV and PNG result files.")


if __name__ == "__main__":
    main()

The implementation uses NumPy for vector operations and Matplotlib for plots. For higher-fidelity studies, SciPy, python-control, pandapower, or a validated co-simulation interface may be added, but the adaptive law should first be tested against explicit actuator and measurement constraints.

Chapter27_Lesson3.cpp

/*
Chapter27_Lesson3.cpp
Conceptual adaptive frequency and voltage regulation.

Build:
    g++ -std=c++17 -O2 Chapter27_Lesson3.cpp -o Chapter27_Lesson3

Run:
    ./Chapter27_Lesson3
*/

#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>

using Vec3 = std::array<double, 3>;

double dot(const Vec3& a, const Vec3& b) {
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

Vec3 projectedEuler(
    const Vec3& theta,
    const Vec3& derivative,
    double dt,
    double lower,
    double upper
) {
    Vec3 next{};
    for (std::size_t i = 0; i < next.size(); ++i) {
        next[i] = std::clamp(theta[i] + dt * derivative[i], lower, upper);
    }
    return next;
}

int main() {
    constexpr double dt = 0.01;
    constexpr double finalTime = 80.0;
    constexpr int steps = static_cast<int>(finalTime / dt) + 1;

    const Vec3 thetaFTrue{-0.55, -0.20, 1.00};
    const Vec3 thetaVTrue{-0.75, -0.28, 0.85};
    Vec3 thetaFHat{-0.10, -0.05, 0.20};
    Vec3 thetaVHat{-0.15, -0.05, 0.10};

    const Vec3 gammaF{2.0, 0.8, 3.0};
    const Vec3 gammaV{1.6, 0.7, 2.4};

    constexpr double bF = 1.20;
    constexpr double bV = 1.05;
    constexpr double kF = 1.80;
    constexpr double kV = 2.10;

    double eF = 0.025;
    double eV = -0.030;
    double xiF = 0.0;
    double xiV = 0.0;

    std::ofstream csv("Chapter27_Lesson3_cpp_results.csv");
    if (!csv) {
        throw std::runtime_error("Cannot create output CSV file.");
    }
    csv << "time_s,frequency_error_pu,voltage_error_pu,"
           "frequency_control_pu,voltage_control_pu,"
           "active_load_step_pu,reactive_load_step_pu\n";
    csv << std::setprecision(12);

    for (int k = 0; k < steps; ++k) {
        const double t = k * dt;
        const double loadP = (t < 8.0) ? 0.0 : ((t < 42.0) ? 0.035 : 0.055);
        const double loadQ = (t < 15.0) ? 0.0 : ((t < 50.0) ? -0.025 : 0.040);

        const Vec3 phiF{eF, xiF, loadP};
        const Vec3 phiV{eV, xiV, loadQ};

        const double uF = (-dot(thetaFHat, phiF) - kF * eF) / bF;
        const double uV = (-dot(thetaVHat, phiV) - kV * eV) / bV;

        csv << t << ',' << eF << ',' << eV << ',' << uF << ',' << uV
            << ',' << loadP << ',' << loadQ << '\n';

        if (k == steps - 1) {
            break;
        }

        const double unmatchedF = 0.0015 * std::sin(0.35 * t);
        const double unmatchedV = 0.0010 * std::cos(0.27 * t);

        const double eFDot = dot(thetaFTrue, phiF) + bF * uF + unmatchedF;
        const double eVDot = dot(thetaVTrue, phiV) + bV * uV + unmatchedV;

        Vec3 thetaFDot{};
        Vec3 thetaVDot{};
        for (std::size_t i = 0; i < 3; ++i) {
            thetaFDot[i] = gammaF[i] * phiF[i] * eF;
            thetaVDot[i] = gammaV[i] * phiV[i] * eV;
        }

        eF += dt * eFDot;
        eV += dt * eVDot;
        xiF += dt * eF;
        xiV += dt * eV;
        thetaFHat = projectedEuler(thetaFHat, thetaFDot, dt, -3.0, 3.0);
        thetaVHat = projectedEuler(thetaVHat, thetaVDot, dt, -3.0, 3.0);
    }

    std::cout << std::scientific
              << "Final frequency error: " << eF << '\n'
              << "Final voltage error:   " << eV << '\n'
              << "Created Chapter27_Lesson3_cpp_results.csv\n";
    return 0;
}

This version uses only the C++17 standard library. For matrix-rich models, Eigen is a common linear-algebra choice; production real-time deployment also requires deterministic scheduling, finite-value checks, and explicit failover behavior.

Chapter27_Lesson3.java

/*
Chapter27_Lesson3.java
Conceptual adaptive frequency and voltage regulation.

Build:
    javac Chapter27_Lesson3.java

Run:
    java Chapter27_Lesson3
*/

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.charset.StandardCharsets;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Locale;

public final class Chapter27_Lesson3 {
    private Chapter27_Lesson3() {}

    private static double dot(double[] a, double[] b) {
        return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
    }

    private static void projectedEuler(
            double[] theta,
            double[] derivative,
            double dt,
            double lower,
            double upper) {
        for (int i = 0; i < theta.length; i++) {
            double candidate = theta[i] + dt * derivative[i];
            theta[i] = Math.max(lower, Math.min(upper, candidate));
        }
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);

        final double dt = 0.01;
        final double finalTime = 80.0;
        final int steps = (int) (finalTime / dt) + 1;

        final double[] thetaFTrue = {-0.55, -0.20, 1.00};
        final double[] thetaVTrue = {-0.75, -0.28, 0.85};
        final double[] thetaFHat = {-0.10, -0.05, 0.20};
        final double[] thetaVHat = {-0.15, -0.05, 0.10};
        final double[] gammaF = {2.0, 0.8, 3.0};
        final double[] gammaV = {1.6, 0.7, 2.4};

        final double bF = 1.20;
        final double bV = 1.05;
        final double kF = 1.80;
        final double kV = 2.10;

        double eF = 0.025;
        double eV = -0.030;
        double xiF = 0.0;
        double xiV = 0.0;

        Path output = Path.of("Chapter27_Lesson3_java_results.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(
                output, StandardCharsets.UTF_8)) {
            writer.write(
                "time_s,frequency_error_pu,voltage_error_pu,"
                + "frequency_control_pu,voltage_control_pu,"
                + "active_load_step_pu,reactive_load_step_pu\n"
            );

            for (int k = 0; k < steps; k++) {
                double t = k * dt;
                double loadP = t < 8.0 ? 0.0 : (t < 42.0 ? 0.035 : 0.055);
                double loadQ = t < 15.0 ? 0.0 : (t < 50.0 ? -0.025 : 0.040);

                double[] phiF = {eF, xiF, loadP};
                double[] phiV = {eV, xiV, loadQ};

                double uF = (-dot(thetaFHat, phiF) - kF * eF) / bF;
                double uV = (-dot(thetaVHat, phiV) - kV * eV) / bV;

                writer.write(String.format(
                    Locale.US,
                    "%.8f,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                    t, eF, eV, uF, uV, loadP, loadQ
                ));

                if (k == steps - 1) {
                    break;
                }

                double unmatchedF = 0.0015 * Math.sin(0.35 * t);
                double unmatchedV = 0.0010 * Math.cos(0.27 * t);

                double eFDot = dot(thetaFTrue, phiF) + bF * uF + unmatchedF;
                double eVDot = dot(thetaVTrue, phiV) + bV * uV + unmatchedV;

                double[] thetaFDot = new double[3];
                double[] thetaVDot = new double[3];
                for (int i = 0; i < 3; i++) {
                    thetaFDot[i] = gammaF[i] * phiF[i] * eF;
                    thetaVDot[i] = gammaV[i] * phiV[i] * eV;
                }

                eF += dt * eFDot;
                eV += dt * eVDot;
                xiF += dt * eF;
                xiV += dt * eV;
                projectedEuler(thetaFHat, thetaFDot, dt, -3.0, 3.0);
                projectedEuler(thetaVHat, thetaVDot, dt, -3.0, 3.0);
            }
        }

        System.out.printf(Locale.US, "Final frequency error: %.6e%n", eF);
        System.out.printf(Locale.US, "Final voltage error:   %.6e%n", eV);
        System.out.println("Created " + output);
    }
}

The Java implementation is dependency-free. EJML or Apache Commons Math can support larger state-space calculations, while industrial communication normally requires a separately validated interface layer.

Chapter27_Lesson3.m

% Chapter27_Lesson3.m
% Conceptual adaptive frequency and voltage regulation.
% The model is intentionally reduced-order and intended for teaching.
%
% Run:
%   Chapter27_Lesson3

clear; clc; close all;

dt = 0.01;
tf = 80;
t = (0:dt:tf)';
N = numel(t);

thetaFTrue = [-0.55; -0.20; 1.00];
thetaVTrue = [-0.75; -0.28; 0.85];
thetaFHat = zeros(3, N);
thetaVHat = zeros(3, N);
thetaFHat(:,1) = [-0.10; -0.05; 0.20];
thetaVHat(:,1) = [-0.15; -0.05; 0.10];

GammaF = diag([2.0, 0.8, 3.0]);
GammaV = diag([1.6, 0.7, 2.4]);
bF = 1.20; bV = 1.05;
kF = 1.80; kV = 2.10;

eF = zeros(N,1); eV = zeros(N,1);
xiF = zeros(N,1); xiV = zeros(N,1);
uF = zeros(N,1); uV = zeros(N,1);
loadP = zeros(N,1); loadQ = zeros(N,1);
eF(1) = 0.025;
eV(1) = -0.030;

lower = -3 * ones(3,1);
upper =  3 * ones(3,1);

for k = 1:N-1
    tk = t(k);
    if tk < 8
        loadP(k) = 0;
    elseif tk < 42
        loadP(k) = 0.035;
    else
        loadP(k) = 0.055;
    end

    if tk < 15
        loadQ(k) = 0;
    elseif tk < 50
        loadQ(k) = -0.025;
    else
        loadQ(k) = 0.040;
    end

    phiF = [eF(k); xiF(k); loadP(k)];
    phiV = [eV(k); xiV(k); loadQ(k)];

    uF(k) = (-thetaFHat(:,k)' * phiF - kF * eF(k)) / bF;
    uV(k) = (-thetaVHat(:,k)' * phiV - kV * eV(k)) / bV;

    dF = 0.0015 * sin(0.35 * tk);
    dV = 0.0010 * cos(0.27 * tk);

    eFdot = thetaFTrue' * phiF + bF * uF(k) + dF;
    eVdot = thetaVTrue' * phiV + bV * uV(k) + dV;

    thetaFdot = GammaF * phiF * eF(k);
    thetaVdot = GammaV * phiV * eV(k);

    eF(k+1) = eF(k) + dt * eFdot;
    eV(k+1) = eV(k) + dt * eVdot;
    xiF(k+1) = xiF(k) + dt * eF(k);
    xiV(k+1) = xiV(k) + dt * eV(k);

    thetaFHat(:,k+1) = min(max(thetaFHat(:,k) + dt*thetaFdot, lower), upper);
    thetaVHat(:,k+1) = min(max(thetaVHat(:,k) + dt*thetaVdot, lower), upper);
end

loadP(end) = loadP(end-1);
loadQ(end) = loadQ(end-1);
uF(end) = uF(end-1);
uV(end) = uV(end-1);

results = table(t, eF, eV, uF, uV, loadP, loadQ);
writetable(results, 'Chapter27_Lesson3_matlab_results.csv');

figure('Name', 'Adaptive Frequency and Voltage Regulation');
tiledlayout(2,1);

nexttile;
plot(t, eF, 'LineWidth', 1.3); hold on;
plot(t, loadP, '--', 'LineWidth', 1.0);
grid on;
ylabel('Per unit');
legend('Frequency error', 'Active-load disturbance', 'Location', 'best');

nexttile;
plot(t, eV, 'LineWidth', 1.3); hold on;
plot(t, loadQ, '--', 'LineWidth', 1.0);
grid on;
xlabel('Time (s)');
ylabel('Per unit');
legend('Voltage error', 'Reactive-load disturbance', 'Location', 'best');

exportgraphics(gcf, 'Chapter27_Lesson3_matlab_results.png', 'Resolution', 180);

fprintf('Final frequency error: %.6e\n', eF(end));
fprintf('Final voltage error:   %.6e\n', eV(end));

% Simulink mapping:
% 1. Use one MATLAB Function block per adaptive channel.
% 2. Inputs: e, integral(e), measured disturbance proxy.
% 3. Outputs: certainty-equivalent control and parameter estimates.
% 4. Place Saturation blocks on control commands and parameter estimates.
% 5. Use Rate Limiter blocks for turbine/governor and exciter limits.
% 6. Log eF, eV, uF, uV, and all adaptive estimates.

In Simulink, place the adaptive update in a discrete MATLAB Function block, keep the physical plant and actuator dynamics in separate subsystems, and add Saturation, Rate Limiter, Enabled Subsystem, and fallback-controller logic. Simscape Electrical can supply richer network components, but its parameter set and solver configuration must be validated for the study.

Chapter27_Lesson3.nb

Notebook[{
 Cell["Chapter 27, Lesson 3: Adaptive Frequency/Voltage Control", "Title"],
 Cell["Conceptual reduced-order simulation with explicit Euler integration and componentwise projection.", "Text"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"dt", "=", "0.01"}], ";",
   RowBox[{"tf", "=", "80.0"}], ";",
   RowBox[{"time", "=", RowBox[{"Range", "[", RowBox[{"0.0", ",", "tf", ",", "dt"}], "]"}]}], ";",
   RowBox[{"n", "=", RowBox[{"Length", "[", "time", "]"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"thetaFTrue", "=", RowBox[{"{", RowBox[{"-0.55", ",", "-0.20", ",", "1.00"}], "}"}]}], ";",
   RowBox[{"thetaVTrue", "=", RowBox[{"{", RowBox[{"-0.75", ",", "-0.28", ",", "0.85"}], "}"}]}], ";",
   RowBox[{"gammaF", "=", RowBox[{"{", RowBox[{"2.0", ",", "0.8", ",", "3.0"}], "}"}]}], ";",
   RowBox[{"gammaV", "=", RowBox[{"{", RowBox[{"1.6", ",", "0.7", ",", "2.4"}], "}"}]}], ";",
   RowBox[{"bF", "=", "1.20"}], ";", RowBox[{"bV", "=", "1.05"}], ";",
   RowBox[{"kF", "=", "1.80"}], ";", RowBox[{"kV", "=", "2.10"}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"eF", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", "n"}], "]"}]}], ";",
   RowBox[{"eV", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", "n"}], "]"}]}], ";",
   RowBox[{"xiF", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", "n"}], "]"}]}], ";",
   RowBox[{"xiV", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", "n"}], "]"}]}], ";",
   RowBox[{"uF", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", "n"}], "]"}]}], ";",
   RowBox[{"uV", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", "n"}], "]"}]}], ";",
   RowBox[{"thetaFHat", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", RowBox[{"{", RowBox[{"n", ",", "3"}], "}"}]}], "]"}]}], ";",
   RowBox[{"thetaVHat", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", RowBox[{"{", RowBox[{"n", ",", "3"}], "}"}]}], "]"}]}], ";",
   RowBox[{RowBox[{"eF", "[", RowBox[{"[", "1", "]"}], "]"}], "=", "0.025"}], ";",
   RowBox[{RowBox[{"eV", "[", RowBox[{"[", "1", "]"}], "]"}], "=", "-0.030"}], ";",
   RowBox[{RowBox[{"thetaFHat", "[", RowBox[{"[", "1", "]"}], "]"}], "=", RowBox[{"{", RowBox[{"-0.10", ",", "-0.05", ",", "0.20"}], "}"}]}], ";",
   RowBox[{RowBox[{"thetaVHat", "[", RowBox[{"[", "1", "]"}], "]"}], "=", RowBox[{"{", RowBox[{"-0.15", ",", "-0.05", ",", "0.10"}], "}"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{"Do", "[",
   RowBox[{
    RowBox[{
     RowBox[{"t", "=", RowBox[{"time", "[", RowBox[{"[", "k", "]"}], "]"}]}], ";",
     RowBox[{"loadP", "=", RowBox[{"Piecewise", "[", RowBox[{"{", RowBox[{
       RowBox[{"{", RowBox[{"0.0", ",", RowBox[{"t", "<", "8.0"}]}], "}"}], ",",
       RowBox[{"{", RowBox[{"0.035", ",", RowBox[{"t", "<", "42.0"}]}], "}"}]
       }], "}", ",", "0.055"}], "]"}]}], ";",
     RowBox[{"loadQ", "=", RowBox[{"Piecewise", "[", RowBox[{"{", RowBox[{
       RowBox[{"{", RowBox[{"0.0", ",", RowBox[{"t", "<", "15.0"}]}], "}"}], ",",
       RowBox[{"{", RowBox[{"-0.025", ",", RowBox[{"t", "<", "50.0"}]}], "}"}]
       }], "}", ",", "0.040"}], "]"}]}], ";",
     RowBox[{"phiF", "=", RowBox[{"{", RowBox[{RowBox[{"eF", "[", RowBox[{"[", "k", "]"}], "]"}], ",", RowBox[{"xiF", "[", RowBox[{"[", "k", "]"}], "]"}], ",", "loadP"}], "}"}]}], ";",
     RowBox[{"phiV", "=", RowBox[{"{", RowBox[{RowBox[{"eV", "[", RowBox[{"[", "k", "]"}], "]"}], ",", RowBox[{"xiV", "[", RowBox[{"[", "k", "]"}], "]"}], ",", "loadQ"}], "}"}]}], ";",
     RowBox[{RowBox[{"uF", "[", RowBox[{"[", "k", "]"}], "]"}], "=", RowBox[{
       RowBox[{"(", RowBox[{RowBox[{"-", RowBox[{RowBox[{"thetaFHat", "[", RowBox[{"[", "k", "]"}], "]"}], ".", "phiF"}]}], "-", RowBox[{"kF", " ", RowBox[{"eF", "[", RowBox[{"[", "k", "]"}], "]"}]}]}], ")"}], "/", "bF"}]}], ";",
     RowBox[{RowBox[{"uV", "[", RowBox[{"[", "k", "]"}], "]"}], "=", RowBox[{
       RowBox[{"(", RowBox[{RowBox[{"-", RowBox[{RowBox[{"thetaVHat", "[", RowBox[{"[", "k", "]"}], "]"}], ".", "phiV"}]}], "-", RowBox[{"kV", " ", RowBox[{"eV", "[", RowBox[{"[", "k", "]"}], "]"}]}]}], ")"}], "/", "bV"}]}], ";",
     RowBox[{"eFDot", "=", RowBox[{"thetaFTrue", ".", "phiF", "+", RowBox[{"bF", " ", RowBox[{"uF", "[", RowBox[{"[", "k", "]"}], "]"}]}], "+", RowBox[{"0.0015", " ", RowBox[{"Sin", "[", RowBox[{"0.35", " ", "t"}], "]"}]}]}]}], ";",
     RowBox[{"eVDot", "=", RowBox[{"thetaVTrue", ".", "phiV", "+", RowBox[{"bV", " ", RowBox[{"uV", "[", RowBox[{"[", "k", "]"}], "]"}]}], "+", RowBox[{"0.0010", " ", RowBox[{"Cos", "[", RowBox[{"0.27", " ", "t"}], "]"}]}]}]}], ";",
     RowBox[{RowBox[{"eF", "[", RowBox[{"[", RowBox[{"k", "+", "1"}], "]"}], "]"}], "=", RowBox[{RowBox[{"eF", "[", RowBox[{"[", "k", "]"}], "]"}], "+", RowBox[{"dt", " ", "eFDot"}]}]}], ";",
     RowBox[{RowBox[{"eV", "[", RowBox[{"[", RowBox[{"k", "+", "1"}], "]"}], "]"}], "=", RowBox[{RowBox[{"eV", "[", RowBox[{"[", "k", "]"}], "]"}], "+", RowBox[{"dt", " ", "eVDot"}]}]}], ";",
     RowBox[{RowBox[{"xiF", "[", RowBox[{"[", RowBox[{"k", "+", "1"}], "]"}], "]"}], "=", RowBox[{RowBox[{"xiF", "[", RowBox[{"[", "k", "]"}], "]"}], "+", RowBox[{"dt", " ", RowBox[{"eF", "[", RowBox[{"[", "k", "]"}], "]"}]}]}]}], ";",
     RowBox[{RowBox[{"xiV", "[", RowBox[{"[", RowBox[{"k", "+", "1"}], "]"}], "]"}], "=", RowBox[{RowBox[{"xiV", "[", RowBox[{"[", "k", "]"}], "]"}], "+", RowBox[{"dt", " ", RowBox[{"eV", "[", RowBox[{"[", "k", "]"}], "]"}]}]}]}], ";",
     RowBox[{RowBox[{"thetaFHat", "[", RowBox[{"[", RowBox[{"k", "+", "1"}], "]"}], "]"}], "=", RowBox[{"Clip", "[", RowBox[{
       RowBox[{RowBox[{"thetaFHat", "[", RowBox[{"[", "k", "]"}], "]"}], "+", RowBox[{"dt", " ", "gammaF", " ", "phiF", " ", RowBox[{"eF", "[", RowBox[{"[", "k", "]"}], "]"}]}]}], ",", RowBox[{"{", RowBox[{"-3.0", ",", "3.0"}], "}"}]}], "]"}]}], ";",
     RowBox[{RowBox[{"thetaVHat", "[", RowBox[{"[", RowBox[{"k", "+", "1"}], "]"}], "]"}], "=", RowBox[{"Clip", "[", RowBox[{
       RowBox[{RowBox[{"thetaVHat", "[", RowBox[{"[", "k", "]"}], "]"}], "+", RowBox[{"dt", " ", "gammaV", " ", "phiV", " ", RowBox[{"eV", "[", RowBox[{"[", "k", "]"}], "]"}]}]}]}], ",", RowBox[{"{", RowBox[{"-3.0", ",", "3.0"}], "}"}]}], "]"}]}], ";"
    }], ",", RowBox[{"{", RowBox[{"k", ",", "1", ",", RowBox[{"n", "-", "1"}]}], "}"}]}], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"GraphicsGrid", "[", RowBox[{"{", RowBox[{"{",
   RowBox[{
    RowBox[{"ListLinePlot", "[", RowBox[{RowBox[{"Transpose", "[", RowBox[{"{", RowBox[{"time", ",", "eF"}], "}"}], "]"}], ",", RowBox[{"PlotLabel", "->", "\"Frequency error\""}], ",", RowBox[{"AxesLabel", "->", RowBox[{"{", RowBox[{"\"t (s)\"", ",", "\"e_f\""}], "}"}]}]}], "]"}], ",",
    RowBox[{"ListLinePlot", "[", RowBox[{RowBox[{"Transpose", "[", RowBox[{"{", RowBox[{"time", ",", "eV"}], "}"}], "]"}], ",", RowBox[{"PlotLabel", "->", "\"Voltage error\""}], ",", RowBox[{"AxesLabel", "->", RowBox[{"{", RowBox[{"\"t (s)\"", ",", "\"e_v\""}], "}"}]}]}], "]"}]
   }], "}"}], "}"}], "]"}]], "Input"]
},
WindowTitle -> "Chapter27_Lesson3",
StyleDefinitions -> "Default.nb"
]

The notebook expression creates an executable notebook using built-in list operations, explicit Euler integration, projection with Clip, and ListLinePlot.

The software examples use two intentionally simple matched error channels:

\[ \dot e_f=\boldsymbol{\theta}_f^{*\mathsf T} [e_f,\xi_f,\Delta P_L]^{\mathsf T}+b_fu_f+d_f, \]

\[ \dot e_v=\boldsymbol{\theta}_v^{*\mathsf T} [e_v,\xi_v,\Delta Q_L]^{\mathsf T}+b_vu_v+d_v. \]

Step changes represent operating-condition changes, while small sinusoidal terms represent unmatched residual dynamics. The plots should be read as a demonstration of adaptive cancellation and bounded parameter updates—not as a prediction of grid frequency nadir, voltage-security margin, relay operation, or inter-area oscillation damping.

Problem 1 — Primary droop error. An isolated area has \( D=1.0 \) pu load damping and \( R=0.05 \) pu droop. A load increase of \( \Delta P_L=0.02 \) pu occurs. Find the quasi-steady frequency deviation with primary control only.

Solution:

\[ \Delta f_{\mathrm{ss}} =-\frac{0.02}{1.0+1/0.05} =-\frac{0.02}{21} =-9.5238\times 10^{-4}\ \mathrm{pu}. \]

On a 50 Hz base, this is approximately \( -0.0476\ \mathrm{Hz} \). The nonzero value demonstrates why secondary restoration is required.

Problem 2 — ACE equilibrium. For a two-area system, show that \( \mathrm{ACE}_1=\mathrm{ACE}_2=0 \) implies nominal steady frequency and scheduled tie-line interchange.

Solution:

\[ B_1\Delta f+\Delta P_{12}=0,\qquad B_2\Delta f-\Delta P_{12}=0. \]

Adding the equations gives

\[ (B_1+B_2)\Delta f=0. \]

Because \( B_1+B_2\neq 0 \), \( \Delta f=0 \). Substitution into either ACE equation gives \( \Delta P_{12}=0 \), meaning that the incremental interchange error has been removed.

Problem 3 — Lyapunov cancellation. For \( \dot e=-ke+\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\phi} \), derive an adaptive law that makes \( \dot V=-ke^2 \) for \( V=\frac12e^2+\frac12\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta}} \).

Solution:

Since \( \tilde{\boldsymbol{\theta}}= \boldsymbol{\theta}^*-\hat{\boldsymbol{\theta}} \), choose

\[ \dot{\hat{\boldsymbol{\theta}}} =\boldsymbol{\Gamma}\boldsymbol{\phi}e, \qquad \dot{\tilde{\boldsymbol{\theta}}} =-\boldsymbol{\Gamma}\boldsymbol{\phi}e. \]

\[ \begin{aligned} \dot V &=e(-ke+\tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\phi}) +\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1} (-\boldsymbol{\Gamma}\boldsymbol{\phi}e)\\ &=-ke^2. \end{aligned} \]

The cancellation proves boundedness and error-energy decay, but not parameter convergence without persistent excitation.

Problem 4 — Effect of leakage. Consider \( \dot{\hat{\boldsymbol{\theta}}} =\boldsymbol{\Gamma}\boldsymbol{\phi}e -\sigma\boldsymbol{\Gamma}\hat{\boldsymbol{\theta}} \). Show qualitatively why the origin is generally replaced by an ultimate bound when \( \boldsymbol{\theta}^*\neq\mathbf{0} \).

Solution:

The parameter-error term in the Lyapunov derivative becomes

\[ -\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1}\dot{\hat{\boldsymbol{\theta}}} = -\tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\phi}e +\sigma\tilde{\boldsymbol{\theta}}^{\mathsf T} \hat{\boldsymbol{\theta}}. \]

The first term cancels the cross term, but the second does not vanish. Using \( \hat{\boldsymbol{\theta}}= \boldsymbol{\theta}^*-\tilde{\boldsymbol{\theta}} \),

\[ \sigma\tilde{\boldsymbol{\theta}}^{\mathsf T} \hat{\boldsymbol{\theta}} = \sigma\tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\theta}^* -\sigma\|\tilde{\boldsymbol{\theta}}\|^2. \]

Young's inequality bounds the first term by a constant plus a quadratic term. Consequently, \( \dot V \) is negative outside a compact set, giving uniform ultimate boundedness rather than exact convergence.

Problem 5 — Voltage command under a reactive-load step. For \( \dot e_v=\theta_1e_v+\theta_2\xi_v+ \theta_3\Delta Q_L+b_vu_v \), write the certainty-equivalent controller and update law.

Solution:

\[ \boldsymbol{\phi}_v= [e_v,\xi_v,\Delta Q_L]^{\mathsf T}, \qquad u_v=\frac{-\hat{\boldsymbol{\theta}}_v^{\mathsf T} \boldsymbol{\phi}_v-k_ve_v}{b_v}, \]

\[ \dot{\hat{\boldsymbol{\theta}}}_v =\boldsymbol{\Gamma}_v\boldsymbol{\phi}_ve_v. \]

In implementation, apply voltage-command saturation, reactive-current limits, projection, signal validation, and a fallback AVR or converter controller. Without these constraints, the mathematical command may be physically infeasible.

Frequency adaptation was formulated around an ACE or frequency-error channel, and voltage adaptation around a voltage/reactive-power channel. In each case, linearly parameterized matched uncertainty permits a certainty-equivalent controller whose Lyapunov update law cancels the state-parameter cross term. Projection, normalization, leakage, dead zones, actuator constraints, and certified fallback logic are essential practical additions. The reduced proofs establish properties of the assumed error models; they do not by themselves certify a real interconnected grid.