After completing this lesson, students should be able to:

• formulate a quadratic tracking-error cost for first-order MRAC;
• derive the MIT update rule from steepest descent in parameter space;
• derive exact sensitivity differential equations for adjustable controller gains;
• construct implementable filtered sensitivity signals;
• select adaptation gains with attention to sign, scaling, and numerical stiffness;
• distinguish gradient descent from a Lyapunov-based stability design; and
• implement the adaptive loop in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The lesson uses only concepts already developed in the course: the first-order uncertain plant, the reference model, the direct MRAC controller structure, tracking error, parameter error, and elementary gradient descent. Lyapunov redesign is intentionally deferred to Lesson 4.

Consider the unknown first-order plant introduced in Lesson 1:

\[ \dot y(t) = -a y(t) + b u(t), \qquad a > 0, \qquad b \neq 0. \]

The desired closed-loop response is specified by the stable reference model:

\[ \dot y_m(t) = -a_m y_m(t) + b_m r(t), \qquad a_m > 0. \]

We use the adjustable direct-control law:

\[ u(t) = \theta_y(t)y(t) + \theta_r(t)r(t). \]

Substitution into the plant produces:

\[ \dot y = \left(-a+b\theta_y\right)y+b\theta_r r. \]

If the plant parameters were known, exact model matching would be obtained with the ideal constant gains:

\[ \theta_y^{\ast}=\frac{a-a_m}{b}, \qquad \theta_r^{\ast}=\frac{b_m}{b}. \]

Define tracking error and parameter errors by:

\[ e = y-y_m, \qquad \widetilde\theta_y=\theta_y-\theta_y^{\ast}, \qquad \widetilde\theta_r=\theta_r-\theta_r^{\ast}. \]

Using the matching identities, the error dynamics may be written:

\[ \dot e = -a_m e +b\widetilde\theta_y y +b\widetilde\theta_r r. \]

This equation shows why adaptation is required: unknown parameter errors act as endogenous forcing terms in the tracking-error dynamics.

The MIT rule begins with the instantaneous quadratic performance index:

\[ J(\boldsymbol\theta,t) =\frac{1}{2}e^2(t), \qquad \boldsymbol\theta= \begin{bmatrix}\theta_y & \theta_r\end{bmatrix}^{T}. \]

The gradient of the cost with respect to an adjustable scalar parameter \( \theta_i \) is:

\[ \frac{\partial J}{\partial\theta_i} = e\frac{\partial e}{\partial\theta_i}. \]

Continuous-time steepest descent moves each parameter opposite to this gradient:

\[ \dot\theta_i = -\gamma_i\frac{\partial J}{\partial\theta_i} = -\gamma_i e \frac{\partial e}{\partial\theta_i}, \qquad \gamma_i > 0. \]

The quantity \( \partial e/\partial\theta_i \) is called a sensitivity derivative. It measures how a small change in a controller parameter changes the tracking error. Therefore, the MIT rule is not merely “error times signal”; it is error multiplied by a model of the error sensitivity.

\[ \boxed{\dot{\boldsymbol\theta} = -\boldsymbol\Gamma e\boldsymbol\psi} \]

where \( \boldsymbol\Gamma =\operatorname{diag}(\gamma_y,\gamma_r) \) and \( \boldsymbol\psi =\partial e/\partial\boldsymbol\theta \).

Because the reference model does not depend on the controller parameters, its parameter derivatives vanish:

\[ \frac{\partial y_m}{\partial\theta_y}=0, \qquad \frac{\partial y_m}{\partial\theta_r}=0. \]

Hence the error sensitivities equal the plant-output sensitivities:

\[ s_y \equiv \frac{\partial e}{\partial\theta_y} = \frac{\partial y}{\partial\theta_y}, \qquad s_r \equiv \frac{\partial e}{\partial\theta_r} = \frac{\partial y}{\partial\theta_r}. \]

Differentiate the closed-loop plant equation with respect to \( \theta_y \). During the local sensitivity calculation, the current parameter values are treated as frozen:

\[ \begin{aligned} \dot s_y &= \frac{\partial}{\partial\theta_y} \left[ (-a+b\theta_y)y+b\theta_r r \right] \\ &= (-a+b\theta_y)s_y+by. \end{aligned} \]

Similarly:

\[ \dot s_r = (-a+b\theta_y)s_r+br. \]

With parameter-independent initial plant state, the natural sensitivity initial conditions are:

\[ s_y(0)=0, \qquad s_r(0)=0. \]

The exact MIT law is therefore:

\[ \dot\theta_y=-\gamma_y e s_y, \qquad \dot\theta_r=-\gamma_r e s_r. \]

This exact form is conceptually important but often impractical because the sensitivity equations contain the unknown plant coefficients \( a \) and \( b \). A practical controller therefore replaces them with a stable approximate sensitivity model.

Introduce the differential operator \( p=d/dt \). For frozen parameters, the exact sensitivities may be represented formally as:

\[ s_y = \frac{b}{p+a-b\theta_y}y, \qquad s_r = \frac{b}{p+a-b\theta_y}r. \]

Near successful model matching:

\[ a-b\theta_y \approx a_m. \]

Thus a practical approximation is:

\[ s_y \approx \frac{b}{p+a_m}y, \qquad s_r \approx \frac{b}{p+a_m}r. \]

When the sign of \( b \) is known, its magnitude can be absorbed into the positive adaptation gains. Define filtered regressors:

\[ \phi_y=\frac{1}{p+a_m}y, \qquad \phi_r=\frac{1}{p+a_m}r. \]

Their realizable differential equations are:

\[ \dot\phi_y=-a_m\phi_y+y, \qquad \dot\phi_r=-a_m\phi_r+r. \]

For known \( \operatorname{sgn}(b) \), one useful convention is:

\[ \dot\theta_y = -\gamma_y\operatorname{sgn}(b)e\phi_y, \qquad \dot\theta_r = -\gamma_r\operatorname{sgn}(b)e\phi_r. \]

In the numerical examples below, \( b > 0 \), so \( \operatorname{sgn}(b)=1 \). The implemented rule is:

\[ \boxed{ \dot\theta_y=-\gamma_y e\phi_y, \qquad \dot\theta_r=-\gamma_r e\phi_r } \]

The complete first-order MIT-rule MRAC contains four coupled subsystems: reference model, adjustable controller, uncertain plant, and adaptation mechanism.

flowchart TD
  R["Command r"] --> RM["Reference model: y_m"]
  R --> C["Controller: u = theta_y*y + theta_r*r"]
  C --> P["Unknown first-order plant"]
  P --> Y["Plant output y"]
  Y --> C
  Y --> E["Tracking error: e = y - y_m"]
  RM --> E
  Y --> FY["Stable filter: phi_y"]
  R --> FR["Stable filter: phi_r"]
  E --> A["Gradient adaptation"]
  FY --> A
  FR --> A
  A --> TY["Update theta_y"]
  A --> TR["Update theta_r"]
  TY --> C
  TR --> C
        

The reference model is not used merely to produce a target output. Its stable denominator also supplies a practical approximation to the unavailable plant sensitivity dynamics.

Define the controller regressor and filtered regressor vectors:

\[ \boldsymbol\omega= \begin{bmatrix}y\\r\end{bmatrix}, \qquad \boldsymbol\phi= \frac{1}{p+a_m}\boldsymbol\omega = \begin{bmatrix}\phi_y\\\phi_r\end{bmatrix}. \]

The control and adaptation laws become:

\[ u=\boldsymbol\theta^T\boldsymbol\omega, \qquad \dot{\boldsymbol\theta} = -\boldsymbol\Gamma e\boldsymbol\phi. \]

With exact sensitivities, the parameter velocity is antiparallel to the local gradient:

\[ \dot{\boldsymbol\theta} = -\boldsymbol\Gamma \nabla_{\boldsymbol\theta}J. \]

For an infinitesimal frozen-time parameter displacement \( d\boldsymbol\theta =-\boldsymbol\Gamma\nabla J\,dt \), the corresponding first-order cost variation is:

\[ dJ = \nabla J^T d\boldsymbol\theta = -\nabla J^T\boldsymbol\Gamma\nabla J\,dt \leq 0. \]

This is a local optimization statement. It is not a proof that the total derivative \( \dot J \) of the coupled dynamic system is always nonpositive, because plant and model state evolution also contribute to \( \dot e \).

Adaptation gains are not dimensionless tuning knobs. From \( \dot\theta_i=-\gamma_i e\phi_i \):

\[ [\gamma_i] = \frac{[\theta_i]}{[e][\phi_i][t]}. \]

Therefore, changing sensor units or reference amplitude changes the effective adaptation rate. A practical scaling procedure is:

\[ \bar y=\frac{y}{Y_0}, \qquad \bar r=\frac{r}{R_0}, \qquad \bar e=\frac{e}{Y_0}. \]

where \(Y_0\) and \(R_0\) are representative amplitudes. Adaptation gains can then be tuned in the normalized coordinates and transformed back consistently.

A common practical normalized gradient law is:

\[ \dot{\boldsymbol\theta} = -\frac{\boldsymbol\Gamma e\boldsymbol\phi} {\varepsilon+\boldsymbol\phi^T\boldsymbol\phi}, \qquad \varepsilon > 0. \]

This normalization limits very large parameter rates when the filtered regressor becomes large. Formal projection and normalization properties are developed later in Chapter 9; here the expression is introduced only as a gain-scaling preview.

The MIT rule is attractive because it is systematic, intuitive, and easy to implement. Nevertheless, several distinctions are essential:

• Gradient descent is local. Exact cost sensitivities are replaced by approximate filtered sensitivities.
• Cost reduction is not automatically stability. The adaptive parameters are part of a nonlinear time-varying feedback loop.
• Fast adaptation can destabilize numerical or physical dynamics. Large gains may excite neglected dynamics, noise, delays, and saturation.
• Tracking does not imply parameter identification. Several gain combinations may produce similar closed-loop behavior for a non-exciting command.
• The sign of the high-frequency gain matters. A wrong assumed sign reverses the descent direction.

The next lesson replaces the heuristic gradient argument with a Lyapunov-based adaptive-law design for the same first-order plant. That redesign will produce an explicit stability certificate under stated assumptions.

The implementations use:

\[ a=1.2, \quad b=0.8, \quad a_m=2, \quad b_m=2, \quad \gamma_y=\gamma_r=2. \]

The command contains constant and sinusoidal components:

\[ r(t)=1+0.35\sin(0.5t). \]

The ideal matching gains are:

\[ \theta_y^{\ast} = \frac{1.2-2}{0.8} =-1, \qquad \theta_r^{\ast} = \frac{2}{0.8} =2.5. \]

The command is deliberately not a pure constant. The sinusoidal component supplies additional excitation and makes the parameter trajectories more informative, although this finite experiment alone does not establish a persistent-excitation result.

flowchart TD
  S["Initialize y, y_m, filters, and parameters"] --> R["Evaluate command r(t)"]
  R --> U["Compute u = theta_y*y + theta_r*r"]
  U --> D["Integrate plant and reference model"]
  D --> E["Compute e = y - y_m"]
  E --> F["Integrate stable sensitivity filters"]
  F --> G["Apply gradient parameter updates"]
  G --> L["Log y, y_m, e, u, theta_y, theta_r"]
  L --> Q["Final time reached?"]
  Q -->|No| R
  Q -->|Yes| O["Plot results and calculate RMS error"]
        

Because the implementation uses approximate sensitivities and a finite command, good tracking should not be interpreted as a guarantee that both parameter estimates converge exactly to their ideal values.

Chapter5_Lesson3.py

This implementation uses NumPy, Matplotlib, and a from-scratch fourth-order Runge–Kutta integrator. It writes a CSV file and displays separate tracking, error, and parameter plots.

"""
Chapter5_Lesson3.py
MIT-rule MRAC simulation for a first-order plant.

Dependencies:
    numpy
    matplotlib
"""

from __future__ import annotations

import csv
from dataclasses import dataclass
from typing import Callable

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Parameters:
    plant_a: float = 1.2
    plant_b: float = 0.8
    model_a: float = 2.0
    model_b: float = 2.0
    gamma_y: float = 2.0
    gamma_r: float = 2.0
    dt: float = 0.002
    final_time: float = 30.0


def reference(t: float) -> float:
    """Bounded command with a constant and sinusoidal component."""
    return 1.0 + 0.35 * np.sin(0.5 * t)


def dynamics(t: float, x: np.ndarray, p: Parameters) -> np.ndarray:
    """
    State ordering:
        x = [y, y_m, phi_y, phi_r, theta_y, theta_r]
    """
    y, y_m, phi_y, phi_r, theta_y, theta_r = x
    r = reference(t)

    u = theta_y * y + theta_r * r
    error = y - y_m

    y_dot = -p.plant_a * y + p.plant_b * u
    y_m_dot = -p.model_a * y_m + p.model_b * r

    # Approximate MIT sensitivities generated by stable reference-model filters.
    phi_y_dot = -p.model_a * phi_y + y
    phi_r_dot = -p.model_a * phi_r + r

    theta_y_dot = -p.gamma_y * error * phi_y
    theta_r_dot = -p.gamma_r * error * phi_r

    return np.array(
        [y_dot, y_m_dot, phi_y_dot, phi_r_dot, theta_y_dot, theta_r_dot],
        dtype=float,
    )


def rk4_step(
    f: Callable[[float, np.ndarray, Parameters], np.ndarray],
    t: float,
    x: np.ndarray,
    dt: float,
    p: Parameters,
) -> np.ndarray:
    k1 = f(t, x, p)
    k2 = f(t + 0.5 * dt, x + 0.5 * dt * k1, p)
    k3 = f(t + 0.5 * dt, x + 0.5 * dt * k2, p)
    k4 = f(t + dt, x + dt * k3, p)
    return x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(p: Parameters) -> dict[str, np.ndarray]:
    steps = int(round(p.final_time / p.dt)) + 1
    time = np.linspace(0.0, p.final_time, steps)
    state = np.zeros((steps, 6), dtype=float)
    control = np.zeros(steps, dtype=float)
    command = np.zeros(steps, dtype=float)

    for k, t in enumerate(time):
        y, _, _, _, theta_y, theta_r = state[k]
        command[k] = reference(t)
        control[k] = theta_y * y + theta_r * command[k]

        if k + 1 < steps:
            state[k + 1] = rk4_step(dynamics, t, state[k], p.dt, p)

            if not np.all(np.isfinite(state[k + 1])):
                raise FloatingPointError(
                    "The adaptive simulation diverged numerically. "
                    "Reduce the adaptation gains or integration step."
                )

    return {
        "time": time,
        "state": state,
        "reference": command,
        "control": control,
    }


def save_csv(result: dict[str, np.ndarray], filename: str) -> None:
    time = result["time"]
    state = result["state"]
    reference_signal = result["reference"]
    control = result["control"]

    with open(filename, "w", newline="", encoding="utf-8") as csv_file:
        writer = csv.writer(csv_file)
        writer.writerow(
            [
                "time",
                "reference",
                "plant_output",
                "model_output",
                "tracking_error",
                "control",
                "theta_y",
                "theta_r",
                "phi_y",
                "phi_r",
            ]
        )
        for k, t in enumerate(time):
            y, y_m, phi_y, phi_r, theta_y, theta_r = state[k]
            writer.writerow(
                [
                    t,
                    reference_signal[k],
                    y,
                    y_m,
                    y - y_m,
                    control[k],
                    theta_y,
                    theta_r,
                    phi_y,
                    phi_r,
                ]
            )


def plot_result(result: dict[str, np.ndarray], p: Parameters) -> None:
    time = result["time"]
    state = result["state"]
    y = state[:, 0]
    y_m = state[:, 1]
    theta_y = state[:, 4]
    theta_r = state[:, 5]
    error = y - y_m

    ideal_theta_y = (p.plant_a - p.model_a) / p.plant_b
    ideal_theta_r = p.model_b / p.plant_b

    plt.figure()
    plt.plot(time, y, label="Plant output y")
    plt.plot(time, y_m, "--", label="Reference-model output y_m")
    plt.xlabel("Time [s]")
    plt.ylabel("Output")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(time, error)
    plt.xlabel("Time [s]")
    plt.ylabel("Tracking error e = y - y_m")
    plt.grid(True)
    plt.tight_layout()

    plt.figure()
    plt.plot(time, theta_y, label="theta_y")
    plt.plot(time, theta_r, label="theta_r")
    plt.axhline(ideal_theta_y, linestyle="--", label="ideal theta_y")
    plt.axhline(ideal_theta_r, linestyle="--", label="ideal theta_r")
    plt.xlabel("Time [s]")
    plt.ylabel("Adaptive parameters")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


def main() -> None:
    parameters = Parameters()
    result = simulate(parameters)
    save_csv(result, "Chapter5_Lesson3_results.csv")

    final_state = result["state"][-1]
    rms_error = float(
        np.sqrt(np.mean((result["state"][:, 0] - result["state"][:, 1]) ** 2))
    )
    print(f"Final theta_y = {final_state[4]:.6f}")
    print(f"Final theta_r = {final_state[5]:.6f}")
    print(f"RMS tracking error = {rms_error:.6f}")
    print("Saved Chapter5_Lesson3_results.csv")

    plot_result(result, parameters)


if __name__ == "__main__":
    main()

Chapter5_Lesson3.cpp

The C++17 version uses only the standard library and writes simulation results to CSV for plotting in any numerical environment.

// Chapter5_Lesson3.cpp
// MIT-rule MRAC simulation for a first-order plant.
// Build: g++ -std=c++17 -O2 Chapter5_Lesson3.cpp -o Chapter5_Lesson3

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

struct Parameters {
    double plantA = 1.2;
    double plantB = 0.8;
    double modelA = 2.0;
    double modelB = 2.0;
    double gammaY = 2.0;
    double gammaR = 2.0;
    double dt = 0.002;
    double finalTime = 30.0;
};

using State = std::array<double, 6>;

double referenceSignal(double t) {
    return 1.0 + 0.35 * std::sin(0.5 * t);
}

State addScaled(const State& x, const State& dx, double scale) {
    State result{};
    for (std::size_t i = 0; i < x.size(); ++i) {
        result[i] = x[i] + scale * dx[i];
    }
    return result;
}

State dynamics(double t, const State& x, const Parameters& p) {
    const double y = x[0];
    const double yModel = x[1];
    const double phiY = x[2];
    const double phiR = x[3];
    const double thetaY = x[4];
    const double thetaR = x[5];

    const double r = referenceSignal(t);
    const double u = thetaY * y + thetaR * r;
    const double error = y - yModel;

    return State{
        -p.plantA * y + p.plantB * u,
        -p.modelA * yModel + p.modelB * r,
        -p.modelA * phiY + y,
        -p.modelA * phiR + r,
        -p.gammaY * error * phiY,
        -p.gammaR * error * phiR
    };
}

State rk4Step(double t, const State& x, const Parameters& p) {
    const State k1 = dynamics(t, x, p);
    const State k2 = dynamics(t + 0.5 * p.dt, addScaled(x, k1, 0.5 * p.dt), p);
    const State k3 = dynamics(t + 0.5 * p.dt, addScaled(x, k2, 0.5 * p.dt), p);
    const State k4 = dynamics(t + p.dt, addScaled(x, k3, p.dt), p);

    State next{};
    for (std::size_t i = 0; i < x.size(); ++i) {
        next[i] = x[i] + (p.dt / 6.0) *
            (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        if (!std::isfinite(next[i])) {
            throw std::runtime_error(
                "Non-finite state detected. Reduce adaptation gains or time step."
            );
        }
    }
    return next;
}

int main() {
    try {
        const Parameters p;
        State x{0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

        std::ofstream csv("Chapter5_Lesson3_results.csv");
        if (!csv) {
            throw std::runtime_error("Cannot create output CSV file.");
        }

        csv << "time,reference,plant_output,model_output,tracking_error,"
               "control,theta_y,theta_r,phi_y,phi_r\n";
        csv << std::setprecision(12);

        double squaredErrorIntegral = 0.0;
        const int steps = static_cast<int>(std::llround(p.finalTime / p.dt));

        for (int k = 0; k <= steps; ++k) {
            const double t = k * p.dt;
            const double r = referenceSignal(t);
            const double u = x[4] * x[0] + x[5] * r;
            const double error = x[0] - x[1];

            csv << t << ',' << r << ',' << x[0] << ',' << x[1] << ','
                << error << ',' << u << ',' << x[4] << ',' << x[5] << ','
                << x[2] << ',' << x[3] << '\n';

            squaredErrorIntegral += error * error * p.dt;
            if (k < steps) {
                x = rk4Step(t, x, p);
            }
        }

        const double rmsError = std::sqrt(squaredErrorIntegral / p.finalTime);
        std::cout << "Final theta_y = " << x[4] << '\n';
        std::cout << "Final theta_r = " << x[5] << '\n';
        std::cout << "RMS tracking error = " << rmsError << '\n';
        std::cout << "Saved Chapter5_Lesson3_results.csv\n";
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

Chapter5_Lesson3.java

The Java implementation uses the standard library, an explicit RK4 solver, finite-value checks, and CSV output.

// Chapter5_Lesson3.java
// MIT-rule MRAC simulation for a first-order plant.
// Build and run:
//   javac Chapter5_Lesson3.java
//   java Chapter5_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 Chapter5_Lesson3 {
    private static final class Parameters {
        final double plantA = 1.2;
        final double plantB = 0.8;
        final double modelA = 2.0;
        final double modelB = 2.0;
        final double gammaY = 2.0;
        final double gammaR = 2.0;
        final double dt = 0.002;
        final double finalTime = 30.0;
    }

    private static double referenceSignal(double t) {
        return 1.0 + 0.35 * Math.sin(0.5 * t);
    }

    private static double[] dynamics(double t, double[] x, Parameters p) {
        double y = x[0];
        double yModel = x[1];
        double phiY = x[2];
        double phiR = x[3];
        double thetaY = x[4];
        double thetaR = x[5];

        double r = referenceSignal(t);
        double u = thetaY * y + thetaR * r;
        double error = y - yModel;

        return new double[] {
            -p.plantA * y + p.plantB * u,
            -p.modelA * yModel + p.modelB * r,
            -p.modelA * phiY + y,
            -p.modelA * phiR + r,
            -p.gammaY * error * phiY,
            -p.gammaR * error * phiR
        };
    }

    private static double[] addScaled(double[] x, double[] dx, double scale) {
        double[] result = new double[x.length];
        for (int i = 0; i < x.length; i++) {
            result[i] = x[i] + scale * dx[i];
        }
        return result;
    }

    private static double[] rk4Step(double t, double[] x, Parameters p) {
        double[] k1 = dynamics(t, x, p);
        double[] k2 = dynamics(t + 0.5 * p.dt, addScaled(x, k1, 0.5 * p.dt), p);
        double[] k3 = dynamics(t + 0.5 * p.dt, addScaled(x, k2, 0.5 * p.dt), p);
        double[] k4 = dynamics(t + p.dt, addScaled(x, k3, p.dt), p);

        double[] next = new double[x.length];
        for (int i = 0; i < x.length; i++) {
            next[i] = x[i] + (p.dt / 6.0)
                * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
            if (!Double.isFinite(next[i])) {
                throw new ArithmeticException(
                    "Non-finite state detected. Reduce adaptation gains or time step."
                );
            }
        }
        return next;
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        Parameters p = new Parameters();
        double[] x = new double[] {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
        int steps = (int) Math.round(p.finalTime / p.dt);
        double squaredErrorIntegral = 0.0;

        Path output = Path.of("Chapter5_Lesson3_results.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(
                output, StandardCharsets.UTF_8)) {
            writer.write(
                "time,reference,plant_output,model_output,tracking_error,"
                + "control,theta_y,theta_r,phi_y,phi_r\n"
            );

            for (int k = 0; k <= steps; k++) {
                double t = k * p.dt;
                double r = referenceSignal(t);
                double u = x[4] * x[0] + x[5] * r;
                double error = x[0] - x[1];

                writer.write(String.format(
                    Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,"
                    + "%.12f,%.12f,%.12f,%.12f%n",
                    t, r, x[0], x[1], error, u, x[4], x[5], x[2], x[3]
                ));

                squaredErrorIntegral += error * error * p.dt;
                if (k < steps) {
                    x = rk4Step(t, x, p);
                }
            }
        } catch (IOException exception) {
            System.err.println("File error: " + exception.getMessage());
            System.exit(1);
        }

        double rmsError = Math.sqrt(squaredErrorIntegral / p.finalTime);
        System.out.printf(Locale.US, "Final theta_y = %.6f%n", x[4]);
        System.out.printf(Locale.US, "Final theta_r = %.6f%n", x[5]);
        System.out.printf(Locale.US, "RMS tracking error = %.6f%n", rmsError);
        System.out.println("Saved Chapter5_Lesson3_results.csv");
    }
}

Chapter5_Lesson3.m

The MATLAB script uses ode45, calculates RMS tracking error, exports a table, and plots the main adaptive signals.

% Chapter5_Lesson3.m
% MIT-rule MRAC simulation for a first-order plant.
% Requires base MATLAB only.

clear; clc; close all;

p.plantA = 1.2;
p.plantB = 0.8;
p.modelA = 2.0;
p.modelB = 2.0;
p.gammaY = 2.0;
p.gammaR = 2.0;
p.finalTime = 30.0;

x0 = zeros(6, 1);
solverOptions = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[t, x] = ode45(@(time, state) mracDynamics(time, state, p), ...
               [0, p.finalTime], x0, solverOptions);

y = x(:, 1);
yModel = x(:, 2);
phiY = x(:, 3);
phiR = x(:, 4);
thetaY = x(:, 5);
thetaR = x(:, 6);
r = 1.0 + 0.35 .* sin(0.5 .* t);
u = thetaY .* y + thetaR .* r;
errorSignal = y - yModel;

idealThetaY = (p.plantA - p.modelA) / p.plantB;
idealThetaR = p.modelB / p.plantB;
rmsError = sqrt(trapz(t, errorSignal.^2) / p.finalTime);

results = table(t, r, y, yModel, errorSignal, u, thetaY, thetaR, phiY, phiR, ...
    'VariableNames', {'time', 'reference', 'plant_output', 'model_output', ...
    'tracking_error', 'control', 'theta_y', 'theta_r', 'phi_y', 'phi_r'});
writetable(results, 'Chapter5_Lesson3_results.csv');

fprintf('Final theta_y = %.6f\n', thetaY(end));
fprintf('Final theta_r = %.6f\n', thetaR(end));
fprintf('RMS tracking error = %.6f\n', rmsError);
fprintf('Saved Chapter5_Lesson3_results.csv\n');

figure;
plot(t, y, 'LineWidth', 1.4); hold on;
plot(t, yModel, '--', 'LineWidth', 1.4);
grid on;
xlabel('Time [s]');
ylabel('Output');
legend('Plant output y', 'Reference-model output y_m', 'Location', 'best');
title('MIT-Rule MRAC Tracking');

figure;
plot(t, errorSignal, 'LineWidth', 1.4);
grid on;
xlabel('Time [s]');
ylabel('e = y - y_m');
title('Tracking Error');

figure;
plot(t, thetaY, 'LineWidth', 1.4); hold on;
plot(t, thetaR, 'LineWidth', 1.4);
yline(idealThetaY, '--');
yline(idealThetaR, '--');
grid on;
xlabel('Time [s]');
ylabel('Adaptive parameters');
legend('\theta_y', '\theta_r', '\theta_y^*', '\theta_r^*', ...
       'Location', 'best');
title('Adaptive Parameter Evolution');

function dx = mracDynamics(time, x, p)
    y = x(1);
    yModel = x(2);
    phiY = x(3);
    phiR = x(4);
    thetaY = x(5);
    thetaR = x(6);

    r = 1.0 + 0.35 * sin(0.5 * time);
    u = thetaY * y + thetaR * r;
    errorSignal = y - yModel;

    dx = zeros(6, 1);
    dx(1) = -p.plantA * y + p.plantB * u;
    dx(2) = -p.modelA * yModel + p.modelB * r;
    dx(3) = -p.modelA * phiY + y;
    dx(4) = -p.modelA * phiR + r;
    dx(5) = -p.gammaY * errorSignal * phiY;
    dx(6) = -p.gammaR * errorSignal * phiR;
end

Chapter5_Lesson3_Simulink.m

This script programmatically creates and simulates a Simulink block diagram. It requires Simulink and saves Chapter5_Lesson3_MIT_MRAC.slx.

% Chapter5_Lesson3_Simulink.m
% Programmatically builds a Simulink implementation of the MIT-rule MRAC.
% Requires Simulink.

clear; clc;

modelName = 'Chapter5_Lesson3_MIT_MRAC';
if bdIsLoaded(modelName)
    close_system(modelName, 0);
end
if exist([modelName '.slx'], 'file')
    delete([modelName '.slx']);
end

load_system('simulink');
new_system(modelName);
open_system(modelName);

% Command r(t) = 1 + 0.35 sin(0.5 t)
add_block('simulink/Sources/Constant', [modelName '/Constant'], ...
    'Value', '1', 'Position', [30 35 70 65]);
add_block('simulink/Sources/Sine Wave', [modelName '/Sine'], ...
    'Amplitude', '0.35', 'Frequency', '0.5', ...
    'Position', [30 85 70 115]);
add_block('simulink/Math Operations/Sum', [modelName '/ReferenceSum'], ...
    'Inputs', '++', 'Position', [110 50 135 100]);

% Reference model y_m = 2/(s+2) r
add_block('simulink/Continuous/Transfer Fcn', [modelName '/ReferenceModel'], ...
    'Numerator', '[2]', 'Denominator', '[1 2]', ...
    'Position', [185 55 285 95]);

% Plant y = 0.8/(s+1.2) u
add_block('simulink/Continuous/Transfer Fcn', [modelName '/Plant'], ...
    'Numerator', '[0.8]', 'Denominator', '[1 1.2]', ...
    'Position', [710 205 810 245]);

% Controller u = theta_y*y + theta_r*r
add_block('simulink/Math Operations/Product', [modelName '/ThetaY_times_y'], ...
    'Position', [525 175 565 205]);
add_block('simulink/Math Operations/Product', [modelName '/ThetaR_times_r'], ...
    'Position', [525 245 565 275]);
add_block('simulink/Math Operations/Sum', [modelName '/ControlSum'], ...
    'Inputs', '++', 'Position', [625 205 650 255]);

% Tracking error e = y - y_m
add_block('simulink/Math Operations/Sum', [modelName '/TrackingError'], ...
    'Inputs', '+-', 'Position', [875 115 900 165]);

% Stable sensitivity filters phi_y = 1/(s+2)y, phi_r = 1/(s+2)r
add_block('simulink/Continuous/Transfer Fcn', [modelName '/FilterY'], ...
    'Numerator', '[1]', 'Denominator', '[1 2]', ...
    'Position', [890 250 980 285]);
add_block('simulink/Continuous/Transfer Fcn', [modelName '/FilterR'], ...
    'Numerator', '[1]', 'Denominator', '[1 2]', ...
    'Position', [330 350 420 385]);

% theta_y_dot = -gamma_y*e*phi_y
add_block('simulink/Math Operations/Product', [modelName '/Error_times_PhiY'], ...
    'Position', [1015 170 1055 205]);
add_block('simulink/Math Operations/Gain', [modelName '/MinusGammaY'], ...
    'Gain', '-2', 'Position', [1090 170 1145 205]);
add_block('simulink/Continuous/Integrator', [modelName '/ThetaY'], ...
    'InitialCondition', '0', 'Position', [1180 170 1210 205]);

% theta_r_dot = -gamma_r*e*phi_r
add_block('simulink/Math Operations/Product', [modelName '/Error_times_PhiR'], ...
    'Position', [1015 335 1055 370]);
add_block('simulink/Math Operations/Gain', [modelName '/MinusGammaR'], ...
    'Gain', '-2', 'Position', [1090 335 1145 370]);
add_block('simulink/Continuous/Integrator', [modelName '/ThetaR'], ...
    'InitialCondition', '0', 'Position', [1180 335 1210 370]);

% Logging
signals = {'PlantOutput', 'ModelOutput', 'Error', 'ThetaYLog', 'ThetaRLog'};
positions = {[980 65 1065 95], [330 55 415 85], [940 115 1000 145], ...
             [1250 170 1335 200], [1250 335 1335 365]};
variables = {'y_sim', 'ym_sim', 'e_sim', 'theta_y_sim', 'theta_r_sim'};
for k = 1:numel(signals)
    add_block('simulink/Sinks/To Workspace', [modelName '/' signals{k}], ...
        'VariableName', variables{k}, 'SaveFormat', 'Structure With Time', ...
        'Position', positions{k});
end

% Connections
add_line(modelName, 'Constant/1', 'ReferenceSum/1', 'autorouting', 'on');
add_line(modelName, 'Sine/1', 'ReferenceSum/2', 'autorouting', 'on');
add_line(modelName, 'ReferenceSum/1', 'ReferenceModel/1', 'autorouting', 'on');
add_line(modelName, 'ReferenceSum/1', 'ThetaR_times_r/2', 'autorouting', 'on');
add_line(modelName, 'ReferenceSum/1', 'FilterR/1', 'autorouting', 'on');

add_line(modelName, 'ThetaY_times_y/1', 'ControlSum/1', 'autorouting', 'on');
add_line(modelName, 'ThetaR_times_r/1', 'ControlSum/2', 'autorouting', 'on');
add_line(modelName, 'ControlSum/1', 'Plant/1', 'autorouting', 'on');

add_line(modelName, 'Plant/1', 'ThetaY_times_y/2', 'autorouting', 'on');
add_line(modelName, 'Plant/1', 'TrackingError/1', 'autorouting', 'on');
add_line(modelName, 'Plant/1', 'FilterY/1', 'autorouting', 'on');
add_line(modelName, 'Plant/1', 'PlantOutput/1', 'autorouting', 'on');

add_line(modelName, 'ReferenceModel/1', 'TrackingError/2', 'autorouting', 'on');
add_line(modelName, 'ReferenceModel/1', 'ModelOutput/1', 'autorouting', 'on');
add_line(modelName, 'TrackingError/1', 'Error/1', 'autorouting', 'on');

add_line(modelName, 'TrackingError/1', 'Error_times_PhiY/1', 'autorouting', 'on');
add_line(modelName, 'FilterY/1', 'Error_times_PhiY/2', 'autorouting', 'on');
add_line(modelName, 'Error_times_PhiY/1', 'MinusGammaY/1', 'autorouting', 'on');
add_line(modelName, 'MinusGammaY/1', 'ThetaY/1', 'autorouting', 'on');
add_line(modelName, 'ThetaY/1', 'ThetaY_times_y/1', 'autorouting', 'on');
add_line(modelName, 'ThetaY/1', 'ThetaYLog/1', 'autorouting', 'on');

add_line(modelName, 'TrackingError/1', 'Error_times_PhiR/1', 'autorouting', 'on');
add_line(modelName, 'FilterR/1', 'Error_times_PhiR/2', 'autorouting', 'on');
add_line(modelName, 'Error_times_PhiR/1', 'MinusGammaR/1', 'autorouting', 'on');
add_line(modelName, 'MinusGammaR/1', 'ThetaR/1', 'autorouting', 'on');
add_line(modelName, 'ThetaR/1', 'ThetaR_times_r/1', 'autorouting', 'on');
add_line(modelName, 'ThetaR/1', 'ThetaRLog/1', 'autorouting', 'on');

set_param(modelName, 'StopTime', '30', 'Solver', 'ode45');
save_system(modelName);
sim(modelName);

figure;
plot(y_sim.time, y_sim.signals.values, 'LineWidth', 1.3); hold on;
plot(ym_sim.time, ym_sim.signals.values, '--', 'LineWidth', 1.3);
grid on;
xlabel('Time [s]');
ylabel('Output');
legend('y', 'y_m', 'Location', 'best');
title('Simulink MIT-Rule MRAC');

figure;
plot(e_sim.time, e_sim.signals.values, 'LineWidth', 1.3);
grid on;
xlabel('Time [s]');
ylabel('e = y - y_m');
title('Tracking Error');

fprintf('Created and simulated %s.slx\n', modelName);

Chapter5_Lesson3.nb

The notebook uses NDSolveValue for the coupled adaptive differential equations and generates tracking, error, and parameter plots.

Notebook[{
 Cell["Chapter 5, Lesson 3: MIT Rule and Gradient Adaptive Law", "Title"],
 Cell["First-order MRAC simulation using NDSolveValue.", "Text"],
 Cell[BoxData[
  ToBoxes[
   ClearAll["Global`*"];
   plantA = 1.2; plantB = 0.8;
   modelA = 2.0; modelB = 2.0;
   gammaY = 2.0; gammaR = 2.0;
   finalTime = 30.0;

   r[t_] := 1.0 + 0.35 Sin[0.5 t];

   solution = NDSolveValue[
     {
      y'[t] == -plantA y[t] + plantB (thetaY[t] y[t] + thetaR[t] r[t]),
      yModel'[t] == -modelA yModel[t] + modelB r[t],
      phiY'[t] == -modelA phiY[t] + y[t],
      phiR'[t] == -modelA phiR[t] + r[t],
      thetaY'[t] == -gammaY (y[t] - yModel[t]) phiY[t],
      thetaR'[t] == -gammaR (y[t] - yModel[t]) phiR[t],
      y[0] == 0, yModel[0] == 0,
      phiY[0] == 0, phiR[0] == 0,
      thetaY[0] == 0, thetaR[0] == 0
     },
     {y, yModel, phiY, phiR, thetaY, thetaR},
     {t, 0, finalTime},
     Method -> {"TimeIntegration" -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5}}
   ];

   ySol[t_] := solution[[1]][t];
   yModelSol[t_] := solution[[2]][t];
   thetaYSol[t_] := solution[[5]][t];
   thetaRSol[t_] := solution[[6]][t];

   idealThetaY = (plantA - modelA)/plantB;
   idealThetaR = modelB/plantB;

   trackingPlot = Plot[
     Evaluate[{ySol[t], yModelSol[t]}],
     {t, 0, finalTime},
     PlotLegends -> {"Plant output y", "Reference-model output y_m"},
     AxesLabel -> {"Time [s]", "Output"},
     GridLines -> Automatic,
     PlotRange -> All
   ];

   errorPlot = Plot[
     Evaluate[ySol[t] - yModelSol[t]],
     {t, 0, finalTime},
     AxesLabel -> {"Time [s]", "Tracking error"},
     GridLines -> Automatic,
     PlotRange -> All
   ];

   parameterPlot = Plot[
     Evaluate[{thetaYSol[t], thetaRSol[t], idealThetaY, idealThetaR}],
     {t, 0, finalTime},
     PlotLegends -> {"theta_y", "theta_r", "ideal theta_y", "ideal theta_r"},
     AxesLabel -> {"Time [s]", "Parameter"},
     GridLines -> Automatic,
     PlotRange -> All
   ];

   rmsError = Sqrt[
     NIntegrate[(ySol[t] - yModelSol[t])^2, {t, 0, finalTime}] / finalTime
   ];

   Print["Final theta_y = ", thetaYSol[finalTime]];
   Print["Final theta_r = ", thetaRSol[finalTime]];
   Print["RMS tracking error = ", rmsError];

   GraphicsRow[{trackingPlot, errorPlot, parameterPlot}, ImageSize -> Large]
  ]
 ], "Input"]
},
WindowSize -> {1280, 780},
StyleDefinitions -> "Default.nb"
]

Students should inspect four features rather than only the final tracking plot:

1. Initial transient: parameters start at zero, so the controller initially differs substantially from the ideal matching law.
2. Error-driven motion: parameter changes are largest when both tracking error and the corresponding filtered regressor are large.
3. Reduced adaptation near tracking: when \(e\) becomes small, the gradient signal weakens.
4. Parameter non-uniqueness under finite excitation: low tracking error can occur before the gains approach their ideal values.

A gain sweep is educational. Repeating the simulation with \(\gamma_y=\gamma_r\) equal to \(0.5\), \(2\), and \(8\) typically reveals the basic speed–oscillation trade-off. The largest value should be treated as an experiment, not as a recommended design.

Problem 1 (Ideal matching gains): For \(\dot y=-1.5y+0.5u\), \(\dot y_m=-3y_m+2r\), and \(u=\theta_y y+\theta_r r\), calculate the ideal gains.

Solution:

\[ \theta_y^{\ast} = \frac{a-a_m}{b} = \frac{1.5-3}{0.5} =-3, \qquad \theta_r^{\ast} = \frac{b_m}{b} = \frac{2}{0.5} =4. \]

Substitution gives \( -a+b\theta_y^{\ast}=-3 \) and \(b\theta_r^{\ast}=2\), which exactly reproduce the reference-model coefficients.

Problem 2 (Exact sensitivity derivation): Derive the sensitivity dynamics for \(s_r=\partial y/\partial\theta_r\) when \(\dot y=(-a+b\theta_y)y+b\theta_r r\).

Solution:

\[ \begin{aligned} \dot s_r &= \frac{\partial \dot y}{\partial\theta_r} \\ &= \frac{\partial}{\partial\theta_r} \left[(-a+b\theta_y)y+b\theta_r r\right] \\ &= (-a+b\theta_y) \frac{\partial y}{\partial\theta_r} +br \\ &= (-a+b\theta_y)s_r+br. \end{aligned} \]

The term involving \(s_r\) appears because the plant output itself changes when \(\theta_r\) changes.

Problem 3 (Instantaneous update): Suppose at one instant \(e=0.4\), \(\phi_y=-0.25\), \(\phi_r=0.6\), \(\gamma_y=3\), and \(\gamma_r=1.5\). Calculate the parameter rates.

Solution:

\[ \dot\theta_y = -(3)(0.4)(-0.25) =0.30, \]

\[ \dot\theta_r = -(1.5)(0.4)(0.6) =-0.36. \]

Thus \(\theta_y\) increases while \(\theta_r\) decreases at this instant. The directions depend jointly on error sign and sensitivity sign.

Problem 4 (Frozen-time descent property): Let \(\dot{\boldsymbol\theta} =-\boldsymbol\Gamma\nabla J\), where \(\boldsymbol\Gamma\) is symmetric positive definite. Show that the parameter-adjustment contribution to the first-order cost variation is nonpositive.

Solution:

\[ \left.\frac{dJ}{dt}\right|_{\mathrm{adapt}} = \nabla J^T\dot{\boldsymbol\theta} = -\nabla J^T\boldsymbol\Gamma\nabla J \leq 0. \]

Positive definiteness implies the quadratic form is positive whenever the gradient is nonzero. This result concerns the adaptation direction at a frozen plant/model state; it does not alone prove stability of the full time-varying closed loop.

Problem 5 (Filter realization): Starting from \(\phi_r=(p+a_m)^{-1}r\), derive a state equation suitable for software or Simulink implementation.

Solution:

\[ (p+a_m)\phi_r=r. \]

Since \(p\phi_r=\dot\phi_r\):

\[ \dot\phi_r+a_m\phi_r=r, \]

and therefore:

\[ \boxed{\dot\phi_r=-a_m\phi_r+r}. \]

Problem 6 (Effect of the wrong control-gain sign): Explain what occurs if the plant has \(b < 0\) but the implementation uses the update law derived for \(b > 0\).

Solution:

The true sensitivities contain the factor \(b\). Replacing them by filtered regressors without including \(\operatorname{sgn}(b)\) reverses the gradient direction when \(b < 0\). The parameters then move approximately uphill rather than downhill in the local cost landscape. Therefore, the sign of the plant input gain must be known for this basic MIT rule. Chapter 16 later treats the unknown-control-direction problem.

The MIT rule chooses adjustable controller parameters by steepest descent on the instantaneous cost \(J=e^2/2\). Exact implementation requires error sensitivities, which obey their own differential equations but depend on unknown plant parameters. For the first-order MRAC problem, stable filters based on the reference-model denominator provide practical sensitivity approximations, yielding \(\dot\theta_y=-\gamma_y e\phi_y\) and \(\dot\theta_r=-\gamma_r e\phi_r\). The rule is intuitive and computationally simple, but it is a gradient heuristic rather than a complete closed-loop stability proof. Lesson 4 develops the Lyapunov-based alternative.