After completing this lesson, students should be able to:

• prove that a projection operator preserves the nominal MRAC Lyapunov inequality when the ideal parameter belongs to the admissible set;
• derive a global bound on the magnitude of a normalized gradient update;
• explain why normalization does not automatically preserve an unnormalized direct-MRAC cancellation proof;
• identify the tracking penalty caused by an incorrectly chosen projection set;
• compare tracking error, control effort, parameter excursions, and update-rate peaks using reproducible simulations.

Three distinct questions must not be confused:

1. Closed-loop stability: Are all internal signals bounded, and does the tracking error converge to zero in the nominal model?
2. Parameter safety: Do the estimated gains remain inside physically or computationally acceptable limits?
3. Transient performance: How large are the tracking error, control input, and parameter-update peaks before steady behavior is reached?

Projection primarily answers the second question while preserving the first under a correct set assumption. Normalization primarily modifies the third question, but its effect on the first question depends on the exact adaptive architecture.

Consider the scalar plant and first-order 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 the nominal stability proof begins with \(d(t)=0\). Use the direct adaptive control law

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

Exact matching is obtained when the ideal parameters satisfy

\[ -a+b\theta_x^{*}=-a_m, \qquad b\theta_r^{*}=b_m. \]

Therefore,

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

Define the tracking and parameter errors as \(e=x-x_m\) and \(\tilde{\boldsymbol{\theta} } =\hat{\boldsymbol{\theta} }-\boldsymbol{\theta}^{*}\). Under the nominal matching assumptions,

\[ \dot{e}=-a_m e+b\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\phi}. \]

This compact error model is the basis for all stability comparisons in this lesson.

Let \(\boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}>0\) be the adaptation-gain matrix. The classical update law is

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

Choose the Lyapunov candidate

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

Because the ideal parameter is constant, \(\dot{\tilde{\boldsymbol{\theta} } } =\dot{\hat{\boldsymbol{\theta} } }\). Differentiation gives

\[ \begin{aligned} \dot{V} &=e\left(-a_m e+ b\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\right) +|b|\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta} } }\\ &=-a_m e^{2} +b e\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} -|b|\operatorname{sgn}(b)e \tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\\ &=-a_m e^{2}\le 0. \end{aligned} \]

Thus \(e\) and \(\tilde{\boldsymbol{\theta} }\) are bounded and \(e\in L_2\). With bounded reference signals and bounded error dynamics, the Barbalat argument introduced earlier in the course yields \(e(t)\rightarrow 0\).

Notice what this theorem does not provide: it does not prevent large but finite parameter excursions, and it does not guarantee that the estimated parameters converge to their ideal values. The additional signal-richness conditions for parameter convergence are the subject of Chapter 10.

Let \(\Omega\subset\mathbb{R}^{p}\) be a known compact, convex parameter set and assume \(\boldsymbol{\theta}^{*}\in\Omega\). Define the unconstrained update direction

\[ \mathbf{Y} =-\boldsymbol{\Gamma}\boldsymbol{\phi}e\, \operatorname{sgn}(b) \]

and the projected law

\[ \dot{\hat{\boldsymbol{\theta} } } =\operatorname{Proj}_{\Omega} \left(\hat{\boldsymbol{\theta} },\mathbf{Y}\right). \]

The essential projection property is

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

Add and subtract \(\mathbf{Y}\) in the derivative of \(V\):

\[ \begin{aligned} \dot{V} &=-a_m e^2+ |b|\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \left[ \operatorname{Proj}_{\Omega} \left(\hat{\boldsymbol{\theta} },\mathbf{Y}\right) -\mathbf{Y} \right]\\ &\le -a_m e^2. \end{aligned} \]

Therefore projection does not weaken the nominal Lyapunov conclusion. It adds the invariant-set guarantee

\[ \hat{\boldsymbol{\theta} }(0)\in\Omega \quad\Longrightarrow\quad \hat{\boldsymbol{\theta} }(t)\in\Omega \quad\text{for all }t\ge 0. \]

For componentwise bounds \(\theta_{i,\min}\le\hat{\theta}_i \le\theta_{i,\max}\), a simple tangent-cone projection sets an outward update component to zero at an active boundary:

\[ \operatorname{Proj}_i(\hat{\theta}_i,Y_i)= \begin{cases} 0, &\hat{\theta}_i=\theta_{i,\max}\text{ and }Y_i>0,\\ 0, &\hat{\theta}_i=\theta_{i,\min}\text{ and }Y_i<0,\\ Y_i,&\text{otherwise}. \end{cases} \]

In sampled software, the projected derivative should be combined with a final numerical clamp because a finite integration step can cross a boundary by a small amount.

Projection is inactive in the interior of \(\Omega\). Consequently, if the unconstrained estimate never reaches the boundary, projected and unprojected controllers have identical trajectories. Performance changes only when a boundary becomes active.

5.1 Useful effects

• Parameter estimates cannot enter physically impossible regions.
• If \(\|\hat{\boldsymbol{\theta} }\|\le\theta_{\max}\), then the adaptive control component satisfies \(|u|\le\theta_{\max}\|\boldsymbol{\phi}\|\).
• Numerical overflow and excessively large adaptive gains become less likely.
• The admissible parameter set provides an explicit engineering contract that can be monitored during implementation.

5.2 Possible penalties

If the ideal parameter is excluded from the selected set, exact model matching is impossible. Let \(\bar{\boldsymbol{\theta} }\) be the closest feasible parameter and define the residual mismatch

\[ \Delta_{\Omega}(t)= b\left( \bar{\boldsymbol{\theta} }-\boldsymbol{\theta}^{*} \right)^{T}\boldsymbol{\phi}(t). \]

The Lyapunov derivative contains the additional term \(e\Delta_{\Omega}\). If \(|\Delta_{\Omega}(t)|\le\delta_{\Omega}\), then

\[ \dot{V} \le -a_m e^2+|e|\delta_{\Omega} \le -\frac{a_m}{2}e^2+ \frac{\delta_{\Omega}^{2} }{2a_m}. \]

The conclusion is practical boundedness rather than exact asymptotic tracking. Therefore bounds should be wide enough to contain credible ideal values, but not so wide that they cease to protect the implementation.

A discontinuous box projection can also switch repeatedly near a boundary when noise is present. The smooth projection variants introduced earlier in this chapter reduce this boundary switching at the cost of a thin transition layer.

A common normalized direction uses

\[ m^2(t)=m_0^2+\nu \boldsymbol{\phi}^{T}(t)\boldsymbol{\phi}(t), \qquad m_0>0,\quad \nu>0, \]

\[ \mathbf{Y}_N= -\frac{\boldsymbol{\Gamma}\boldsymbol{\phi}e} {m^2(t)}\operatorname{sgn}(b). \]

Let \(\rho=\|\boldsymbol{\phi}\|\). The scalar function

\[ f(\rho)=\frac{\rho}{m_0^2+\nu\rho^2} \]

has derivative

\[ f'(\rho)= \frac{m_0^2-\nu\rho^2} {(m_0^2+\nu\rho^2)^2}. \]

Its maximum occurs at \(\rho=m_0/\sqrt{\nu}\), giving

\[ \max_{\rho\ge 0}f(\rho) =\frac{1}{2m_0\sqrt{\nu} }. \]

Hence the normalized update satisfies the global bound

\[ \|\mathbf{Y}_N\| \le \frac{\lambda_{\max}(\boldsymbol{\Gamma})} {2m_0\sqrt{\nu} }|e|. \]

This result explains the practical value of normalization: the parameter update cannot grow proportionally without limit as the regressor amplitude increases. Large commands or measured-state spikes reduce the effective adaptation gain

\[ \boldsymbol{\Gamma}_{\mathrm{eff} }(t) =\frac{\boldsymbol{\Gamma} }{m^2(t)}. \]

The trade-off is immediate. A larger \(\nu\) reduces update peaks and noise sensitivity, but also slows adaptation whenever the regressor is large.

It is tempting to replace the classical update by the normalized update and assume that the original proof remains unchanged. For the direct MRAC error model, this assumption is generally false.

Take \(m_0=1\) and use \(\dot{\hat{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}e \operatorname{sgn}(b)/m^2\) in the original Lyapunov function. The derivative becomes

\[ \dot{V} =-a_m e^2+ b e\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} \left(1-\frac{1}{m^2}\right). \]

The remaining cross term has no fixed sign. Normalization may still work well in simulation, but the unnormalized cancellation theorem can no longer be quoted as its proof.

One proof-consistent approach is to introduce normalization in an auxiliary prediction-error channel. Define

\[ \bar{\boldsymbol{\phi} }= \frac{\boldsymbol{\phi} }{m}, \qquad \dot{\varepsilon} =-\lambda\varepsilon+ \tilde{\boldsymbol{\theta} }^{T} \bar{\boldsymbol{\phi} }, \]

and use

\[ \dot{\hat{\boldsymbol{\theta} } } =\operatorname{Proj}_{\Omega} \left( \hat{\boldsymbol{\theta} }, -\boldsymbol{\Gamma} \bar{\boldsymbol{\phi} }\varepsilon \right). \]

With \(V_{\varepsilon}= \varepsilon^2/2+ \tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta} }/2\),

\[ \dot{V}_{\varepsilon} \le -\lambda\varepsilon^2. \]

This is a rigorous normalized-estimator result. A complete adaptive controller must then connect the prediction channel to the tracking channel using a compatible control design. The general design rule is: modify the adaptive law and its proof together.

flowchart TD
  A["Start from a proven adaptive error model"] --> B["Choose credible parameter set Omega"]
  B --> C["Verify ideal parameter is expected inside Omega"]
  C --> D["Apply projection with the required Lyapunov inequality"]
  D --> E["Inspect regressor amplitude and update-rate peaks"]
  E --> F{"Are large-signal update \nspikes unacceptable?"}
  F -->|No| G["Retain projected Lyapunov update"]
  F -->|Yes| H["Introduce a normalized predictor or \nre-derive the tracking proof"]
  H --> I["Tune normalization strength \nand adaptation gain together"]
  G --> J["Validate tracking, control effort, \nbounds, and numerical step"]
  I --> J
        

Projection is usually a theorem-preserving wrapper around a valid update law. Normalization is a gain-scheduling operation on the update direction, so it must be evaluated as part of the dynamic design rather than as a purely numerical afterthought.

The implementations compare three modes using the same plant, model, reference, disturbance, measurement perturbation, initial conditions, and explicit-Euler step:

1. plain: unnormalized gradient MRAC;
2. projection: gradient MRAC with box projection;
3. proj_norm: projected update divided by \(1+0.7\|\boldsymbol{\phi}\|^2\).

\[ a=0.8,\quad b=1.2,\quad a_m=1.5,\quad b_m=1.5, \]

\[ \boldsymbol{\Gamma}=12\mathbf{I}, \qquad -1.5\le\hat{\theta}_x\le 0.2, \qquad 0.2\le\hat{\theta}_r\le 2.0. \]

The corresponding ideal values are

\[ \theta_x^{*}= \frac{0.8-1.5}{1.2} =-0.5833, \qquad \theta_r^{*}= \frac{1.5}{1.2} =1.25, \]

so both ideal values lie inside the selected projection box.

flowchart TD
  R["Reference r"] --> RM["Reference model"]
  R --> PHI["Regressor: measured state and reference"]
  P["Plant"] --> Y["Measured state"]
  Y --> PHI
  RM --> E["Measured tracking error"]
  Y --> E
  E --> AD["Adaptive update"]
  PHI --> AD
  AD --> PN["Optional normalization and projection"]
  PN --> TH["Parameter integrators"]
  TH --> U["Adaptive control"]
  PHI --> U
  U --> P
        

The reference contains two short high-amplitude intervals. These intervals make the effect of normalization visible without introducing concepts beyond the current chapter.

The programs calculate

\[ \operatorname{IAE}=\int_0^T |e(t)|\,dt, \qquad e_{\mathrm{RMS} }= \sqrt{\frac{1}{T}\int_0^T e^2(t)\,dt}, \]

\[ u_{\mathrm{RMS} }= \sqrt{\frac{1}{T}\int_0^T u^2(t)\,dt}, \qquad \max_t\|\dot{\hat{\boldsymbol{\theta} } }(t)\|. \]

With \(T=20\) s and \(\Delta t=10^{-3}\) s, the supplied implementations produce approximately:

The plain controller is not unstable in this experiment; its bound-violation count only records time samples outside the chosen engineering box. Projection removes those excursions with almost no tracking penalty because the ideal parameters are feasible. Normalization reduces the peak update norm by about one half, but the slower effective adaptation increases the transient tracking measures. This is the central performance trade-off.

Useful libraries include NumPy for vector operations, SciPy for higher-order ODE solvers, python-control for linear reference-model analysis, and Matplotlib for visualization. The supplied program implements the adaptive loop from scratch and uses Matplotlib only when available.

Chapter9_Lesson5.py

"""
Chapter9_Lesson5.py
Impact of projection and normalization on a first-order MRAC simulation.

Modes:
  plain      - classical gradient update
  projection - componentwise box projection
  proj_norm  - box projection plus regressor normalization

The normalized direct-MRAC update is included as a performance experiment.
As explained in the lesson, it does not inherit the unnormalized Lyapunov
cancellation automatically; a proof-consistent normalized estimator must be
derived as part of the complete adaptive architecture.
"""

from __future__ import annotations

import csv
import math
from dataclasses import dataclass
from pathlib import Path
from typing import Dict, List, Tuple


@dataclass
class Metrics:
    mode: str
    iae: float
    rmse: float
    max_error: float
    control_rms: float
    max_theta_dot: float
    theta_x_final: float
    theta_r_final: float
    bound_violations: int


def reference(t: float) -> float:
    """Bounded command with two short high-amplitude intervals."""
    r = math.sin(0.7 * t) + 0.5 * math.sin(1.9 * t)
    if 6.0 <= t < 6.5:
        r += 4.0
    if 12.0 <= t < 12.5:
        r -= 3.0
    return r


def projected_rate(
    theta: List[float],
    rate: List[float],
    lower: List[float],
    upper: List[float],
) -> List[float]:
    """Project an update direction onto the tangent cone of a box."""
    result = rate.copy()
    for i in range(len(theta)):
        outward_at_lower = theta[i] <= lower[i] and result[i] < 0.0
        outward_at_upper = theta[i] >= upper[i] and result[i] > 0.0
        if outward_at_lower or outward_at_upper:
            result[i] = 0.0
    return result


def simulate(
    mode: str,
    dt: float = 1.0e-3,
    final_time: float = 20.0,
) -> Tuple[Metrics, List[Dict[str, float]]]:
    if mode not in {"plain", "projection", "proj_norm"}:
        raise ValueError(f"Unsupported mode: {mode}")

    # Plant: x_dot = -a*x + b*u + d(t)
    # Model: xm_dot = -am*xm + bm*r
    a, b = 0.8, 1.2
    am, bm = 1.5, 1.5

    gamma = [12.0, 12.0]
    lower = [-1.5, 0.2]
    upper = [0.2, 2.0]
    nu = 0.7

    x = 0.0
    xm = 0.0
    theta = [0.0, 0.5]

    iae = 0.0
    ise = 0.0
    control_energy = 0.0
    max_error = 0.0
    max_theta_dot = 0.0
    bound_violations = 0
    samples: List[Dict[str, float]] = []

    steps = int(round(final_time / dt))
    for k in range(steps + 1):
        t = k * dt
        r = reference(t)

        # Deterministic high-frequency measurement perturbation.
        y = x + 0.01 * math.sin(47.0 * t)
        e = y - xm
        phi = [y, r]
        u = theta[0] * phi[0] + theta[1] * phi[1]

        raw_rate = [
            -gamma[0] * phi[0] * e,
            -gamma[1] * phi[1] * e,
        ]

        if mode == "proj_norm":
            m2 = 1.0 + nu * (phi[0] ** 2 + phi[1] ** 2)
            theta_dot = [raw_rate[0] / m2, raw_rate[1] / m2]
        else:
            theta_dot = raw_rate

        if mode in {"projection", "proj_norm"}:
            theta_dot = projected_rate(theta, theta_dot, lower, upper)

        theta_dot_norm = math.hypot(theta_dot[0], theta_dot[1])
        max_theta_dot = max(max_theta_dot, theta_dot_norm)

        disturbance = 0.03 * math.sin(3.1 * t)
        x_dot = -a * x + b * u + disturbance
        xm_dot = -am * xm + bm * r

        # Explicit Euler is used deliberately so every language reproduces
        # the same experiment. Decrease dt for a convergence check.
        x += dt * x_dot
        xm += dt * xm_dot
        theta[0] += dt * theta_dot[0]
        theta[1] += dt * theta_dot[1]

        if mode in {"projection", "proj_norm"}:
            theta[0] = min(max(theta[0], lower[0]), upper[0])
            theta[1] = min(max(theta[1], lower[1]), upper[1])

        outside = (
            theta[0] < lower[0] - 1.0e-12
            or theta[0] > upper[0] + 1.0e-12
            or theta[1] < lower[1] - 1.0e-12
            or theta[1] > upper[1] + 1.0e-12
        )
        if outside:
            bound_violations += 1

        iae += abs(e) * dt
        ise += e * e * dt
        control_energy += u * u * dt
        max_error = max(max_error, abs(e))

        if k % 20 == 0:
            samples.append(
                {
                    "time": t,
                    "reference": r,
                    "plant": x,
                    "model": xm,
                    "error": e,
                    "control": u,
                    "theta_x": theta[0],
                    "theta_r": theta[1],
                    "theta_dot_norm": theta_dot_norm,
                }
            )

    metrics = Metrics(
        mode=mode,
        iae=iae,
        rmse=math.sqrt(ise / final_time),
        max_error=max_error,
        control_rms=math.sqrt(control_energy / final_time),
        max_theta_dot=max_theta_dot,
        theta_x_final=theta[0],
        theta_r_final=theta[1],
        bound_violations=bound_violations,
    )
    return metrics, samples


def write_csv(all_samples: Dict[str, List[Dict[str, float]]]) -> Path:
    output = Path("Chapter9_Lesson5_results.csv")
    fieldnames = [
        "mode",
        "time",
        "reference",
        "plant",
        "model",
        "error",
        "control",
        "theta_x",
        "theta_r",
        "theta_dot_norm",
    ]
    with output.open("w", newline="", encoding="utf-8") as file:
        writer = csv.DictWriter(file, fieldnames=fieldnames)
        writer.writeheader()
        for mode, rows in all_samples.items():
            for row in rows:
                writer.writerow({"mode": mode, **row})
    return output


def optional_plot(all_samples: Dict[str, List[Dict[str, float]]]) -> None:
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        print("matplotlib is not installed; skipping plots.")
        return

    for mode, rows in all_samples.items():
        time = [row["time"] for row in rows]
        error = [row["error"] for row in rows]
        plt.plot(time, error, label=mode)
    plt.xlabel("Time (s)")
    plt.ylabel("Tracking error")
    plt.title("Projection and normalization comparison")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter9_Lesson5_error.png", dpi=180)
    plt.show()


def main() -> None:
    all_samples: Dict[str, List[Dict[str, float]]] = {}
    for mode in ("plain", "projection", "proj_norm"):
        metrics, samples = simulate(mode)
        all_samples[mode] = samples
        print(
            f"{metrics.mode:10s} "
            f"IAE={metrics.iae:.6f} "
            f"RMSE={metrics.rmse:.6f} "
            f"max|e|={metrics.max_error:.6f} "
            f"u_RMS={metrics.control_rms:.6f} "
            f"max||theta_dot||={metrics.max_theta_dot:.6f} "
            f"theta_f=({metrics.theta_x_final:.6f}, "
            f"{metrics.theta_r_final:.6f}) "
            f"bound_violations={metrics.bound_violations}"
        )

    csv_path = write_csv(all_samples)
    print(f"Wrote {csv_path.resolve()}")
    optional_plot(all_samples)


if __name__ == "__main__":
    main()

For larger adaptive systems, Eigen is appropriate for matrix calculations and Boost.Odeint provides higher-order numerical integration. The implementation below uses only the C++17 standard library so that the projection and normalization logic remains explicit.

Chapter9_Lesson5.cpp

// Chapter9_Lesson5.cpp
// Projection and normalization effects in a first-order MRAC simulation.
//
// Build:
//   g++ -std=c++17 -O2 Chapter9_Lesson5.cpp -o Chapter9_Lesson5
//
// The normalized direct-MRAC update is a performance experiment. It does not
// automatically inherit the classical unnormalized Lyapunov cancellation.

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

struct Sample {
    double time{};
    double reference{};
    double plant{};
    double model{};
    double error{};
    double control{};
    double theta_x{};
    double theta_r{};
    double theta_dot_norm{};
};

struct Metrics {
    std::string mode;
    double iae{};
    double rmse{};
    double max_error{};
    double control_rms{};
    double max_theta_dot{};
    double theta_x_final{};
    double theta_r_final{};
    int bound_violations{};
};

double reference_signal(double t) {
    double r = std::sin(0.7 * t) + 0.5 * std::sin(1.9 * t);
    if (t >= 6.0 && t < 6.5) {
        r += 4.0;
    }
    if (t >= 12.0 && t < 12.5) {
        r -= 3.0;
    }
    return r;
}

std::array<double, 2> projected_rate(
    const std::array<double, 2>& theta,
    std::array<double, 2> rate,
    const std::array<double, 2>& lower,
    const std::array<double, 2>& upper
) {
    for (std::size_t i = 0; i < theta.size(); ++i) {
        const bool outward_at_lower = theta[i] <= lower[i] && rate[i] < 0.0;
        const bool outward_at_upper = theta[i] >= upper[i] && rate[i] > 0.0;
        if (outward_at_lower || outward_at_upper) {
            rate[i] = 0.0;
        }
    }
    return rate;
}

std::pair<Metrics, std::vector<Sample>> simulate(
    const std::string& mode,
    double dt = 1.0e-3,
    double final_time = 20.0
) {
    if (mode != "plain" && mode != "projection" && mode != "proj_norm") {
        throw std::invalid_argument("Unsupported mode: " + mode);
    }

    const double a = 0.8;
    const double b = 1.2;
    const double am = 1.5;
    const double bm = 1.5;
    const std::array<double, 2> gamma{12.0, 12.0};
    const std::array<double, 2> lower{-1.5, 0.2};
    const std::array<double, 2> upper{0.2, 2.0};
    const double nu = 0.7;

    double x = 0.0;
    double xm = 0.0;
    std::array<double, 2> theta{0.0, 0.5};

    double iae = 0.0;
    double ise = 0.0;
    double control_energy = 0.0;
    double max_error = 0.0;
    double max_theta_dot = 0.0;
    int bound_violations = 0;
    std::vector<Sample> samples;

    const int steps = static_cast<int>(std::llround(final_time / dt));
    for (int k = 0; k <= steps; ++k) {
        const double t = k * dt;
        const double r = reference_signal(t);
        const double y = x + 0.01 * std::sin(47.0 * t);
        const double e = y - xm;
        const std::array<double, 2> phi{y, r};
        const double u = theta[0] * phi[0] + theta[1] * phi[1];

        std::array<double, 2> theta_dot{
            -gamma[0] * phi[0] * e,
            -gamma[1] * phi[1] * e
        };

        if (mode == "proj_norm") {
            const double m2 = 1.0 + nu * (phi[0] * phi[0] + phi[1] * phi[1]);
            theta_dot[0] /= m2;
            theta_dot[1] /= m2;
        }

        if (mode == "projection" || mode == "proj_norm") {
            theta_dot = projected_rate(theta, theta_dot, lower, upper);
        }

        const double theta_dot_norm = std::hypot(theta_dot[0], theta_dot[1]);
        max_theta_dot = std::max(max_theta_dot, theta_dot_norm);

        const double disturbance = 0.03 * std::sin(3.1 * t);
        const double x_dot = -a * x + b * u + disturbance;
        const double xm_dot = -am * xm + bm * r;

        x += dt * x_dot;
        xm += dt * xm_dot;
        theta[0] += dt * theta_dot[0];
        theta[1] += dt * theta_dot[1];

        if (mode == "projection" || mode == "proj_norm") {
            theta[0] = std::clamp(theta[0], lower[0], upper[0]);
            theta[1] = std::clamp(theta[1], lower[1], upper[1]);
        }

        const bool outside =
            theta[0] < lower[0] - 1.0e-12 ||
            theta[0] > upper[0] + 1.0e-12 ||
            theta[1] < lower[1] - 1.0e-12 ||
            theta[1] > upper[1] + 1.0e-12;
        if (outside) {
            ++bound_violations;
        }

        iae += std::abs(e) * dt;
        ise += e * e * dt;
        control_energy += u * u * dt;
        max_error = std::max(max_error, std::abs(e));

        if (k % 20 == 0) {
            samples.push_back(
                Sample{t, r, x, xm, e, u, theta[0], theta[1], theta_dot_norm}
            );
        }
    }

    Metrics metrics{
        mode,
        iae,
        std::sqrt(ise / final_time),
        max_error,
        std::sqrt(control_energy / final_time),
        max_theta_dot,
        theta[0],
        theta[1],
        bound_violations
    };
    return {metrics, samples};
}

int main() {
    try {
        std::ofstream csv("Chapter9_Lesson5_results.csv");
        if (!csv) {
            throw std::runtime_error("Cannot open output CSV file.");
        }
        csv << "mode,time,reference,plant,model,error,control,"
               "theta_x,theta_r,theta_dot_norm\n";
        csv << std::setprecision(12);

        for (const std::string mode : {"plain", "projection", "proj_norm"}) {
            const auto [metrics, samples] = simulate(mode);
            for (const Sample& row : samples) {
                csv << mode << ','
                    << row.time << ','
                    << row.reference << ','
                    << row.plant << ','
                    << row.model << ','
                    << row.error << ','
                    << row.control << ','
                    << row.theta_x << ','
                    << row.theta_r << ','
                    << row.theta_dot_norm << '\n';
            }

            std::cout << std::fixed << std::setprecision(6)
                      << std::setw(10) << metrics.mode
                      << " IAE=" << metrics.iae
                      << " RMSE=" << metrics.rmse
                      << " max|e|=" << metrics.max_error
                      << " u_RMS=" << metrics.control_rms
                      << " max||theta_dot||=" << metrics.max_theta_dot
                      << " theta_f=(" << metrics.theta_x_final
                      << ", " << metrics.theta_r_final << ')'
                      << " bound_violations=" << metrics.bound_violations
                      << '\n';
        }

        std::cout << "Wrote Chapter9_Lesson5_results.csv\n";
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

Apache Commons Math supplies matrix classes and ODE integrators when a production implementation needs them. This Java version uses the JDK alone, writes a CSV file, and reproduces the same numerical experiment.

Chapter9_Lesson5.java

// Chapter9_Lesson5.java
// Projection and normalization effects in a first-order MRAC simulation.
//
// Build and run:
//   javac Chapter9_Lesson5.java
//   java Chapter9_Lesson5
//
// The normalized direct-MRAC update is a performance experiment. It does not
// automatically inherit the classical unnormalized Lyapunov cancellation.

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.ArrayList;
import java.util.List;
import java.util.Locale;

public final class Chapter9_Lesson5 {
    private Chapter9_Lesson5() {}

    private static final class Sample {
        final double time;
        final double reference;
        final double plant;
        final double model;
        final double error;
        final double control;
        final double thetaX;
        final double thetaR;
        final double thetaDotNorm;

        Sample(
            double time,
            double reference,
            double plant,
            double model,
            double error,
            double control,
            double thetaX,
            double thetaR,
            double thetaDotNorm
        ) {
            this.time = time;
            this.reference = reference;
            this.plant = plant;
            this.model = model;
            this.error = error;
            this.control = control;
            this.thetaX = thetaX;
            this.thetaR = thetaR;
            this.thetaDotNorm = thetaDotNorm;
        }
    }

    private static final class Metrics {
        final String mode;
        final double iae;
        final double rmse;
        final double maxError;
        final double controlRms;
        final double maxThetaDot;
        final double thetaXFinal;
        final double thetaRFinal;
        final int boundViolations;

        Metrics(
            String mode,
            double iae,
            double rmse,
            double maxError,
            double controlRms,
            double maxThetaDot,
            double thetaXFinal,
            double thetaRFinal,
            int boundViolations
        ) {
            this.mode = mode;
            this.iae = iae;
            this.rmse = rmse;
            this.maxError = maxError;
            this.controlRms = controlRms;
            this.maxThetaDot = maxThetaDot;
            this.thetaXFinal = thetaXFinal;
            this.thetaRFinal = thetaRFinal;
            this.boundViolations = boundViolations;
        }
    }

    private static final class Result {
        final Metrics metrics;
        final List<Sample> samples;

        Result(Metrics metrics, List<Sample> samples) {
            this.metrics = metrics;
            this.samples = samples;
        }
    }

    private static double referenceSignal(double t) {
        double r = Math.sin(0.7 * t) + 0.5 * Math.sin(1.9 * t);
        if (t >= 6.0 && t < 6.5) {
            r += 4.0;
        }
        if (t >= 12.0 && t < 12.5) {
            r -= 3.0;
        }
        return r;
    }

    private static double[] projectedRate(
        double[] theta,
        double[] rate,
        double[] lower,
        double[] upper
    ) {
        double[] result = rate.clone();
        for (int i = 0; i < theta.length; i++) {
            boolean outwardAtLower = theta[i] <= lower[i] && result[i] < 0.0;
            boolean outwardAtUpper = theta[i] >= upper[i] && result[i] > 0.0;
            if (outwardAtLower || outwardAtUpper) {
                result[i] = 0.0;
            }
        }
        return result;
    }

    private static Result simulate(
        String mode,
        double dt,
        double finalTime
    ) {
        if (
            !mode.equals("plain")
            && !mode.equals("projection")
            && !mode.equals("proj_norm")
        ) {
            throw new IllegalArgumentException("Unsupported mode: " + mode);
        }

        double a = 0.8;
        double b = 1.2;
        double am = 1.5;
        double bm = 1.5;
        double[] gamma = {12.0, 12.0};
        double[] lower = {-1.5, 0.2};
        double[] upper = {0.2, 2.0};
        double nu = 0.7;

        double x = 0.0;
        double xm = 0.0;
        double[] theta = {0.0, 0.5};

        double iae = 0.0;
        double ise = 0.0;
        double controlEnergy = 0.0;
        double maxError = 0.0;
        double maxThetaDot = 0.0;
        int boundViolations = 0;
        List<Sample> samples = new ArrayList<>();

        int steps = (int) Math.round(finalTime / dt);
        for (int k = 0; k <= steps; k++) {
            double t = k * dt;
            double r = referenceSignal(t);
            double y = x + 0.01 * Math.sin(47.0 * t);
            double e = y - xm;
            double[] phi = {y, r};
            double u = theta[0] * phi[0] + theta[1] * phi[1];

            double[] thetaDot = {
                -gamma[0] * phi[0] * e,
                -gamma[1] * phi[1] * e
            };

            if (mode.equals("proj_norm")) {
                double m2 = 1.0 + nu * (
                    phi[0] * phi[0] + phi[1] * phi[1]
                );
                thetaDot[0] /= m2;
                thetaDot[1] /= m2;
            }

            if (mode.equals("projection") || mode.equals("proj_norm")) {
                thetaDot = projectedRate(theta, thetaDot, lower, upper);
            }

            double thetaDotNorm = Math.hypot(thetaDot[0], thetaDot[1]);
            maxThetaDot = Math.max(maxThetaDot, thetaDotNorm);

            double disturbance = 0.03 * Math.sin(3.1 * t);
            double xDot = -a * x + b * u + disturbance;
            double xmDot = -am * xm + bm * r;

            x += dt * xDot;
            xm += dt * xmDot;
            theta[0] += dt * thetaDot[0];
            theta[1] += dt * thetaDot[1];

            if (mode.equals("projection") || mode.equals("proj_norm")) {
                theta[0] = Math.min(Math.max(theta[0], lower[0]), upper[0]);
                theta[1] = Math.min(Math.max(theta[1], lower[1]), upper[1]);
            }

            boolean outside =
                theta[0] < lower[0] - 1.0e-12
                || theta[0] > upper[0] + 1.0e-12
                || theta[1] < lower[1] - 1.0e-12
                || theta[1] > upper[1] + 1.0e-12;
            if (outside) {
                boundViolations++;
            }

            iae += Math.abs(e) * dt;
            ise += e * e * dt;
            controlEnergy += u * u * dt;
            maxError = Math.max(maxError, Math.abs(e));

            if (k % 20 == 0) {
                samples.add(
                    new Sample(
                        t,
                        r,
                        x,
                        xm,
                        e,
                        u,
                        theta[0],
                        theta[1],
                        thetaDotNorm
                    )
                );
            }
        }

        Metrics metrics = new Metrics(
            mode,
            iae,
            Math.sqrt(ise / finalTime),
            maxError,
            Math.sqrt(controlEnergy / finalTime),
            maxThetaDot,
            theta[0],
            theta[1],
            boundViolations
        );
        return new Result(metrics, samples);
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        String[] modes = {"plain", "projection", "proj_norm"};
        Path output = Path.of("Chapter9_Lesson5_results.csv");

        try (
            BufferedWriter writer = Files.newBufferedWriter(
                output,
                StandardCharsets.UTF_8
            )
        ) {
            writer.write(
                "mode,time,reference,plant,model,error,control,"
                + "theta_x,theta_r,theta_dot_norm\n"
            );

            for (String mode : modes) {
                Result result = simulate(mode, 1.0e-3, 20.0);
                for (Sample row : result.samples) {
                    writer.write(
                        String.format(
                            Locale.US,
                            "%s,%.12g,%.12g,%.12g,%.12g,%.12g,"
                            + "%.12g,%.12g,%.12g,%.12g%n",
                            mode,
                            row.time,
                            row.reference,
                            row.plant,
                            row.model,
                            row.error,
                            row.control,
                            row.thetaX,
                            row.thetaR,
                            row.thetaDotNorm
                        )
                    );
                }

                Metrics m = result.metrics;
                System.out.printf(
                    Locale.US,
                    "%-10s IAE=%.6f RMSE=%.6f max|e|=%.6f "
                    + "u_RMS=%.6f max||theta_dot||=%.6f "
                    + "theta_f=(%.6f, %.6f) bound_violations=%d%n",
                    m.mode,
                    m.iae,
                    m.rmse,
                    m.maxError,
                    m.controlRms,
                    m.maxThetaDot,
                    m.thetaXFinal,
                    m.thetaRFinal,
                    m.boundViolations
                );
            }
            System.out.println("Wrote " + output.toAbsolutePath());
        } catch (IOException error) {
            System.err.println("I/O error: " + error.getMessage());
            System.exit(1);
        }
    }
}

MATLAB provides convenient matrix notation, while Control System Toolbox supports reference-model analysis. In Simulink, use Integrator blocks for the plant, model, and parameter states; implement projection and normalization in a MATLAB Function block; and log the tracking and parameter signals. The script below performs the complete numerical experiment and documents the corresponding Simulink mapping.

Chapter9_Lesson5.m

% Chapter9_Lesson5.m
% Projection and normalization effects in a first-order MRAC simulation.
%
% Modes:
%   plain      - classical gradient update
%   projection - componentwise box projection
%   proj_norm  - box projection plus regressor normalization
%
% The normalized direct-MRAC update is a performance experiment. It does not
% automatically inherit the classical unnormalized Lyapunov cancellation.

clear; clc; close all;

modes = ["plain", "projection", "proj_norm"];
allRows = table();

for mode = modes
    [metrics, rows] = simulateMode(mode, 1.0e-3, 20.0);
    allRows = [allRows; rows]; %#ok<AGROW>

    fprintf(['%-10s IAE=%.6f RMSE=%.6f max|e|=%.6f ', ...
             'u_RMS=%.6f max||theta_dot||=%.6f ', ...
             'theta_f=(%.6f, %.6f) bound_violations=%d\n'], ...
        metrics.mode, metrics.iae, metrics.rmse, metrics.maxError, ...
        metrics.controlRms, metrics.maxThetaDot, ...
        metrics.thetaXFinal, metrics.thetaRFinal, ...
        metrics.boundViolations);
end

writetable(allRows, 'Chapter9_Lesson5_results.csv');
disp('Wrote Chapter9_Lesson5_results.csv');

figure('Name', 'Chapter 9 Lesson 5');
hold on;
for mode = modes
    selector = allRows.mode == mode;
    plot(allRows.time(selector), allRows.error(selector), ...
        'DisplayName', mode);
end
xlabel('Time (s)');
ylabel('Tracking error');
title('Projection and normalization comparison');
grid on;
legend('Location', 'best');
hold off;

% Simulink mapping:
% 1. Use Integrator blocks for x and x_m.
% 2. Use Sum/Gain blocks for x_dot=-a*x+b*u+d and
%    x_m_dot=-a_m*x_m+b_m*r.
% 3. Use a MATLAB Function block for projectedNormalizedRate().
% 4. Feed theta_dot into two Integrator blocks and use theta in u=theta''*phi.
% 5. Log e, u, theta, and theta_dot_norm with To Workspace blocks.
%
% Control System Toolbox is useful for reference-model analysis. Simulink
% implements the block architecture. This script itself uses only base MATLAB.

function [metrics, rows] = simulateMode(mode, dt, finalTime)
    arguments
        mode (1,1) string
        dt (1,1) double {mustBePositive}
        finalTime (1,1) double {mustBePositive}
    end

    validModes = ["plain", "projection", "proj_norm"];
    if ~any(mode == validModes)
        error('Unsupported mode: %s', mode);
    end

    a = 0.8;
    b = 1.2;
    am = 1.5;
    bm = 1.5;
    gamma = [12.0; 12.0];
    lower = [-1.5; 0.2];
    upper = [0.2; 2.0];
    nu = 0.7;

    x = 0.0;
    xm = 0.0;
    theta = [0.0; 0.5];

    iae = 0.0;
    ise = 0.0;
    controlEnergy = 0.0;
    maxError = 0.0;
    maxThetaDot = 0.0;
    boundViolations = 0;

    steps = round(finalTime / dt);
    sampleCount = floor(steps / 20) + 1;
    data = zeros(sampleCount, 9);
    rowIndex = 1;

    for k = 0:steps
        t = k * dt;
        r = referenceSignal(t);
        y = x + 0.01 * sin(47.0 * t);
        e = y - xm;
        phi = [y; r];
        u = theta.' * phi;

        rawRate = -gamma .* phi .* e;
        if mode == "proj_norm"
            m2 = 1.0 + nu * (phi.' * phi);
            thetaDot = rawRate / m2;
        else
            thetaDot = rawRate;
        end

        if mode == "projection" || mode == "proj_norm"
            thetaDot = projectedRate(theta, thetaDot, lower, upper);
        end

        thetaDotNorm = norm(thetaDot);
        maxThetaDot = max(maxThetaDot, thetaDotNorm);

        disturbance = 0.03 * sin(3.1 * t);
        xDot = -a * x + b * u + disturbance;
        xmDot = -am * xm + bm * r;

        x = x + dt * xDot;
        xm = xm + dt * xmDot;
        theta = theta + dt * thetaDot;

        if mode == "projection" || mode == "proj_norm"
            theta = min(max(theta, lower), upper);
        end

        outside = any(theta < lower - 1.0e-12) || ...
                  any(theta > upper + 1.0e-12);
        if outside
            boundViolations = boundViolations + 1;
        end

        iae = iae + abs(e) * dt;
        ise = ise + e^2 * dt;
        controlEnergy = controlEnergy + u^2 * dt;
        maxError = max(maxError, abs(e));

        if mod(k, 20) == 0
            data(rowIndex, :) = [t, r, x, xm, e, u, ...
                                 theta(1), theta(2), thetaDotNorm];
            rowIndex = rowIndex + 1;
        end
    end

    n = size(data, 1);
    rows = table( ...
        repmat(mode, n, 1), ...
        data(:, 1), data(:, 2), data(:, 3), data(:, 4), ...
        data(:, 5), data(:, 6), data(:, 7), data(:, 8), data(:, 9), ...
        'VariableNames', {'mode', 'time', 'reference', 'plant', 'model', ...
                          'error', 'control', 'theta_x', 'theta_r', ...
                          'theta_dot_norm'} ...
    );

    metrics = struct( ...
        'mode', mode, ...
        'iae', iae, ...
        'rmse', sqrt(ise / finalTime), ...
        'maxError', maxError, ...
        'controlRms', sqrt(controlEnergy / finalTime), ...
        'maxThetaDot', maxThetaDot, ...
        'thetaXFinal', theta(1), ...
        'thetaRFinal', theta(2), ...
        'boundViolations', boundViolations ...
    );
end

function r = referenceSignal(t)
    r = sin(0.7 * t) + 0.5 * sin(1.9 * t);
    if t >= 6.0 && t < 6.5
        r = r + 4.0;
    end
    if t >= 12.0 && t < 12.5
        r = r - 3.0;
    end
end

function result = projectedRate(theta, rate, lower, upper)
    result = rate;
    for i = 1:numel(theta)
        outwardAtLower = theta(i) <= lower(i) && result(i) < 0.0;
        outwardAtUpper = theta(i) >= upper(i) && result(i) > 0.0;
        if outwardAtLower || outwardAtUpper
            result(i) = 0.0;
        end
    end
end

Wolfram Language functions useful for more advanced versions include NDSolve for ODE integration, WhenEvent for boundary events, Clip for box constraints, and LyapunovSolve for linear Lyapunov equations. The downloadable notebook wraps the following source in an input cell.

Chapter9_Lesson5.nb

ClearAll["Global`*"];

referenceSignal[t_] := Module[{r},
  r = Sin[0.7 t] + 0.5 Sin[1.9 t];
  If[6.0 <= t < 6.5, r += 4.0];
  If[12.0 <= t < 12.5, r -= 3.0];
  r
];

projectedRate[theta_, rate_, lower_, upper_] := MapThread[
  Function[{th, dr, lo, hi},
    If[(th <= lo && dr < 0.0) || (th >= hi && dr > 0.0), 0.0, dr]
  ],
  {theta, rate, lower, upper}
];

simulate[mode_String, dt_: 0.001, finalTime_: 20.0] := Module[
  {
    a = 0.8, b = 1.2, am = 1.5, bm = 1.5,
    gamma = {12.0, 12.0}, lower = {-1.5, 0.2},
    upper = {0.2, 2.0}, nu = 0.7,
    x = 0.0, xm = 0.0, theta = {0.0, 0.5},
    iae = 0.0, ise = 0.0, controlEnergy = 0.0,
    maxError = 0.0, maxThetaDot = 0.0,
    boundViolations = 0, samples = {},
    steps, t, r, y, e, phi, u, rawRate, thetaDot,
    m2, disturbance, xDot, xmDot, outside
  },

  If[! MemberQ[{"plain", "projection", "proj_norm"}, mode],
    Return[$Failed]
  ];

  steps = Round[finalTime/dt];

  Do[
    t = k dt;
    r = referenceSignal[t];
    y = x + 0.01 Sin[47.0 t];
    e = y - xm;
    phi = {y, r};
    u = theta . phi;

    rawRate = -gamma phi e;
    thetaDot = If[
      mode === "proj_norm",
      m2 = 1.0 + nu (phi . phi);
      rawRate/m2,
      rawRate
    ];

    If[MemberQ[{"projection", "proj_norm"}, mode],
      thetaDot = projectedRate[theta, thetaDot, lower, upper]
    ];

    maxThetaDot = Max[maxThetaDot, Norm[thetaDot]];

    disturbance = 0.03 Sin[3.1 t];
    xDot = -a x + b u + disturbance;
    xmDot = -am xm + bm r;

    x += dt xDot;
    xm += dt xmDot;
    theta += dt thetaDot;

    If[MemberQ[{"projection", "proj_norm"}, mode],
      theta = Clip[theta, Transpose[{lower, upper}]]
    ];

    outside = Or[
      theta[[1]] < lower[[1]] - 10^-12,
      theta[[1]] > upper[[1]] + 10^-12,
      theta[[2]] < lower[[2]] - 10^-12,
      theta[[2]] > upper[[2]] + 10^-12
    ];
    If[outside, boundViolations++];

    iae += Abs[e] dt;
    ise += e^2 dt;
    controlEnergy += u^2 dt;
    maxError = Max[maxError, Abs[e]];

    If[Mod[k, 20] == 0,
      AppendTo[
        samples,
        <|
          "mode" -> mode, "time" -> t, "reference" -> r,
          "plant" -> x, "model" -> xm, "error" -> e,
          "control" -> u, "theta_x" -> theta[[1]],
          "theta_r" -> theta[[2]],
          "theta_dot_norm" -> Norm[thetaDot]
        |>
      ]
    ],
    {k, 0, steps}
  ];

  <|
    "Metrics" -> <|
      "mode" -> mode,
      "IAE" -> iae,
      "RMSE" -> Sqrt[ise/finalTime],
      "MaxError" -> maxError,
      "ControlRMS" -> Sqrt[controlEnergy/finalTime],
      "MaxThetaDot" -> maxThetaDot,
      "ThetaFinal" -> theta,
      "BoundViolations" -> boundViolations
    |>,
    "Samples" -> samples
  |>
];

modes = {"plain", "projection", "proj_norm"};
results = AssociationMap[simulate, modes];

Dataset[AssociationMap[results[#]["Metrics"] &, modes]]

allSamples = Flatten[results[#]["Samples"] & /@ modes];
Export["Chapter9_Lesson5_results.csv", allSamples];

ListLinePlot[
  Table[
    ({ #time, #error } &) /@ results[mode]["Samples"],
    {mode, modes}
  ],
  PlotLegends -> modes,
  Frame -> True,
  FrameLabel -> {"Time (s)", "Tracking error"},
  PlotLabel -> "Projection and normalization comparison",
  ImageSize -> Large
]

(* Useful Wolfram Language facilities:
   NDSolve for higher-accuracy ODE integration,
   WhenEvent for boundary events,
   Clip for box constraints, and
   LyapunovSolve for linear reference-model analysis.
*)

Problem 1 (Projection and the Lyapunov Derivative): For the scalar MRAC error model, prove that the projected update preserves \(\dot{V}\le-a_m e^2\) when \(\boldsymbol{\theta}^{*}\in\Omega\).

Solution:

Write the projected update as \(\dot{\hat{\boldsymbol{\theta} } } =\mathbf{Y}+ [\operatorname{Proj}_{\Omega}(\hat{\boldsymbol{\theta} },\mathbf{Y}) -\mathbf{Y}]\). The term containing \(\mathbf{Y}\) cancels the state–parameter cross term exactly, as in the baseline proof. The remaining term is

\[ |b|\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \left[ \operatorname{Proj}_{\Omega} (\hat{\boldsymbol{\theta} },\mathbf{Y}) -\mathbf{Y} \right], \]

which is nonpositive by the defining projection property. Hence

\[ \dot{V}\le-a_m e^2. \]

Problem 2 (Maximum Normalized Regressor Gain): Find the maximum of \(f(\rho)=\rho/(m_0^2+\nu\rho^2)\) for \(\rho\ge0\).

Solution:

Set the derivative to zero:

\[ f'(\rho)= \frac{m_0^2-\nu\rho^2} {(m_0^2+\nu\rho^2)^2}=0. \]

Therefore \(\rho=m_0/\sqrt{\nu}\). Substitution gives

\[ f_{\max}= \frac{m_0/\sqrt{\nu} } {m_0^2+\nu(m_0^2/\nu)} =\frac{1}{2m_0\sqrt{\nu} }. \]

Problem 3 (Checking the Projection Box): For \(a=0.8\), \(b=1.2\), \(a_m=1.5\), and \(b_m=1.5\), calculate the ideal controller parameters and determine whether they lie in \([-1.5,0.2]\times[0.2,2.0]\).

Solution:

\[ \theta_x^{*}=\frac{0.8-1.5}{1.2} =-0.5833, \qquad \theta_r^{*}=\frac{1.5}{1.2}=1.25. \]

Both values are strictly inside their respective intervals. Therefore the projection box does not exclude exact matching for the nominal plant.

Problem 4 (Incorrect Projection Set): Suppose the upper bound for \(\hat{\theta}_r\) is incorrectly chosen as \(1.0\). Explain the expected stability and performance effect.

Solution:

The ideal value \(\theta_r^{*}=1.25\) is no longer feasible. Projection still prevents parameter escape, but the controller cannot reproduce the required reference gain. A residual forcing term remains in the error dynamics. The tracking error is therefore expected to approach a nonzero neighborhood whose size increases with the command amplitude and the distance between the feasible boundary and the ideal value.

Problem 5 (Interpreting the Simulation): Why does the projected-normalized mode have a smaller maximum update norm but a larger IAE than the projection-only mode?

Solution:

Normalization divides the update by \(1+\nu\|\boldsymbol{\phi}\|^2\). During large-command intervals, this denominator substantially reduces the effective adaptation gain. The parameters change more smoothly and the peak update is smaller, but the controller takes longer to approach the gains needed for tracking. The accumulated transient error therefore increases.

Problem 6 (Proof Audit): A designer applies the normalized direct-MRAC update and writes \(\dot{V}=-a_m e^2\) without further derivation. Identify the missing term.

Solution:

The parameter part of the derivative cancels only \(1/m^2\) of the original cross term. The omitted term is

\[ b e\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\phi} \left(1-\frac{1}{m^2}\right). \]

Because this term has no fixed sign, the unnormalized proof cannot be reused unchanged.

Projection and normalization solve different problems. Projection enforces admissible parameter bounds and, when its geometric inequality is satisfied and the ideal parameter lies in the set, preserves the classical MRAC Lyapunov result. Incorrectly tight bounds introduce model mismatch and a nonzero tracking neighborhood.

Normalization bounds the sensitivity of the adaptive update to regressor amplitude. It reduces update spikes and often improves numerical and noise behavior, but it also lowers the effective adaptation rate. Most importantly, inserting normalization into a direct MRAC law changes the Lyapunov derivative. A rigorous design must normalize a compatible error channel or re-derive the complete tracking proof.

The supplied implementations demonstrate the practical trade-off: projection eliminates parameter-bound violations with almost no penalty when the bounds are correct, whereas normalization substantially reduces peak update rates but can increase transient tracking error.