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

1. derive a discrete-time MRAC tracking-error equation from model-matching conditions;
2. distinguish tracking error from one-step equation or prediction error;
3. explain why continuous-time update laws cannot be copied directly into sampled systems;
4. derive a normalized adaptive law and prove monotonicity of a quadratic parameter-error function;
5. state the distinction between tracking convergence and parameter convergence; and
6. implement the controller with the correct sample-by-sample execution order.

We assume a fixed sampling interval, a deterministic plant, measurable plant output, bounded command input, a known calibrated nonzero input gain, and exact plant order. Quantization, noise, delay, projection, saturation, and robust modifications are addressed in later lessons.

At sample \( k \), measured signals form a regressor. The controller uses the parameter estimate stored at that sample to calculate the plant input. After the plant advances and \( y(k+1) \) is measured, the equation error becomes available and updates the parameters for the next sample.

flowchart TD
  R["Command r(k)"] --> RM["Stable reference model"]
  R --> REG["Regressor phi(k)"]
  Y["Measured plant output y(k)"] --> REG
  REG --> C["Adaptive controller u(k) = theta_hat(k)^T phi(k)"]
  C --> P["Unknown discrete-time plant"]
  P --> YN["Next measurement y(k+1)"]
  RM --> YM["Reference output ym(k+1)"]
  YN --> ERR["Tracking and equation errors"]
  YM --> ERR
  ERR --> A["Normalized parameter update"]
  A --> C
        

The use of \( e(k+1) \) is causal because the update is executed after the transition from sample \( k \) to sample \( k+1 \). The resulting estimate is used in \( u(k+1) \).

Consider the first-order discrete-time plant

\[ y(k+1)=a\,y(k)+b\,u(k), \qquad b\neq 0. \]

The coefficient \( a \) is uncertain, while the nonzero input coefficient \( b \) is known or calibrated for this introductory design. Select the stable reference model

\[ y_m(k+1)=a_m y_m(k)+b_m r(k), \qquad |a_m| < 1. \]

Use the direct adaptive controller

\[ u(k)=\hat\theta_y(k)y(k)+\hat\theta_r(k)r(k) =\hat{\boldsymbol\theta}^{T}(k)\boldsymbol\phi(k), \]

\[ \hat{\boldsymbol\theta}(k)= \begin{bmatrix} \hat\theta_y(k)\\ \hat\theta_r(k) \end{bmatrix}, \qquad \boldsymbol\phi(k)= \begin{bmatrix} y(k)\\ r(k) \end{bmatrix}. \]

This is a direct adaptive controller because its adjustable quantities are controller gains rather than explicit estimates of \( a \) and \( b \).

Suppose constant ideal gains exist such that the plant under the ideal controller has exactly the reference-model dynamics. Substituting \( u^*(k)=\theta_y^* y(k)+\theta_r^* r(k) \) gives

\[ y(k+1)=\left(a+b\theta_y^*\right)y(k) +b\theta_r^*r(k). \]

Coefficient matching requires

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

Hence

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

These gains are used only for analysis. The adaptive controller does not calculate them from plant parameters.

Define the tracking error and parameter error as

\[ e(k)=y(k)-y_m(k), \qquad \widetilde{\boldsymbol\theta}(k) =\hat{\boldsymbol\theta}(k)-\boldsymbol\theta^*. \]

Substitution of the adaptive controller into the plant gives

\[ y(k+1) = \left(a+b\hat\theta_y(k)\right)y(k) +b\hat\theta_r(k)r(k). \]

Add and subtract the ideal terms:

\[ \begin{aligned} y(k+1) &= \left(a+b\theta_y^*\right)y(k) +b\theta_r^*r(k)\\ &\quad+ b\widetilde{\boldsymbol\theta}^{T}(k)\boldsymbol\phi(k)\\ &= a_m y(k)+b_m r(k) +b\widetilde{\boldsymbol\theta}^{T}(k)\boldsymbol\phi(k). \end{aligned} \]

Subtract the reference-model equation:

\[ \boxed{ e(k+1)=a_m e(k) +b\widetilde{\boldsymbol\theta}^{T}(k)\boldsymbol\phi(k) }. \]

The reference pole defines the stable homogeneous part, while parameter mismatch acts as a forcing input.

Define the one-step equation error

\[ \varepsilon(k+1) =e(k+1)-a_m e(k). \]

Using the tracking-error equation,

\[ \boxed{ \varepsilon(k+1) =b\widetilde{\boldsymbol\theta}^{T}(k)\boldsymbol\phi(k) }. \]

This is an algebraic regression in the parameter error. It can be evaluated after \( y(k+1) \) is measured. The tracking error is the output of the stable filter driven by the sequence of equation errors:

\[ e(k) = a_m^k e(0) + \sum_{i=0}^{k-1} a_m^{k-1-i}\varepsilon(i+1). \]

Thus the two errors are related but not identical. Adaptation acts directly on the regression error, whereas tracking reflects the filtered history of that error.

Let \( V_\theta(k)=\widetilde{\boldsymbol\theta}^T(k) \Gamma^{-1}\widetilde{\boldsymbol\theta}(k) \), where \( \Gamma=\Gamma^T\succ 0 \). For a finite parameter increment \( \Delta\hat{\boldsymbol\theta}(k) \),

\[ \begin{aligned} V_\theta(k+1)-V_\theta(k) &= 2\widetilde{\boldsymbol\theta}^{T}(k) \Gamma^{-1}\Delta\hat{\boldsymbol\theta}(k)\\ &\quad+ \Delta\hat{\boldsymbol\theta}^{T}(k) \Gamma^{-1}\Delta\hat{\boldsymbol\theta}(k). \end{aligned} \]

The last term is nonnegative and has no direct counterpart in the derivative of a continuous-time quadratic Lyapunov function. Therefore, cancellation of only the linear cross term is insufficient. A normalization, an implicit update, or an augmented-error construction is needed.

Define the normalizer

\[ m^2(k) = 1+b^2\boldsymbol\phi^T(k) \Gamma\boldsymbol\phi(k). \]

The update law is

\[ \boxed{ \hat{\boldsymbol\theta}(k+1) = \hat{\boldsymbol\theta}(k) - \frac{ \Gamma b\boldsymbol\phi(k)\varepsilon(k+1) }{ 1+b^2\boldsymbol\phi^T(k) \Gamma\boldsymbol\phi(k) } }. \]

The denominator scales the correction inversely with regressor energy. A factor \( 0 < \eta < 2 \) can multiply the correction, but \( \eta=1 \) is used here to keep the proof transparent.

Let

\[ s(k)=b^2\boldsymbol\phi^T(k) \Gamma\boldsymbol\phi(k)\geq 0. \]

Since the ideal parameter is constant,

\[ \widetilde{\boldsymbol\theta}(k+1) = \widetilde{\boldsymbol\theta}(k) - \frac{ \Gamma b\boldsymbol\phi(k)\varepsilon(k+1) }{ 1+s(k) }. \]

The regression identity gives

\[ b\boldsymbol\phi^T(k) \widetilde{\boldsymbol\theta}(k) = \varepsilon(k+1). \]

Expanding the Lyapunov difference yields

\[ \begin{aligned} \Delta V_\theta(k) &= -\frac{2\varepsilon^2(k+1)}{1+s(k)} + \frac{s(k)\varepsilon^2(k+1)}{(1+s(k))^2}\\ &= -\frac{2+s(k)}{(1+s(k))^2} \varepsilon^2(k+1) \leq 0. \end{aligned} \]

Therefore \( V_\theta(k) \) is nonincreasing and the parameter estimate is bounded. Summing the Lyapunov differences gives

\[ \sum_{k=0}^{\infty} \frac{2+s(k)}{(1+s(k))^2} \varepsilon^2(k+1) < \infty. \]

If the regressor is bounded, the weighting factor is bounded away from zero. It follows that \( \varepsilon\in\ell_2 \) and \( \varepsilon(k) \) tends to zero.

The tracking dynamics are

\[ e(k+1)=a_m e(k)+\varepsilon(k+1). \]

Since \( |a_m| < 1 \), this is an exponentially stable linear system driven by the equation error. If \( \varepsilon(k) \) tends to zero, then

\[ \lim_{k\to\infty} e(k)=0. \]

One proof splits the convolution into an old finite part, whose contribution decays geometrically, and a recent tail, whose terms are uniformly small because the equation error tends to zero.

This result establishes output tracking under the stated boundedness conditions. It does not yet establish convergence of every controller parameter to its ideal value.

For a two-parameter controller, persistent excitation requires constants \( N > 0 \) and \( \alpha > 0 \) such that

\[ \sum_{i=k}^{k+N-1} \boldsymbol\phi(i)\boldsymbol\phi^T(i) \succeq \alpha I \qquad \text{for every } k. \]

This condition prevents a nonzero parameter-error direction from remaining invisible to the regression. Under boundedness and standard regularity assumptions, it strengthens equation-error convergence into convergence of \( \widetilde{\boldsymbol\theta}(k) \) to zero.

A constant command can produce excellent tracking without identifying both parameters. In steady state, the regressor approaches one fixed direction, so only one linear combination of the parameter errors is constrained.

For a fully measured plant and reference model,

\[ x(k+1)=A x(k)+B u(k), \qquad x_m(k+1)=A_m x_m(k)+B_m r(k). \]

Use

\[ u(k)=\hat K_x^T(k)x(k)+\hat K_r(k)r(k). \]

If ideal matrices satisfy

\[ A+B K_x^{*T}=A_m, \qquad B K_r^*=B_m, \]

then the state tracking error satisfies

\[ e_x(k+1) = A_m e_x(k) + B\left( \widetilde K_x^T(k)x(k) + \widetilde K_r(k)r(k) \right). \]

For higher-order systems, the measurable error generally does not immediately provide the scalar regression used in the first-order proof. Classical discrete-time MRAC designs therefore employ filtered regressors, augmented errors, positive-real transformations, or prediction-error identifiers.

flowchart TD
  A["Read y(k), r(k), and stored theta_hat(k)"] --> B["Form phi(k) = [y(k), r(k)]"]
  B --> C["Compute u(k)"]
  C --> D["Apply u(k) for one sample interval"]
  D --> E["Measure y(k+1) and update the reference model"]
  E --> F["Compute e(k), e(k+1), and epsilon(k+1)"]
  F --> G["Compute normalizer"]
  G --> H["Update theta_hat(k+1)"]
  H --> I["Store variables and advance k"]
  I --> A
        

An implementation that updates the parameters before \( y(k+1) \) is measured is noncausal. An extra one-sample delay is causal but defines a different adaptive system and requires a separate stability analysis.

The NumPy implementation simulates the plant, reference model, normalized update, tracking error, equation error, and parameter trajectories.

Chapter14_Lesson2.py

"""
Chapter14_Lesson2.py
Discrete-Time MRAC: error equations and normalized adaptation.
"""

from __future__ import annotations

import numpy as np
import matplotlib.pyplot as plt


def reference_signal(k: int) -> float:
    """Bounded command with multiple frequencies."""
    square = 1.0 if (k // 75) % 2 == 0 else -0.6
    return square + 0.2 * np.sin(0.05 * k)


def main() -> None:
    samples = 600

    # The simulator uses the true plant values. The controller assumes b is
    # known/calibrated, while a is not used by the adaptive law.
    a = 0.82
    b = 0.55

    # Stable reference model.
    am = 0.45
    bm = 0.55

    gamma = np.diag([0.50, 0.25])
    theta = np.array([0.0, 0.0], dtype=float)
    theta_star = np.array([(am - a) / b, bm / b], dtype=float)

    y = 0.0
    ym = 0.0

    r_hist = np.zeros(samples)
    y_hist = np.zeros(samples)
    ym_hist = np.zeros(samples)
    e_hist = np.zeros(samples)
    eps_hist = np.zeros(samples)
    u_hist = np.zeros(samples)
    theta_hist = np.zeros((samples, 2))

    for k in range(samples):
        r = reference_signal(k)
        phi = np.array([y, r], dtype=float)

        # Use theta(k) to compute u(k).
        u = float(theta @ phi)

        # Advance plant and reference model.
        y_next = a * y + b * u
        ym_next = am * ym + bm * r

        e = y - ym
        e_next = y_next - ym_next

        # Available after y(k+1) is measured.
        epsilon = e_next - am * e

        normalizer = 1.0 + b * b * float(phi @ gamma @ phi)
        theta = theta - (gamma @ (b * phi * epsilon)) / normalizer

        r_hist[k] = r
        y_hist[k] = y_next
        ym_hist[k] = ym_next
        e_hist[k] = e_next
        eps_hist[k] = epsilon
        u_hist[k] = u
        theta_hist[k, :] = theta

        y = y_next
        ym = ym_next

    print("Ideal parameters:", theta_star)
    print("Final estimates:", theta)
    print("Final tracking error:", e_hist[-1])
    print(
        "RMS tracking error over last 100 samples:",
        np.sqrt(np.mean(e_hist[-100:] ** 2)),
    )

    k_axis = np.arange(samples)

    plt.figure()
    plt.plot(k_axis, y_hist, label="Plant output y")
    plt.plot(k_axis, ym_hist, "--", label="Reference output ym")
    plt.plot(k_axis, r_hist, ":", label="Command r")
    plt.xlabel("Sample k")
    plt.ylabel("Signal")
    plt.title("Discrete-Time MRAC Tracking")
    plt.grid(True)
    plt.legend()

    plt.figure()
    plt.plot(k_axis, e_hist, label="Tracking error e")
    plt.plot(k_axis, eps_hist, label="Equation error epsilon")
    plt.xlabel("Sample k")
    plt.ylabel("Error")
    plt.title("Tracking and Equation Errors")
    plt.grid(True)
    plt.legend()

    plt.figure()
    plt.plot(k_axis, theta_hist[:, 0], label="theta_y")
    plt.plot(k_axis, theta_hist[:, 1], label="theta_r")
    plt.axhline(theta_star[0], linestyle="--", label="theta_y ideal")
    plt.axhline(theta_star[1], linestyle="--", label="theta_r ideal")
    plt.xlabel("Sample k")
    plt.ylabel("Parameter")
    plt.title("Adaptive Controller Parameters")
    plt.grid(True)
    plt.legend()

    plt.show()


if __name__ == "__main__":
    main()

This C++17 implementation uses only the standard library and exports all sample histories to CSV.

Chapter14_Lesson2.cpp

/*
Chapter14_Lesson2.cpp
Discrete-time first-order MRAC with normalized adaptation.

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

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

double referenceSignal(int k) {
    const double square = ((k / 75) % 2 == 0) ? 1.0 : -0.6;
    return square + 0.2 * std::sin(0.05 * static_cast<double>(k));
}

int main() {
    constexpr int samples = 600;

    const double a = 0.82;
    const double b = 0.55;
    const double am = 0.45;
    const double bm = 0.55;

    const double gammaY = 0.50;
    const double gammaR = 0.25;

    std::array<double, 2> theta{0.0, 0.0};
    const std::array<double, 2> thetaStar{(am - a) / b, bm / b};

    double y = 0.0;
    double ym = 0.0;
    double squaredTailError = 0.0;

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

    csv << "k,r,y,ym,e,epsilon,u,theta_y,theta_r\n";
    csv << std::setprecision(12);

    for (int k = 0; k < samples; ++k) {
        const double r = referenceSignal(k);
        const double phiY = y;
        const double phiR = r;

        const double u = theta[0] * phiY + theta[1] * phiR;

        const double yNext = a * y + b * u;
        const double ymNext = am * ym + bm * r;

        const double e = y - ym;
        const double eNext = yNext - ymNext;
        const double epsilon = eNext - am * e;

        const double quadratic =
            gammaY * phiY * phiY + gammaR * phiR * phiR;
        const double normalizer = 1.0 + b * b * quadratic;

        theta[0] -= gammaY * b * phiY * epsilon / normalizer;
        theta[1] -= gammaR * b * phiR * epsilon / normalizer;

        csv << k << ',' << r << ',' << yNext << ',' << ymNext << ','
            << eNext << ',' << epsilon << ',' << u << ','
            << theta[0] << ',' << theta[1] << '\n';

        if (k >= samples - 100) {
            squaredTailError += eNext * eNext;
        }

        y = yNext;
        ym = ymNext;
    }

    std::cout << "Ideal theta_y = " << thetaStar[0]
              << ", ideal theta_r = " << thetaStar[1] << '\n';
    std::cout << "Final theta_y = " << theta[0]
              << ", final theta_r = " << theta[1] << '\n';
    std::cout << "Tail RMS tracking error = "
              << std::sqrt(squaredTailError / 100.0) << '\n';

    return 0;
}

The Java program follows the same causal sample ordering and writes a CSV trajectory for plotting or verification.

Chapter14_Lesson2.java

/*
Chapter14_Lesson2.java
Discrete-time first-order MRAC with normalized adaptation.

Compile:
    javac Chapter14_Lesson2.java
Run:
    java Chapter14_Lesson2
*/

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Locale;

public final class Chapter14_Lesson2 {
    private Chapter14_Lesson2() {
    }

    private static double referenceSignal(int k) {
        double square = ((k / 75) % 2 == 0) ? 1.0 : -0.6;
        return square + 0.2 * Math.sin(0.05 * k);
    }

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

        final int samples = 600;

        final double a = 0.82;
        final double b = 0.55;
        final double am = 0.45;
        final double bm = 0.55;

        final double gammaY = 0.50;
        final double gammaR = 0.25;

        double thetaY = 0.0;
        double thetaR = 0.0;

        final double thetaYStar = (am - a) / b;
        final double thetaRStar = bm / b;

        double y = 0.0;
        double ym = 0.0;
        double squaredTailError = 0.0;

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter14_Lesson2_results_java.csv"))) {

            writer.write("k,r,y,ym,e,epsilon,u,theta_y,theta_r");
            writer.newLine();

            for (int k = 0; k < samples; k++) {
                double r = referenceSignal(k);
                double phiY = y;
                double phiR = r;

                double u = thetaY * phiY + thetaR * phiR;

                double yNext = a * y + b * u;
                double ymNext = am * ym + bm * r;

                double e = y - ym;
                double eNext = yNext - ymNext;
                double epsilon = eNext - am * e;

                double quadratic =
                        gammaY * phiY * phiY + gammaR * phiR * phiR;
                double normalizer = 1.0 + b * b * quadratic;

                thetaY -= gammaY * b * phiY * epsilon / normalizer;
                thetaR -= gammaR * b * phiR * epsilon / normalizer;

                writer.write(String.format(
                        "%d,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f",
                        k, r, yNext, ymNext, eNext, epsilon, u,
                        thetaY, thetaR));
                writer.newLine();

                if (k >= samples - 100) {
                    squaredTailError += eNext * eNext;
                }

                y = yNext;
                ym = ymNext;
            }
        } catch (IOException exception) {
            System.err.println("I/O failure: " + exception.getMessage());
            System.exit(1);
        }

        System.out.printf("Ideal theta_y = %.8f, ideal theta_r = %.8f%n",
                thetaYStar, thetaRStar);
        System.out.printf("Final theta_y = %.8f, final theta_r = %.8f%n",
                thetaY, thetaR);
        System.out.printf("Tail RMS tracking error = %.10f%n",
                Math.sqrt(squaredTailError / 100.0));
    }
}

The MATLAB script performs an array-based simulation. The second script programmatically creates a fixed-step Simulink model with a MATLAB Function block containing the plant, model, and adaptive update.

Chapter14_Lesson2.m

% Chapter14_Lesson2.m
% Discrete-Time MRAC: scalar error equation and normalized adaptation.

clear;
clc;
close all;

samples = 600;

a = 0.82;
b = 0.55;

am = 0.45;
bm = 0.55;

Gamma = diag([0.50, 0.25]);

theta = [0; 0];
theta_star = [(am - a) / b; bm / b];

y = 0;
ym = 0;

r_hist = zeros(samples, 1);
y_hist = zeros(samples, 1);
ym_hist = zeros(samples, 1);
e_hist = zeros(samples, 1);
epsilon_hist = zeros(samples, 1);
u_hist = zeros(samples, 1);
theta_hist = zeros(samples, 2);

for index = 1:samples
    k = index - 1;

    if mod(floor(k / 75), 2) == 0
        square = 1.0;
    else
        square = -0.6;
    end
    r = square + 0.2 * sin(0.05 * k);

    phi = [y; r];
    u = theta.' * phi;

    y_next = a * y + b * u;
    ym_next = am * ym + bm * r;

    e = y - ym;
    e_next = y_next - ym_next;
    epsilon = e_next - am * e;

    normalizer = 1 + b^2 * (phi.' * Gamma * phi);
    theta = theta - Gamma * (b * phi * epsilon) / normalizer;

    r_hist(index) = r;
    y_hist(index) = y_next;
    ym_hist(index) = ym_next;
    e_hist(index) = e_next;
    epsilon_hist(index) = epsilon;
    u_hist(index) = u;
    theta_hist(index, :) = theta.';

    y = y_next;
    ym = ym_next;
end

fprintf('Ideal parameters: theta_y = %.8f, theta_r = %.8f\n', ...
    theta_star(1), theta_star(2));
fprintf('Final estimates: theta_y = %.8f, theta_r = %.8f\n', ...
    theta(1), theta(2));
fprintf('RMS error over last 100 samples = %.10f\n', ...
    sqrt(mean(e_hist(end-99:end).^2)));

k_axis = 0:(samples - 1);

figure;
plot(k_axis, y_hist, 'LineWidth', 1.2);
hold on;
plot(k_axis, ym_hist, '--', 'LineWidth', 1.2);
plot(k_axis, r_hist, ':', 'LineWidth', 1.0);
grid on;
xlabel('Sample k');
ylabel('Signal');
title('Discrete-Time MRAC Tracking');
legend('Plant output y', 'Reference output y_m', 'Command r');

figure;
plot(k_axis, e_hist, 'LineWidth', 1.2);
hold on;
plot(k_axis, epsilon_hist, 'LineWidth', 1.0);
grid on;
xlabel('Sample k');
ylabel('Error');
title('Tracking Error and Equation Error');
legend('e', '\epsilon');

figure;
plot(k_axis, theta_hist(:, 1), 'LineWidth', 1.2);
hold on;
plot(k_axis, theta_hist(:, 2), 'LineWidth', 1.2);
yline(theta_star(1), '--');
yline(theta_star(2), '--');
grid on;
xlabel('Sample k');
ylabel('Parameter');
title('Adaptive Controller Parameters');
legend('\theta_y', '\theta_r', '\theta_y^*', '\theta_r^*');

Chapter14_Lesson2_Simulink.m

% Chapter14_Lesson2_Simulink.m
% Creates a fixed-step Simulink model.
% Requires Simulink.

clear;
clc;

model = 'Chapter14_Lesson2_Simulink_Model';

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

if exist([model '.slx'], 'file')
    delete([model '.slx']);
end

new_system(model);
open_system(model);

add_block('simulink/Sources/In1', [model '/Reference'], ...
    'Position', [40 110 70 130]);

add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [model '/Discrete MRAC'], ...
    'Position', [150 50 390 210]);

outputNames = {'Plant Output', 'Model Output', 'Tracking Error', ...
               'Theta Y', 'Theta R'};

for i = 1:numel(outputNames)
    add_block('simulink/Sinks/Out1', [model '/' outputNames{i}], ...
        'Position', [500 30 + 45*i 530 50 + 45*i]);
end

root = sfroot;
chart = root.find('-isa', 'Stateflow.EMChart', ...
    'Path', [model '/Discrete MRAC']);

chart.Script = sprintf([ ...
'function [yOut, ymOut, eOut, thetaYOut, thetaROut] = fcn(r)\n' ...
'%%#codegen\n' ...
'persistent y ym thetaY thetaR\n' ...
'if isempty(y)\n' ...
'    y = 0.0;\n' ...
'    ym = 0.0;\n' ...
'    thetaY = 0.0;\n' ...
'    thetaR = 0.0;\n' ...
'end\n' ...
'a = 0.82;\n' ...
'b = 0.55;\n' ...
'am = 0.45;\n' ...
'bm = 0.55;\n' ...
'gammaY = 0.50;\n' ...
'gammaR = 0.25;\n' ...
'u = thetaY*y + thetaR*r;\n' ...
'yNext = a*y + b*u;\n' ...
'ymNext = am*ym + bm*r;\n' ...
'e = y - ym;\n' ...
'eNext = yNext - ymNext;\n' ...
'epsilon = eNext - am*e;\n' ...
'normalizer = 1.0 + b*b*(gammaY*y*y + gammaR*r*r);\n' ...
'thetaY = thetaY - gammaY*b*y*epsilon/normalizer;\n' ...
'thetaR = thetaR - gammaR*b*r*epsilon/normalizer;\n' ...
'y = yNext;\n' ...
'ym = ymNext;\n' ...
'yOut = y;\n' ...
'ymOut = ym;\n' ...
'eOut = eNext;\n' ...
'thetaYOut = thetaY;\n' ...
'thetaROut = thetaR;\n' ...
'end\n']);

set_param(model, 'SolverType', 'Fixed-step', ...
    'Solver', 'FixedStepDiscrete', ...
    'FixedStep', '1', ...
    'StopTime', '599');

set_param(model, 'SimulationCommand', 'update');

add_line(model, 'Reference/1', 'Discrete MRAC/1', 'autorouting', 'on');
for i = 1:numel(outputNames)
    add_line(model, sprintf('Discrete MRAC/%d', i), ...
        sprintf('%s/1', outputNames{i}), 'autorouting', 'on');
end

save_system(model);
open_system(model);

fprintf('Created %s.slx\n', model);
fprintf('Drive the Reference input with a bounded discrete signal.\n');

The notebook executes the same recursion and produces tracking, error, and parameter plots.

Chapter14_Lesson2.nb

Notebook[{
 Cell["Chapter 14, Lesson 2: Discrete-Time MRAC", "Title"],
 Cell["Normalized one-step equation-error adaptation.", "Text"],
 Cell[BoxData[
"ClearAll[\"Global`*\"];\n\
samples = 600;\n\
a = 0.82; b = 0.55;\n\
am = 0.45; bm = 0.55;\n\
gammaY = 0.50; gammaR = 0.25;\n\
thetaY = 0.; thetaR = 0.;\n\
thetaYStar = (am - a)/b;\n\
thetaRStar = bm/b;\n\
y = 0.; ym = 0.;\n\
reference[k_] := If[EvenQ[Floor[k/75]], 1.0, -0.6] + 0.2 Sin[0.05 k];\n\
history = Reap[\n\
   Do[\n\
    r = reference[k];\n\
    u = thetaY y + thetaR r;\n\
    yNext = a y + b u;\n\
    ymNext = am ym + bm r;\n\
    e = y - ym;\n\
    eNext = yNext - ymNext;\n\
    epsilon = eNext - am e;\n\
    normalizer = 1 + b^2 (gammaY y^2 + gammaR r^2);\n\
    thetaY = thetaY - gammaY b y epsilon/normalizer;\n\
    thetaR = thetaR - gammaR b r epsilon/normalizer;\n\
    Sow[{k, r, yNext, ymNext, eNext, epsilon, u, thetaY, thetaR}];\n\
    y = yNext;\n\
    ym = ymNext,\n\
    {k, 0, samples - 1}\n\
    ]\n\
   ][[2, 1]];\n\
Print[\"Ideal parameters: \", {thetaYStar, thetaRStar}];\n\
Print[\"Final estimates: \", {thetaY, thetaR}];\n\
Print[\"Tail RMS error: \", Sqrt[Mean[history[[-100 ;;, 5]]^2]]];\n\
trackingPlot = ListLinePlot[\n\
   {history[[All, {1, 3}]], history[[All, {1, 4}]], history[[All, {1, 2}]]},\n\
   PlotLegends -> {\"Plant output y\", \"Reference output ym\", \"Command r\"},\n\
   AxesLabel -> {\"k\", \"signal\"}, PlotLabel -> \"Discrete-Time MRAC Tracking\"];\n\
errorPlot = ListLinePlot[\n\
   {history[[All, {1, 5}]], history[[All, {1, 6}]]},\n\
   PlotLegends -> {\"Tracking error e\", \"Equation error epsilon\"},\n\
   AxesLabel -> {\"k\", \"error\"}, PlotLabel -> \"Error Signals\"];\n\
parameterPlot = ListLinePlot[\n\
   {history[[All, {1, 8}]], history[[All, {1, 9}]]},\n\
   PlotLegends -> {\"thetaY\", \"thetaR\"},\n\
   AxesLabel -> {\"k\", \"parameter\"}, PlotLabel -> \"Adaptive Parameters\"];\n\
Column[{trackingPlot, errorPlot, parameterPlot}]"
 ], "Input"]
},
WindowTitle -> "Chapter14_Lesson2",
StyleDefinitions -> "Default.nb"
]

Reference pole: Smaller \( |a_m| \) requests faster model dynamics. The ideal feedback gain may therefore have larger magnitude and demand greater control effort.

Adaptation matrix: Larger diagonal entries of \( \Gamma \) accelerate motion in the corresponding parameter directions. The normalizer controls the effective step size when the regressor is large, but aggressive gains can still create undesirable transients.

Error plots: The Lyapunov proof acts directly on \( \varepsilon(k) \). The tracking error may decay more slowly because it is the output of the stable reference-model error filter.

Parameter plots: Estimates that stop away from the ideal values do not necessarily indicate failure. Without persistent excitation, several gain vectors can be indistinguishable on the observed trajectory.

Input-gain information: The proved update uses the calibrated magnitude of \( b \). A sign-only design requires a gain bound, a rescaled normalization, or an augmented-error construction. Completely unknown control direction is treated separately with Nussbaum-type methods.

Problem 1 — Matching gains: Consider \( y(k+1)=0.9y(k)+0.4u(k) \) and \( y_m(k+1)=0.5y_m(k)+0.4r(k) \). Determine the ideal controller gains.

Solution:

\[ \theta_y^* = \frac{0.5-0.9}{0.4} = -1, \qquad \theta_r^* = \frac{0.4}{0.4} = 1. \]

The ideal controller is \( u(k)=-y(k)+r(k) \). Substitution produces the desired reference-model coefficients.

Problem 2 — Measurable equation error: Derive a one-step regression from the tracking-error equation.

Solution:

\[ \varepsilon(k+1) \equiv e(k+1)-a_m e(k) = b\widetilde{\boldsymbol\theta}^{T}(k) \boldsymbol\phi(k). \]

Both tracking errors on the left are known after \( y(k+1) \) is measured, so the regression is causal for computing the next parameter estimate.

Problem 3 — Quadratic update term: Verify the nonnegative quadratic term in the Lyapunov difference.

Solution: With

\[ \Delta\hat{\boldsymbol\theta} = - \frac{ \Gamma b\boldsymbol\phi\varepsilon }{ 1+s }, \qquad s=b^2\boldsymbol\phi^T\Gamma\boldsymbol\phi, \]

we obtain

\[ \begin{aligned} \Delta\hat{\boldsymbol\theta}^T \Gamma^{-1} \Delta\hat{\boldsymbol\theta} &= \frac{ b^2\varepsilon^2 \boldsymbol\phi^T\Gamma\boldsymbol\phi }{ (1+s)^2 }\\ &= \frac{s\varepsilon^2}{(1+s)^2}. \end{aligned} \]

Problem 4 — Constant command: Explain why exact tracking for a constant reference does not guarantee convergence of both controller parameters.

Solution:

At steady state the regressor approaches a fixed vector. Vanishing equation error then imposes only

\[ \widetilde{\boldsymbol\theta}^{T} \boldsymbol\phi_{ss}=0. \]

This is one scalar constraint on two unknown parameter-error coordinates. A time-varying persistently exciting regressor is needed for unique identification.

Problem 5 — Uniform tracking bound: Suppose \( |a_m| < 1 \) and \( |\varepsilon(k)|\leq\bar\varepsilon \). Derive a uniform bound on the tracking error.

Solution:

\[ \begin{aligned} |e(k)| &\leq |a_m|^k|e(0)| + \bar\varepsilon \sum_{j=0}^{k-1}|a_m|^j\\ &\leq |e(0)| + \frac{\bar\varepsilon}{1-|a_m|}. \end{aligned} \]

The stable reference-model pole therefore gives a finite input-to-state bound from equation error to tracking error.

Problem 6 — Programming experiment: Replace the multi-level command in the supplied programs with a constant step and compare terminal tracking and parameter errors.

Solution:

Tracking should become small, but the estimates generally stop away from the ideal vector. The constant command fails to excite two independent regressor directions, demonstrating that successful control and exact parameter identification are different objectives.

The discrete-time MRAC error model separates stable reference-model dynamics from forcing caused by controller-parameter mismatch. Subtracting the stable error recursion gives a causal one-step equation error that is linear in the parameter error. A normalized gradient update controls the finite quadratic parameter increment and makes a quadratic parameter-error function nonincreasing. With a bounded regressor, equation error converges to zero; stability of the reference-model error filter then gives tracking convergence. Persistent excitation is additionally required for unique convergence of the controller parameters.