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

• distinguish tracking-error convergence from parameter boundedness and parameter convergence;
• derive the integral mechanism by which a small update residual produces drift;
• explain why bounded or vanishing signals need not have a bounded time integral;
• identify physical and numerical consequences of implausible adaptive parameters;
• state the invariance and Lyapunov-compatibility requirements that motivate projection;
• compare projection conceptually with the robust modifications introduced in Chapter 8.

Recall from the earlier MRAC chapters that the adaptive parameter vector is not generally an independently measured physical quantity. It is an internal controller state generated by integration. Therefore, any persistent component entering its update law can accumulate over time.

Consider the uncertain first-order plant and stable reference model

\[ \dot{x}=-a x+b u+d(t), \qquad \dot{x}_m=-a_m x_m+b_m r, \]

where \(a_m>0\), the sign of \(b\) is known, and \(d(t)\) collects disturbances and modeling mismatch. Use the direct adaptive control law

\[ u=\hat{k}_x x+\hat{k}_r r =\hat{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}, \qquad \boldsymbol{\phi}= \begin{bmatrix}x\\r\end{bmatrix}. \]

If ideal parameters exist, they satisfy

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

Define the tracking error \(e=x-x_m\) and parameter error \(\tilde{\boldsymbol{\theta}} =\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\). The error dynamics become

\[ \dot{e}=-a_m e+ b\,\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}+d(t). \]

A standard gradient law with known control direction is

\[ \dot{\hat{\boldsymbol{\theta}}} =-\boldsymbol{\Gamma}\boldsymbol{\phi}e\,\operatorname{sgn}(b), \qquad \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}>0. \]

Parameter drift is the slow or rapid migration of an adaptive estimate away from a physically meaningful or design-valid region, usually because the update law continues integrating a residual that does not represent true parametric error. Drift does not require a large instantaneous error. A very small bias can be dangerous because adaptation is an accumulation process.

\[ \hat{\boldsymbol{\theta}}(t) = \hat{\boldsymbol{\theta}}(0) - \int_{0}^{t} \boldsymbol{\Gamma} \boldsymbol{\phi}(\tau)e(\tau) \operatorname{sgn}(b)\,d\tau. \]

Let the update-channel residual be \(\mathbf{w}(t)= \boldsymbol{\phi}(t)e(t)\operatorname{sgn}(b)\). Then

\[ \hat{\boldsymbol{\theta}}(t) = \hat{\boldsymbol{\theta}}(0) - \boldsymbol{\Gamma}\int_{0}^{t}\mathbf{w}(\tau)\,d\tau. \]

Thus, boundedness of \(\mathbf{w}\) is not sufficient. What directly prevents drift in this integral representation is a stronger property such as absolute integrability:

\[ \int_{0}^{\infty}\|\mathbf{w}(\tau)\|\,d\tau < \infty. \]

flowchart TD
  A["Small disturbance, sensor offset, or model mismatch"] --> B["Nonzero residual in tracking or prediction error"]
  B --> C["Residual multiplied by regressor"]
  C --> D["Adaptive law integrates the update signal"]
  D --> E["Parameter estimate moves outside a credible region"]
  E --> F["Large gain, wrong sign, saturation, noise amplification, or instability"]
  F --> G["Need an admissible parameter set and inward boundary logic"]
        

4.1 Nonzero-Mean Residual: Linear Drift

In a scalar update channel, suppose the residual has decomposition \(w(t)=\bar{w}+\tilde{w}(t)\), where \(\bar{w}\neq 0\) and the oscillatory component has a bounded integral. Then

\[ \hat{\theta}(t) = \hat{\theta}(0)-\gamma\bar{w}t -\gamma\int_{0}^{t}\tilde{w}(\tau)\,d\tau. \]

The first term after the initial condition grows linearly. Therefore, even an arbitrarily small constant sensor offset produces unbounded drift if it enters an unconstrained pure-integral update law.

4.2 Vanishing and Square-Integrable Residual: Logarithmic Drift

A residual may converge to zero and still have an unbounded integral. Take \(w(t)=1/(1+t)\). It satisfies

\[ \lim_{t\to\infty}w(t)=0, \qquad \int_{0}^{\infty}w^{2}(t)\,dt=1, \]

but

\[ \int_{0}^{t}w(\tau)\,d\tau=\log(1+t). \]

Hence \(\hat{\theta}(t)=\hat{\theta}(0)-\gamma\log(1+t)\) drifts without bound. This counterexample is important: convergence of the tracking error to zero, or membership of an update signal in \(\mathcal{L}_{2}\), does not by itself make the parameter update absolutely integrable.

4.3 Zero-Mean Discrete Noise: Mean-Square Drift

Consider the sampled update

\[ \hat{\theta}_{k+1} = \hat{\theta}_{k} - h\gamma n_k, \]

where \(n_k\) are independent, zero-mean samples with variance \(\sigma_n^2\). Then

\[ \mathbb{E}[\hat{\theta}_{N}] = \hat{\theta}_{0}, \qquad \operatorname{Var}(\hat{\theta}_{N}) = N h^{2}\gamma^{2}\sigma_n^{2}. \]

The expected estimate does not move, but its variance increases linearly with time. The estimate behaves as a random walk, so zero-mean noise does not imply a bounded realization.

For the MRAC error model, consider the familiar candidate

\[ V= \frac{1}{2}e^{2} + \frac{|b|}{2} \tilde{\boldsymbol{\theta}}^{T} \boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta}}. \]

Substitution of the nominal adaptive law cancels the cross term and gives

\[ \dot{V}=-a_m e^{2}+e\,d(t). \]

With \(|d(t)|\le\bar{d}\), Young's inequality yields

\[ e\,d \le \frac{a_m}{2}e^{2} + \frac{\bar{d}^{2}}{2a_m}, \]

and therefore

\[ \dot{V} \le -\frac{a_m}{2}e^{2} + \frac{\bar{d}^{2}}{2a_m}. \]

This inequality can support an ultimate-bound argument for the tracking error under additional conditions, but it does not provide a negative term in \(\tilde{\boldsymbol{\theta}}\). In particular, it permits the conservative estimate

\[ V(t) \le V(0)+\frac{\bar{d}^{2}}{2a_m}t, \]

which does not rule out parameter growth. The cancellation that is elegant in the ideal model also removes the direct damping of parameter error. Persistent excitation, robust modifications, or explicit constraints are needed for stronger conclusions.

In direct adaptive control, the estimates are controller gains. When they leave the design-valid region, several failures can occur:

• loss of control direction: an estimated gain may acquire a sign incompatible with the plant;
• actuator saturation: large gains create commands that cannot be realized;
• noise amplification: high feedback gains transmit measurement noise into the actuator;
• excitation of neglected dynamics: aggressive control can excite high-frequency modes omitted from the nominal model;
• invalid certainty-equivalent redesign: an indirect controller may be recomputed from an estimate for which the algebraic design equations are singular or non-stabilizing;
• poor restart behavior: stored implausible estimates can cause a dangerous transient after reboot or mode switching.

Consequently, bounded tracking over a short simulation interval is not sufficient evidence that the adaptive controller is safe for prolonged operation.

1. Measurement bias: a constant offset can create a nonzero-mean adaptation signal.
2. Unmodeled dynamics: the adaptive law may interpret neglected poles, zeros, delays, or flexible modes as parameter error.
3. External disturbances: load torque, aerodynamic gusts, or process disturbances may be absorbed incorrectly into the estimated parameter vector.
4. Insufficient excitation: when regressors do not contain enough independent information, many parameter vectors can yield similar tracking behavior.
5. Actuator saturation: the controller assumes that the requested input was applied, while the plant receives a clipped input.
6. Sampling and numerical integration: finite step size, quantization, asynchronous sensing, and roundoff can introduce a systematic update bias.

In many systems, exact parameters are unknown but credible bounds are available from physics, identification data, manufacturing tolerances, or sign information. Represent this knowledge by a closed convex set

\[ \Omega= \left\{ \boldsymbol{\theta}\in\mathbb{R}^{p} : f(\boldsymbol{\theta})\le 0 \right\}. \]

A common special case is the componentwise box

\[ \Omega= [\theta_{1,\min},\theta_{1,\max}] \times\cdots\times [\theta_{p,\min},\theta_{p,\max}]. \]

The purpose of projection is not to force the estimate to equal the true parameter. Its immediate purpose is to ensure that the estimate remains in a region where the controller structure, sign assumptions, gain limits, and stability analysis remain meaningful.

flowchart TD
  A["Estimate inside admissible set"] --> B["Nominal update points \ninward or tangent"]
  B --> C["Apply nominal update"]
  A --> D["Estimate at boundary"]
  D --> E["Update points \ninward or tangent"]
  E --> C
  D --> F["Update points outward"]
  F --> G["Remove outward component"]
  G --> H["Estimate remains admissible"]
        

Before the multidimensional definition in Lesson 2, consider a scalar estimate constrained to \(\theta_{\min}\le\hat{\theta}\le\theta_{\max}\). Let \(y(t)\) denote the nominal update direction. The intuitive boundary rule is

\[ \dot{\hat{\theta}}= \begin{cases} 0, & \hat{\theta}=\theta_{\min} \text{ and }y<0,\\ 0, & \hat{\theta}=\theta_{\max} \text{ and }y>0,\\ y, & \text{otherwise}. \end{cases} \]

The nominal law is unchanged in the interior. At a boundary, only motion that would leave the admissible interval is suppressed. This produces forward invariance:

\[ \hat{\theta}(0)\in[\theta_{\min},\theta_{\max}] \quad\Longrightarrow\quad \hat{\theta}(t)\in[\theta_{\min},\theta_{\max}] \quad\text{for all }t\ge 0. \]

For a vector estimate, simply clipping every component is useful for a rectangular set but does not capture the geometry of a general convex boundary or the exact Lyapunov inequality. Lesson 2 develops the proper operator.

Suppose the ideal parameter lies in the admissible set: \(\boldsymbol{\theta}^{*}\in\Omega\). If \(\mathbf{y}\) is the nominal update direction, a projection mechanism should satisfy an inequality of the form

\[ (\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*})^{T} \boldsymbol{\Gamma}^{-1} \left[ \operatorname{Proj}(\hat{\boldsymbol{\theta}},\mathbf{y}) - \mathbf{y} \right] \le 0. \]

This property says that replacing the nominal update by its projected version does not add a harmful positive term to the standard parameter-error part of the Lyapunov derivative. The complete construction and proof are deferred to Lesson 2; here the inequality explains why arbitrary post-processing or ad hoc clipping is not always theoretically interchangeable with a projection operator.

These methods can be combined. Projection answers a geometric question: where may the estimate go? Leakage and dead zones answer a dynamical question: when and how strongly should adaptation move?

To expose drift without hiding it inside a larger closed-loop simulation, the implementations below isolate the parameter-update channel

\[ \dot{\hat{\theta}}=-\gamma w(t), \]

with the bounded residual

\[ w(t)=0.030 +0.018\sin(3t) +0.012\sin(17t) +\frac{0.008}{1+t}. \]

The positive mean component causes the unconstrained estimate to drift downward approximately linearly. The constrained experiment uses the scalar boundary logic of Section 9 with \(\Omega=[-2,2]\). This is intentionally a diagnostic experiment, not a complete MRAC design.

Chapter9_Lesson1.py

"""
Chapter9_Lesson1.py
Parameter-drift diagnostic and a projection-inspired scalar boundary guard.

Dependencies:
    pip install numpy matplotlib
"""

from __future__ import annotations

import numpy as np
import matplotlib.pyplot as plt


def boundary_guard(theta: float, raw_rate: float, lower: float, upper: float) -> float:
    """Suppress only an outward update at an active scalar boundary."""
    if theta <= lower and raw_rate < 0.0:
        return 0.0
    if theta >= upper and raw_rate > 0.0:
        return 0.0
    return raw_rate


def simulate(use_guard: bool, dt: float = 0.002, final_time: float = 120.0):
    gamma = 1.5
    theta0 = 0.7
    lower, upper = -2.0, 2.0

    time = np.arange(0.0, final_time + dt, dt)
    theta = np.empty_like(time)
    residual = np.empty_like(time)
    theta[0] = theta0

    for k, t in enumerate(time[:-1]):
        # A small nonzero bias plus bounded high-frequency contamination.
        # This signal models the residual entering a gradient update channel.
        residual[k] = (
            0.030
            + 0.018 * np.sin(3.0 * t)
            + 0.012 * np.sin(17.0 * t)
            + 0.008 / (1.0 + t)
        )

        raw_rate = -gamma * residual[k]
        applied_rate = (
            boundary_guard(theta[k], raw_rate, lower, upper)
            if use_guard
            else raw_rate
        )

        theta[k + 1] = theta[k] + dt * applied_rate

        # Numerical safeguard for the constrained experiment only.
        if use_guard:
            theta[k + 1] = np.clip(theta[k + 1], lower, upper)

    residual[-1] = residual[-2]
    return time, theta, residual


def main() -> None:
    time, theta_free, residual = simulate(use_guard=False)
    _, theta_guarded, _ = simulate(use_guard=True)

    print(f"Unconstrained final estimate: {theta_free[-1]: .6f}")
    print(f"Guarded final estimate:       {theta_guarded[-1]: .6f}")
    print(f"Residual time average:        {np.trapz(residual, time) / time[-1]: .6f}")

    plt.figure(figsize=(10, 5))
    plt.plot(time, theta_free, label="unconstrained update")
    plt.plot(time, theta_guarded, label="boundary-guarded update")
    plt.axhline(-2.0, linestyle="--", label="lower admissible bound")
    plt.axhline(2.0, linestyle="--", label="upper admissible bound")
    plt.xlabel("time (s)")
    plt.ylabel("parameter estimate")
    plt.title("Parameter drift caused by a small persistent residual")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

Chapter9_Lesson1.cpp

/*
Chapter9_Lesson1.cpp
Parameter-drift diagnostic and a projection-inspired scalar boundary guard.

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

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

double boundary_guard(
    const double theta,
    const double raw_rate,
    const double lower,
    const double upper
) {
    if (theta <= lower && raw_rate < 0.0) {
        return 0.0;
    }
    if (theta >= upper && raw_rate > 0.0) {
        return 0.0;
    }
    return raw_rate;
}

double residual_signal(const double t) {
    return 0.030
        + 0.018 * std::sin(3.0 * t)
        + 0.012 * std::sin(17.0 * t)
        + 0.008 / (1.0 + t);
}

int main() {
    try {
        constexpr double dt = 0.002;
        constexpr double final_time = 120.0;
        constexpr double gamma = 1.5;
        constexpr double lower = -2.0;
        constexpr double upper = 2.0;

        double theta_free = 0.7;
        double theta_guarded = 0.7;
        double residual_integral = 0.0;

        std::ofstream csv("Chapter9_Lesson1.csv");
        if (!csv) {
            throw std::runtime_error("Cannot create Chapter9_Lesson1.csv");
        }

        csv << "time,residual,theta_free,theta_guarded\n";
        csv << std::setprecision(12);

        const int steps = static_cast<int>(final_time / dt);
        for (int k = 0; k <= steps; ++k) {
            const double t = k * dt;
            const double residual = residual_signal(t);
            const double raw_rate = -gamma * residual;

            csv << t << ',' << residual << ','
                << theta_free << ',' << theta_guarded << '\n';

            if (k == steps) {
                break;
            }

            theta_free += dt * raw_rate;

            const double guarded_rate =
                boundary_guard(theta_guarded, raw_rate, lower, upper);
            theta_guarded += dt * guarded_rate;
            theta_guarded = std::fmax(lower, std::fmin(upper, theta_guarded));

            residual_integral += dt * residual;
        }

        std::cout << std::fixed << std::setprecision(6);
        std::cout << "Unconstrained final estimate: " << theta_free << '\n';
        std::cout << "Guarded final estimate:       " << theta_guarded << '\n';
        std::cout << "Residual time average:        "
                  << residual_integral / final_time << '\n';
        std::cout << "CSV written to Chapter9_Lesson1.csv\n";
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

Chapter9_Lesson1.java

/*
Chapter9_Lesson1.java
Parameter-drift diagnostic and a projection-inspired scalar boundary guard.

Build and run:
    javac Chapter9_Lesson1.java
    java Chapter9_Lesson1
*/

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 Chapter9_Lesson1 {
    private Chapter9_Lesson1() {
    }

    private static double boundaryGuard(
            double theta,
            double rawRate,
            double lower,
            double upper) {
        if (theta <= lower && rawRate < 0.0) {
            return 0.0;
        }
        if (theta >= upper && rawRate > 0.0) {
            return 0.0;
        }
        return rawRate;
    }

    private static double residualSignal(double t) {
        return 0.030
                + 0.018 * Math.sin(3.0 * t)
                + 0.012 * Math.sin(17.0 * t)
                + 0.008 / (1.0 + t);
    }

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

        final double dt = 0.002;
        final double finalTime = 120.0;
        final double gamma = 1.5;
        final double lower = -2.0;
        final double upper = 2.0;

        double thetaFree = 0.7;
        double thetaGuarded = 0.7;
        double residualIntegral = 0.0;

        Path output = Path.of("Chapter9_Lesson1.csv");

        try (BufferedWriter writer = Files.newBufferedWriter(
                output, StandardCharsets.UTF_8)) {
            writer.write("time,residual,theta_free,theta_guarded");
            writer.newLine();

            int steps = (int) Math.round(finalTime / dt);
            for (int k = 0; k <= steps; k++) {
                double t = k * dt;
                double residual = residualSignal(t);
                double rawRate = -gamma * residual;

                writer.write(String.format(
                        Locale.US,
                        "%.12f,%.12f,%.12f,%.12f%n",
                        t, residual, thetaFree, thetaGuarded));

                if (k == steps) {
                    break;
                }

                thetaFree += dt * rawRate;

                double guardedRate =
                        boundaryGuard(thetaGuarded, rawRate, lower, upper);
                thetaGuarded += dt * guardedRate;
                thetaGuarded = Math.max(lower, Math.min(upper, thetaGuarded));

                residualIntegral += dt * residual;
            }
        } catch (IOException error) {
            System.err.println("I/O error: " + error.getMessage());
            System.exit(1);
        }

        System.out.printf("Unconstrained final estimate: %.6f%n", thetaFree);
        System.out.printf("Guarded final estimate:       %.6f%n", thetaGuarded);
        System.out.printf(
                "Residual time average:        %.6f%n",
                residualIntegral / finalTime);
        System.out.println("CSV written to " + output);
    }
}

The MATLAB script performs the numerical experiment and writes a CSV file. When Simulink is licensed, it also constructs Chapter9_Lesson1_Simulink.slx. The model uses a free integrator and a limited integrator to visualize unconstrained and boundary-limited adaptation.

Chapter9_Lesson1.m

% Chapter9_Lesson1.m
% Parameter-drift diagnostic and a projection-inspired scalar boundary guard.
% The optional final section creates a comparable Simulink model when
% Simulink is installed.

clear; clc; close all;

dt = 0.002;
finalTime = 120.0;
gamma = 1.5;
theta0 = 0.7;
lower = -2.0;
upper = 2.0;

time = (0:dt:finalTime).';
n = numel(time);
thetaFree = zeros(n, 1);
thetaGuarded = zeros(n, 1);
residual = zeros(n, 1);

thetaFree(1) = theta0;
thetaGuarded(1) = theta0;

for k = 1:n-1
    t = time(k);
    residual(k) = 0.030 ...
        + 0.018 * sin(3.0 * t) ...
        + 0.012 * sin(17.0 * t) ...
        + 0.008 / (1.0 + t);

    rawRate = -gamma * residual(k);

    thetaFree(k + 1) = thetaFree(k) + dt * rawRate;

    guardedRate = boundaryGuard( ...
        thetaGuarded(k), rawRate, lower, upper);
    thetaGuarded(k + 1) = thetaGuarded(k) + dt * guardedRate;
    thetaGuarded(k + 1) = min(upper, max(lower, thetaGuarded(k + 1)));
end

residual(end) = residual(end - 1);

fprintf('Unconstrained final estimate: %.6f\n', thetaFree(end));
fprintf('Guarded final estimate:       %.6f\n', thetaGuarded(end));
fprintf('Residual time average:        %.6f\n', trapz(time, residual) / finalTime);

figure('Name', 'Parameter Drift');
plot(time, thetaFree, 'LineWidth', 1.2); hold on;
plot(time, thetaGuarded, 'LineWidth', 1.2);
yline(lower, '--');
yline(upper, '--');
grid on;
xlabel('time (s)');
ylabel('parameter estimate');
title('Parameter drift caused by a small persistent residual');
legend('unconstrained update', 'boundary-guarded update', ...
       'lower admissible bound', 'upper admissible bound', ...
       'Location', 'best');

resultTable = table(time, residual, thetaFree, thetaGuarded);
writetable(resultTable, 'Chapter9_Lesson1.csv');

% Optional Simulink construction. A Limited Integrator gives the scalar
% boundary behavior used in this introductory lesson.
if license('test', 'Simulink')
    buildProjectionPreviewModel(gamma, theta0, lower, upper, finalTime);
end

function rate = boundaryGuard(theta, rawRate, lower, upper)
    if theta <= lower && rawRate < 0
        rate = 0;
    elseif theta >= upper && rawRate > 0
        rate = 0;
    else
        rate = rawRate;
    end
end

function buildProjectionPreviewModel(gamma, theta0, lower, upper, finalTime)
    model = 'Chapter9_Lesson1_Simulink';
    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/Constant', [model '/Bias'], ...
        'Value', '0.030', 'Position', [40 40 90 70]);
    add_block('simulink/Sources/Sine Wave', [model '/LowFreq'], ...
        'Amplitude', '0.018', 'Frequency', '3.0', ...
        'Position', [40 100 90 130]);
    add_block('simulink/Sources/Sine Wave', [model '/HighFreq'], ...
        'Amplitude', '0.012', 'Frequency', '17.0', ...
        'Position', [40 160 90 190]);
    add_block('simulink/Math Operations/Sum', [model '/ResidualSum'], ...
        'Inputs', '+++', 'Position', [150 85 180 165]);
    add_block('simulink/Math Operations/Gain', [model '/MinusGamma'], ...
        'Gain', num2str(-gamma), 'Position', [230 105 290 145]);

    add_block('simulink/Continuous/Integrator', [model '/FreeIntegrator'], ...
        'InitialCondition', num2str(theta0), ...
        'Position', [350 50 390 80]);
    add_block('simulink/Continuous/Limited Integrator', ...
        [model '/GuardedIntegrator'], ...
        'InitialCondition', num2str(theta0), ...
        'LowerSaturationLimit', num2str(lower), ...
        'UpperSaturationLimit', num2str(upper), ...
        'Position', [350 150 410 190]);

    add_block('simulink/Sinks/To Workspace', [model '/ThetaFree'], ...
        'VariableName', 'thetaFreeSim', 'SaveFormat', 'Structure With Time', ...
        'Position', [480 45 560 85]);
    add_block('simulink/Sinks/To Workspace', [model '/ThetaGuarded'], ...
        'VariableName', 'thetaGuardedSim', ...
        'SaveFormat', 'Structure With Time', ...
        'Position', [480 145 570 195]);

    add_line(model, 'Bias/1', 'ResidualSum/1');
    add_line(model, 'LowFreq/1', 'ResidualSum/2');
    add_line(model, 'HighFreq/1', 'ResidualSum/3');
    add_line(model, 'ResidualSum/1', 'MinusGamma/1');
    add_line(model, 'MinusGamma/1', 'FreeIntegrator/1');
    add_line(model, 'MinusGamma/1', 'GuardedIntegrator/1');
    add_line(model, 'FreeIntegrator/1', 'ThetaFree/1');
    add_line(model, 'GuardedIntegrator/1', 'ThetaGuarded/1');

    set_param(model, 'StopTime', num2str(finalTime), ...
        'Solver', 'ode4', 'FixedStep', '0.002');
    save_system(model);
    fprintf('Created %s.slx\n', model);
end

Chapter9_Lesson1.nb

ClearAll["Global`*"];

dt = 0.002;
finalTime = 120.0;
gamma = 1.5;
theta0 = 0.7;
lower = -2.0;
upper = 2.0;

residual[t_] := 0.030
  + 0.018 Sin[3.0 t]
  + 0.012 Sin[17.0 t]
  + 0.008/(1.0 + t);

boundaryGuard[theta_, rawRate_, lo_, hi_] := Which[
  theta <= lo && rawRate < 0.0, 0.0,
  theta >= hi && rawRate > 0.0, 0.0,
  True, rawRate
];

steps = Round[finalTime/dt];
time = N[Range[0, steps] dt];

thetaFree = ConstantArray[0.0, steps + 1];
thetaGuarded = ConstantArray[0.0, steps + 1];
thetaFree[[1]] = theta0;
thetaGuarded[[1]] = theta0;

Do[
  t = time[[k]];
  rawRate = -gamma residual[t];

  thetaFree[[k + 1]] = thetaFree[[k]] + dt rawRate;

  guardedRate = boundaryGuard[
    thetaGuarded[[k]], rawRate, lower, upper
  ];
  thetaGuarded[[k + 1]] = Clip[
    thetaGuarded[[k]] + dt guardedRate,
    {lower, upper}
  ],
  {k, 1, steps}
];

Print[
  "Unconstrained final estimate: ",
  NumberForm[Last[thetaFree], {10, 6}]
];
Print[
  "Guarded final estimate:       ",
  NumberForm[Last[thetaGuarded], {10, 6}]
];

data = Transpose[{time, residual /@ time, thetaFree, thetaGuarded}];
Export[
  "Chapter9_Lesson1.csv",
  Prepend[data, {"time", "residual", "theta_free", "theta_guarded"}]
];

ListLinePlot[
  {
    Transpose[{time, thetaFree}],
    Transpose[{time, thetaGuarded}]
  },
  PlotLegends -> {"unconstrained update", "boundary-guarded update"},
  PlotRange -> All,
  AxesLabel -> {"time (s)", "parameter estimate"},
  PlotLabel -> "Parameter drift caused by a small persistent residual",
  GridLines -> {None, {lower, upper}},
  ImageSize -> Large
]

Let \(\bar{w}\) denote the long-time average of the residual. The oscillatory terms contribute bounded integrals, whereas the bias contributes \(\bar{w}t\). Consequently,

\[ \hat{\theta}_{\mathrm{free}}(t) = \hat{\theta}(0)-\gamma\bar{w}t+O(1). \]

With \(\gamma=1.5\) and \(\bar{w}\approx0.030\), the asymptotic drift rate is approximately \(-0.045\) parameter units per second. The boundary-guarded estimate follows the nominal law until reaching the lower bound and then rejects only the outward component. Thus the experiment separates two facts:

1. a bounded residual can generate an unbounded unconstrained estimate;
2. a forward-invariant parameter set can prevent that estimate from leaving a credible region.

1. Derive parameter bounds from plant physics or controller-validity conditions rather than from an arbitrary plotting range.
2. Include the ideal parameter vector in the interior of the admissible set whenever the uncertainty description permits it.
3. Leave margin for modeling error; bounds that are too tight create persistent boundary activity and bias.
4. Monitor time spent on the boundary. Frequent activation may indicate incorrect bounds, actuator saturation, inadequate modeling, or insufficient excitation.
5. Treat projection as one robustness layer, not as a substitute for correct sampling, anti-saturation logic, filtering, and validation against neglected dynamics.

Problem 1 (Lyapunov derivative with disturbance): Starting from the first-order MRAC error model in Section 2, derive \(\dot{V}\) for the nominal gradient law and explain why bounded disturbance does not automatically prove bounded parameter estimates.

Solution:

Differentiate the candidate \(V=\frac12e^2+\frac{|b|}{2} \tilde{\boldsymbol{\theta}}^T\boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta}}\). Since the ideal parameters are constant, \(\dot{\tilde{\boldsymbol{\theta}}} =\dot{\hat{\boldsymbol{\theta}}}\). Therefore,

\[ \dot{V} = e\left( -a_m e+ b\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}+d \right) + |b|\tilde{\boldsymbol{\theta}}^{T} \boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta}}}. \]

Substituting \(\dot{\hat{\boldsymbol{\theta}}} =-\boldsymbol{\Gamma}\boldsymbol{\phi}e\operatorname{sgn}(b)\) cancels the parameter cross term, leaving

\[ \dot{V}=-a_m e^2+ed. \]

Bounding \(ed\) by Young's inequality leaves a positive constant term and no negative term proportional to \(\|\tilde{\boldsymbol{\theta}}\|^2\). Hence the result does not exclude growth of the parameter part of \(V\).

Problem 2 (A vanishing drift signal): Let \(\dot{\hat{\theta}}=-\gamma/(1+t)\) with \(\gamma>0\). Show that the update signal tends to zero and is square integrable, but the estimate is unbounded.

Solution:

\[ \lim_{t\to\infty}\frac{1}{1+t}=0, \qquad \int_0^\infty\frac{1}{(1+t)^2}\,dt=1. \]

Integration of the adaptive law gives

\[ \hat{\theta}(t) = \hat{\theta}(0)-\gamma\log(1+t). \]

Since \(\log(1+t)\) is unbounded, the estimate drifts even though its instantaneous update tends to zero.

Problem 3 (Noise-driven random walk): For \(\hat{\theta}_{k+1}=\hat{\theta}_k-h\gamma n_k\), assume independent samples with \(\mathbb{E}[n_k]=0\) and \(\operatorname{Var}(n_k)=\sigma_n^2\). Derive the mean and variance after \(N\) updates.

Solution:

\[ \hat{\theta}_N = \hat{\theta}_0-h\gamma\sum_{k=0}^{N-1}n_k. \]

Linearity of expectation gives \(\mathbb{E}[\hat{\theta}_N]=\hat{\theta}_0\). Because the samples are independent, variances add:

\[ \operatorname{Var}(\hat{\theta}_N) = h^2\gamma^2 \sum_{k=0}^{N-1}\sigma_n^2 = Nh^2\gamma^2\sigma_n^2. \]

The root-mean-square deviation grows proportionally to \(\sqrt{N}\).

Problem 4 (Forward invariance of a scalar interval): Prove that the boundary rule in Section 9 keeps \(\hat{\theta}(t)\) in \([\theta_{\min},\theta_{\max}]\) when the initial estimate is inside the interval.

Solution:

At the lower boundary, every outward direction has \(y<0\) and is replaced by zero; hence the vector field cannot point below the interval. At the upper boundary, every outward direction has \(y>0\) and is replaced by zero; hence the vector field cannot point above the interval. Inward and tangent directions remain unchanged. Therefore no first exit from the closed interval is possible.

Problem 5 (Projection versus leakage): For the scalar biased update \(\dot{\hat{\theta}}=-\gamma b_0\), compare the long-time behavior of: (a) the unconstrained law; (b) leakage \(\dot{\hat{\theta}}=-\gamma b_0-\sigma\hat{\theta}\); and (c) boundary-limited adaptation on \([\theta_{\min},\theta_{\max}]\).

Solution:

The unconstrained estimate is \(\hat{\theta}(0)-\gamma b_0t\) and is unbounded for \(b_0\neq0\). With \(\sigma>0\), leakage gives

\[ \hat{\theta}(t) = e^{-\sigma t}\hat{\theta}(0) - \frac{\gamma b_0}{\sigma} \left(1-e^{-\sigma t}\right), \]

so the estimate approaches \(-\gamma b_0/\sigma\), a soft bound determined by tuning. Boundary-limited adaptation follows the nominal update until it reaches the relevant boundary, after which outward motion is blocked. It enforces the prescribed hard interval but may remain active on the boundary.

Adaptive parameters are integrator states. Small persistent residuals can therefore cause linear, logarithmic, or noise-driven mean-square drift. The nominal Lyapunov cancellation that establishes ideal-model tracking does not automatically damp parameter error under disturbances and unmodeled dynamics. Physical prior knowledge can be encoded as an admissible set, and a projection mechanism can preserve that set while retaining a Lyapunov-compatible update. Lesson 2 formalizes the multidimensional projection operator and proves its principal properties.