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

• identify the adjustable quantities in direct and indirect MRAC;
• write the two nested maps used by an indirect adaptive controller;
• derive the ideal controller gains for a simple matched plant;
• prove structural equivalence when plant estimates are exact;
• explain why equivalent ideal controllers can have different transients;
• implement both structures without confusing tracking and prediction errors.

The required background is limited to linear state-space models, reference models, tracking-error dynamics, quadratic Lyapunov functions, gradient adaptation, and parameter constraints. No persistent-excitation theorem, robust modification, or higher-order MRAC construction is required here.

Both architectures begin with the same plant, reference model, and tracking objective. Consider a linearly parameterized plant

\[ \dot{\mathbf{x}} = \mathbf{A}(\boldsymbol{\theta}_p)\mathbf{x} + \mathbf{B}(\boldsymbol{\theta}_p)u, \qquad y=\mathbf{C}\mathbf{x}, \]

and a stable reference model

\[ \dot{\mathbf{x}}_m = \mathbf{A}_m\mathbf{x}_m+\mathbf{B}_m r, \qquad y_m=\mathbf{C}_m\mathbf{x}_m. \]

The tracking error introduced in the previous lesson is \( e=y-y_m \), or \( \mathbf{e}=\mathbf{x}-\mathbf{x}_m \) when the full state is used. The model-following objective is

\[ \lim_{t\rightarrow\infty} e(t)=0, \]

subject to bounded closed-loop signals. The distinction between direct and indirect MRAC is therefore not the desired output behavior; it is the internal route by which the controller parameters are changed.

Let the controller be parameterized by an adjustable vector \( \hat{\boldsymbol{\theta}}_c \):

\[ u = \hat{\boldsymbol{\theta}}_c^{T}\boldsymbol{\phi}_c(\mathbf{x},r). \]

A direct adaptive law changes this controller vector without first constructing an explicit estimate of the physical plant vector \( \boldsymbol{\theta}_p \):

\[ \dot{\hat{\boldsymbol{\theta}}}_c = \mathcal{A}_d \left( e,\boldsymbol{\phi}_c,\hat{\boldsymbol{\theta}}_c \right). \]

Thus, the principal adaptive state is the controller parameter itself. The update is usually driven by tracking information because the controller gains are being chosen to make the plant follow the reference model. The adjustable gains need not have a unique physical interpretation as mass, damping, resistance, or time constant.

In compact functional notation, direct MRAC is the one-stage adaptive map

\[ (e,\mathbf{x},r) \longmapsto \dot{\hat{\boldsymbol{\theta}}}_c \longmapsto \hat{\boldsymbol{\theta}}_c \longmapsto u. \]

An indirect architecture introduces an estimated plant model

\[ \dot{\hat{\mathbf{x}}} = \mathbf{A}(\hat{\boldsymbol{\theta}}_p)\hat{\mathbf{x}} + \mathbf{B}(\hat{\boldsymbol{\theta}}_p)u + \mathbf{L}(y-\hat y), \qquad \hat y=\mathbf{C}\hat{\mathbf{x}}. \]

The identifier uses the prediction error \( \varepsilon=y-\hat y \) to update \( \hat{\boldsymbol{\theta}}_p \):

\[ \dot{\hat{\boldsymbol{\theta}}}_p = \mathcal{A}_i \left( \varepsilon,\boldsymbol{\phi}_p, \hat{\boldsymbol{\theta}}_p \right). \]

A controller-design map then converts the estimated plant parameters into controller parameters:

\[ \hat{\boldsymbol{\theta}}_c = \mathcal{D} \left( \hat{\boldsymbol{\theta}}_p, \mathbf{A}_m,\mathbf{B}_m \right). \]

Indirect MRAC is therefore a two-stage adaptive map:

\[ (y,u) \longmapsto \hat{\boldsymbol{\theta}}_p \longmapsto \mathcal{D}(\hat{\boldsymbol{\theta}}_p) = \hat{\boldsymbol{\theta}}_c \longmapsto u. \]

This separation makes the plant estimate inspectable, but it also means that estimation error is transformed through the possibly nonlinear redesign map \( \mathcal{D} \).

flowchart TD
  R["Command r"] --> RM["Reference model"]
  RM --> YM["Desired output ym"]

  R --> DC["Direct controller: \ngains ky_hat, kr_hat"]
  YD["Plant output y"] --> DC
  DC --> UDP["Control u"]
  UDP --> PD["Unknown plant"]
  PD --> YD
  YD --> ED["Tracking error e = y - ym"]
  YM --> ED
  ED --> DA["Direct adaptation"]
  DA -->|"updates controller gains"| DC

  R --> IC["Indirect controller: \ngains from estimates"]
  YI["Plant output y"] --> IC
  IC --> UIP["Control u"]
  UIP --> PI["Unknown plant"]
  PI --> YI
  YI --> ID["Plant identifier"]
  UIP --> ID
  ID --> PE["Prediction error epsilon = y - yhat"]
  PE --> PA["Parameter estimator"]
  PA -->|"updates a_hat, b_hat"| MAP["Controller-design map"]
  MAP -->|"computes ky_hat, kr_hat"| IC
        

The direct loop closes adaptation around the tracking error. The indirect loop contains an additional internal model and a redesign operation. In practical designs, both loops may use filtered regressors and additional measured signals, but the structural distinction remains.

To expose the structural difference with minimal algebra, consider the unknown first-order plant

\[ \dot y=-a y+b u, \qquad a>0, \qquad b>0, \]

and the desired reference model

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

Use the controller structure

\[ u=-k_y y+k_r r. \]

Substitution into the plant gives

\[ \dot y = -(a+b k_y)y+b k_r r. \]

Exact coefficient matching with the reference model requires

\[ a+b k_y^{*}=a_m, \qquad b k_r^{*}=b_m. \]

Hence the ideal controller parameters are

\[ k_y^{*}=\frac{a_m-a}{b}, \qquad k_r^{*}=\frac{b_m}{b}. \]

These equations define the controller-design map for this example:

\[ \mathcal{D}(a,b) = \begin{bmatrix} (a_m-a)/b\\ b_m/b \end{bmatrix}. \]

In direct MRAC, \( \hat k_y \) and \( \hat k_r \) are adaptive states:

\[ u=-\hat k_y y+\hat k_r r. \]

Define controller-parameter errors \( \tilde k_y=\hat k_y-k_y^{*} \) and \( \tilde k_r=\hat k_r-k_r^{*} \). With \( e=y-y_m \), subtraction of the reference-model equation from the plant equation yields

\[ \dot e = -a_m e - b\tilde k_y y + b\tilde k_r r. \]

For known positive control direction, one illustrative Lyapunov candidate is

\[ V_d = \frac{1}{2}e^2 + \frac{b}{2\gamma_y}\tilde k_y^2 + \frac{b}{2\gamma_r}\tilde k_r^2. \]

Choosing the direct updates

\[ \dot{\hat k}_y=\gamma_y e y, \qquad \dot{\hat k}_r=-\gamma_r e r \]

cancels the cross terms:

\[ \begin{aligned} \dot V_d &= e\dot e + \frac{b}{\gamma_y}\tilde k_y\dot{\hat k}_y + \frac{b}{\gamma_r}\tilde k_r\dot{\hat k}_r\\ &= -a_m e^2. \end{aligned} \]

This calculation is included to reveal what “direct” means: the Lyapunov function contains controller-gain errors rather than explicit physical-plant estimation errors.

In indirect MRAC, the adaptive states are \( \hat a \) and \( \hat b \). The controller is recomputed algebraically:

\[ \hat k_y = \frac{a_m-\hat a}{\hat b}, \qquad \hat k_r = \frac{b_m}{\hat b}. \]

Introduce the series-parallel identifier

\[ \dot{\hat y} = -\hat a y+\hat b u+\ell(y-\hat y), \qquad \ell>0. \]

With prediction error \( \varepsilon=y-\hat y \) and plant-parameter errors \( \tilde a=\hat a-a \), \( \tilde b=\hat b-b \), the identifier error satisfies

\[ \dot\varepsilon = -\ell\varepsilon + \tilde a y - \tilde b u. \]

Consider

\[ V_i = \frac{1}{2}\varepsilon^2 + \frac{1}{2\gamma_a}\tilde a^2 + \frac{1}{2\gamma_b}\tilde b^2. \]

The identifier updates

\[ \dot{\hat a} = -\gamma_a\varepsilon y, \qquad \dot{\hat b} = \gamma_b\varepsilon u \]

give

\[ \dot V_i=-\ell\varepsilon^2. \]

This is an identifier result, not by itself a complete closed-loop MRAC proof. The indirect controller couples the identifier to the redesign map, so boundedness and admissibility of that map must also be established.

Proposition. Assume the scalar plant has \( b\neq 0 \). If an indirect architecture reaches \( \hat a=a \) and \( \hat b=b \), then its certainty-equivalent controller equals the ideal direct controller and produces the reference-model dynamics.

Proof. Exact estimates imply

\[ \hat k_y = \frac{a_m-\hat a}{\hat b} = \frac{a_m-a}{b} = k_y^{*}, \qquad \hat k_r = \frac{b_m}{\hat b} = \frac{b_m}{b} = k_r^{*}. \]

Substituting these gains into the plant gives

\[ \begin{aligned} \dot y &= -a y + b\left( -\frac{a_m-a}{b}y+\frac{b_m}{b}r \right)\\ &= -a_m y+b_m r. \end{aligned} \]

Thus, if \( y(0)=y_m(0) \), uniqueness of solutions gives \( y(t)=y_m(t) \) for all \( t\geq 0 \). For unequal initial conditions, \( \dot e=-a_m e \), so \( e(t)=e(0)e^{-a_m t} \). Therefore the ideal direct gains and the exact-estimate indirect gains define the same nominal closed-loop system. \( \square \)

Structural equivalence at the ideal point does not imply identical adaptive trajectories. The direct parameters obey differential equations in controller space, while the indirect controller parameters are the image of estimator dynamics under \( \mathcal{D} \).

\[ \dot{\hat{\boldsymbol{\theta}}}_c = \frac{\partial\mathcal{D}} {\partial\hat{\boldsymbol{\theta}}_p} \dot{\hat{\boldsymbol{\theta}}}_p. \]

For the scalar map, its Jacobian is

\[ \frac{\partial(\hat k_y,\hat k_r)} {\partial(\hat a,\hat b)} = \begin{bmatrix} -1/\hat b & -(a_m-\hat a)/\hat b^2\\ 0 & -b_m/\hat b^2 \end{bmatrix}. \]

Therefore, small errors in \( \hat b \) can be strongly amplified when \( |\hat b| \) is small. A practical indirect design must keep the estimate inside a region such as

\[ \hat b(t)\geq b_{\min}>0. \]

Direct MRAC does not divide by a plant-gain estimate, but its controller gains can still drift or react strongly to tracking-error noise. Consequently, neither architecture is automatically superior; the dominant issue is which internal variables, constraints, and error signals best fit the available knowledge.

flowchart TD
  S["Measured signals: r, u, y"] --> Q{"Which quantities \nare adaptive?"}
  Q -->|"Controller gains"| D1["Direct MRAC"]
  D1 --> D2["Form tracking error e = y - ym"]
  D2 --> D3["Update controller gains"]
  D3 --> D4["Apply new control immediately"]

  Q -->|"Plant parameters"| I1["Indirect MRAC"]
  I1 --> I2["Run estimated plant model"]
  I2 --> I3["Form prediction error epsilon = y - yhat"]
  I3 --> I4["Update plant estimates"]
  I4 --> I5["Redesign controller from estimates"]
  I5 --> I6["Apply recomputed control"]
        

The tracking error and prediction error answer different questions:

• \( e=y-y_m \) asks whether closed-loop behavior matches the desired model.
• \( \varepsilon=y-\hat y \) asks whether the estimated plant predicts the measured plant.

A small prediction error does not automatically imply small tracking error if the redesign map or reference model is inappropriate. Conversely, small tracking error does not guarantee that estimated physical parameters equal their true values.

The programs simulate two parallel controllers on identical unknown plants. Before \( t=10 \) seconds, \( a=1.4 \) and \( b=0.8 \); after that instant, \( a=1.0 \) and \( b=1.1 \). The reference model is

\[ \dot y_m=-2y_m+2r, \]

with command

\[ r(t) = 1+0.5\sin(0.8t)+0.25\sin(1.7t). \]

The direct controller updates \( \hat k_y,\hat k_r \) from tracking error. The indirect controller updates \( \hat a,\hat b \) from prediction error and then computes controller gains from the matching equations. Simple clipping keeps \( \hat b \) away from zero; it is a numerical parameter constraint, not a substitute for a complete projection-based stability analysis.

The parameter change is included to reveal architectural response. The experiment should not be interpreted as a theorem that either adaptive law identifies arbitrary time-varying parameters.

The Python version uses NumPy for arrays and Matplotlib for visualization. The python-control package is useful for linear-system analysis, but the adaptive loop is integrated explicitly here so that the update equations remain visible.

Chapter4_Lesson3.py

"""
Chapter4_Lesson3.py
Structural comparison of direct and indirect MRAC for a first-order plant.

Dependencies:
    numpy
    matplotlib

The numerical integration uses forward Euler deliberately so that every update
equation remains visible. This is an educational architecture comparison, not a
production-grade adaptive-control implementation.
"""

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
import csv
import math

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Configuration:
    dt: float = 0.002
    final_time: float = 20.0
    model_a: float = 2.0
    model_b: float = 2.0
    direct_gamma_y: float = 4.0
    direct_gamma_r: float = 4.0
    identifier_gamma_a: float = 8.0
    identifier_gamma_b: float = 8.0
    observer_injection: float = 8.0
    b_hat_min: float = 0.20
    b_hat_max: float = 2.00
    a_hat_min: float = 0.10
    a_hat_max: float = 3.00


def command(t: float) -> float:
    """Bounded command containing more than one frequency."""
    return 1.0 + 0.50 * math.sin(0.8 * t) + 0.25 * math.sin(1.7 * t)


def plant_parameters(t: float) -> tuple[float, float]:
    """Unknown plant parameters; they change once to expose adaptation."""
    if t < 10.0:
        return 1.4, 0.8
    return 1.0, 1.1


def simulate(cfg: Configuration) -> dict[str, np.ndarray]:
    time = np.arange(0.0, cfg.final_time + cfg.dt, cfg.dt)
    n = time.size

    # Direct MRAC states and controller parameters.
    y_direct = np.zeros(n)
    ym_direct = np.zeros(n)
    ky_direct = np.zeros(n)
    kr_direct = np.zeros(n)
    kr_direct[0] = 1.0
    u_direct = np.zeros(n)

    # Indirect MRAC plant, reference model, identifier, and estimates.
    y_indirect = np.zeros(n)
    ym_indirect = np.zeros(n)
    y_hat = np.zeros(n)
    a_hat = np.zeros(n)
    b_hat = np.zeros(n)
    a_hat[0] = 0.7
    b_hat[0] = 1.2
    ky_indirect = np.zeros(n)
    kr_indirect = np.zeros(n)
    u_indirect = np.zeros(n)

    reference = np.zeros(n)
    actual_a = np.zeros(n)
    actual_b = np.zeros(n)

    for k in range(n - 1):
        t = time[k]
        r = command(t)
        a, b = plant_parameters(t)
        reference[k] = r
        actual_a[k] = a
        actual_b[k] = b

        # ---------------- Direct MRAC ----------------
        e_d = y_direct[k] - ym_direct[k]
        u_d = -ky_direct[k] * y_direct[k] + kr_direct[k] * r
        u_direct[k] = u_d

        y_dot_d = -a * y_direct[k] + b * u_d
        ym_dot_d = -cfg.model_a * ym_direct[k] + cfg.model_b * r

        # Direct adaptation changes controller gains themselves.
        ky_dot = cfg.direct_gamma_y * e_d * y_direct[k]
        kr_dot = -cfg.direct_gamma_r * e_d * r

        y_direct[k + 1] = y_direct[k] + cfg.dt * y_dot_d
        ym_direct[k + 1] = ym_direct[k] + cfg.dt * ym_dot_d
        ky_direct[k + 1] = ky_direct[k] + cfg.dt * ky_dot
        kr_direct[k + 1] = kr_direct[k] + cfg.dt * kr_dot

        # ---------------- Indirect MRAC ----------------
        # Controller gains are recomputed from the estimated plant parameters.
        b_effective = max(b_hat[k], cfg.b_hat_min)
        ky_i = (cfg.model_a - a_hat[k]) / b_effective
        kr_i = cfg.model_b / b_effective
        ky_indirect[k] = ky_i
        kr_indirect[k] = kr_i

        u_i = -ky_i * y_indirect[k] + kr_i * r
        u_indirect[k] = u_i

        y_dot_i = -a * y_indirect[k] + b * u_i
        ym_dot_i = -cfg.model_a * ym_indirect[k] + cfg.model_b * r

        prediction_error = y_indirect[k] - y_hat[k]
        y_hat_dot = (
            -a_hat[k] * y_indirect[k]
            + b_hat[k] * u_i
            + cfg.observer_injection * prediction_error
        )

        # Identifier adaptation uses prediction error, not tracking error.
        a_hat_dot = -cfg.identifier_gamma_a * prediction_error * y_indirect[k]
        b_hat_dot = cfg.identifier_gamma_b * prediction_error * u_i

        y_indirect[k + 1] = y_indirect[k] + cfg.dt * y_dot_i
        ym_indirect[k + 1] = ym_indirect[k] + cfg.dt * ym_dot_i
        y_hat[k + 1] = y_hat[k] + cfg.dt * y_hat_dot
        a_hat[k + 1] = np.clip(
            a_hat[k] + cfg.dt * a_hat_dot,
            cfg.a_hat_min,
            cfg.a_hat_max,
        )
        b_hat[k + 1] = np.clip(
            b_hat[k] + cfg.dt * b_hat_dot,
            cfg.b_hat_min,
            cfg.b_hat_max,
        )

    reference[-1] = command(time[-1])
    actual_a[-1], actual_b[-1] = plant_parameters(time[-1])
    ky_indirect[-1] = (cfg.model_a - a_hat[-1]) / max(
        b_hat[-1], cfg.b_hat_min
    )
    kr_indirect[-1] = cfg.model_b / max(b_hat[-1], cfg.b_hat_min)
    u_direct[-1] = (
        -ky_direct[-1] * y_direct[-1] + kr_direct[-1] * reference[-1]
    )
    u_indirect[-1] = (
        -ky_indirect[-1] * y_indirect[-1]
        + kr_indirect[-1] * reference[-1]
    )

    return {
        "time": time,
        "reference": reference,
        "actual_a": actual_a,
        "actual_b": actual_b,
        "y_direct": y_direct,
        "ym_direct": ym_direct,
        "ky_direct": ky_direct,
        "kr_direct": kr_direct,
        "u_direct": u_direct,
        "y_indirect": y_indirect,
        "ym_indirect": ym_indirect,
        "y_hat": y_hat,
        "a_hat": a_hat,
        "b_hat": b_hat,
        "ky_indirect": ky_indirect,
        "kr_indirect": kr_indirect,
        "u_indirect": u_indirect,
    }


def rms(values: np.ndarray) -> float:
    return float(np.sqrt(np.mean(np.square(values))))


def save_csv(data: dict[str, np.ndarray], target: Path) -> None:
    keys = list(data.keys())
    with target.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(keys)
        for row in zip(*(data[key] for key in keys), strict=True):
            writer.writerow(row)


def save_figures(data: dict[str, np.ndarray], target: Path) -> None:
    t = data["time"]

    plt.figure(figsize=(10, 6))
    plt.plot(t, data["ym_direct"], label="Reference model")
    plt.plot(t, data["y_direct"], label="Direct MRAC")
    plt.plot(t, data["y_indirect"], label="Indirect MRAC")
    plt.xlabel("Time (s)")
    plt.ylabel("Output")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(target, dpi=180)
    plt.close()


def main() -> None:
    cfg = Configuration()
    data = simulate(cfg)

    direct_error = data["y_direct"] - data["ym_direct"]
    indirect_error = data["y_indirect"] - data["ym_indirect"]

    print(f"Direct MRAC RMS tracking error:   {rms(direct_error):.6f}")
    print(f"Indirect MRAC RMS tracking error: {rms(indirect_error):.6f}")
    print(
        "Final direct gains: "
        f"ky={data['ky_direct'][-1]:.6f}, "
        f"kr={data['kr_direct'][-1]:.6f}"
    )
    print(
        "Final indirect estimates: "
        f"a_hat={data['a_hat'][-1]:.6f}, "
        f"b_hat={data['b_hat'][-1]:.6f}"
    )

    base = Path(__file__).resolve().parent
    save_csv(data, base / "Chapter4_Lesson3_results.csv")
    save_figures(data, base / "Chapter4_Lesson3_tracking.png")


if __name__ == "__main__":
    main()

The C++17 implementation uses only the standard library. For higher-dimensional MRAC, Eigen is a common choice for matrix and vector operations.

Chapter4_Lesson3.cpp

/*
Chapter4_Lesson3.cpp
Structural comparison of direct and indirect MRAC for a first-order plant.

Build:
    g++ -std=c++17 -O2 Chapter4_Lesson3.cpp -o Chapter4_Lesson3
*/

#include <algorithm>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>

struct Configuration {
    double dt = 0.002;
    double finalTime = 20.0;
    double modelA = 2.0;
    double modelB = 2.0;
    double directGammaY = 4.0;
    double directGammaR = 4.0;
    double identifierGammaA = 8.0;
    double identifierGammaB = 8.0;
    double observerInjection = 8.0;
    double bHatMin = 0.20;
    double bHatMax = 2.00;
    double aHatMin = 0.10;
    double aHatMax = 3.00;
};

double command(double t) {
    return 1.0 + 0.50 * std::sin(0.8 * t) + 0.25 * std::sin(1.7 * t);
}

void plantParameters(double t, double& a, double& b) {
    if (t < 10.0) {
        a = 1.4;
        b = 0.8;
    } else {
        a = 1.0;
        b = 1.1;
    }
}

double rms(const std::vector<double>& values) {
    double sumSquares = 0.0;
    for (double value : values) {
        sumSquares += value * value;
    }
    return std::sqrt(sumSquares / static_cast<double>(values.size()));
}

int main() {
    try {
        const Configuration cfg;
        const std::size_t n =
            static_cast<std::size_t>(cfg.finalTime / cfg.dt) + 1U;

        std::vector<double> time(n), reference(n), actualA(n), actualB(n);
        std::vector<double> yDirect(n), ymDirect(n), kyDirect(n), krDirect(n);
        std::vector<double> uDirect(n);
        std::vector<double> yIndirect(n), ymIndirect(n), yHat(n);
        std::vector<double> aHat(n), bHat(n), kyIndirect(n), krIndirect(n);
        std::vector<double> uIndirect(n);

        krDirect[0] = 1.0;
        aHat[0] = 0.7;
        bHat[0] = 1.2;

        for (std::size_t k = 0; k + 1U < n; ++k) {
            const double t = static_cast<double>(k) * cfg.dt;
            time[k] = t;
            reference[k] = command(t);
            plantParameters(t, actualA[k], actualB[k]);

            const double r = reference[k];
            const double a = actualA[k];
            const double b = actualB[k];

            // ---------------- Direct MRAC ----------------
            const double eDirect = yDirect[k] - ymDirect[k];
            uDirect[k] = -kyDirect[k] * yDirect[k] + krDirect[k] * r;

            const double yDirectDot = -a * yDirect[k] + b * uDirect[k];
            const double ymDirectDot =
                -cfg.modelA * ymDirect[k] + cfg.modelB * r;
            const double kyDot =
                cfg.directGammaY * eDirect * yDirect[k];
            const double krDot =
                -cfg.directGammaR * eDirect * r;

            yDirect[k + 1U] = yDirect[k] + cfg.dt * yDirectDot;
            ymDirect[k + 1U] = ymDirect[k] + cfg.dt * ymDirectDot;
            kyDirect[k + 1U] = kyDirect[k] + cfg.dt * kyDot;
            krDirect[k + 1U] = krDirect[k] + cfg.dt * krDot;

            // ---------------- Indirect MRAC ----------------
            const double bEffective = std::max(bHat[k], cfg.bHatMin);
            kyIndirect[k] = (cfg.modelA - aHat[k]) / bEffective;
            krIndirect[k] = cfg.modelB / bEffective;
            uIndirect[k] =
                -kyIndirect[k] * yIndirect[k] + krIndirect[k] * r;

            const double yIndirectDot =
                -a * yIndirect[k] + b * uIndirect[k];
            const double ymIndirectDot =
                -cfg.modelA * ymIndirect[k] + cfg.modelB * r;
            const double predictionError = yIndirect[k] - yHat[k];
            const double yHatDot =
                -aHat[k] * yIndirect[k]
                + bHat[k] * uIndirect[k]
                + cfg.observerInjection * predictionError;

            const double aHatDot =
                -cfg.identifierGammaA * predictionError * yIndirect[k];
            const double bHatDot =
                cfg.identifierGammaB * predictionError * uIndirect[k];

            yIndirect[k + 1U] = yIndirect[k] + cfg.dt * yIndirectDot;
            ymIndirect[k + 1U] = ymIndirect[k] + cfg.dt * ymIndirectDot;
            yHat[k + 1U] = yHat[k] + cfg.dt * yHatDot;
            aHat[k + 1U] = std::clamp(
                aHat[k] + cfg.dt * aHatDot,
                cfg.aHatMin,
                cfg.aHatMax
            );
            bHat[k + 1U] = std::clamp(
                bHat[k] + cfg.dt * bHatDot,
                cfg.bHatMin,
                cfg.bHatMax
            );
        }

        time[n - 1U] = cfg.finalTime;
        reference[n - 1U] = command(cfg.finalTime);
        plantParameters(
            cfg.finalTime,
            actualA[n - 1U],
            actualB[n - 1U]
        );
        const double bEffective =
            std::max(bHat[n - 1U], cfg.bHatMin);
        kyIndirect[n - 1U] =
            (cfg.modelA - aHat[n - 1U]) / bEffective;
        krIndirect[n - 1U] = cfg.modelB / bEffective;
        uDirect[n - 1U] =
            -kyDirect[n - 1U] * yDirect[n - 1U]
            + krDirect[n - 1U] * reference[n - 1U];
        uIndirect[n - 1U] =
            -kyIndirect[n - 1U] * yIndirect[n - 1U]
            + krIndirect[n - 1U] * reference[n - 1U];

        std::vector<double> directError(n), indirectError(n);
        for (std::size_t k = 0; k < n; ++k) {
            directError[k] = yDirect[k] - ymDirect[k];
            indirectError[k] = yIndirect[k] - ymIndirect[k];
        }

        std::cout << std::fixed << std::setprecision(6);
        std::cout << "Direct MRAC RMS tracking error:   "
                  << rms(directError) << '\n';
        std::cout << "Indirect MRAC RMS tracking error: "
                  << rms(indirectError) << '\n';
        std::cout << "Final direct gains: ky="
                  << kyDirect.back()
                  << ", kr=" << krDirect.back() << '\n';
        std::cout << "Final indirect estimates: a_hat="
                  << aHat.back()
                  << ", b_hat=" << bHat.back() << '\n';

        std::ofstream csv("Chapter4_Lesson3_results_cpp.csv");
        if (!csv) {
            throw std::runtime_error("Could not create CSV output.");
        }
        csv << "time,reference,actual_a,actual_b,"
               "y_direct,ym_direct,ky_direct,kr_direct,u_direct,"
               "y_indirect,ym_indirect,y_hat,a_hat,b_hat,"
               "ky_indirect,kr_indirect,u_indirect\n";
        csv << std::setprecision(12);
        for (std::size_t k = 0; k < n; ++k) {
            csv << time[k] << ','
                << reference[k] << ','
                << actualA[k] << ','
                << actualB[k] << ','
                << yDirect[k] << ','
                << ymDirect[k] << ','
                << kyDirect[k] << ','
                << krDirect[k] << ','
                << uDirect[k] << ','
                << yIndirect[k] << ','
                << ymIndirect[k] << ','
                << yHat[k] << ','
                << aHat[k] << ','
                << bHat[k] << ','
                << kyIndirect[k] << ','
                << krIndirect[k] << ','
                << uIndirect[k] << '\n';
        }
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }

    return 0;
}

The Java implementation uses primitive arrays and writes a CSV file. Libraries such as EJML or Apache Commons Math can support larger matrix-based adaptive systems.

Chapter4_Lesson3.java

/*
Chapter4_Lesson3.java
Structural comparison of direct and indirect MRAC for a first-order plant.

Build and run:
    javac Chapter4_Lesson3.java
    java Chapter4_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 Chapter4_Lesson3 {
    private static final double DT = 0.002;
    private static final double FINAL_TIME = 20.0;
    private static final double MODEL_A = 2.0;
    private static final double MODEL_B = 2.0;
    private static final double DIRECT_GAMMA_Y = 4.0;
    private static final double DIRECT_GAMMA_R = 4.0;
    private static final double IDENTIFIER_GAMMA_A = 8.0;
    private static final double IDENTIFIER_GAMMA_B = 8.0;
    private static final double OBSERVER_INJECTION = 8.0;
    private static final double B_HAT_MIN = 0.20;
    private static final double B_HAT_MAX = 2.00;
    private static final double A_HAT_MIN = 0.10;
    private static final double A_HAT_MAX = 3.00;

    private Chapter4_Lesson3() {
    }

    private static double command(double t) {
        return 1.0
            + 0.50 * Math.sin(0.8 * t)
            + 0.25 * Math.sin(1.7 * t);
    }

    private static double plantA(double t) {
        return t < 10.0 ? 1.4 : 1.0;
    }

    private static double plantB(double t) {
        return t < 10.0 ? 0.8 : 1.1;
    }

    private static double clamp(double value, double lower, double upper) {
        return Math.max(lower, Math.min(value, upper));
    }

    private static double rms(double[] values) {
        double sumSquares = 0.0;
        for (double value : values) {
            sumSquares += value * value;
        }
        return Math.sqrt(sumSquares / values.length);
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);

        final int n = (int) Math.round(FINAL_TIME / DT) + 1;

        double[] time = new double[n];
        double[] reference = new double[n];
        double[] actualA = new double[n];
        double[] actualB = new double[n];

        double[] yDirect = new double[n];
        double[] ymDirect = new double[n];
        double[] kyDirect = new double[n];
        double[] krDirect = new double[n];
        double[] uDirect = new double[n];

        double[] yIndirect = new double[n];
        double[] ymIndirect = new double[n];
        double[] yHat = new double[n];
        double[] aHat = new double[n];
        double[] bHat = new double[n];
        double[] kyIndirect = new double[n];
        double[] krIndirect = new double[n];
        double[] uIndirect = new double[n];

        krDirect[0] = 1.0;
        aHat[0] = 0.7;
        bHat[0] = 1.2;

        for (int k = 0; k < n - 1; k++) {
            final double t = k * DT;
            final double r = command(t);
            final double a = plantA(t);
            final double b = plantB(t);

            time[k] = t;
            reference[k] = r;
            actualA[k] = a;
            actualB[k] = b;

            // ---------------- Direct MRAC ----------------
            final double eDirect = yDirect[k] - ymDirect[k];
            uDirect[k] =
                -kyDirect[k] * yDirect[k] + krDirect[k] * r;

            final double yDirectDot =
                -a * yDirect[k] + b * uDirect[k];
            final double ymDirectDot =
                -MODEL_A * ymDirect[k] + MODEL_B * r;
            final double kyDot =
                DIRECT_GAMMA_Y * eDirect * yDirect[k];
            final double krDot =
                -DIRECT_GAMMA_R * eDirect * r;

            yDirect[k + 1] = yDirect[k] + DT * yDirectDot;
            ymDirect[k + 1] = ymDirect[k] + DT * ymDirectDot;
            kyDirect[k + 1] = kyDirect[k] + DT * kyDot;
            krDirect[k + 1] = krDirect[k] + DT * krDot;

            // ---------------- Indirect MRAC ----------------
            final double bEffective = Math.max(bHat[k], B_HAT_MIN);
            kyIndirect[k] = (MODEL_A - aHat[k]) / bEffective;
            krIndirect[k] = MODEL_B / bEffective;
            uIndirect[k] =
                -kyIndirect[k] * yIndirect[k]
                + krIndirect[k] * r;

            final double yIndirectDot =
                -a * yIndirect[k] + b * uIndirect[k];
            final double ymIndirectDot =
                -MODEL_A * ymIndirect[k] + MODEL_B * r;
            final double predictionError = yIndirect[k] - yHat[k];
            final double yHatDot =
                -aHat[k] * yIndirect[k]
                + bHat[k] * uIndirect[k]
                + OBSERVER_INJECTION * predictionError;
            final double aHatDot =
                -IDENTIFIER_GAMMA_A
                * predictionError
                * yIndirect[k];
            final double bHatDot =
                IDENTIFIER_GAMMA_B
                * predictionError
                * uIndirect[k];

            yIndirect[k + 1] =
                yIndirect[k] + DT * yIndirectDot;
            ymIndirect[k + 1] =
                ymIndirect[k] + DT * ymIndirectDot;
            yHat[k + 1] = yHat[k] + DT * yHatDot;
            aHat[k + 1] = clamp(
                aHat[k] + DT * aHatDot,
                A_HAT_MIN,
                A_HAT_MAX
            );
            bHat[k + 1] = clamp(
                bHat[k] + DT * bHatDot,
                B_HAT_MIN,
                B_HAT_MAX
            );
        }

        time[n - 1] = FINAL_TIME;
        reference[n - 1] = command(FINAL_TIME);
        actualA[n - 1] = plantA(FINAL_TIME);
        actualB[n - 1] = plantB(FINAL_TIME);
        final double bEffective = Math.max(bHat[n - 1], B_HAT_MIN);
        kyIndirect[n - 1] = (MODEL_A - aHat[n - 1]) / bEffective;
        krIndirect[n - 1] = MODEL_B / bEffective;
        uDirect[n - 1] =
            -kyDirect[n - 1] * yDirect[n - 1]
            + krDirect[n - 1] * reference[n - 1];
        uIndirect[n - 1] =
            -kyIndirect[n - 1] * yIndirect[n - 1]
            + krIndirect[n - 1] * reference[n - 1];

        double[] directError = new double[n];
        double[] indirectError = new double[n];
        for (int k = 0; k < n; k++) {
            directError[k] = yDirect[k] - ymDirect[k];
            indirectError[k] = yIndirect[k] - ymIndirect[k];
        }

        System.out.printf(
            "Direct MRAC RMS tracking error:   %.6f%n",
            rms(directError)
        );
        System.out.printf(
            "Indirect MRAC RMS tracking error: %.6f%n",
            rms(indirectError)
        );
        System.out.printf(
            "Final direct gains: ky=%.6f, kr=%.6f%n",
            kyDirect[n - 1],
            krDirect[n - 1]
        );
        System.out.printf(
            "Final indirect estimates: a_hat=%.6f, b_hat=%.6f%n",
            aHat[n - 1],
            bHat[n - 1]
        );

        Path csvPath = Path.of("Chapter4_Lesson3_results_java.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(
                csvPath,
                StandardCharsets.UTF_8
        )) {
            writer.write(
                "time,reference,actual_a,actual_b,"
                + "y_direct,ym_direct,ky_direct,kr_direct,u_direct,"
                + "y_indirect,ym_indirect,y_hat,a_hat,b_hat,"
                + "ky_indirect,kr_indirect,u_indirect"
            );
            writer.newLine();

            for (int k = 0; k < n; k++) {
                writer.write(String.format(
                    Locale.ROOT,
                    "%.12f,%.12f,%.12f,%.12f,"
                    + "%.12f,%.12f,%.12f,%.12f,%.12f,"
                    + "%.12f,%.12f,%.12f,%.12f,%.12f,"
                    + "%.12f,%.12f,%.12f",
                    time[k],
                    reference[k],
                    actualA[k],
                    actualB[k],
                    yDirect[k],
                    ymDirect[k],
                    kyDirect[k],
                    krDirect[k],
                    uDirect[k],
                    yIndirect[k],
                    ymIndirect[k],
                    yHat[k],
                    aHat[k],
                    bHat[k],
                    kyIndirect[k],
                    krIndirect[k],
                    uIndirect[k]
                ));
                writer.newLine();
            }
        } catch (IOException error) {
            System.err.println("Could not write CSV: " + error.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB script performs the same numerical simulation and, when Simulink is licensed, creates a structural model named Chapter4_Lesson3_Structure.slx. The generated Simulink model is intentionally a scaffold: its parallel subsystems emphasize direct versus indirect signal flow, while the script contains the complete numerical update equations.

Chapter4_Lesson3.m

% Chapter4_Lesson3.m
% Structural comparison of direct and indirect MRAC for a first-order plant.
% The script also creates a conceptual Simulink scaffold when Simulink exists.

clear; close all; clc;

dt = 0.002;
finalTime = 20.0;
time = (0:dt:finalTime).';
n = numel(time);

modelA = 2.0;
modelB = 2.0;

directGammaY = 4.0;
directGammaR = 4.0;

identifierGammaA = 8.0;
identifierGammaB = 8.0;
observerInjection = 8.0;

aHatMin = 0.10;
aHatMax = 3.00;
bHatMin = 0.20;
bHatMax = 2.00;

reference = zeros(n, 1);
actualA = zeros(n, 1);
actualB = zeros(n, 1);

yDirect = zeros(n, 1);
ymDirect = zeros(n, 1);
kyDirect = zeros(n, 1);
krDirect = zeros(n, 1);
uDirect = zeros(n, 1);
krDirect(1) = 1.0;

yIndirect = zeros(n, 1);
ymIndirect = zeros(n, 1);
yHat = zeros(n, 1);
aHat = zeros(n, 1);
bHat = zeros(n, 1);
kyIndirect = zeros(n, 1);
krIndirect = zeros(n, 1);
uIndirect = zeros(n, 1);
aHat(1) = 0.7;
bHat(1) = 1.2;

for k = 1:n-1
    t = time(k);
    reference(k) = commandSignal(t);

    if t < 10.0
        actualA(k) = 1.4;
        actualB(k) = 0.8;
    else
        actualA(k) = 1.0;
        actualB(k) = 1.1;
    end

    r = reference(k);
    a = actualA(k);
    b = actualB(k);

    % ---------------- Direct MRAC ----------------
    eDirect = yDirect(k) - ymDirect(k);
    uDirect(k) = -kyDirect(k) * yDirect(k) + krDirect(k) * r;

    yDirectDot = -a * yDirect(k) + b * uDirect(k);
    ymDirectDot = -modelA * ymDirect(k) + modelB * r;

    kyDot = directGammaY * eDirect * yDirect(k);
    krDot = -directGammaR * eDirect * r;

    yDirect(k+1) = yDirect(k) + dt * yDirectDot;
    ymDirect(k+1) = ymDirect(k) + dt * ymDirectDot;
    kyDirect(k+1) = kyDirect(k) + dt * kyDot;
    krDirect(k+1) = krDirect(k) + dt * krDot;

    % ---------------- Indirect MRAC ----------------
    bEffective = max(bHat(k), bHatMin);
    kyIndirect(k) = (modelA - aHat(k)) / bEffective;
    krIndirect(k) = modelB / bEffective;
    uIndirect(k) = ...
        -kyIndirect(k) * yIndirect(k) + krIndirect(k) * r;

    yIndirectDot = -a * yIndirect(k) + b * uIndirect(k);
    ymIndirectDot = -modelA * ymIndirect(k) + modelB * r;

    predictionError = yIndirect(k) - yHat(k);
    yHatDot = ...
        -aHat(k) * yIndirect(k) ...
        + bHat(k) * uIndirect(k) ...
        + observerInjection * predictionError;

    aHatDot = ...
        -identifierGammaA * predictionError * yIndirect(k);
    bHatDot = ...
        identifierGammaB * predictionError * uIndirect(k);

    yIndirect(k+1) = yIndirect(k) + dt * yIndirectDot;
    ymIndirect(k+1) = ymIndirect(k) + dt * ymIndirectDot;
    yHat(k+1) = yHat(k) + dt * yHatDot;

    aHat(k+1) = min( ...
        max(aHat(k) + dt * aHatDot, aHatMin), ...
        aHatMax ...
    );
    bHat(k+1) = min( ...
        max(bHat(k) + dt * bHatDot, bHatMin), ...
        bHatMax ...
    );
end

reference(end) = commandSignal(time(end));
actualA(end) = 1.0;
actualB(end) = 1.1;
bEffective = max(bHat(end), bHatMin);
kyIndirect(end) = (modelA - aHat(end)) / bEffective;
krIndirect(end) = modelB / bEffective;
uDirect(end) = ...
    -kyDirect(end) * yDirect(end) + krDirect(end) * reference(end);
uIndirect(end) = ...
    -kyIndirect(end) * yIndirect(end) ...
    + krIndirect(end) * reference(end);

directError = yDirect - ymDirect;
indirectError = yIndirect - ymIndirect;

fprintf( ...
    'Direct MRAC RMS tracking error:   %.6f\n', ...
    sqrt(mean(directError.^2)) ...
);
fprintf( ...
    'Indirect MRAC RMS tracking error: %.6f\n', ...
    sqrt(mean(indirectError.^2)) ...
);
fprintf( ...
    'Final direct gains: ky=%.6f, kr=%.6f\n', ...
    kyDirect(end), ...
    krDirect(end) ...
);
fprintf( ...
    'Final indirect estimates: a_hat=%.6f, b_hat=%.6f\n', ...
    aHat(end), ...
    bHat(end) ...
);

results = table( ...
    time, reference, actualA, actualB, ...
    yDirect, ymDirect, kyDirect, krDirect, uDirect, ...
    yIndirect, ymIndirect, yHat, aHat, bHat, ...
    kyIndirect, krIndirect, uIndirect ...
);
writetable(results, 'Chapter4_Lesson3_results_matlab.csv');

figure('Name', 'Direct vs Indirect MRAC');
plot(time, ymDirect, 'LineWidth', 1.4); hold on;
plot(time, yDirect, 'LineWidth', 1.1);
plot(time, yIndirect, 'LineWidth', 1.1);
grid on;
xlabel('Time (s)');
ylabel('Output');
legend('Reference model', 'Direct MRAC', 'Indirect MRAC', ...
    'Location', 'best');
title('Structural comparison of direct and indirect MRAC');

if license('test', 'Simulink')
    buildStructuralSimulinkModel();
end

function r = commandSignal(t)
    r = 1.0 + 0.50 * sin(0.8 * t) + 0.25 * sin(1.7 * t);
end

function buildStructuralSimulinkModel()
    % This model is a visual scaffold, while the script above performs the
    % numerical adaptive simulation. Later lessons can replace each
    % pass-through subsystem with detailed controller and identifier blocks.

    model = 'Chapter4_Lesson3_Structure';

    if bdIsLoaded(model)
        close_system(model, 0);
    end

    new_system(model);
    open_system(model);

    add_block( ...
        'simulink/Sources/Step', ...
        [model '/Command'], ...
        'Position', [35 115 65 145] ...
    );
    add_block( ...
        'simulink/Continuous/Transfer Fcn', ...
        [model '/Reference Model'], ...
        'Numerator', '[2]', ...
        'Denominator', '[1 2]', ...
        'Position', [120 25 245 75] ...
    );
    add_block( ...
        'simulink/Ports & Subsystems/Subsystem', ...
        [model '/Direct MRAC'], ...
        'Position', [120 100 245 155] ...
    );
    add_block( ...
        'simulink/Ports & Subsystems/Subsystem', ...
        [model '/Indirect MRAC'], ...
        'Position', [120 190 245 245] ...
    );
    add_block( ...
        'simulink/Signal Routing/Mux', ...
        [model '/Mux'], ...
        'Inputs', '3', ...
        'Position', [315 70 320 215] ...
    );
    add_block( ...
        'simulink/Sinks/Scope', ...
        [model '/Scope'], ...
        'Position', [370 120 405 155] ...
    );

    makePassThroughSubsystem( ...
        [model '/Direct MRAC'], ...
        'tracking error updates controller gains' ...
    );
    makePassThroughSubsystem( ...
        [model '/Indirect MRAC'], ...
        'prediction error updates plant estimates' ...
    );

    add_line(model, 'Command/1', 'Reference Model/1', 'autorouting', 'on');
    add_line(model, 'Command/1', 'Direct MRAC/1', 'autorouting', 'on');
    add_line(model, 'Command/1', 'Indirect MRAC/1', 'autorouting', 'on');
    add_line(model, 'Reference Model/1', 'Mux/1', 'autorouting', 'on');
    add_line(model, 'Direct MRAC/1', 'Mux/2', 'autorouting', 'on');
    add_line(model, 'Indirect MRAC/1', 'Mux/3', 'autorouting', 'on');
    add_line(model, 'Mux/1', 'Scope/1', 'autorouting', 'on');

    set_param(model, 'StopTime', '20');
    save_system(model);
end

function makePassThroughSubsystem(path, annotationText)
    Simulink.SubSystem.deleteContents(path);

    add_block( ...
        'simulink/Ports & Subsystems/In1', ...
        [path '/r'], ...
        'Position', [30 40 60 60] ...
    );
    add_block( ...
        'simulink/Ports & Subsystems/Out1', ...
        [path '/y'], ...
        'Position', [190 40 220 60] ...
    );
    add_block( ...
        'simulink/Signal Routing/Goto', ...
        [path '/Architecture Note'], ...
        'GotoTag', matlab.lang.makeValidName(annotationText), ...
        'TagVisibility', 'local', ...
        'Position', [90 80 155 105] ...
    );
    add_line(path, 'r/1', 'y/1', 'autorouting', 'on');
end

The notebook uses native list operations, Clip, ListLinePlot, and Export. The downloadable notebook contains this complete input cell.

Chapter4_Lesson3.nb

(* Chapter4_Lesson3.nb
   Structural comparison of direct and indirect MRAC.
   Evaluate the notebook input cell to run the simulation. *)

ClearAll["Global`*"];

dt = 0.002;
finalTime = 20.0;
time = Range[0.0, finalTime, dt];
n = Length[time];

modelA = 2.0;
modelB = 2.0;

directGammaY = 4.0;
directGammaR = 4.0;
identifierGammaA = 8.0;
identifierGammaB = 8.0;
observerInjection = 8.0;

aHatMin = 0.10;
aHatMax = 3.00;
bHatMin = 0.20;
bHatMax = 2.00;

command[t_] := 1.0 + 0.50 Sin[0.8 t] + 0.25 Sin[1.7 t];
plantA[t_] := If[t < 10.0, 1.4, 1.0];
plantB[t_] := If[t < 10.0, 0.8, 1.1];

reference = ConstantArray[0.0, n];
actualA = ConstantArray[0.0, n];
actualB = ConstantArray[0.0, n];

yDirect = ConstantArray[0.0, n];
ymDirect = ConstantArray[0.0, n];
kyDirect = ConstantArray[0.0, n];
krDirect = ConstantArray[0.0, n];
uDirect = ConstantArray[0.0, n];
krDirect[[1]] = 1.0;

yIndirect = ConstantArray[0.0, n];
ymIndirect = ConstantArray[0.0, n];
yHat = ConstantArray[0.0, n];
aHat = ConstantArray[0.0, n];
bHat = ConstantArray[0.0, n];
kyIndirect = ConstantArray[0.0, n];
krIndirect = ConstantArray[0.0, n];
uIndirect = ConstantArray[0.0, n];
aHat[[1]] = 0.7;
bHat[[1]] = 1.2;

Do[
  t = time[[k]];
  r = command[t];
  a = plantA[t];
  b = plantB[t];

  reference[[k]] = r;
  actualA[[k]] = a;
  actualB[[k]] = b;

  (* Direct MRAC *)
  eDirect = yDirect[[k]] - ymDirect[[k]];
  uDirect[[k]] = -kyDirect[[k]] yDirect[[k]] + krDirect[[k]] r;

  yDirectDot = -a yDirect[[k]] + b uDirect[[k]];
  ymDirectDot = -modelA ymDirect[[k]] + modelB r;
  kyDot = directGammaY eDirect yDirect[[k]];
  krDot = -directGammaR eDirect r;

  yDirect[[k + 1]] = yDirect[[k]] + dt yDirectDot;
  ymDirect[[k + 1]] = ymDirect[[k]] + dt ymDirectDot;
  kyDirect[[k + 1]] = kyDirect[[k]] + dt kyDot;
  krDirect[[k + 1]] = krDirect[[k]] + dt krDot;

  (* Indirect MRAC *)
  bEffective = Max[bHat[[k]], bHatMin];
  kyIndirect[[k]] = (modelA - aHat[[k]])/bEffective;
  krIndirect[[k]] = modelB/bEffective;
  uIndirect[[k]] =
    -kyIndirect[[k]] yIndirect[[k]] + krIndirect[[k]] r;

  yIndirectDot = -a yIndirect[[k]] + b uIndirect[[k]];
  ymIndirectDot = -modelA ymIndirect[[k]] + modelB r;

  predictionError = yIndirect[[k]] - yHat[[k]];
  yHatDot =
    -aHat[[k]] yIndirect[[k]]
    + bHat[[k]] uIndirect[[k]]
    + observerInjection predictionError;

  aHatDot =
    -identifierGammaA predictionError yIndirect[[k]];
  bHatDot =
    identifierGammaB predictionError uIndirect[[k]];

  yIndirect[[k + 1]] = yIndirect[[k]] + dt yIndirectDot;
  ymIndirect[[k + 1]] = ymIndirect[[k]] + dt ymIndirectDot;
  yHat[[k + 1]] = yHat[[k]] + dt yHatDot;
  aHat[[k + 1]] =
    Clip[aHat[[k]] + dt aHatDot, {aHatMin, aHatMax}];
  bHat[[k + 1]] =
    Clip[bHat[[k]] + dt bHatDot, {bHatMin, bHatMax}],
  {k, 1, n - 1}
];

reference[[-1]] = command[time[[-1]]];
actualA[[-1]] = plantA[time[[-1]]];
actualB[[-1]] = plantB[time[[-1]]];

bEffective = Max[bHat[[-1]], bHatMin];
kyIndirect[[-1]] = (modelA - aHat[[-1]])/bEffective;
krIndirect[[-1]] = modelB/bEffective;
uDirect[[-1]] =
  -kyDirect[[-1]] yDirect[[-1]] + krDirect[[-1]] reference[[-1]];
uIndirect[[-1]] =
  -kyIndirect[[-1]] yIndirect[[-1]]
  + krIndirect[[-1]] reference[[-1]];

directError = yDirect - ymDirect;
indirectError = yIndirect - ymIndirect;

rms[x_] := Sqrt[Mean[x^2]];

Print[
  "Direct MRAC RMS tracking error: ",
  NumberForm[rms[directError], {8, 6}]
];
Print[
  "Indirect MRAC RMS tracking error: ",
  NumberForm[rms[indirectError], {8, 6}]
];
Print[
  "Final direct gains: ky=",
  NumberForm[kyDirect[[-1]], {8, 6}],
  ", kr=",
  NumberForm[krDirect[[-1]], {8, 6}]
];
Print[
  "Final indirect estimates: aHat=",
  NumberForm[aHat[[-1]], {8, 6}],
  ", bHat=",
  NumberForm[bHat[[-1]], {8, 6}]
];

trackingPlot = ListLinePlot[
  {
    Transpose[{time, ymDirect}],
    Transpose[{time, yDirect}],
    Transpose[{time, yIndirect}]
  },
  PlotLegends -> {
    "Reference model",
    "Direct MRAC",
    "Indirect MRAC"
  },
  Frame -> True,
  FrameLabel -> {"Time (s)", "Output"},
  PlotRange -> All,
  ImageSize -> Large
];

Export["Chapter4_Lesson3_tracking_mathematica.png", trackingPlot];

rows = Transpose[{
  time,
  reference,
  actualA,
  actualB,
  yDirect,
  ymDirect,
  kyDirect,
  krDirect,
  uDirect,
  yIndirect,
  ymIndirect,
  yHat,
  aHat,
  bHat,
  kyIndirect,
  krIndirect,
  uIndirect
}];

header = {
  "time",
  "reference",
  "actual_a",
  "actual_b",
  "y_direct",
  "ym_direct",
  "ky_direct",
  "kr_direct",
  "u_direct",
  "y_indirect",
  "ym_indirect",
  "y_hat",
  "a_hat",
  "b_hat",
  "ky_indirect",
  "kr_indirect",
  "u_indirect"
};

Export[
  "Chapter4_Lesson3_results_mathematica.csv",
  Prepend[rows, header]
];

trackingPlot

The most important outputs are not merely the RMS errors. Students should inspect the internal trajectories:

• \( \hat k_y,\hat k_r \) show how direct adaptation moves in controller space.
• \( \hat a,\hat b \) show how the indirect identifier moves in plant-model space.
• \( \mathcal{D}(\hat a,\hat b) \) shows how estimator transients are transformed into control-gain transients.
• \( y-\hat y \) and \( y-y_m \) should be plotted separately because they measure different objectives.

Numerical closeness of \( y \) and \( y_m \) does not establish convergence of physical estimates. Similarly, accurate prediction alone does not establish the desired model-following transient.

Problem 1 (Matching Map): For \( \dot y=-a y+b u \), \( \dot y_m=-a_m y_m+b_m r \), and \( u=-k_y y+k_r r \), derive the ideal gains.

Solution: Substitution gives

\[ \dot y=-(a+b k_y)y+b k_r r. \]

Matching the coefficients of \( y \) and \( r \) with the reference model yields \( a+b k_y=a_m \) and \( b k_r=b_m \). Therefore,

\[ k_y^{*}=\frac{a_m-a}{b}, \qquad k_r^{*}=\frac{b_m}{b}. \]

Problem 2 (Direct Lyapunov Cancellation): Starting from \( \dot e=-a_m e-b\tilde k_y y+b\tilde k_r r \), verify that the direct laws in Section 7 make \( \dot V_d=-a_m e^2 \).

Solution: Differentiate

\[ V_d = \frac{1}{2}e^2 + \frac{b}{2\gamma_y}\tilde k_y^2 + \frac{b}{2\gamma_r}\tilde k_r^2. \]

Then

\[ \begin{aligned} \dot V_d &= -a_m e^2 - b e\tilde k_y y + b e\tilde k_r r\\ &\quad + \frac{b}{\gamma_y}\tilde k_y(\gamma_y e y) + \frac{b}{\gamma_r}\tilde k_r(-\gamma_r e r)\\ &= -a_m e^2. \end{aligned} \]

The cancellation occurs because the adaptive laws are paired with the parameter-error terms in the tracking-error dynamics.

Problem 3 (Indirect Sensitivity): Compute the first order effect of small estimate errors \( \delta a,\delta b \) on \( \hat k_r=b_m/\hat b \).

Solution: Linearization about \( \hat b=b \) gives

\[ \delta k_r \approx \frac{\partial k_r}{\partial b}\delta b = -\frac{b_m}{b^2}\delta b. \]

The sensitivity magnitude is \( b_m/b^2 \), which becomes large as \( |b| \) becomes small. This explains why an indirect design must constrain the estimated input gain away from zero.

Problem 4 (Error-Signal Diagnosis): Suppose \( \varepsilon=y-\hat y=0 \) but \( e=y-y_m\neq 0 \). What can be concluded?

Solution: The estimated model predicts the measured plant along the current trajectory, but the closed-loop plant does not match the chosen reference model. Possible causes include an inappropriate redesign map, inaccurate controller implementation, a reference model incompatible with achievable matching, or estimates that reproduce the observed trajectory without producing the required controller gains. Zero prediction error is therefore not equivalent to zero tracking error.

Problem 5 (Adaptive-State Count): For the scalar example, count the additional adaptive states in each architecture, excluding the plant and reference-model states.

Solution: Direct MRAC has two adaptive controller states, \( \hat k_y,\hat k_r \). The indirect example has two parameter-estimate states, \( \hat a,\hat b \), plus the identifier state \( \hat y \). Its controller gains are algebraic outputs of the redesign map rather than independent adaptive states.

Direct and indirect MRAC share a reference model and model-following objective, but differ in their internal coordinates. Direct MRAC adapts controller gains from tracking information. Indirect MRAC adapts a plant model from prediction information and then redesigns the controller. Exact plant estimates recover the same ideal matching gains as a direct design, yet the adaptive transients generally differ because the two parameter trajectories obey different dynamics and because the indirect redesign map can amplify estimator errors. This structural distinction is essential before deriving detailed first- and second-order MRAC laws in later chapters.

1. Parks, P.C. (1966). Lyapunov redesign of model reference adaptive control systems. IEEE Transactions on Automatic Control, 11(3), 362–367.
2. Landau, I.D. (1974). A survey of model reference adaptive techniques—Theory and applications. Automatica, 10(4), 353–379.
3. Monopoli, R.V. (1974). Model reference adaptive control with an augmented error signal. IEEE Transactions on Automatic Control, 19(5), 474–484.
4. Feuer, A., & Morse, A.S. (1978). Adaptive control of single-input, single-output linear systems. IEEE Transactions on Automatic Control, 23(4), 557–569.
5. Narendra, K.S., & Valavani, L.S. (1978). Stable adaptive controller design—Direct control. IEEE Transactions on Automatic Control, 23(4), 570–583.
6. Narendra, K.S., & Valavani, L.S. (1979). Direct and indirect adaptive control. Automatica, 15(6), 653–664.
7. Morse, A.S. (1980). Global stability of parameter-adaptive control systems. IEEE Transactions on Automatic Control, 25(3), 433–439.
8. Narendra, K.S., Lin, Y.-H., & Valavani, L.S. (1980). Stable adaptive controller design, Part II: Proof of stability. IEEE Transactions on Automatic Control, 25(3), 440–448.