In fixed-parameter feedback, the controller parameters are selected once and then kept constant. In adaptive feedback, at least one controller parameter or model parameter is updated online from measured signals. At this introductory level, write a broad adaptive controller as

\[ \dot{x}=f(x,u,\theta),\qquad u=C(x,r,\hat{\theta},\rho),\qquad \dot{\hat{\theta} }=\Gamma\,\Psi(x,r,e,\hat{\theta}) \]

where \( x \) is the plant state or measured output vector, \( r \) is the command/reference, \( \theta \) is the unknown plant parameter vector, \( \hat{\theta} \) is its online estimate or an adapted controller parameter, and \( \rho \) is a scheduling variable such as speed, altitude, operating point, load, or temperature.

The major classes differ by the answer to one question: What object is adjusted online? MRAC adjusts controller parameters to make a plant imitate a reference model; STR estimates plant parameters and recomputes a controller; adaptive backstepping recursively designs stabilizing virtual controls for nonlinear strict-feedback systems with unknown parameters; gain scheduling switches or interpolates fixed controllers using a measured scheduling variable.

flowchart TD
  A["Adaptive / parameter-varying
control problem"] --> B["Reference behavior
specified?"]
  B -->|yes| C["MRAC:
adapt controller gains
to track reference model"]
  B -->|no| D["Plant model estimated
online?"]
  D -->|yes| E["STR:
identify plant parameters
then recompute controller"]
  D -->|no| F["Nonlinear strict-feedback
structure?"]
  F -->|yes| G["Adaptive backstepping:
recursive virtual-control
design"]
  F -->|no| H["Known operating variable
available?"]
  H -->|yes| I["Gain scheduling:
interpolate predesigned
fixed controllers"]
  H -->|no| J["Need more modeling
before selecting an
adaptive architecture"]
        

MRAC begins with a desired closed-loop model. For a first-order scalar teaching example, let the uncertain plant and reference model be

\[ \dot{x}=a x+b u,\qquad \dot{x}_m=a_m x_m+b_m r,\qquad a_m < 0 \]

The adaptive controller is chosen as

\[ u=\hat{\theta}_x x+\hat{\theta}_r r. \]

If there exist ideal constants \( \theta_x^*,\theta_r^* \) satisfying the matching equations

\[ a+b\theta_x^*=a_m,\qquad b\theta_r^*=b_m, \]

then the plant can behave like the reference model when \( \hat{\theta}_x=\theta_x^* \) and \( \hat{\theta}_r=\theta_r^* \). Define tracking error and parameter errors by

\[ e=x-x_m,\qquad \tilde{\theta}_x=\hat{\theta}_x-\theta_x^*,\qquad \tilde{\theta}_r=\hat{\theta}_r-\theta_r^*. \]

Substitution gives the MRAC error model

\[ \dot{e}=a_m e+b\tilde{\theta}_x x+b\tilde{\theta}_r r. \]

For known positive input-gain sign, a standard Lyapunov-motivated update law is

\[ \dot{\hat{\theta} }_x=-\gamma_x bPex,\qquad \dot{\hat{\theta} }_r=-\gamma_r bPer, \qquad 2a_mP=-1,\quad P > 0. \]

Use the candidate

\[ V=Pe^2+\frac{b}{\gamma_x}\tilde{\theta}_x^2+ \frac{b}{\gamma_r}\tilde{\theta}_r^2. \]

Differentiating along the trajectories and using the above adaptive laws gives

\[ \dot{V}=2Pa_m e^2=-e^2\le 0. \]

This small derivation shows the central MRAC idea: the adaptation law is not arbitrary tuning; it is chosen to cancel cross-terms in the derivative of a Lyapunov function. Later chapters will generalize this to higher-order systems, unknown gain signs, normalization, projection, and robustness modifications.

A self-tuning regulator is usually an indirect adaptive controller. Instead of adapting the controller gains directly, STR first estimates the plant model and then recomputes the controller from that estimate. In a discrete-time first-order teaching example,

\[ y_{k+1}=a y_k+b u_k=\phi_k^T\theta,\qquad \phi_k=\begin{bmatrix}y_k\\u_k\end{bmatrix},\quad \theta=\begin{bmatrix}a\\b\end{bmatrix}. \]

A common online estimator is recursive least squares (RLS):

\[ K_k=\frac{P_{k-1}\phi_k}{\lambda+\phi_k^T P_{k-1}\phi_k}, \qquad \hat{\theta}_k=\hat{\theta}_{k-1} +K_k\left(y_{k+1}-\phi_k^T\hat{\theta}_{k-1}\right) \]

\[ P_k=\frac{1}{\lambda} \left(P_{k-1}-K_k\phi_k^T P_{k-1}\right),\qquad 0 < \lambda\le 1. \]

Under certainty-equivalence control, the controller is recomputed as if \( \hat{a}_k,\hat{b}_k \) were the true plant parameters. A one-step tracking law is

\[ u_k=\frac{r_{k+1}-\hat{a}_k y_k}{\hat{b}_k}. \]

STR is therefore a two-loop architecture: an estimation loop plus a control-design loop. Its practical strength is that it can reuse familiar non-adaptive design formulas after identifying a plant model. Its weakness is that poor online estimates can temporarily produce poor controllers unless safeguards are used.

Adaptive backstepping applies Lyapunov design recursively to nonlinear systems with a triangular or strict-feedback structure. To keep this introduction within Chapter 1 scope, consider the second-order uncertain plant

\[ \dot{x}_1=x_2,\qquad \dot{x}_2=\theta x_1+u, \]

where \( \theta \) is unknown. Choose the first error coordinate and virtual control as

\[ z_1=x_1,\qquad \alpha=-c_1z_1,\qquad z_2=x_2-\alpha=x_2+c_1x_1, \qquad c_1 > 0. \]

The adaptive backstepping controller

\[ u=-\hat{\theta}z_1-z_1-c_2z_2-c_1z_2+c_1^2z_1,\qquad c_2 > 0 \]

with update law

\[ \dot{\hat{\theta} }=\gamma z_1z_2,\qquad \gamma > 0 \]

yields the transformed dynamics

\[ \dot{z}_1=z_2-c_1z_1,\qquad \dot{z}_2=\tilde{\theta}z_1-z_1-c_2z_2, \qquad \tilde{\theta}=\theta-\hat{\theta}. \]

With

\[ V=\frac{1}{2}z_1^2+\frac{1}{2}z_2^2+ \frac{1}{2\gamma}\tilde{\theta}^2, \]

the derivative becomes

\[ \dot{V}=-c_1z_1^2-c_2z_2^2. \]

The important lesson is not the final formula itself; it is the recursive construction: stabilize the first subsystem with a virtual control, define an error between the real state and the virtual control, then design the actual input and parameter update so that the Lyapunov derivative is negative semidefinite.

Gain scheduling is often discussed beside adaptive control because controller gains change during operation. However, gain scheduling does not necessarily identify unknown parameters and does not necessarily contain a Lyapunov adaptation law. A typical linear parameter-varying description is

\[ \dot{x}=A(\rho)x+B(\rho)u,\qquad u=-K(\rho)x+N(\rho)r. \]

The gain is commonly interpolated from fixed controllers designed at frozen operating points:

\[ K(\rho)=\sum_{i=1}^{M}w_i(\rho)K_i,\qquad \sum_{i=1}^{M}w_i(\rho)=1,\qquad w_i(\rho)\ge 0. \]

If each frozen closed-loop matrix \( A_i-B_iK_i \) is stable, this alone does not automatically prove stability for fast scheduling changes. A common sufficient check is to seek a common quadratic Lyapunov matrix:

\[ (A_i-B_iK_i)^TP+P(A_i-B_iK_i)=-Q_i,\qquad P=P^T > 0,\qquad Q_i=Q_i^T > 0. \]

This emphasizes the key distinction: MRAC and adaptive backstepping adapt based on error signals, STR adapts by identification, while gain scheduling varies according to a measured operating condition.

The following table summarizes the four classes at a high level. Detailed design will be developed gradually in later chapters.

flowchart LR
  M["MRAC"] --> M1["Reference model"]
  M --> M2["Tracking error drives gain update"]
  S["STR"] --> S1["Online identifier"]
  S --> S2["Controller recomputed from estimated model"]
  B["Adaptive Backstepping"] --> B1["Virtual controls"]
  B --> B2["Recursive Lyapunov design"]
  G["Gain Scheduling"] --> G1["Measured operating variable"]
  G --> G2["Interpolate fixed controllers"]
        

The examples below are deliberately small and mostly implemented from scratch. For larger projects, useful numerical and control libraries include:

• Python: numpy, scipy, matplotlib, and control for state-space models, simulation, and plotting.
• C++: Eigen for matrix algebra, Boost.Odeint for ODE integration, and qpOASES or OSQP when constrained optimization is later needed.
• Java: Apache Commons Math or EJML for linear algebra and numerical routines.
• MATLAB/Simulink: Control System Toolbox, System Identification Toolbox, Simulink, and MATLAB Function blocks for adaptive laws.
• Wolfram Mathematica: NDSolve, StateSpaceModel, symbolic manipulation, and interactive plotting.

In Simulink, the MRAC structure can be built from four subsystems: Plant, Reference Model, Adaptive Law, and Control Law. The adaptive law subsystem integrates \( \dot{\hat{\theta} } \); the control subsystem computes \( u \) from the current parameter estimates.

File: Chapter1_Lesson3.py

# Chapter1_Lesson3.py
# Classes of Adaptive Controllers: MRAC, STR, Adaptive Backstepping, Gain Scheduling
# This script compares four simplified controller classes on small teaching examples.
#
# Required packages:
#   pip install numpy matplotlib
#
# Optional control libraries for larger projects:
#   python-control: pip install control
#   scipy: pip install scipy

import math
import numpy as np
import matplotlib.pyplot as plt


def simulate_mrac(dt=0.001, tf=8.0):
    """
    Direct MRAC for a scalar first-order plant:
        x_dot = a x + b u
        x_m_dot = a_m x_m + b_m r
        u = theta_x x + theta_r r
    Matching requires:
        a + b theta_x* = a_m,  b theta_r* = b_m
    A Lyapunov-style update is used for demonstration.
    """
    n = int(tf / dt)
    t = np.linspace(0.0, tf, n)
    a = -0.6
    b = 1.8
    a_m = -2.0
    b_m = 2.0
    gamma_x = 4.0
    gamma_r = 4.0

    x = 0.0
    xm = 0.0
    theta_x = 0.0
    theta_r = 0.0

    X = np.zeros(n)
    XM = np.zeros(n)
    U = np.zeros(n)
    THX = np.zeros(n)
    THR = np.zeros(n)

    # P solves 2 a_m P = -Q for scalar reference error model; choose Q=1.
    P = 1.0 / (-2.0 * a_m)

    for k in range(n):
        r = 1.0 if t[k] > 0.5 else 0.0
        e = x - xm
        u = theta_x * x + theta_r * r

        # Lyapunov-inspired direct MRAC laws when sign(b) is known positive.
        theta_x_dot = -gamma_x * x * e * P * b
        theta_r_dot = -gamma_r * r * e * P * b

        x_dot = a * x + b * u
        xm_dot = a_m * xm + b_m * r

        x += dt * x_dot
        xm += dt * xm_dot
        theta_x += dt * theta_x_dot
        theta_r += dt * theta_r_dot

        X[k] = x
        XM[k] = xm
        U[k] = u
        THX[k] = theta_x
        THR[k] = theta_r

    return t, X, XM, U, THX, THR


def simulate_str(dt=0.02, tf=10.0):
    """
    Indirect self-tuning regulator (STR) for a discrete ARX-like plant:
        y[k+1] = a y[k] + b u[k]
    The parameter vector [a, b]^T is estimated by recursive least squares (RLS),
    then the controller is recomputed online:
        u[k] = (r[k+1] - a_hat y[k]) / b_hat
    """
    n = int(tf / dt)
    t = np.linspace(0.0, tf, n)

    # True discrete plant parameters
    a_true = 0.92
    b_true = 0.18

    theta = np.array([0.4, 0.4], dtype=float)     # [a_hat, b_hat]
    P = 100.0 * np.eye(2)
    lam = 0.995

    y = 0.0
    u = 0.0

    Y = np.zeros(n)
    U = np.zeros(n)
    AH = np.zeros(n)
    BH = np.zeros(n)

    for k in range(n):
        r_next = 1.0 if t[k] > 1.0 else 0.0

        # Measurement generated by plant
        noise = 0.0
        y_next = a_true * y + b_true * u + noise

        # RLS update using phi[k] = [y[k], u[k]]
        phi = np.array([y, u], dtype=float)
        denom = lam + phi @ P @ phi
        K = (P @ phi) / denom
        prediction_error = y_next - phi @ theta
        theta = theta + K * prediction_error
        P = (P - np.outer(K, phi) @ P) / lam

        # Certainty-equivalence controller update
        a_hat = theta[0]
        b_hat = theta[1]
        if abs(b_hat) < 0.05:
            b_hat = 0.05 if b_hat >= 0.0 else -0.05
        u = (r_next - a_hat * y_next) / b_hat
        u = float(np.clip(u, -8.0, 8.0))

        y = y_next
        Y[k] = y
        U[k] = u
        AH[k] = theta[0]
        BH[k] = theta[1]

    return t, Y, U, AH, BH


def simulate_adaptive_backstepping(dt=0.001, tf=8.0):
    """
    Adaptive backstepping regulation for:
        x1_dot = x2
        x2_dot = theta x1 + u
    Define:
        z1 = x1
        alpha = -c1 z1
        z2 = x2 - alpha = x2 + c1 x1
    Control:
        u = -theta_hat z1 - z1 - c2 z2 - c1 z2 + c1^2 z1
    Adaptive law:
        theta_hat_dot = gamma z1 z2
    Then:
        V_dot = -c1 z1^2 - c2 z2^2
    for the ideal continuous-time equations.
    """
    n = int(tf / dt)
    t = np.linspace(0.0, tf, n)
    theta_true = 1.5
    c1 = 2.0
    c2 = 2.5
    gamma = 5.0

    x1 = 1.5
    x2 = -0.5
    theta_hat = 0.0

    X1 = np.zeros(n)
    X2 = np.zeros(n)
    U = np.zeros(n)
    TH = np.zeros(n)

    for k in range(n):
        z1 = x1
        z2 = x2 + c1 * x1
        u = -theta_hat * z1 - z1 - c2 * z2 - c1 * z2 + c1 * c1 * z1
        theta_hat_dot = gamma * z1 * z2

        x1_dot = x2
        x2_dot = theta_true * x1 + u

        x1 += dt * x1_dot
        x2 += dt * x2_dot
        theta_hat += dt * theta_hat_dot

        X1[k] = x1
        X2[k] = x2
        U[k] = u
        TH[k] = theta_hat

    return t, X1, X2, U, TH


def simulate_gain_scheduling(dt=0.001, tf=8.0):
    """
    Gain scheduling is shown as interpolation between fixed linear gains.
    Plant:
        x_dot = a(rho) x + b u
    Controller:
        u = -K(rho) x + r
    No online identification or parameter adaptation occurs here.
    """
    n = int(tf / dt)
    t = np.linspace(0.0, tf, n)
    b = 1.0

    # Scheduling grid rho and stabilizing gains K for frozen operating points
    rho_grid = np.array([0.0, 0.5, 1.0])
    k_grid = np.array([1.5, 2.5, 3.5])

    x = 0.0
    X = np.zeros(n)
    RHO = np.zeros(n)
    KG = np.zeros(n)

    for k in range(n):
        rho = 0.5 + 0.5 * math.sin(0.7 * t[k])
        a_rho = -0.5 + 1.2 * rho
        K = float(np.interp(rho, rho_grid, k_grid))
        r = 1.0 if t[k] > 0.5 else 0.0
        u = -K * x + r
        x_dot = a_rho * x + b * u
        x += dt * x_dot

        X[k] = x
        RHO[k] = rho
        KG[k] = K

    return t, X, RHO, KG


def main():
    tm, x, xm, u_m, thx, thr = simulate_mrac()
    ts, y, u_s, ah, bh = simulate_str()
    tb, x1, x2, u_b, th = simulate_adaptive_backstepping()
    tg, xg, rho, kg = simulate_gain_scheduling()

    plt.figure()
    plt.plot(tm, x, label="MRAC plant x")
    plt.plot(tm, xm, "--", label="Reference model x_m")
    plt.xlabel("time [s]")
    plt.ylabel("state")
    plt.title("Direct MRAC Tracking")
    plt.legend()
    plt.grid(True)

    plt.figure()
    plt.plot(ts, y, label="STR output y")
    plt.plot(ts, ah, label="a_hat")
    plt.plot(ts, bh, label="b_hat")
    plt.xlabel("time [s]")
    plt.title("Self-Tuning Regulator with RLS")
    plt.legend()
    plt.grid(True)

    plt.figure()
    plt.plot(tb, x1, label="x1")
    plt.plot(tb, x2, label="x2")
    plt.plot(tb, th, label="theta_hat")
    plt.xlabel("time [s]")
    plt.title("Adaptive Backstepping Regulation")
    plt.legend()
    plt.grid(True)

    plt.figure()
    plt.plot(tg, xg, label="scheduled state x")
    plt.plot(tg, rho, label="rho")
    plt.plot(tg, kg, label="K(rho)")
    plt.xlabel("time [s]")
    plt.title("Gain Scheduling")
    plt.legend()
    plt.grid(True)

    plt.show()


if __name__ == "__main__":
    main()

      

File: Chapter1_Lesson3.cpp

// Chapter1_Lesson3.cpp
// Classes of Adaptive Controllers: MRAC, STR, Adaptive Backstepping, Gain Scheduling
// Compile:
//   g++ -std=c++17 -O2 Chapter1_Lesson3.cpp -o Chapter1_Lesson3
// Run:
//   ./Chapter1_Lesson3
//
// The program writes CSV files that can be plotted in Python, MATLAB, Excel, or GNUplot.

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

void simulate_mrac() {
    const double dt = 0.001, tf = 8.0;
    const int n = static_cast<int>(tf / dt);
    const double a = -0.6, b = 1.8, am = -2.0, bm = 2.0;
    const double gammax = 4.0, gammar = 4.0;
    const double P = 1.0 / (-2.0 * am);

    double x = 0.0, xm = 0.0, theta_x = 0.0, theta_r = 0.0;

    std::ofstream file("Chapter1_Lesson3_MRAC.csv");
    file << "t,x,xm,u,theta_x,theta_r\n";

    for (int k = 0; k < n; ++k) {
        double t = k * dt;
        double r = (t > 0.5) ? 1.0 : 0.0;
        double e = x - xm;
        double u = theta_x * x + theta_r * r;

        double theta_x_dot = -gammax * x * e * P * b;
        double theta_r_dot = -gammar * r * e * P * b;

        double x_dot = a * x + b * u;
        double xm_dot = am * xm + bm * r;

        x += dt * x_dot;
        xm += dt * xm_dot;
        theta_x += dt * theta_x_dot;
        theta_r += dt * theta_r_dot;

        file << t << "," << x << "," << xm << "," << u << ","
             << theta_x << "," << theta_r << "\n";
    }
}

void simulate_str() {
    const double dt = 0.02, tf = 10.0;
    const int n = static_cast<int>(tf / dt);
    const double a_true = 0.92, b_true = 0.18, lambda = 0.995;

    std::array<double, 2> theta{0.4, 0.4}; // [a_hat, b_hat]
    double P00 = 100.0, P01 = 0.0, P10 = 0.0, P11 = 100.0;
    double y = 0.0, u = 0.0;

    std::ofstream file("Chapter1_Lesson3_STR.csv");
    file << "t,y,u,a_hat,b_hat\n";

    for (int k = 0; k < n; ++k) {
        double t = k * dt;
        double r_next = (t > 1.0) ? 1.0 : 0.0;
        double y_next = a_true * y + b_true * u;

        double phi0 = y;
        double phi1 = u;
        double Pphi0 = P00 * phi0 + P01 * phi1;
        double Pphi1 = P10 * phi0 + P11 * phi1;
        double denom = lambda + phi0 * Pphi0 + phi1 * Pphi1;

        double K0 = Pphi0 / denom;
        double K1 = Pphi1 / denom;

        double pred = phi0 * theta[0] + phi1 * theta[1];
        double eps = y_next - pred;
        theta[0] += K0 * eps;
        theta[1] += K1 * eps;

        double oldP00 = P00, oldP01 = P01, oldP10 = P10, oldP11 = P11;
        P00 = (oldP00 - K0 * (phi0 * oldP00 + phi1 * oldP10)) / lambda;
        P01 = (oldP01 - K0 * (phi0 * oldP01 + phi1 * oldP11)) / lambda;
        P10 = (oldP10 - K1 * (phi0 * oldP00 + phi1 * oldP10)) / lambda;
        P11 = (oldP11 - K1 * (phi0 * oldP01 + phi1 * oldP11)) / lambda;

        double bhat_safe = theta[1];
        if (std::abs(bhat_safe) < 0.05) {
            bhat_safe = (bhat_safe >= 0.0) ? 0.05 : -0.05;
        }

        u = (r_next - theta[0] * y_next) / bhat_safe;
        u = std::clamp(u, -8.0, 8.0);
        y = y_next;

        file << t << "," << y << "," << u << "," << theta[0] << "," << theta[1] << "\n";
    }
}

void simulate_adaptive_backstepping() {
    const double dt = 0.001, tf = 8.0;
    const int n = static_cast<int>(tf / dt);
    const double theta_true = 1.5, c1 = 2.0, c2 = 2.5, gamma = 5.0;

    double x1 = 1.5, x2 = -0.5, theta_hat = 0.0;

    std::ofstream file("Chapter1_Lesson3_Backstepping.csv");
    file << "t,x1,x2,u,theta_hat\n";

    for (int k = 0; k < n; ++k) {
        double t = k * dt;
        double z1 = x1;
        double z2 = x2 + c1 * x1;
        double u = -theta_hat * z1 - z1 - c2 * z2 - c1 * z2 + c1 * c1 * z1;
        double theta_hat_dot = gamma * z1 * z2;

        double x1_dot = x2;
        double x2_dot = theta_true * x1 + u;

        x1 += dt * x1_dot;
        x2 += dt * x2_dot;
        theta_hat += dt * theta_hat_dot;

        file << t << "," << x1 << "," << x2 << "," << u << "," << theta_hat << "\n";
    }
}

void simulate_gain_scheduling() {
    const double dt = 0.001, tf = 8.0;
    const int n = static_cast<int>(tf / dt);
    double x = 0.0;

    std::ofstream file("Chapter1_Lesson3_GainScheduling.csv");
    file << "t,x,rho,K\n";

    for (int k = 0; k < n; ++k) {
        double t = k * dt;
        double rho = 0.5 + 0.5 * std::sin(0.7 * t);
        double a_rho = -0.5 + 1.2 * rho;

        double K = 0.0;
        if (rho <= 0.5) {
            K = 1.5 + (2.5 - 1.5) * (rho - 0.0) / 0.5;
        } else {
            K = 2.5 + (3.5 - 2.5) * (rho - 0.5) / 0.5;
        }

        double r = (t > 0.5) ? 1.0 : 0.0;
        double u = -K * x + r;
        double x_dot = a_rho * x + u;
        x += dt * x_dot;

        file << t << "," << x << "," << rho << "," << K << "\n";
    }
}

int main() {
    simulate_mrac();
    simulate_str();
    simulate_adaptive_backstepping();
    simulate_gain_scheduling();

    std::cout << "Generated CSV files:\n";
    std::cout << "  Chapter1_Lesson3_MRAC.csv\n";
    std::cout << "  Chapter1_Lesson3_STR.csv\n";
    std::cout << "  Chapter1_Lesson3_Backstepping.csv\n";
    std::cout << "  Chapter1_Lesson3_GainScheduling.csv\n";
    return 0;
}

      

File: Chapter1_Lesson3.java

// Chapter1_Lesson3.java
// Classes of Adaptive Controllers: MRAC, STR, Adaptive Backstepping, Gain Scheduling
// Compile:
//   javac Chapter1_Lesson3.java
// Run:
//   java Chapter1_Lesson3
//
// This program writes CSV files for plotting.

import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;

public class Chapter1_Lesson3 {

    static double clip(double x, double lo, double hi) {
        return Math.max(lo, Math.min(hi, x));
    }

    static void simulateMRAC() throws IOException {
        double dt = 0.001, tf = 8.0;
        int n = (int)(tf / dt);
        double a = -0.6, b = 1.8, am = -2.0, bm = 2.0;
        double gammaX = 4.0, gammaR = 4.0;
        double P = 1.0 / (-2.0 * am);

        double x = 0.0, xm = 0.0, thetaX = 0.0, thetaR = 0.0;

        try (PrintWriter out = new PrintWriter(new FileWriter("Chapter1_Lesson3_MRAC.csv"))) {
            out.println("t,x,xm,u,theta_x,theta_r");

            for (int k = 0; k < n; k++) {
                double t = k * dt;
                double r = (t > 0.5) ? 1.0 : 0.0;
                double e = x - xm;
                double u = thetaX * x + thetaR * r;

                double thetaXDot = -gammaX * x * e * P * b;
                double thetaRDot = -gammaR * r * e * P * b;

                double xDot = a * x + b * u;
                double xmDot = am * xm + bm * r;

                x += dt * xDot;
                xm += dt * xmDot;
                thetaX += dt * thetaXDot;
                thetaR += dt * thetaRDot;

                out.printf("%.6f,%.10f,%.10f,%.10f,%.10f,%.10f%n",
                        t, x, xm, u, thetaX, thetaR);
            }
        }
    }

    static void simulateSTR() throws IOException {
        double dt = 0.02, tf = 10.0;
        int n = (int)(tf / dt);
        double aTrue = 0.92, bTrue = 0.18, lambda = 0.995;

        double aHat = 0.4, bHat = 0.4;
        double P00 = 100.0, P01 = 0.0, P10 = 0.0, P11 = 100.0;
        double y = 0.0, u = 0.0;

        try (PrintWriter out = new PrintWriter(new FileWriter("Chapter1_Lesson3_STR.csv"))) {
            out.println("t,y,u,a_hat,b_hat");

            for (int k = 0; k < n; k++) {
                double t = k * dt;
                double rNext = (t > 1.0) ? 1.0 : 0.0;
                double yNext = aTrue * y + bTrue * u;

                double phi0 = y;
                double phi1 = u;
                double Pphi0 = P00 * phi0 + P01 * phi1;
                double Pphi1 = P10 * phi0 + P11 * phi1;
                double denom = lambda + phi0 * Pphi0 + phi1 * Pphi1;

                double K0 = Pphi0 / denom;
                double K1 = Pphi1 / denom;

                double pred = phi0 * aHat + phi1 * bHat;
                double err = yNext - pred;
                aHat += K0 * err;
                bHat += K1 * err;

                double oldP00 = P00, oldP01 = P01, oldP10 = P10, oldP11 = P11;
                P00 = (oldP00 - K0 * (phi0 * oldP00 + phi1 * oldP10)) / lambda;
                P01 = (oldP01 - K0 * (phi0 * oldP01 + phi1 * oldP11)) / lambda;
                P10 = (oldP10 - K1 * (phi0 * oldP00 + phi1 * oldP10)) / lambda;
                P11 = (oldP11 - K1 * (phi0 * oldP01 + phi1 * oldP11)) / lambda;

                double bSafe = bHat;
                if (Math.abs(bSafe) < 0.05) {
                    bSafe = (bSafe >= 0.0) ? 0.05 : -0.05;
                }
                u = clip((rNext - aHat * yNext) / bSafe, -8.0, 8.0);
                y = yNext;

                out.printf("%.6f,%.10f,%.10f,%.10f,%.10f%n", t, y, u, aHat, bHat);
            }
        }
    }

    static void simulateAdaptiveBackstepping() throws IOException {
        double dt = 0.001, tf = 8.0;
        int n = (int)(tf / dt);
        double thetaTrue = 1.5, c1 = 2.0, c2 = 2.5, gamma = 5.0;

        double x1 = 1.5, x2 = -0.5, thetaHat = 0.0;

        try (PrintWriter out = new PrintWriter(new FileWriter("Chapter1_Lesson3_Backstepping.csv"))) {
            out.println("t,x1,x2,u,theta_hat");

            for (int k = 0; k < n; k++) {
                double t = k * dt;
                double z1 = x1;
                double z2 = x2 + c1 * x1;
                double u = -thetaHat * z1 - z1 - c2 * z2 - c1 * z2 + c1 * c1 * z1;
                double thetaHatDot = gamma * z1 * z2;

                double x1Dot = x2;
                double x2Dot = thetaTrue * x1 + u;

                x1 += dt * x1Dot;
                x2 += dt * x2Dot;
                thetaHat += dt * thetaHatDot;

                out.printf("%.6f,%.10f,%.10f,%.10f,%.10f%n", t, x1, x2, u, thetaHat);
            }
        }
    }

    static void simulateGainScheduling() throws IOException {
        double dt = 0.001, tf = 8.0;
        int n = (int)(tf / dt);
        double x = 0.0;

        try (PrintWriter out = new PrintWriter(new FileWriter("Chapter1_Lesson3_GainScheduling.csv"))) {
            out.println("t,x,rho,K");

            for (int k = 0; k < n; k++) {
                double t = k * dt;
                double rho = 0.5 + 0.5 * Math.sin(0.7 * t);
                double aRho = -0.5 + 1.2 * rho;

                double K;
                if (rho <= 0.5) {
                    K = 1.5 + (2.5 - 1.5) * rho / 0.5;
                } else {
                    K = 2.5 + (3.5 - 2.5) * (rho - 0.5) / 0.5;
                }

                double r = (t > 0.5) ? 1.0 : 0.0;
                double u = -K * x + r;
                double xDot = aRho * x + u;
                x += dt * xDot;

                out.printf("%.6f,%.10f,%.10f,%.10f%n", t, x, rho, K);
            }
        }
    }

    public static void main(String[] args) throws IOException {
        simulateMRAC();
        simulateSTR();
        simulateAdaptiveBackstepping();
        simulateGainScheduling();
        System.out.println("Generated CSV files for MRAC, STR, adaptive backstepping, and gain scheduling.");
    }
}

      

File: Chapter1_Lesson3.m

% Chapter1_Lesson3.m
% Classes of Adaptive Controllers: MRAC, STR, Adaptive Backstepping, Gain Scheduling
% Run in MATLAB or GNU Octave:
%   Chapter1_Lesson3

clear; close all; clc;

%% Direct MRAC
dt = 0.001; tf = 8.0; n = floor(tf/dt);
t = linspace(0, tf, n);
a = -0.6; b = 1.8; am = -2.0; bm = 2.0;
gamma_x = 4.0; gamma_r = 4.0;
P = 1.0 / (-2.0 * am);

x = 0; xm = 0; theta_x = 0; theta_r = 0;
X = zeros(1,n); XM = zeros(1,n); THX = zeros(1,n); THR = zeros(1,n);

for k = 1:n
    r = double(t(k) > 0.5);
    e = x - xm;
    u = theta_x * x + theta_r * r;

    theta_x_dot = -gamma_x * x * e * P * b;
    theta_r_dot = -gamma_r * r * e * P * b;

    x_dot = a * x + b * u;
    xm_dot = am * xm + bm * r;

    x = x + dt * x_dot;
    xm = xm + dt * xm_dot;
    theta_x = theta_x + dt * theta_x_dot;
    theta_r = theta_r + dt * theta_r_dot;

    X(k) = x; XM(k) = xm; THX(k) = theta_x; THR(k) = theta_r;
end

figure; plot(t, X, t, XM, '--', t, THX, t, THR);
grid on; xlabel('time [s]'); title('Direct MRAC');
legend('x', 'x_m', 'theta_x', 'theta_r');

%% Self-Tuning Regulator with RLS
dt = 0.02; tf = 10.0; n = floor(tf/dt);
t = linspace(0, tf, n);
a_true = 0.92; b_true = 0.18;
theta = [0.4; 0.4];
Pmat = 100 * eye(2);
lambda = 0.995;
y = 0; u = 0;
Y = zeros(1,n); U = zeros(1,n); AH = zeros(1,n); BH = zeros(1,n);

for k = 1:n
    r_next = double(t(k) > 1.0);
    y_next = a_true * y + b_true * u;

    phi = [y; u];
    K = (Pmat * phi) / (lambda + phi' * Pmat * phi);
    eps = y_next - phi' * theta;
    theta = theta + K * eps;
    Pmat = (Pmat - K * phi' * Pmat) / lambda;

    b_safe = theta(2);
    if abs(b_safe) < 0.05
        b_safe = 0.05 * sign(b_safe + 1e-12);
    end
    u = (r_next - theta(1) * y_next) / b_safe;
    u = min(max(u, -8), 8);

    y = y_next;
    Y(k) = y; U(k) = u; AH(k) = theta(1); BH(k) = theta(2);
end

figure; plot(t, Y, t, AH, t, BH);
grid on; xlabel('time [s]'); title('STR with RLS');
legend('y', 'a hat', 'b hat');

%% Adaptive Backstepping
dt = 0.001; tf = 8.0; n = floor(tf/dt);
t = linspace(0, tf, n);
theta_true = 1.5; c1 = 2.0; c2 = 2.5; gamma = 5.0;

x1 = 1.5; x2 = -0.5; theta_hat = 0.0;
X1 = zeros(1,n); X2 = zeros(1,n); TH = zeros(1,n);

for k = 1:n
    z1 = x1;
    z2 = x2 + c1 * x1;
    u = -theta_hat * z1 - z1 - c2 * z2 - c1 * z2 + c1^2 * z1;
    theta_hat_dot = gamma * z1 * z2;

    x1_dot = x2;
    x2_dot = theta_true * x1 + u;

    x1 = x1 + dt * x1_dot;
    x2 = x2 + dt * x2_dot;
    theta_hat = theta_hat + dt * theta_hat_dot;

    X1(k) = x1; X2(k) = x2; TH(k) = theta_hat;
end

figure; plot(t, X1, t, X2, t, TH);
grid on; xlabel('time [s]'); title('Adaptive Backstepping');
legend('x1', 'x2', 'theta hat');

%% Gain Scheduling
dt = 0.001; tf = 8.0; n = floor(tf/dt);
t = linspace(0, tf, n);
x = 0; XG = zeros(1,n); RHO = zeros(1,n); KG = zeros(1,n);
rho_grid = [0.0 0.5 1.0];
K_grid = [1.5 2.5 3.5];

for k = 1:n
    rho = 0.5 + 0.5 * sin(0.7 * t(k));
    a_rho = -0.5 + 1.2 * rho;
    K = interp1(rho_grid, K_grid, rho, 'linear');
    r = double(t(k) > 0.5);
    u = -K * x + r;
    x_dot = a_rho * x + u;
    x = x + dt * x_dot;

    XG(k) = x; RHO(k) = rho; KG(k) = K;
end

figure; plot(t, XG, t, RHO, t, KG);
grid on; xlabel('time [s]'); title('Gain Scheduling');
legend('x', 'rho', 'K(rho)');

      

For a Simulink version, implement the same equations with Integrator blocks. For MRAC, use one integrator for \( x \), one for \( x_m \), and one integrator for each adaptive parameter. For STR, implement the RLS matrix update in a MATLAB Function block with persistent variables.

File: Chapter1_Lesson3.nb

(* Chapter1_Lesson3.nb *)
(* Classes of Adaptive Controllers: MRAC, STR, Adaptive Backstepping, Gain Scheduling *)
(* This is Wolfram Language code saved with .nb extension for direct copying into a Mathematica notebook. *)

ClearAll["Global`*"];

(* Direct MRAC *)
mracSol = NDSolve[
  {
    x'[t] == -0.6 x[t] + 1.8 (thetax[t] x[t] + thetar[t] UnitStep[t - 0.5]),
    xm'[t] == -2.0 xm[t] + 2.0 UnitStep[t - 0.5],
    thetax'[t] == -4.0 x[t] (x[t] - xm[t]) 0.25 1.8,
    thetar'[t] == -4.0 UnitStep[t - 0.5] (x[t] - xm[t]) 0.25 1.8,
    x[0] == 0, xm[0] == 0, thetax[0] == 0, thetar[0] == 0
  },
  {x, xm, thetax, thetar}, {t, 0, 8}
];

Plot[
  Evaluate[{x[t], xm[t], thetax[t], thetar[t]} /. mracSol],
  {t, 0, 8},
  PlotLegends -> {"x", "xm", "thetax", "thetar"},
  PlotLabel -> "Direct MRAC"
]

(* Self-Tuning Regulator with RLS *)
strData = Module[
  {dt = 0.02, n = 500, aTrue = 0.92, bTrue = 0.18, theta = {0.4, 0.4},
   p = 100 IdentityMatrix[2], lambda = 0.995, y = 0., u = 0.,
   rows = {}, phi, gain, eps, ynext, rnext, bsafe},
  Do[
    rnext = If[k dt > 1, 1., 0.];
    ynext = aTrue y + bTrue u;
    phi = {y, u};
    gain = (p.phi)/(lambda + phi.p.phi);
    eps = ynext - phi.theta;
    theta = theta + gain eps;
    p = (p - Outer[Times, gain, phi].p)/lambda;
    bsafe = If[Abs[theta[[2]]] < 0.05, 0.05, theta[[2]]];
    u = Clip[(rnext - theta[[1]] ynext)/bsafe, {-8, 8}];
    y = ynext;
    AppendTo[rows, {k dt, y, theta[[1]], theta[[2]]}],
    {k, 1, n}
  ];
  rows
];

ListLinePlot[
  {strData[[All, {1, 2}]], strData[[All, {1, 3}]], strData[[All, {1, 4}]]},
  PlotLegends -> {"y", "aHat", "bHat"},
  PlotLabel -> "STR with RLS"
]

(* Adaptive Backstepping *)
backSol = NDSolve[
  {
    x1'[t] == x2[t],
    x2'[t] == 1.5 x1[t] +
      (-thetahat[t] x1[t] - x1[t] - 2.5 (x2[t] + 2 x1[t]) -
        2 (x2[t] + 2 x1[t]) + 4 x1[t]),
    thetahat'[t] == 5 x1[t] (x2[t] + 2 x1[t]),
    x1[0] == 1.5, x2[0] == -0.5, thetahat[0] == 0
  },
  {x1, x2, thetahat}, {t, 0, 8}
];

Plot[
  Evaluate[{x1[t], x2[t], thetahat[t]} /. backSol],
  {t, 0, 8},
  PlotLegends -> {"x1", "x2", "thetaHat"},
  PlotLabel -> "Adaptive Backstepping"
]

(* Gain Scheduling *)
Kfun = Interpolation[{ {0, 1.5}, {0.5, 2.5}, {1, 3.5} }, InterpolationOrder -> 1];

gainSol = NDSolve[
  {
    rho[t] == 0.5 + 0.5 Sin[0.7 t],
    xg'[t] == (-0.5 + 1.2 rho[t]) xg[t] + (-Kfun[rho[t]] xg[t] + UnitStep[t - 0.5]),
    xg[0] == 0
  },
  xg, {t, 0, 8}
];

Plot[
  Evaluate[xg[t] /. gainSol],
  {t, 0, 8},
  PlotLabel -> "Gain-Scheduled Response"
]

      

Problem 1 (MRAC Error Equation): For the scalar plant \( \dot{x}=ax+bu \), reference model \( \dot{x}_m=a_mx_m+b_mr \), and controller \( u=\hat{\theta}_x x+\hat{\theta}_r r \), derive the tracking-error dynamics.

Solution: Define \( e=x-x_m \). Then

\[ \dot{e}=\dot{x}-\dot{x}_m =ax+b\hat{\theta}_x x+b\hat{\theta}_r r-a_mx_m-b_mr. \]

Add and subtract the matched ideal terms \( b\theta_x^*x \) and \( b\theta_r^*r \), using \( a+b\theta_x^*=a_m \) and \( b\theta_r^*=b_m \). This gives

\[ \dot{e}=a_m(x-x_m)+b(\hat{\theta}_x-\theta_x^*)x+ b(\hat{\theta}_r-\theta_r^*)r =a_me+b\tilde{\theta}_x x+b\tilde{\theta}_r r. \]

Problem 2 (Lyapunov Cancellation in MRAC): For \( b > 0 \) and \( a_m < 0 \), show why \( \dot{\hat{\theta} }_x=-\gamma_x bPex \) and \( \dot{\hat{\theta} }_r=-\gamma_r bPer \) are natural.

Solution: Use

\[ V=Pe^2+\frac{b}{\gamma_x}\tilde{\theta}_x^2+ \frac{b}{\gamma_r}\tilde{\theta}_r^2,\qquad 2a_mP=-1. \]

Since \( \dot{\tilde{\theta} }=\dot{\hat{\theta} } \) for constant ideal parameters,

\[ \dot{V}=2Pe(a_me+b\tilde{\theta}_x x+b\tilde{\theta}_r r) +2\frac{b}{\gamma_x}\tilde{\theta}_x\dot{\hat{\theta} }_x +2\frac{b}{\gamma_r}\tilde{\theta}_r\dot{\hat{\theta} }_r. \]

The selected adaptive laws cancel the mixed terms involving \( e\tilde{\theta}_x x \) and \( e\tilde{\theta}_r r \), leaving

\[ \dot{V}=2Pa_me^2=-e^2\le 0. \]

Problem 3 (STR Certainty Equivalence): For \( y_{k+1}=ay_k+bu_k \), derive the certainty-equivalence control law that enforces \( y_{k+1}=r_{k+1} \) using estimates.

Solution: Replace unknown parameters by estimates:

\[ r_{k+1}=\hat{a}_ky_k+\hat{b}_ku_k. \]

Solving for \( u_k \) gives

\[ u_k=\frac{r_{k+1}-\hat{a}_ky_k}{\hat{b}_k}. \]

In real implementations, one must protect against division by small \( \hat{b}_k \), actuator saturation, and noise.

Problem 4 (Adaptive Backstepping Derivative): For the system in Section 4, verify that the proposed controller and adaptive law give \( \dot{V}=-c_1z_1^2-c_2z_2^2 \).

Solution: From \( z_1=x_1 \) and \( z_2=x_2+c_1x_1 \),

\[ \dot{z}_1=z_2-c_1z_1. \]

Substitution of the controller into \( \dot{z}_2=\dot{x}_2+c_1\dot{x}_1 \) yields

\[ \dot{z}_2=\tilde{\theta}z_1-z_1-c_2z_2. \]

Therefore,

\[ \dot{V}=z_1(z_2-c_1z_1)+z_2(\tilde{\theta}z_1-z_1-c_2z_2) -\frac{1}{\gamma}\tilde{\theta}\dot{\hat{\theta} }. \]

Using \( \dot{\hat{\theta} }=\gamma z_1z_2 \), the parameter cross-term cancels, so

\[ \dot{V}=-c_1z_1^2-c_2z_2^2. \]

Problem 5 (Class Selection): A process has slowly changing operating temperature \( \rho \), and linear controllers have already been designed at several temperature points. No online identification is available. Which class is the most natural first design?

Solution: Gain scheduling is the most natural first design. The controller can interpolate among fixed gains

\[ K(\rho)=\sum_{i=1}^{M}w_i(\rho)K_i. \]

However, the designer should still analyze transitions between operating points. Frozen-point stability alone is not a full proof of scheduled closed-loop stability.

MRAC, STR, adaptive backstepping, and gain scheduling represent four different ways of dealing with changing or uncertain dynamics. MRAC adapts controller parameters to match a reference model. STR identifies a model online and recomputes the controller. Adaptive backstepping recursively constructs stabilizing controllers for nonlinear strict-feedback systems. Gain scheduling changes controller gains according to measured operating variables but is not necessarily true adaptation. These distinctions will guide the detailed mathematical development in later chapters.