Students are assumed to know the scalar and vector MRAC constructions from Chapters 5–7, the robust-modification motivation from Chapter 8, and the definition of a projection operator from Lessons 1–2 of this chapter. By the end of this lesson, students should be able to:

• insert a projection operator into a continuous-time gradient or Lyapunov adaptive law;
• implement box constraints and smooth convex constraints;
• prove that the parameter estimate remains in its admissible set;
• use the projection inequality in a Lyapunov derivative;
• handle boundary tolerances, ODE-solver overshoot, and numerical retraction; and
• verify the implementation in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The nominal update direction will be denoted by \( \mathbf{y}(t) \in \mathbb{R}^{p} \). In a gradient MRAC law it commonly has the form \( \mathbf{y}=-\Gamma\boldsymbol{\phi}e \). The direction is intentionally left unnormalized here; signal normalization is the subject of Lesson 4.

Let \( \hat{\boldsymbol{\theta} }(t) \) be the adaptive parameter vector and let \( \Omega \subset \mathbb{R}^{p} \) be a closed, convex set known to contain the ideal parameter vector \( \boldsymbol{\theta}^{\ast} \). A nominal law is

\[ \dot{\hat{\boldsymbol{\theta} } }=\mathbf{y} \]

and its projected replacement is

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

Projection changes the update only when the estimate is on, or sufficiently close to, the boundary and the nominal direction points outward. Tangential and inward directions are preserved.

flowchart TD
  R["Reference r"] --> C["Adaptive controller"]
  X["Measured state x"] --> C
  C --> P["Plant"]
  P --> X
  X --> E["Tracking error e"]
  R --> M["Reference model"]
  M --> E
  E --> N["Nominal update direction y"]
  X --> N
  N --> B["Boundary and outward-direction test"]
  H["Current parameter estimate"] --> B
  B --> D["Projected parameter derivative"]
  D --> I["Continuous-time integrator"]
  I --> H
  H --> C
        

This architecture separates three objects that should not be confused: the nominal adaptation direction, the projection map, and the numerical integrator. Correct continuous-time mathematics does not automatically imply a correct finite-step implementation.

The most direct implementation uses independently specified lower and upper bounds:

\[ \Omega_{\mathrm{box} } = \left\{ \boldsymbol{\theta}\in\mathbb{R}^{p}: \ell_i\leq\theta_i\leq u_i,\; i=1,\ldots,p \right\}. \]

For the nominal component \( y_i \), define

\[ \operatorname{Proj}_{[\ell_i,u_i]}(\hat{\theta}_i,y_i) = \begin{cases} 0, & \hat{\theta}_i=\ell_i \text{ and } y_i < 0, \\ 0, & \hat{\theta}_i=u_i \text{ and } y_i > 0, \\ y_i, & \text{otherwise}. \end{cases} \]

The vector projection is obtained component by component:

\[ \operatorname{Proj}_{\Omega_{\mathrm{box} } } (\hat{\boldsymbol{\theta} },\mathbf{y}) = \begin{bmatrix} \operatorname{Proj}_{[\ell_1,u_1]}(\hat{\theta}_1,y_1) \\ \vdots \\ \operatorname{Proj}_{[\ell_p,u_p]}(\hat{\theta}_p,y_p) \end{bmatrix}. \]

In floating-point software, exact equality with a boundary is unreliable. A practical test uses a tolerance \( \varepsilon_b > 0 \):

\[ \hat{\theta}_i\leq\ell_i+\varepsilon_b \quad\text{or}\quad \hat{\theta}_i\geq u_i-\varepsilon_b. \]

A tolerance must be small relative to the width \( u_i-\ell_i \). It is a numerical device, not an additional theoretical boundary layer.

Box constraints are simple but create corners. For an admissible set

\[ \Omega = \left\{ \boldsymbol{\theta}:h(\boldsymbol{\theta})\leq0 \right\}, \]

where \( h \) is differentiable and convex, the outward normal at a regular boundary point is \( \mathbf{n}=\nabla h(\hat{\boldsymbol{\theta} }) \). The tangent-cone projection is

\[ \operatorname{Proj}_{\Omega} (\hat{\boldsymbol{\theta} },\mathbf{y}) = \begin{cases} \mathbf{y}, & h(\hat{\boldsymbol{\theta} }) < 0 \text{ or } \nabla h(\hat{\boldsymbol{\theta} })^{T}\mathbf{y}\leq0, \\ \mathbf{y} - \dfrac{ \nabla h(\hat{\boldsymbol{\theta} }) \nabla h(\hat{\boldsymbol{\theta} })^{T} }{ \left\| \nabla h(\hat{\boldsymbol{\theta} }) \right\|^{2} }\mathbf{y}, & h(\hat{\boldsymbol{\theta} })=0 \text{ and } \nabla h(\hat{\boldsymbol{\theta} })^{T}\mathbf{y} > 0. \end{cases} \]

The second branch removes only the outward normal component. The tangential component remains unchanged.

A common ellipsoid is

\[ h(\boldsymbol{\theta}) = (\boldsymbol{\theta}-\mathbf{c})^{T} W (\boldsymbol{\theta}-\mathbf{c}) -\rho^{2}, \qquad W=W^{T}\succ0, \]

with normal

\[ \nabla h(\boldsymbol{\theta}) = 2W(\boldsymbol{\theta}-\mathbf{c}). \]

Smooth projection operators may introduce an inner set and a thin outer shell so the correction turns on continuously. Such operators are useful when later controller derivations differentiate the adaptive law. The box implementation in this lesson is deliberately transparent and directly exposes the boundary logic.

Theorem 1 (box invariance). Suppose \( \hat{\boldsymbol{\theta} }(0)\in\Omega_{\mathrm{box} } \) and

\[ \dot{\hat{\boldsymbol{\theta} } } = \operatorname{Proj}_{\Omega_{\mathrm{box} } } (\hat{\boldsymbol{\theta} },\mathbf{y}(t)), \]

where \( \mathbf{y}(t) \) is locally bounded. Then \( \hat{\boldsymbol{\theta} }(t)\in\Omega_{\mathrm{box} } \) for every time for which the solution exists.

Proof. Consider component \( i \).

• At \( \hat{\theta}_i=\ell_i \), every outward direction has \( y_i < 0 \) and is replaced by zero. Hence \( \dot{\hat{\theta} }_i\geq0 \) at the lower boundary.
• At \( \hat{\theta}_i=u_i \), every outward direction has \( y_i > 0 \) and is replaced by zero. Hence \( \dot{\hat{\theta} }_i\leq0 \) at the upper boundary.

Therefore the vector field belongs to the tangent cone of the interval at both endpoints. No component can cross its interval boundary, so the Cartesian product of the intervals is forward invariant. \( \square \)

For a smooth convex set, the same conclusion follows because the active branch enforces

\[ \nabla h(\hat{\boldsymbol{\theta} })^{T} \dot{\hat{\boldsymbol{\theta} } } = 0 \]

whenever the nominal direction has a positive outward normal component.

The key analytical property is not merely boundedness. For every ideal parameter vector \( \boldsymbol{\theta}^{\ast}\in\Omega \), a standard projection is selected to satisfy

\[ \left( \hat{\boldsymbol{\theta} } - \boldsymbol{\theta}^{\ast} \right)^{T} \left[ \operatorname{Proj}_{\Omega} (\hat{\boldsymbol{\theta} },\mathbf{y}) - \mathbf{y} \right] \leq0. \]

Proof for the box. If component \( i \) is not projected, its contribution is zero. At the upper boundary with \( y_i > 0 \),

\[ (\hat{\theta}_i-\theta_i^{\ast}) \left[ \operatorname{Proj}(\hat{\theta}_i,y_i)-y_i \right] = -(u_i-\theta_i^{\ast})y_i \leq0. \]

At the lower boundary with \( y_i < 0 \),

\[ (\hat{\theta}_i-\theta_i^{\ast}) \left[ \operatorname{Proj}(\hat{\theta}_i,y_i)-y_i \right] = -(\ell_i-\theta_i^{\ast})y_i \leq0. \]

Summing the component contributions proves the inequality. \( \square \)

Proof for a smooth convex boundary. On an active boundary, let \( \mathbf{n}=\nabla h(\hat{\boldsymbol{\theta} }) \). Then

\[ \operatorname{Proj}_{\Omega} (\hat{\boldsymbol{\theta} },\mathbf{y}) - \mathbf{y} = - \mathbf{n} \dfrac{\mathbf{n}^{T}\mathbf{y} }{\|\mathbf{n}\|^{2} }. \]

Convexity and \( h(\boldsymbol{\theta}^{\ast})\leq0=h(\hat{\boldsymbol{\theta} }) \) imply

\[ \mathbf{n}^{T} (\hat{\boldsymbol{\theta} }-\boldsymbol{\theta}^{\ast}) \geq0. \]

Because the branch is active only when \( \mathbf{n}^{T}\mathbf{y} > 0 \), their product appears with a minus sign and the projection inequality follows.

Reuse the first-order plant and reference-model structure introduced in Chapter 5:

\[ \dot{x}=ax+bu, \qquad \dot{x}_m=a_mx_m+b_mr, \qquad a_m < 0. \]

Assume the sign of \( b \) is known and, for this derivation, take \( b > 0 \). Use

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

The matching parameters satisfy

\[ a+bk_x^{\ast}=a_m, \qquad bk_r^{\ast}=b_m. \]

With tracking error \( e=x-x_m \) and parameter error \( \tilde{\boldsymbol{\theta} } =\hat{\boldsymbol{\theta} }-\boldsymbol{\theta}^{\ast} \),

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

Select the projected adaptive law

\[ \dot{\hat{\boldsymbol{\theta} } } = \operatorname{Proj}_{\Omega} \left( \hat{\boldsymbol{\theta} }, -\gamma e\boldsymbol{\omega} \right), \qquad \gamma > 0. \]

Consider

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

Its derivative is

\[ \begin{aligned} \dot{V} &= a_m e^{2} + be\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\omega} + \frac{b}{\gamma} \tilde{\boldsymbol{\theta} }^{T} \operatorname{Proj}_{\Omega} \left( \hat{\boldsymbol{\theta} }, -\gamma e\boldsymbol{\omega} \right) \\ &= a_m e^{2} + \frac{b}{\gamma} \tilde{\boldsymbol{\theta} }^{T} \left[ \operatorname{Proj}_{\Omega} \left( \hat{\boldsymbol{\theta} }, -\gamma e\boldsymbol{\omega} \right) + \gamma e\boldsymbol{\omega} \right]. \end{aligned} \]

Apply the projection inequality with \( \mathbf{y}=-\gamma e\boldsymbol{\omega} \):

\[ \dot{V}\leq a_m e^{2}\leq0. \]

Thus projection does not destroy the Lyapunov cancellation. It replaces exact cancellation by an inequality in the favorable direction. As in the earlier MRAC lessons, boundedness and convergence of the tracking error follow under the standard regularity assumptions. Parameter convergence still requires excitation and is not guaranteed by projection alone.

The ideal projected differential equation prevents outward motion at the boundary. A finite-step method can nevertheless evaluate the vector field inside the set and jump beyond the boundary before the next evaluation. Therefore a robust software implementation normally combines:

1. projection of the derivative at every ODE stage;
2. a boundary tolerance for floating-point comparisons;
3. small integration steps or event detection; and
4. a final numerical retraction after each accepted step.

\[ \hat{\theta}_{i,k+1} \leftarrow \min \left\{ u_i, \max \left\{ \ell_i, \hat{\theta}_{i,k+1}^{\mathrm{raw} } \right\} \right\}. \]

Retraction corrects integration error. It should not replace the continuous-time projected derivative, because clipping alone hides the active-boundary dynamics and may introduce step-size-dependent behavior.

flowchart TD
  S["Start ODE stage"] --> R["Read theta and nominal direction y"]
  R --> L["Near lower bound and y points lower?"]
  L -->|yes| Z1["Set that derivative \ncomponent to zero"]
  L -->|no| U["Near upper bound and y points upper?"]
  U -->|yes| Z2["Set that derivative \ncomponent to zero"]
  U -->|no| K["Keep nominal \nderivative component"]
  Z1 --> A["Assemble projected derivative"]
  Z2 --> A
  K --> A
  A --> O["Advance one numerical step"]
  O --> C["Retract tiny overshoot into parameter box"]
  C --> N["Log active constraint and continue"]
        

Implementation Checklist

• Verify \( \ell_i < u_i \) for every component.
• Initialize inside the set or explicitly project the initial estimate.
• Choose bounds from credible physical or design knowledge, not merely to improve a plot.
• Record which constraints are active; frequent activation can indicate poor bounds or excessive adaptation gain.
• Keep actuator saturation logic separate from parameter projection.
• Test inward, outward, and tangential directions at every boundary face.

The multilingual implementations use

\[ a=-0.5,\quad b=1,\quad a_m=-2,\quad b_m=2,\quad\gamma=8, \]

so the matching parameters are

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

The admissible box is

\[ -2\leq\hat{k}_x\leq-0.2, \qquad 0.5\leq\hat{k}_r\leq2.5. \]

A bounded sinusoidal disturbance is included as a stress test:

\[ d(t)=0.15\sin(4t). \]

In the ideal disturbance-free proof, \( \dot V\leq a_m e^2 \). With a disturbance, an additional term \( ed(t) \) appears and the basic result becomes a boundedness or ultimate-bound argument unless another robust modification is added. The experiment is intended to expose parameter drift pressure, not to replace the Chapter 8 robustness analysis.

With the supplied step size, the projected run keeps both parameters inside the box, while the identical unprojected law leaves it. The implementations also write a CSV file so the trajectories can be checked independently.

This implementation uses NumPy for arrays and Matplotlib for plots. The RK4 stages call the projection function separately, and the accepted state is retracted to the box.

Chapter9_Lesson3.py

"""Chapter 9, Lesson 3: continuous-time MRAC with box projection.

The script compares a projected adaptive law with the same unprojected law.
Dependencies: NumPy and Matplotlib.
"""

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
from typing import Callable

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    a: float = -0.5
    b: float = 1.0
    a_m: float = -2.0
    b_m: float = 2.0
    gamma: float = 8.0
    dt: float = 1.0e-3
    final_time: float = 20.0
    lower: tuple[float, float] = (-2.0, 0.5)
    upper: tuple[float, float] = (-0.2, 2.5)


CFG = Config()
LOWER = np.asarray(CFG.lower, dtype=float)
UPPER = np.asarray(CFG.upper, dtype=float)


def reference(t: float) -> float:
    return 1.0 + 0.5 * np.sin(0.7 * t) + 0.25 * np.sin(1.9 * t)


def disturbance(t: float) -> float:
    return 0.15 * np.sin(4.0 * t)


def box_projection(theta: np.ndarray, direction: np.ndarray, tol: float = 1.0e-10) -> tuple[np.ndarray, bool]:
    """Project a velocity onto the tangent cone of a parameter box."""
    projected = direction.copy()
    active = False
    for i in range(theta.size):
        outward_at_lower = theta[i] <= LOWER[i] + tol and direction[i] < 0.0
        outward_at_upper = theta[i] >= UPPER[i] - tol and direction[i] > 0.0
        if outward_at_lower or outward_at_upper:
            projected[i] = 0.0
            active = True
    return projected, active


def dynamics(t: float, state: np.ndarray, use_projection: bool) -> tuple[np.ndarray, bool]:
    x, x_m, k_x, k_r = state
    r = reference(t)
    u = k_x * x + k_r * r
    e = x - x_m

    nominal_direction = -CFG.gamma * e * np.array([x, r], dtype=float)
    if use_projection:
        theta_dot, active = box_projection(np.array([k_x, k_r]), nominal_direction)
    else:
        theta_dot, active = nominal_direction, False

    return np.array(
        [
            CFG.a * x + CFG.b * u + disturbance(t),
            CFG.a_m * x_m + CFG.b_m * r,
            theta_dot[0],
            theta_dot[1],
        ],
        dtype=float,
    ), active


def rk4_step(t: float, state: np.ndarray, step: float, use_projection: bool) -> tuple[np.ndarray, bool]:
    k1, a1 = dynamics(t, state, use_projection)
    k2, a2 = dynamics(t + 0.5 * step, state + 0.5 * step * k1, use_projection)
    k3, a3 = dynamics(t + 0.5 * step, state + 0.5 * step * k2, use_projection)
    k4, a4 = dynamics(t + step, state + step * k3, use_projection)
    next_state = state + (step / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)

    # Numerical retraction: the continuous-time projection guarantees invariance,
    # but a finite integration step can overshoot a boundary slightly.
    if use_projection:
        next_state[2:4] = np.clip(next_state[2:4], LOWER, UPPER)

    return next_state, bool(a1 or a2 or a3 or a4)


def simulate(use_projection: bool) -> dict[str, np.ndarray]:
    steps = int(round(CFG.final_time / CFG.dt))
    time = np.linspace(0.0, CFG.final_time, steps + 1)
    state = np.zeros((steps + 1, 4), dtype=float)
    active = np.zeros(steps + 1, dtype=bool)
    state[0] = np.array([1.0, 0.0, LOWER[0], LOWER[1]], dtype=float)

    for k in range(steps):
        state[k + 1], active[k] = rk4_step(time[k], state[k], CFG.dt, use_projection)

    return {
        "time": time,
        "state": state,
        "error": state[:, 0] - state[:, 1],
        "active": active,
    }


def save_csv(path: Path, projected: dict[str, np.ndarray], plain: dict[str, np.ndarray]) -> None:
    table = np.column_stack(
        [
            projected["time"],
            projected["state"],
            projected["error"],
            projected["active"].astype(int),
            plain["state"][:, 2],
            plain["state"][:, 3],
            plain["error"],
        ]
    )
    header = "t,x_projected,xm_projected,kx_projected,kr_projected,e_projected,projection_active,kx_plain,kr_plain,e_plain"
    np.savetxt(path, table, delimiter=",", header=header, comments="")


def plot_results(projected: dict[str, np.ndarray], plain: dict[str, np.ndarray]) -> None:
    t = projected["time"]

    plt.figure()
    plt.plot(t, projected["state"][:, 0], label="x projected")
    plt.plot(t, projected["state"][:, 1], "--", label="x_m")
    plt.xlabel("Time (s)")
    plt.ylabel("State")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(t, projected["state"][:, 2], label="k_x projected")
    plt.plot(t, projected["state"][:, 3], label="k_r projected")
    plt.plot(t, plain["state"][:, 2], "--", label="k_x plain")
    plt.plot(t, plain["state"][:, 3], "--", label="k_r plain")
    plt.axhline(LOWER[0], linestyle=":")
    plt.axhline(UPPER[0], linestyle=":")
    plt.axhline(LOWER[1], linestyle=":")
    plt.axhline(UPPER[1], linestyle=":")
    plt.xlabel("Time (s)")
    plt.ylabel("Adaptive parameters")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


def main() -> None:
    projected = simulate(use_projection=True)
    plain = simulate(use_projection=False)
    save_csv(Path("Chapter9_Lesson3_results.csv"), projected, plain)

    theta_projected = projected["state"][:, 2:4]
    theta_plain = plain["state"][:, 2:4]
    print("Projected parameter minima:", theta_projected.min(axis=0))
    print("Projected parameter maxima:", theta_projected.max(axis=0))
    print("Plain parameter minima:", theta_plain.min(axis=0))
    print("Plain parameter maxima:", theta_plain.max(axis=0))
    print("RK4 steps with active projection:", int(projected["active"].sum()))

    plot_results(projected, plain)


if __name__ == "__main__":
    main()

The C++17 version uses only the standard library. It writes the same CSV columns as the other implementations and uses std::clamp for numerical retraction.

Chapter9_Lesson3.cpp

// Chapter 9, Lesson 3: continuous-time MRAC with box projection.
// Build: g++ -std=c++17 -O2 Chapter9_Lesson3.cpp -o Chapter9_Lesson3

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

namespace {
constexpr double A = -0.5;
constexpr double B = 1.0;
constexpr double AM = -2.0;
constexpr double BM = 2.0;
constexpr double GAMMA = 8.0;
constexpr double DT = 1.0e-3;
constexpr double FINAL_TIME = 20.0;
constexpr double TOL = 1.0e-10;

using State = std::array<double, 4>;  // x, x_m, k_x, k_r
using Pair = std::array<double, 2>;

constexpr Pair LOWER{-2.0, 0.5};
constexpr Pair UPPER{-0.2, 2.5};

struct Derivative {
    State value{};
    bool projection_active{false};
};

double reference(double t) {
    return 1.0 + 0.5 * std::sin(0.7 * t) + 0.25 * std::sin(1.9 * t);
}

double disturbance(double t) {
    return 0.15 * std::sin(4.0 * t);
}

std::pair<Pair, bool> box_projection(const Pair& theta, const Pair& direction) {
    Pair projected = direction;
    bool active = false;
    for (std::size_t i = 0; i < theta.size(); ++i) {
        const bool outward_at_lower = theta[i] <= LOWER[i] + TOL && direction[i] < 0.0;
        const bool outward_at_upper = theta[i] >= UPPER[i] - TOL && direction[i] > 0.0;
        if (outward_at_lower || outward_at_upper) {
            projected[i] = 0.0;
            active = true;
        }
    }
    return {projected, active};
}

Derivative dynamics(double t, const State& z, bool use_projection) {
    const double x = z[0];
    const double x_m = z[1];
    const double k_x = z[2];
    const double k_r = z[3];
    const double r = reference(t);
    const double u = k_x * x + k_r * r;
    const double e = x - x_m;

    Pair direction{-GAMMA * e * x, -GAMMA * e * r};
    bool active = false;
    if (use_projection) {
        const auto result = box_projection(Pair{k_x, k_r}, direction);
        direction = result.first;
        active = result.second;
    }

    return { {A * x + B * u + disturbance(t), AM * x_m + BM * r,
             direction[0], direction[1]}, active};
}

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

std::pair<State, bool> rk4_step(double t, const State& z, double h, bool use_projection) {
    const Derivative k1 = dynamics(t, z, use_projection);
    const Derivative k2 = dynamics(t + 0.5 * h, add_scaled(z, k1.value, 0.5 * h), use_projection);
    const Derivative k3 = dynamics(t + 0.5 * h, add_scaled(z, k2.value, 0.5 * h), use_projection);
    const Derivative k4 = dynamics(t + h, add_scaled(z, k3.value, h), use_projection);

    State next{};
    for (std::size_t i = 0; i < z.size(); ++i) {
        next[i] = z[i] + (h / 6.0) *
            (k1.value[i] + 2.0 * k2.value[i] + 2.0 * k3.value[i] + k4.value[i]);
    }

    if (use_projection) {
        next[2] = std::clamp(next[2], LOWER[0], UPPER[0]);
        next[3] = std::clamp(next[3], LOWER[1], UPPER[1]);
    }

    return {next, k1.projection_active || k2.projection_active ||
                  k3.projection_active || k4.projection_active};
}

}  // namespace

int main() {
    try {
        const int steps = static_cast<int>(std::llround(FINAL_TIME / DT));
        State projected{1.0, 0.0, LOWER[0], LOWER[1]};
        State plain = projected;

        std::ofstream csv("Chapter9_Lesson3_results.csv");
        if (!csv) {
            throw std::runtime_error("Cannot open Chapter9_Lesson3_results.csv");
        }
        csv << "t,x_projected,xm_projected,kx_projected,kr_projected,e_projected,"
               "projection_active,kx_plain,kr_plain,e_plain\n";
        csv << std::setprecision(12);

        Pair min_projected{projected[2], projected[3]};
        Pair max_projected = min_projected;
        Pair min_plain{plain[2], plain[3]};
        Pair max_plain = min_plain;
        int active_steps = 0;

        for (int k = 0; k <= steps; ++k) {
            const double t = k * DT;
            bool active = false;
            if (k < steps) {
                const auto projected_step = rk4_step(t, projected, DT, true);
                const auto plain_step = rk4_step(t, plain, DT, false);
                active = projected_step.second;

                csv << t << ',' << projected[0] << ',' << projected[1] << ','
                    << projected[2] << ',' << projected[3] << ','
                    << projected[0] - projected[1] << ',' << (active ? 1 : 0) << ','
                    << plain[2] << ',' << plain[3] << ',' << plain[0] - plain[1] << '\n';

                projected = projected_step.first;
                plain = plain_step.first;
                active_steps += active ? 1 : 0;
            }

            for (std::size_t i = 0; i < 2; ++i) {
                min_projected[i] = std::min(min_projected[i], projected[i + 2]);
                max_projected[i] = std::max(max_projected[i], projected[i + 2]);
                min_plain[i] = std::min(min_plain[i], plain[i + 2]);
                max_plain[i] = std::max(max_plain[i], plain[i + 2]);
            }
        }

        std::cout << "Projected parameter ranges: k_x=[" << min_projected[0] << ", "
                  << max_projected[0] << "], k_r=[" << min_projected[1] << ", "
                  << max_projected[1] << "]\n";
        std::cout << "Plain parameter ranges: k_x=[" << min_plain[0] << ", "
                  << max_plain[0] << "], k_r=[" << min_plain[1] << ", "
                  << max_plain[1] << "]\n";
        std::cout << "RK4 steps with active projection: " << active_steps << '\n';
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

The Java version uses the standard numerical and file APIs. Records are used to return a derivative together with the active-constraint flag.

Chapter9_Lesson3.java

// Chapter 9, Lesson 3: continuous-time MRAC with box projection.
// Build: javac Chapter9_Lesson3.java
// Run:   java Chapter9_Lesson3

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.charset.StandardCharsets;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Locale;

public final class Chapter9_Lesson3 {
    private static final double A = -0.5;
    private static final double B = 1.0;
    private static final double AM = -2.0;
    private static final double BM = 2.0;
    private static final double GAMMA = 8.0;
    private static final double DT = 1.0e-3;
    private static final double FINAL_TIME = 20.0;
    private static final double TOL = 1.0e-10;
    private static final double[] LOWER = {-2.0, 0.5};
    private static final double[] UPPER = {-0.2, 2.5};

    private record Derivative(double[] value, boolean active) {}
    private record StepResult(double[] state, boolean active) {}

    private Chapter9_Lesson3() {}

    private static double reference(double t) {
        return 1.0 + 0.5 * Math.sin(0.7 * t) + 0.25 * Math.sin(1.9 * t);
    }

    private static double disturbance(double t) {
        return 0.15 * Math.sin(4.0 * t);
    }

    private static Derivative dynamics(double t, double[] z, boolean useProjection) {
        double x = z[0];
        double xm = z[1];
        double kx = z[2];
        double kr = z[3];
        double r = reference(t);
        double u = kx * x + kr * r;
        double e = x - xm;

        double[] direction = {-GAMMA * e * x, -GAMMA * e * r};
        boolean active = false;
        if (useProjection) {
            double[] theta = {kx, kr};
            for (int i = 0; i < 2; i++) {
                boolean outwardAtLower = theta[i] <= LOWER[i] + TOL && direction[i] < 0.0;
                boolean outwardAtUpper = theta[i] >= UPPER[i] - TOL && direction[i] > 0.0;
                if (outwardAtLower || outwardAtUpper) {
                    direction[i] = 0.0;
                    active = true;
                }
            }
        }

        return new Derivative(new double[] {
            A * x + B * u + disturbance(t),
            AM * xm + BM * r,
            direction[0],
            direction[1]
        }, active);
    }

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

    private static StepResult rk4Step(double t, double[] z, boolean useProjection) {
        Derivative k1 = dynamics(t, z, useProjection);
        Derivative k2 = dynamics(t + 0.5 * DT, addScaled(z, k1.value(), 0.5 * DT), useProjection);
        Derivative k3 = dynamics(t + 0.5 * DT, addScaled(z, k2.value(), 0.5 * DT), useProjection);
        Derivative k4 = dynamics(t + DT, addScaled(z, k3.value(), DT), useProjection);

        double[] next = new double[z.length];
        for (int i = 0; i < z.length; i++) {
            next[i] = z[i] + (DT / 6.0) *
                (k1.value()[i] + 2.0 * k2.value()[i] + 2.0 * k3.value()[i] + k4.value()[i]);
        }

        if (useProjection) {
            next[2] = Math.max(LOWER[0], Math.min(UPPER[0], next[2]));
            next[3] = Math.max(LOWER[1], Math.min(UPPER[1], next[3]));
        }

        boolean active = k1.active() || k2.active() || k3.active() || k4.active();
        return new StepResult(next, active);
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        int steps = (int) Math.round(FINAL_TIME / DT);
        double[] projected = {1.0, 0.0, LOWER[0], LOWER[1]};
        double[] plain = projected.clone();
        double[] minProjected = {projected[2], projected[3]};
        double[] maxProjected = minProjected.clone();
        double[] minPlain = {plain[2], plain[3]};
        double[] maxPlain = minPlain.clone();
        int activeSteps = 0;

        Path output = Path.of("Chapter9_Lesson3_results.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(output, StandardCharsets.UTF_8)) {
            writer.write("t,x_projected,xm_projected,kx_projected,kr_projected,e_projected," +
                         "projection_active,kx_plain,kr_plain,e_plain\n");

            for (int k = 0; k <= steps; k++) {
                double t = k * DT;
                if (k < steps) {
                    StepResult projectedStep = rk4Step(t, projected, true);
                    StepResult plainStep = rk4Step(t, plain, false);

                    writer.write(String.format(Locale.US,
                        "%.9f,%.12f,%.12f,%.12f,%.12f,%.12f,%d,%.12f,%.12f,%.12f%n",
                        t, projected[0], projected[1], projected[2], projected[3],
                        projected[0] - projected[1], projectedStep.active() ? 1 : 0,
                        plain[2], plain[3], plain[0] - plain[1]));

                    projected = projectedStep.state();
                    plain = plainStep.state();
                    if (projectedStep.active()) {
                        activeSteps++;
                    }
                }

                for (int i = 0; i < 2; i++) {
                    minProjected[i] = Math.min(minProjected[i], projected[i + 2]);
                    maxProjected[i] = Math.max(maxProjected[i], projected[i + 2]);
                    minPlain[i] = Math.min(minPlain[i], plain[i + 2]);
                    maxPlain[i] = Math.max(maxPlain[i], plain[i + 2]);
                }
            }
        } catch (IOException ex) {
            System.err.println("Cannot write results: " + ex.getMessage());
            System.exit(1);
        }

        System.out.printf("Projected parameter ranges: k_x=[%.6f, %.6f], k_r=[%.6f, %.6f]%n",
            minProjected[0], maxProjected[0], minProjected[1], maxProjected[1]);
        System.out.printf("Plain parameter ranges: k_x=[%.6f, %.6f], k_r=[%.6f, %.6f]%n",
            minPlain[0], maxPlain[0], minPlain[1], maxPlain[1]);
        System.out.println("RK4 steps with active projection: " + activeSteps);
    }
}

The MATLAB script uses local functions, a fixed-step RK4 method, CSV export, and plots. It requires base MATLAB.

Chapter9_Lesson3.m

% Chapter 9, Lesson 3: continuous-time MRAC with box projection.
% Uses base MATLAB only and compares projected and unprojected adaptive laws.

clear; clc; close all;

cfg.a = -0.5;
cfg.b = 1.0;
cfg.am = -2.0;
cfg.bm = 2.0;
cfg.gamma = 8.0;
cfg.dt = 1.0e-3;
cfg.finalTime = 20.0;
cfg.lower = [-2.0; 0.5];
cfg.upper = [-0.2; 2.5];

[t, zProjected, active] = simulateMRAC(cfg, true);
[~, zPlain, ~] = simulateMRAC(cfg, false);

eProjected = zProjected(:,1) - zProjected(:,2);
ePlain = zPlain(:,1) - zPlain(:,2);

results = table(t, zProjected(:,1), zProjected(:,2), zProjected(:,3), ...
    zProjected(:,4), eProjected, active, zPlain(:,3), zPlain(:,4), ePlain, ...
    'VariableNames', {'t','x_projected','xm_projected','kx_projected', ...
    'kr_projected','e_projected','projection_active','kx_plain','kr_plain','e_plain'});
writetable(results, 'Chapter9_Lesson3_results.csv');

fprintf('Projected parameter minima: [%.6f, %.6f]\n', min(zProjected(:,3)), min(zProjected(:,4)));
fprintf('Projected parameter maxima: [%.6f, %.6f]\n', max(zProjected(:,3)), max(zProjected(:,4)));
fprintf('Plain parameter minima: [%.6f, %.6f]\n', min(zPlain(:,3)), min(zPlain(:,4)));
fprintf('Plain parameter maxima: [%.6f, %.6f]\n', max(zPlain(:,3)), max(zPlain(:,4)));
fprintf('RK4 steps with active projection: %d\n', nnz(active));

figure;
plot(t, zProjected(:,1), 'LineWidth', 1.2); hold on;
plot(t, zProjected(:,2), '--', 'LineWidth', 1.2);
grid on; xlabel('Time (s)'); ylabel('State');
legend('x projected','x_m','Location','best');
title('Projected MRAC tracking');

figure;
plot(t, zProjected(:,3), 'LineWidth', 1.2); hold on;
plot(t, zProjected(:,4), 'LineWidth', 1.2);
plot(t, zPlain(:,3), '--', 'LineWidth', 1.0);
plot(t, zPlain(:,4), '--', 'LineWidth', 1.0);
yline(cfg.lower(1), ':'); yline(cfg.upper(1), ':');
yline(cfg.lower(2), ':'); yline(cfg.upper(2), ':');
grid on; xlabel('Time (s)'); ylabel('Adaptive parameters');
legend('k_x projected','k_r projected','k_x plain','k_r plain','Location','best');
title('Projection enforces the parameter box');

function [t, z, active] = simulateMRAC(cfg, useProjection)
    steps = round(cfg.finalTime / cfg.dt);
    t = (0:steps)' * cfg.dt;
    z = zeros(steps + 1, 4);
    active = false(steps + 1, 1);
    z(1,:) = [1.0, 0.0, cfg.lower(1), cfg.lower(2)];

    for k = 1:steps
        [z(k+1,:), active(k)] = rk4Step(t(k), z(k,:), cfg, useProjection);
    end
end

function [next, active] = rk4Step(t, z, cfg, useProjection)
    [k1, a1] = rhs(t, z, cfg, useProjection);
    [k2, a2] = rhs(t + 0.5*cfg.dt, z + 0.5*cfg.dt*k1, cfg, useProjection);
    [k3, a3] = rhs(t + 0.5*cfg.dt, z + 0.5*cfg.dt*k2, cfg, useProjection);
    [k4, a4] = rhs(t + cfg.dt, z + cfg.dt*k3, cfg, useProjection);
    next = z + (cfg.dt/6.0) * (k1 + 2*k2 + 2*k3 + k4);

    if useProjection
        next(3:4) = min(max(next(3:4)', cfg.lower), cfg.upper)';
    end
    active = a1 || a2 || a3 || a4;
end

function [dz, active] = rhs(t, z, cfg, useProjection)
    x = z(1); xm = z(2); kx = z(3); kr = z(4);
    r = referenceSignal(t);
    u = kx*x + kr*r;
    e = x - xm;
    direction = -cfg.gamma * e * [x; r];

    if useProjection
        [thetaDot, active] = boxProjection([kx; kr], direction, cfg.lower, cfg.upper);
    else
        thetaDot = direction;
        active = false;
    end

    dz = [cfg.a*x + cfg.b*u + disturbanceSignal(t), ...
          cfg.am*xm + cfg.bm*r, thetaDot(1), thetaDot(2)];
end

function [projected, active] = boxProjection(theta, direction, lower, upper)
    tol = 1.0e-10;
    projected = direction;
    active = false;
    for i = 1:numel(theta)
        outwardLower = theta(i) <= lower(i) + tol && direction(i) < 0.0;
        outwardUpper = theta(i) >= upper(i) - tol && direction(i) > 0.0;
        if outwardLower || outwardUpper
            projected(i) = 0.0;
            active = true;
        end
    end
end

function r = referenceSignal(t)
    r = 1.0 + 0.5*sin(0.7*t) + 0.25*sin(1.9*t);
end

function d = disturbanceSignal(t)
    d = 0.15*sin(4.0*t);
end

The following MATLAB script programmatically creates a Simulink model containing a reference generator, a MATLAB Function block for the projected MRAC dynamics, a continuous integrator, and workspace logging. It requires Simulink.

Chapter9_Lesson3_Simulink.m

% Chapter 9, Lesson 3: programmatically build a Simulink projection subsystem.
% Requires Simulink. The generated model is Chapter9_Lesson3_ProjectionModel.slx.

clear; clc;

model = 'Chapter9_Lesson3_ProjectionModel';
if bdIsLoaded(model)
    close_system(model, 0);
end
new_system(model);
open_system(model);

set_param(model, 'Solver', 'ode4', 'FixedStep', '0.001', 'StopTime', '20');

add_block('simulink/Sources/Clock', [model '/Clock'], 'Position', [30 40 60 70]);
add_block('simulink/User-Defined Functions/MATLAB Function', [model '/Reference'], ...
    'Position', [110 25 260 85]);
add_block('simulink/User-Defined Functions/MATLAB Function', [model '/Projected MRAC Dynamics'], ...
    'Position', [330 20 560 150]);
add_block('simulink/Continuous/Integrator', [model '/State Integrator'], ...
    'InitialCondition', '[1;0;-2;0.5]', 'Position', [630 45 660 125]);
add_block('simulink/Sinks/To Workspace', [model '/Log State'], ...
    'VariableName', 'zLog', 'SaveFormat', 'Structure With Time', ...
    'Position', [730 55 830 95]);

referenceCode = sprintf([ ...
    'function r = fcn(t)\n' ...
    'r = 1.0 + 0.5*sin(0.7*t) + 0.25*sin(1.9*t);\n' ...
    'end\n']);

mracCode = sprintf([ ...
    'function dz = fcn(t,r,z)\n' ...
    '%%#codegen\n' ...
    'a=-0.5; b=1.0; am=-2.0; bm=2.0; gamma=8.0;\n' ...
    'lower=[-2.0;0.5]; upper=[-0.2;2.5]; tol=1e-10;\n' ...
    'x=z(1); xm=z(2); theta=z(3:4);\n' ...
    'u=theta(1)*x + theta(2)*r; e=x-xm;\n' ...
    'v=-gamma*e*[x;r]; thetaDot=v;\n' ...
    'for i=1:2\n' ...
    '  if (theta(i)<=lower(i)+tol && v(i)<0) || ...\n' ...
    '     (theta(i)>=upper(i)-tol && v(i)>0)\n' ...
    '    thetaDot(i)=0;\n' ...
    '  end\n' ...
    'end\n' ...
    'd=0.15*sin(4.0*t);\n' ...
    'dz=[a*x+b*u+d; am*xm+bm*r; thetaDot];\n' ...
    'end\n']);

setFunctionScript([model '/Reference'], referenceCode);
setFunctionScript([model '/Projected MRAC Dynamics'], mracCode);

add_line(model, 'Clock/1', 'Reference/1', 'autorouting', 'on');
add_line(model, 'Clock/1', 'Projected MRAC Dynamics/1', 'autorouting', 'on');
add_line(model, 'Reference/1', 'Projected MRAC Dynamics/2', 'autorouting', 'on');
add_line(model, 'State Integrator/1', 'Projected MRAC Dynamics/3', 'autorouting', 'on');
add_line(model, 'Projected MRAC Dynamics/1', 'State Integrator/1', 'autorouting', 'on');
add_line(model, 'State Integrator/1', 'Log State/1', 'autorouting', 'on');

save_system(model, [model '.slx']);
fprintf('Created %s.slx\n', model);
fprintf('Run the model, then inspect zLog.signals.values.\n');

function setFunctionScript(blockPath, scriptText)
    chart = find(sfroot, '-isa', 'Stateflow.EMChart', 'Path', blockPath);
    if isempty(chart)
        error('Could not access MATLAB Function block: %s', blockPath);
    end
    chart.Script = scriptText;
end

The generated integrator contains the full state \( [x,\;x_m,\;\hat{k}_x,\;\hat{k}_r]^T \). For production use, add explicit parameter-bound monitors and consider zero-crossing detection at each active face.

The notebook uses NDSolveValue for the projected and plain systems and exports sampled trajectories. The projection is written as a reusable Wolfram Language function.

Chapter9_Lesson3.nb

Notebook[
 {
  Cell[
   "Chapter 9, Lesson 3: Continuous-Time MRAC with Box Projection",
   "Title"
  ],

  Cell[
   "The notebook compares projected and unprojected adaptive laws and verifies the parameter bounds.",
   "Text"
  ],

  Cell[
   BoxData[
"ClearAll[\"Global`*\"];

a = -0.5;
b = 1.0;
am = -2.0;
bm = 2.0;
gamma = 8.0;

lower = {-2.0, 0.5};
upper = {-0.2, 2.5};

reference[t_] :=
  1.0 + 0.5 Sin[0.7 t] + 0.25 Sin[1.9 t];

disturbance[t_] :=
  0.15 Sin[4.0 t];

boxProjection[theta_List, direction_List] :=
 Module[
  {p = direction, tol = 10^-10},

  Do[
   If[
    (theta[[i]] <= lower[[i]] + tol &&
       direction[[i]] < 0) ||
     (theta[[i]] >= upper[[i]] - tol &&
       direction[[i]] > 0),

    p[[i]] = 0
   ],
   {i, Length[theta]}
  ];

  p
 ];

projectedSolution =
 NDSolveValue[
  {
   x'[t] ==
    a x[t] +
     b (kx[t] x[t] + kr[t] reference[t]) +
     disturbance[t],

   xm'[t] ==
    am xm[t] + bm reference[t],

   kx'[t] ==
    boxProjection[
      {kx[t], kr[t]},
      -gamma (x[t] - xm[t]) {x[t], reference[t]}
     ][[1]],

   kr'[t] ==
    boxProjection[
      {kx[t], kr[t]},
      -gamma (x[t] - xm[t]) {x[t], reference[t]}
     ][[2]],

   x[0] == 1.0,
   xm[0] == 0.0,
   kx[0] == lower[[1]],
   kr[0] == lower[[2]]
  },

  {x, xm, kx, kr},
  {t, 0, 20},

  Method -> {
   \"TimeIntegration\" -> {
    \"EventLocator\",
    \"EventAction\" -> None
   }
  },

  MaxStepFraction -> 1/20000
 ];

plainSolution =
 NDSolveValue[
  {
   xp'[t] ==
    a xp[t] +
     b (kxp[t] xp[t] + krp[t] reference[t]) +
     disturbance[t],

   xmp'[t] ==
    am xmp[t] + bm reference[t],

   kxp'[t] ==
    -gamma (xp[t] - xmp[t]) xp[t],

   krp'[t] ==
    -gamma (xp[t] - xmp[t]) reference[t],

   xp[0] == 1.0,
   xmp[0] == 0.0,
   kxp[0] == lower[[1]],
   krp[0] == lower[[2]]
  },

  {xp, xmp, kxp, krp},
  {t, 0, 20},

  MaxStepFraction -> 1/20000
 ];

trackingPlot =
 Plot[
  Evaluate[
   {
    projectedSolution[[1]][t],
    projectedSolution[[2]][t]
   }
  ],

  {t, 0, 20},

  PlotLegends -> {
   \"x projected\",
   \"x_m\"
  },

  AxesLabel -> {
   \"Time (s)\",
   \"State\"
  },

  GridLines -> Automatic
 ];

parameterPlot =
 Plot[
  Evaluate[
   {
    projectedSolution[[3]][t],
    projectedSolution[[4]][t],
    plainSolution[[3]][t],
    plainSolution[[4]][t]
   }
  ],

  {t, 0, 20},

  PlotLegends -> {
   \"k_x projected\",
   \"k_r projected\",
   \"k_x plain\",
   \"k_r plain\"
  },

  AxesLabel -> {
   \"Time (s)\",
   \"Adaptive parameters\"
  },

  GridLines -> Automatic
 ];

Show[trackingPlot]
Show[parameterPlot]

samples =
 Table[
  {
   tt,
   Sequence @@ Through[projectedSolution[tt]],
   Sequence @@ Through[plainSolution[tt]]
  },

  {tt, 0, 20, 0.01}
 ];

Export[
 \"Chapter9_Lesson3_results.csv\",
 samples
];

projectedRange =
 MinMax /@
  Transpose[
   samples[[All, {4, 5}]]
  ];

plainRange =
 MinMax /@
  Transpose[
   samples[[All, {8, 9}]]
  ];

Print[
 \"Projected parameter ranges: \",
 projectedRange
];

Print[
 \"Plain parameter ranges: \",
 plainRange
];"
   ],
   "Input"
  ]
 },

 WindowTitle -> "Chapter9_Lesson3",
 StyleDefinitions -> "Default.nb"
]        

A projection function should be unit-tested before it is inserted into a closed-loop simulation. For each parameter component, test:

1. Interior: the output equals the nominal direction.
2. Lower boundary, inward: a positive direction is unchanged.
3. Lower boundary, outward: a negative direction becomes zero.
4. Upper boundary, inward: a negative direction is unchanged.
5. Upper boundary, outward: a positive direction becomes zero.
6. Near-boundary tolerance: behavior matches the chosen numerical policy.
7. Long simulation: every logged estimate satisfies its bounds within roundoff.

A useful automated invariant is

\[ \max_{k,i} \left\{ \ell_i-\hat{\theta}_{i,k}, \hat{\theta}_{i,k}-u_i, 0 \right\} \leq\varepsilon_{\mathrm{test} }. \]

A second test numerically checks the Lyapunov-compatible inequality for randomly sampled ideal parameters inside the box:

\[ (\hat{\boldsymbol{\theta} }-\boldsymbol{\theta}^{\ast})^{T} \left[ \operatorname{Proj} (\hat{\boldsymbol{\theta} },\mathbf{y})-\mathbf{y} \right] \leq\varepsilon_{\mathrm{test} }. \]

Problem 1 (component-wise evaluation). Let \( \Omega=[-2,1]\times[0,3] \), \( \hat{\boldsymbol{\theta} }=[1,\;0]^T \), and \( \mathbf{y}=[4,\;-5]^T \). Evaluate the box projection.

Solution. The first estimate is at its upper bound and its direction is positive, so the first component is outward and is removed. The second estimate is at its lower bound and its direction is negative, so the second component is also outward and is removed:

\[ \operatorname{Proj}_{\Omega} (\hat{\boldsymbol{\theta} },\mathbf{y}) = \begin{bmatrix} 0\\ 0 \end{bmatrix}. \]

Problem 2 (projection inequality for one interval). Let \( \theta^{\ast}\in[\ell,u] \). Prove

\[ (\hat{\theta}-\theta^{\ast}) \left[ \operatorname{Proj}_{[\ell,u]}(\hat{\theta},y)-y \right] \leq0. \]

Solution. There are three cases.

1. If projection is inactive, the bracket is zero.
2. At \( \hat{\theta}=u \) with \( y > 0 \), the expression is \( -(u-\theta^{\ast})y\leq0 \).
3. At \( \hat{\theta}=\ell \) with \( y < 0 \), the expression is \( -(\ell-\theta^{\ast})y\leq0 \), because both \( \ell-\theta^{\ast} \) and \( y \) are nonpositive.

These cases exhaust the definition.

Problem 3 (ellipsoidal tangent projection). Let

\[ h(\boldsymbol{\theta}) = \boldsymbol{\theta}^{T} \begin{bmatrix} 4&0\\ 0&1 \end{bmatrix} \boldsymbol{\theta} -1. \]

At \( \hat{\boldsymbol{\theta} }=[0.5,\;0]^T \) and \( \mathbf{y}=[2,\;3]^T \), compute the projected direction.

Solution. The point is on the boundary because \( 4(0.5)^2=1 \). The outward normal is

\[ \mathbf{n} = 2 \begin{bmatrix} 4&0\\ 0&1 \end{bmatrix} \begin{bmatrix} 0.5\\ 0 \end{bmatrix} = \begin{bmatrix} 4\\ 0 \end{bmatrix}. \]

Since \( \mathbf{n}^{T}\mathbf{y}=8 > 0 \), remove the outward normal component:

\[ \mathbf{y}_{p} = \begin{bmatrix} 2\\ 3 \end{bmatrix} - \begin{bmatrix} 4\\ 0 \end{bmatrix} \frac{8}{16} = \begin{bmatrix} 0\\ 3 \end{bmatrix}. \]

Problem 4 (Lyapunov compatibility). For the MRAC error model

\[ \dot e=a_m e+b\tilde{\boldsymbol{\theta} }^T\boldsymbol{\omega}, \qquad a_m < 0,\quad b > 0, \]

show that the projected law from Section 7 yields \( \dot V\leq a_m e^2 \).

Solution. Differentiate

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

Add and subtract the nominal direction \( -\gamma e\boldsymbol{\omega} \). The nominal cross-term cancels the error-model cross-term, leaving

\[ \dot V = a_m e^2 + \frac{b}{\gamma} \tilde{\boldsymbol{\theta} }^T \left[ \operatorname{Proj} (\hat{\boldsymbol{\theta} },-\gamma e\boldsymbol{\omega}) - (-\gamma e\boldsymbol{\omega}) \right]. \]

The bracketed contribution is nonpositive by the projection inequality. Since \( b/\gamma > 0 \), \( \dot V\leq a_m e^2\leq0 \).

Problem 5 (finite-step overshoot). A scalar estimate has upper bound \( u=1 \), current value \( \hat\theta_k=0.99 \), nominal direction \( y=4 \), and forward-Euler step \( \Delta t=0.01 \). Explain why a boundary-only derivative test can produce an invalid numerical state and give the corrected update.

Solution. Because the current value is strictly inside the interval, the continuous-time projection initially leaves the derivative unchanged. A single Euler step gives

\[ \hat\theta_{k+1}^{\mathrm{raw} } = 0.99+0.01(4) = 1.03, \]

which crosses the upper bound. The true continuous-time trajectory would reach the boundary and then have its outward derivative removed. A finite-step implementation should use event detection, a smaller step, or retraction:

\[ \hat\theta_{k+1} = \min \left\{ 1, \max \left\{ \ell, 1.03 \right\} \right\} = 1. \]

Problem 6 (choosing a valid box). For the numerical plant in Section 9, calculate the matching parameters and verify that the selected box contains them.

Solution.

\[ k_x^{\ast} = \frac{-2-(-0.5)}{1} = -1.5, \qquad k_r^{\ast} = \frac{2}{1} = 2. \]

Since \( -2\leq-1.5\leq-0.2 \) and \( 0.5\leq2\leq2.5 \), the ideal parameter vector lies inside the box. If the chosen set excluded the ideal parameter, the projection inequality relative to that ideal vector would no longer be available in the stated form, and exact model matching could be impossible.

A continuous-time projection law modifies only outward parameter motion at the admissible-set boundary. The box implementation is a component-wise tangent-cone rule, while a smooth convex-set implementation removes the outward normal component. Two properties drive the analysis: forward invariance of the parameter set and the projection inequality that makes the correction term nonpositive in the Lyapunov derivative. In software, derivative projection should be combined with stage-wise evaluation, boundary tolerances, suitable integration steps, and a small numerical retraction. Projection bounds parameter estimates; it does not by itself create persistent excitation, identify the true parameters, normalize regressors, or solve actuator saturation.