After completing this lesson, students should be able to:

• define the Euclidean projection of a point onto a closed convex set;
• define the projected direction associated with a parameter-update vector;
• construct nested admissible sets from a continuously differentiable convex function;
• derive the smooth adaptive projection operator and interpret its normal and tangential components;
• prove the key inequality used in Lyapunov analysis;
• prove that the outer admissible set is forward invariant;
• specialize the operator to spherical, ellipsoidal, and box constraints;
• implement and numerically test the operator in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The required background is convex sets, gradients, quadratic Lyapunov functions, parameter error, and the reason parameter drift motivates bounded parameter estimates. Those ideas were introduced in Chapters 2, 3, 8, and Lesson 1 of this chapter.

The word projection is used for two related but distinct operations. Confusing them is a common source of incorrect adaptive-law implementations.

2.1 Static nearest-point projection

Let \( \boldsymbol{\theta} \) be a candidate parameter vector constrained to a nonempty, closed, convex set \( \Omega \subset \mathbb{R}^{p} \). The Euclidean projection of a point \( \mathbf{z} \) is the unique nearest point in that set:

\[ P_{\Omega}(\mathbf{z}) = \underset{\boldsymbol{\xi}\in\Omega}{\operatorname{arg\,min}} \;\frac{1}{2}\|\mathbf{z}-\boldsymbol{\xi}\|_{2}^{2}. \]

For a closed ball \( \Omega=\{\boldsymbol{\xi}:\|\boldsymbol{\xi}\|_{2}\leq r\} \), the solution is

\[ P_{\Omega}(\mathbf{z})= \begin{cases} \mathbf{z}, & \|\mathbf{z}\|_{2}\leq r,\\ r\dfrac{\mathbf{z}}{\|\mathbf{z}\|_{2}}, & \|\mathbf{z}\|_{2}>r. \end{cases} \]

This is a point correction. It answers: “Where should an already computed point be placed so that it belongs to the feasible set?”

2.2 Projection of an update direction

In a continuous-time adaptive law, the unconstrained algorithm proposes a direction \( \mathbf{y} \) at the current estimate \( \widehat{\boldsymbol{\theta}} \). The relevant object is therefore a projected velocity:

\[ \Pi_{\Omega}(\widehat{\boldsymbol{\theta}},\mathbf{y}) = \lim_{h\downarrow 0} \frac{P_{\Omega}(\widehat{\boldsymbol{\theta}}+h\mathbf{y}) -\widehat{\boldsymbol{\theta}}}{h}. \]

This is a direction correction. It answers: “Which part of the proposed update can be retained without allowing the estimate to leave the feasible region?” In the interior, the complete direction is retained. At a smooth boundary, only an outward normal component must be removed.

flowchart TD
  A["Candidate point z"] --> B["Nearest-point map P_Omega(z)"]
  B --> C["Corrected point inside the set"]
  D["Current estimate theta_hat plus direction y"] --> E["Direction map Pi_Omega(theta_hat, y)"]
  E --> F["Feasible velocity for parameter evolution"]
        

The smooth adaptive projection operator is conveniently constructed from a continuously differentiable convex function \( f:\mathbb{R}^{p}→\mathbb{R} \). Define two nested sublevel sets:

\[ \Omega_{0}=\{\boldsymbol{\theta}:f(\boldsymbol{\theta})\leq 0\}, \qquad \Omega_{1}=\{\boldsymbol{\theta}:f(\boldsymbol{\theta})\leq 1\}. \]

Because every sublevel set of a convex function is convex, \( \Omega_{0} \) and \( \Omega_{1} \) are convex, and \( \Omega_{0}\subseteq\Omega_{1} \). Their roles are different:

• \( \Omega_{0} \) is the inner region where adaptation is completely unmodified;
• \( \Omega_{1} \) is the outer invariant region that the estimate must never leave;
• \( \Omega_{A}=\Omega_{1}\setminus\Omega_{0} \) is a boundary layer in which outward motion is gradually reduced.

3.1 Convex sublevel-set proof

Let \( \boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2}\in\Omega_{c} \), where \( \Omega_{c}=\{\boldsymbol{\theta}:f(\boldsymbol{\theta})\leq c\} \), and choose \( \lambda\in[0,1] \). Convexity of \( f \) gives

\[ f\!\left(\lambda\boldsymbol{\theta}_{1} +(1-\lambda)\boldsymbol{\theta}_{2}\right) \leq \lambda f(\boldsymbol{\theta}_{1}) +(1-\lambda)f(\boldsymbol{\theta}_{2}) \leq \lambda c+(1-\lambda)c=c. \]

Hence every line segment joining two points of the sublevel set remains in that set.

3.2 Supporting-gradient inequality

Differentiable convexity also gives the first-order inequality

\[ f(\boldsymbol{\theta}^{*})\geq f(\boldsymbol{\theta})+ \nabla f(\boldsymbol{\theta})^{T} (\boldsymbol{\theta}^{*}-\boldsymbol{\theta}). \]

Rearranging, for any \( \boldsymbol{\theta}^{*}\in\Omega_{0} \) and any point with \( f(\boldsymbol{\theta})>0 \),

\[ (\boldsymbol{\theta}-\boldsymbol{\theta}^{*})^{T} \nabla f(\boldsymbol{\theta}) \geq f(\boldsymbol{\theta})-f(\boldsymbol{\theta}^{*}) \geq f(\boldsymbol{\theta})>0. \]

This inequality is the geometric fact needed in the Lyapunov proof: from an admissible interior point, the gradient points outward through the boundary layer.

Let \( \widehat{\boldsymbol{\theta}}\in\mathbb{R}^{p} \) be the current parameter estimate and let \( \mathbf{y}\in\mathbb{R}^{p} \) be the unconstrained update direction. The smooth projection operator is

\[ \operatorname{Proj} (\widehat{\boldsymbol{\theta}},\mathbf{y},f) = \begin{cases} \mathbf{y}- \dfrac{\nabla f(\widehat{\boldsymbol{\theta}}) \nabla f(\widehat{\boldsymbol{\theta}})^{T}} {\|\nabla f(\widehat{\boldsymbol{\theta}})\|_{2}^{2}} \mathbf{y}\,f(\widehat{\boldsymbol{\theta}}), & f(\widehat{\boldsymbol{\theta}})>0 \;\wedge\; \mathbf{y}^{T}\nabla f(\widehat{\boldsymbol{\theta}})>0,\\[8pt] \mathbf{y}, & \text{otherwise.} \end{cases} \]

The activation conditions have precise meanings:

• \( f(\widehat{\boldsymbol{\theta}})>0 \) means that the estimate lies in the boundary layer rather than the protected inner set;
• \( \mathbf{y}^{T}\nabla f(\widehat{\boldsymbol{\theta}})>0 \) means that the proposed direction has an outward component;
• if either condition is false, the update is left unchanged.

4.1 Tangential-normal decomposition

Introduce the unit outward normal

\[ \mathbf{n}= \frac{\nabla f(\widehat{\boldsymbol{\theta}})} {\|\nabla f(\widehat{\boldsymbol{\theta}})\|_{2}}, \qquad \mathbf{y}_{N}=\mathbf{n}\mathbf{n}^{T}\mathbf{y}, \qquad \mathbf{y}_{T}=\mathbf{y}-\mathbf{y}_{N}. \]

When projection is active,

\[ \operatorname{Proj} (\widehat{\boldsymbol{\theta}},\mathbf{y},f) =\mathbf{y}_{T}+ \left(1-f(\widehat{\boldsymbol{\theta}})\right)\mathbf{y}_{N}. \]

Thus the tangential component is never changed. The outward normal component is retained in full at \( f=0 \), gradually attenuated for \( 0<f<1 \), and completely removed at \( f=1 \).

flowchart TD
  A["Inputs: current estimate theta_hat and raw direction y"] --> B["Compute f(theta_hat) and grad_f"]
  B --> C{"Is f greater than 0?"}
  C -->|no| D["Return y unchanged"]
  C -->|yes| E{"Is dot(y, grad_f) \ngreater than 0?"}
  E -->|no| D
  E -->|yes| F["Remove f times the outward normal component"]
  F --> G["Return tangential part plus reduced normal part"]
        

5.1 Property 1 — Inactivity in the inner set

If \( \widehat{\boldsymbol{\theta}}\in\Omega_{0} \), then \( f(\widehat{\boldsymbol{\theta}})\leq 0 \), so the first activation condition is false and

\[ \operatorname{Proj} (\widehat{\boldsymbol{\theta}},\mathbf{y},f)=\mathbf{y}. \]

Therefore projection does not alter nominal adaptation while the estimate is comfortably inside its admissible range.

5.2 Property 2 — Inward and tangential directions are unchanged

If \( \mathbf{y}^{T}\nabla f(\widehat{\boldsymbol{\theta}})\leq 0 \), then the directional derivative of the constraint function satisfies

\[ Df(\widehat{\boldsymbol{\theta}})[\mathbf{y}] =\nabla f(\widehat{\boldsymbol{\theta}})^{T}\mathbf{y}\leq 0. \]

The proposed update is tangential or moves toward smaller values of \( f \), so no correction is needed.

5.3 Property 3 — The correction cannot increase update magnitude

During active projection, decompose \( \mathbf{y}=\mathbf{y}_{T}+\mathbf{y}_{N} \). Because these components are orthogonal,

\[ \|\operatorname{Proj}(\widehat{\boldsymbol{\theta}},\mathbf{y},f)\|_{2}^{2} =\|\mathbf{y}_{T}\|_{2}^{2} +\left(1-f(\widehat{\boldsymbol{\theta}})\right)^{2} \|\mathbf{y}_{N}\|_{2}^{2} \leq \|\mathbf{y}_{T}\|_{2}^{2}+\|\mathbf{y}_{N}\|_{2}^{2} =\|\mathbf{y}\|_{2}^{2}, \]

where the estimate is assumed to be in \( \Omega_{1} \), so \( 0<f(\widehat{\boldsymbol{\theta}})\leq 1 \) whenever the correction is active.

5.4 Property 4 — Lyapunov-compatible inequality

Let the unknown ideal parameter satisfy \( \boldsymbol{\theta}^{*}\in\Omega_{0} \). The central projection property is

\[ (\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*})^{T} \left[ \operatorname{Proj}(\widehat{\boldsymbol{\theta}},\mathbf{y},f) -\mathbf{y} \right]\leq 0. \]

Proof. If projection is inactive, the bracket is zero and the result is immediate. If it is active, define

\[ \alpha= \frac{\nabla f(\widehat{\boldsymbol{\theta}})^{T}\mathbf{y}} {\|\nabla f(\widehat{\boldsymbol{\theta}})\|_{2}^{2}} f(\widehat{\boldsymbol{\theta}}). \]

Every factor in \( \alpha \) is positive. The correction is \( -\alpha\nabla f \), so

\[ \begin{aligned} &(\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*})^{T} \left[ \operatorname{Proj}(\widehat{\boldsymbol{\theta}},\mathbf{y},f)-\mathbf{y} \right]\\ &\qquad = -\alpha (\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*})^{T} \nabla f(\widehat{\boldsymbol{\theta}}). \end{aligned} \]

From the supporting-gradient inequality in Section 3,

\[ (\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*})^{T} \nabla f(\widehat{\boldsymbol{\theta}}) \geq f(\widehat{\boldsymbol{\theta}})-f(\boldsymbol{\theta}^{*}) \geq f(\widehat{\boldsymbol{\theta}})>0. \]

Multiplication by \( -\alpha<0 \) proves the result.

5.5 Why this inequality matters in adaptive control

For parameter error \( \widetilde{\boldsymbol{\theta}} =\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*} \), a typical Lyapunov term is

\[ V_{\theta}= \frac{1}{2}\widetilde{\boldsymbol{\theta}}^{T} \Gamma^{-1}\widetilde{\boldsymbol{\theta}}, \qquad \Gamma=\Gamma^{T}>0. \]

Replacing an unconstrained update by a correctly matched projection adds a term to \( \dot V_{\theta} \) that is nonpositive rather than destabilizing. Projection therefore preserves the sign structure used in the original Lyapunov argument.

5.6 Property 5 — Forward invariance of the outer set

Consider the projected differential equation

\[ \dot{\widehat{\boldsymbol{\theta}}} =\operatorname{Proj} (\widehat{\boldsymbol{\theta}},\mathbf{y}(t),f), \qquad \widehat{\boldsymbol{\theta}}(0)\in\Omega_{1}. \]

Along a differentiable trajectory,

\[ \dot f= \nabla f(\widehat{\boldsymbol{\theta}})^{T} \dot{\widehat{\boldsymbol{\theta}}}. \]

When projection is active, substitution gives

\[ \dot f= \nabla f(\widehat{\boldsymbol{\theta}})^{T}\mathbf{y} \left(1-f(\widehat{\boldsymbol{\theta}})\right). \]

At the outer boundary \( f=1 \), every outward direction gives \( \dot f=0 \). Every inward or tangential direction gives \( \dot f=\nabla f^{T}\mathbf{y}\leq 0 \). Hence no trajectory starting in \( \Omega_{1} \) can cross to \( f>1 \):

\[ \widehat{\boldsymbol{\theta}}(0)\in\Omega_{1} \quad\Longrightarrow\quad \widehat{\boldsymbol{\theta}}(t)\in\Omega_{1}, \qquad \forall t\geq 0. \]

5.7 Property 6 — Continuity

At \( f=0 \), the correction is multiplied by zero. At \( \mathbf{y}^{T}\nabla f=0 \), the normal correction is also zero. Therefore the two branches agree on both switching surfaces and the operator is continuous. It is not necessarily continuously differentiable across those surfaces; numerical solvers should therefore avoid assuming global smoothness of higher derivatives.

6.1 Spherical set with a smooth boundary layer

Suppose prior knowledge gives a preferred parameter bound \( \theta_{\max}>0 \), and choose a boundary-layer width \( \varepsilon>0 \). Define

\[ f(\boldsymbol{\theta})= \frac{\|\boldsymbol{\theta}\|_{2}^{2}-\theta_{\max}^{2}} {2\varepsilon\theta_{\max}+\varepsilon^{2}}. \]

Then

\[ f(\boldsymbol{\theta})=0 \quad\Longleftrightarrow\quad \|\boldsymbol{\theta}\|_{2}=\theta_{\max}, \]

\[ f(\boldsymbol{\theta})=1 \quad\Longleftrightarrow\quad \|\boldsymbol{\theta}\|_{2}=\theta_{\max}+\varepsilon. \]

The gradient is

\[ \nabla f(\boldsymbol{\theta})= \frac{2\boldsymbol{\theta}} {2\varepsilon\theta_{\max}+\varepsilon^{2}}. \]

Consequently, the outer invariant bound is explicit:

\[ \|\widehat{\boldsymbol{\theta}}(t)\|_{2} \leq \theta_{\max}+\varepsilon, \qquad \forall t\geq 0, \]

provided the initial estimate is inside that outer ball.

6.2 Ellipsoidal set

When different parameter combinations have different scales or physical units, use an ellipsoid. Let \( W=W^{T}>0 \), center \( \mathbf{c} \), and radii \( 0<r_{0}<r_{1} \). Define

\[ f(\boldsymbol{\theta})= \frac{(\boldsymbol{\theta}-\mathbf{c})^{T} W(\boldsymbol{\theta}-\mathbf{c})-r_{0}^{2}} {r_{1}^{2}-r_{0}^{2}}. \]

Its gradient is

\[ \nabla f(\boldsymbol{\theta})= \frac{2W(\boldsymbol{\theta}-\mathbf{c})} {r_{1}^{2}-r_{0}^{2}}. \]

The inner ellipsoid corresponds to \( f\leq 0 \), while the larger ellipsoid corresponds to \( f\leq 1 \). Choosing \( W \) after nondimensionalizing the parameters prevents one large-magnitude parameter from dominating the geometry.

Physical knowledge is often componentwise:

\[ \underline{\theta}_{i}\leq \widehat{\theta}_{i}\leq\overline{\theta}_{i}, \qquad i=1,\ldots,p. \]

A hard projected direction for the box is

\[ [\Pi_{\mathrm{box}}(\widehat{\boldsymbol{\theta}},\mathbf{y})]_{i} =\begin{cases} 0, & \widehat{\theta}_{i}=\overline{\theta}_{i} \;\wedge\;y_{i}>0,\\ 0, & \widehat{\theta}_{i}=\underline{\theta}_{i} \;\wedge\;y_{i}<0,\\ y_{i}, & \text{otherwise.} \end{cases} \]

This operator is easy to implement but is nonsmooth at active faces and corners. Elementwise clipping of the parameter after integration is the static nearest-point projection onto a box; it is not automatically the same mathematical object as projecting the update direction. For a Lyapunov-based controller, the selected projection rule must be the one used in the stability proof.

Adaptive laws frequently use a symmetric positive-definite gain matrix \( \Gamma=\Gamma^{T}>0 \). If \( \mathbf{y} \) denotes the pre-gain update direction, a metric-compatible projection is

\[ \operatorname{Proj}_{\Gamma} (\widehat{\boldsymbol{\theta}},\mathbf{y},f) =\begin{cases} \Gamma\mathbf{y} -\Gamma\nabla f \dfrac{\nabla f^{T}\Gamma\mathbf{y}} {\nabla f^{T}\Gamma\nabla f} f(\widehat{\boldsymbol{\theta}}), & f(\widehat{\boldsymbol{\theta}})>0 \;\wedge\; \mathbf{y}^{T}\Gamma\nabla f>0,\\[8pt] \Gamma\mathbf{y}, & \text{otherwise.} \end{cases} \]

Here \( \nabla f \) is evaluated at \( \widehat{\boldsymbol{\theta}} \). The corresponding Lyapunov-compatible property is

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

A common implementation error is to apply the scalar-gain formula after multiplying by \( \Gamma \) without changing the normal direction and denominator. Unless \( \Gamma \) is a scalar multiple of the identity, that generally does not match the metric in the Lyapunov function.

9.1 Initial-condition requirement

The forward-invariance theorem assumes \( \widehat{\boldsymbol{\theta}}(0)\in\Omega_{1} \). If an estimate is initialized outside the outer set, first apply a static nearest-point projection or select a valid initial estimate.

9.2 Gradient denominator

The denominator \( \|\nabla f\|_{2}^{2} \) must not be evaluated as zero in the active branch. Well-designed spherical and ellipsoidal constraints have nonzero gradients on their boundary layers. Implementations should nevertheless include a small numerical tolerance and fail explicitly if an invalid constraint function is supplied.

9.3 Discrete integration is only approximately invariant

The theorem concerns a continuous-time differential equation. Explicit Euler integration can overshoot the outer boundary by a small amount proportional to the step size. Reducing the step size, using a higher-order solver, or applying a final static feasibility correction can control numerical drift. Such a correction should be documented separately from the mathematical projection operator.

9.4 Scaling and units

Projection geometry should not compare quantities with incompatible units without normalization. Define dimensionless coordinates \( \boldsymbol{\eta}=S^{-1} (\boldsymbol{\theta}-\mathbf{c}) \), or use a physically meaningful ellipsoidal weight matrix. Signal normalization itself is developed in Lesson 4.

9.5 Useful libraries

The program implements both nearest-point projection onto a ball and the smooth direction projection. It then checks the Lyapunov inequality over random samples and integrates a prescribed projected direction to visualize the invariant outer set.

Chapter9_Lesson2.py

"""Chapter9_Lesson2.py

Projection-operator definitions and numerical property checks for
Chapter 9, Lesson 2 of an Adaptive Control course.

Dependencies:
    numpy
    matplotlib
"""

from __future__ import annotations

from dataclasses import dataclass
from typing import Tuple

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class BallProjectionConfig:
    """Inner radius and boundary-layer thickness for smooth projection."""

    theta_max: float = 1.0
    epsilon: float = 0.2

    def __post_init__(self) -> None:
        if self.theta_max <= 0.0:
            raise ValueError("theta_max must be positive")
        if self.epsilon <= 0.0:
            raise ValueError("epsilon must be positive")

    @property
    def outer_radius(self) -> float:
        return self.theta_max + self.epsilon

    @property
    def denominator(self) -> float:
        return 2.0 * self.epsilon * self.theta_max + self.epsilon**2


def convex_function(theta: np.ndarray, cfg: BallProjectionConfig) -> float:
    """Return f(theta): f=0 at theta_max and f=1 at theta_max+epsilon."""
    theta = np.asarray(theta, dtype=float)
    return float((theta @ theta - cfg.theta_max**2) / cfg.denominator)


def convex_gradient(theta: np.ndarray, cfg: BallProjectionConfig) -> np.ndarray:
    """Gradient of the convex boundary function f."""
    theta = np.asarray(theta, dtype=float)
    return 2.0 * theta / cfg.denominator


def smooth_projection(
    theta: np.ndarray,
    y: np.ndarray,
    cfg: BallProjectionConfig,
    tolerance: float = 1.0e-14,
) -> np.ndarray:
    """Project an unconstrained update direction y at parameter estimate theta.

    The correction is active only in the boundary layer f(theta)>0 and only
    when y has an outward component along grad f(theta).
    """
    theta = np.asarray(theta, dtype=float)
    y = np.asarray(y, dtype=float)
    if theta.shape != y.shape:
        raise ValueError("theta and y must have the same shape")

    f_value = convex_function(theta, cfg)
    gradient = convex_gradient(theta, cfg)
    gradient_norm_sq = float(gradient @ gradient)
    outward_measure = float(y @ gradient)

    if f_value > 0.0 and outward_measure > 0.0:
        if gradient_norm_sq <= tolerance:
            raise FloatingPointError("active projection encountered a zero gradient")
        return y - gradient * (outward_measure / gradient_norm_sq) * f_value
    return y.copy()


def euclidean_project_ball(z: np.ndarray, radius: float) -> np.ndarray:
    """Nearest-point projection of z onto a closed Euclidean ball."""
    z = np.asarray(z, dtype=float)
    norm_z = float(np.linalg.norm(z))
    if norm_z <= radius or norm_z == 0.0:
        return z.copy()
    return (radius / norm_z) * z


def verify_lyapunov_inequality(
    cfg: BallProjectionConfig,
    samples: int = 5000,
    seed: int = 7,
) -> Tuple[float, int]:
    """Numerically check (theta-theta_star)^T(Proj(theta,y)-y) <= 0."""
    rng = np.random.default_rng(seed)
    worst_left_side = -np.inf
    violations = 0

    for _ in range(samples):
        theta_direction = rng.normal(size=2)
        theta_direction /= np.linalg.norm(theta_direction)
        theta = theta_direction * rng.uniform(0.0, cfg.outer_radius)

        star_direction = rng.normal(size=2)
        star_direction /= np.linalg.norm(star_direction)
        theta_star = star_direction * rng.uniform(0.0, cfg.theta_max)

        y = rng.normal(size=2)
        projected_y = smooth_projection(theta, y, cfg)
        left_side = float((theta - theta_star) @ (projected_y - y))
        worst_left_side = max(worst_left_side, left_side)
        if left_side > 1.0e-11:
            violations += 1

    return worst_left_side, violations


def simulate_projected_direction(
    cfg: BallProjectionConfig,
    dt: float = 0.001,
    final_time: float = 16.0,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Integrate theta_dot = Proj(theta,y(t)) for a prescribed direction signal."""
    steps = int(round(final_time / dt))
    time = np.linspace(0.0, final_time, steps + 1)
    theta_history = np.zeros((steps + 1, 2))
    raw_direction_history = np.zeros((steps + 1, 2))
    theta_history[0] = np.array([0.15, -0.10])

    for k, t in enumerate(time[:-1]):
        if t < 7.0:
            y = np.array([0.55, 0.30])
        elif t < 10.0:
            y = 0.18 * np.array([np.cos(2.0 * t), np.sin(2.0 * t)])
        else:
            y = np.array([-0.65, -0.20])

        raw_direction_history[k] = y
        theta_dot = smooth_projection(theta_history[k], y, cfg)
        theta_history[k + 1] = theta_history[k] + dt * theta_dot

    raw_direction_history[-1] = raw_direction_history[-2]
    return time, theta_history, raw_direction_history


def plot_results(
    time: np.ndarray,
    theta_history: np.ndarray,
    cfg: BallProjectionConfig,
) -> None:
    """Create trajectory and parameter-norm figures."""
    angle = np.linspace(0.0, 2.0 * np.pi, 500)

    plt.figure(figsize=(7.0, 6.0))
    plt.plot(
        cfg.theta_max * np.cos(angle),
        cfg.theta_max * np.sin(angle),
        "--",
        label="inner set boundary: f=0",
    )
    plt.plot(
        cfg.outer_radius * np.cos(angle),
        cfg.outer_radius * np.sin(angle),
        "-",
        label="outer invariant boundary: f=1",
    )
    plt.plot(theta_history[:, 0], theta_history[:, 1], label="parameter trajectory")
    plt.xlabel("theta_1")
    plt.ylabel("theta_2")
    plt.title("Smooth adaptive projection in parameter space")
    plt.axis("equal")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure(figsize=(7.0, 4.5))
    norms = np.linalg.norm(theta_history, axis=1)
    plt.plot(time, norms, label="||theta(t)||")
    plt.axhline(cfg.theta_max, linestyle="--", label="theta_max")
    plt.axhline(cfg.outer_radius, linestyle="-", label="theta_max + epsilon")
    plt.xlabel("time")
    plt.ylabel("parameter norm")
    plt.title("Forward-invariance numerical check")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


def main() -> None:
    cfg = BallProjectionConfig(theta_max=1.0, epsilon=0.2)

    # Static nearest-point projection example.
    z = np.array([1.4, -0.8])
    print("Euclidean projection:", euclidean_project_ball(z, cfg.theta_max))

    # Direction-projection examples.
    theta = np.array([1.05, 0.20])
    outward_y = np.array([0.8, 0.3])
    inward_y = np.array([-0.8, -0.3])
    print("f(theta):", convex_function(theta, cfg))
    print("outward direction before/after:", outward_y, smooth_projection(theta, outward_y, cfg))
    print("inward direction before/after:", inward_y, smooth_projection(theta, inward_y, cfg))

    worst_left_side, violations = verify_lyapunov_inequality(cfg)
    print("largest sampled Lyapunov left side:", worst_left_side)
    print("sampled inequality violations:", violations)

    time, theta_history, _ = simulate_projected_direction(cfg)
    maximum_norm = float(np.max(np.linalg.norm(theta_history, axis=1)))
    print("maximum simulated norm:", maximum_norm)
    print("outer radius:", cfg.outer_radius)
    plot_results(time, theta_history, cfg)


if __name__ == "__main__":
    main()

This portable C++17 version uses std::vector so that the projection logic is visible without an external dependency. In a larger controller, the same operations can be expressed with Eigen vectors and matrices.

Chapter9_Lesson2.cpp

// Chapter9_Lesson2.cpp
// Smooth projection-operator definitions and numerical property checks.
// Build: g++ -std=c++17 -O2 Chapter9_Lesson2.cpp -o Chapter9_Lesson2

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <random>
#include <stdexcept>
#include <vector>

using Vector = std::vector<double>;

struct BallProjectionConfig {
    double thetaMax{1.0};
    double epsilon{0.2};

    double outerRadius() const { return thetaMax + epsilon; }
    double denominator() const {
        return 2.0 * epsilon * thetaMax + epsilon * epsilon;
    }
    void validate() const {
        if (thetaMax <= 0.0 || epsilon <= 0.0) {
            throw std::invalid_argument("thetaMax and epsilon must be positive");
        }
    }
};

double dot(const Vector& a, const Vector& b) {
    if (a.size() != b.size()) {
        throw std::invalid_argument("vector dimensions must agree");
    }
    double value = 0.0;
    for (std::size_t i = 0; i < a.size(); ++i) {
        value += a[i] * b[i];
    }
    return value;
}

double norm(const Vector& a) { return std::sqrt(dot(a, a)); }

Vector add(const Vector& a, const Vector& b) {
    if (a.size() != b.size()) {
        throw std::invalid_argument("vector dimensions must agree");
    }
    Vector result(a.size());
    for (std::size_t i = 0; i < a.size(); ++i) {
        result[i] = a[i] + b[i];
    }
    return result;
}

Vector subtract(const Vector& a, const Vector& b) {
    if (a.size() != b.size()) {
        throw std::invalid_argument("vector dimensions must agree");
    }
    Vector result(a.size());
    for (std::size_t i = 0; i < a.size(); ++i) {
        result[i] = a[i] - b[i];
    }
    return result;
}

Vector scale(const Vector& a, double scalar) {
    Vector result(a.size());
    for (std::size_t i = 0; i < a.size(); ++i) {
        result[i] = scalar * a[i];
    }
    return result;
}

double convexFunction(const Vector& theta, const BallProjectionConfig& cfg) {
    return (dot(theta, theta) - cfg.thetaMax * cfg.thetaMax) / cfg.denominator();
}

Vector convexGradient(const Vector& theta, const BallProjectionConfig& cfg) {
    return scale(theta, 2.0 / cfg.denominator());
}

Vector smoothProjection(
    const Vector& theta,
    const Vector& y,
    const BallProjectionConfig& cfg,
    double tolerance = 1.0e-14) {
    if (theta.size() != y.size()) {
        throw std::invalid_argument("theta and y must have the same dimension");
    }

    const double fValue = convexFunction(theta, cfg);
    const Vector gradient = convexGradient(theta, cfg);
    const double gradientNormSquared = dot(gradient, gradient);
    const double outwardMeasure = dot(y, gradient);

    if (fValue > 0.0 && outwardMeasure > 0.0) {
        if (gradientNormSquared <= tolerance) {
            throw std::runtime_error("active projection encountered a zero gradient");
        }
        const double multiplier = outwardMeasure / gradientNormSquared * fValue;
        return subtract(y, scale(gradient, multiplier));
    }
    return y;
}

Vector euclideanProjectBall(const Vector& z, double radius) {
    const double zNorm = norm(z);
    if (zNorm <= radius || zNorm == 0.0) {
        return z;
    }
    return scale(z, radius / zNorm);
}

void printVector(const Vector& vector) {
    std::cout << "[";
    for (std::size_t i = 0; i < vector.size(); ++i) {
        std::cout << vector[i] << (i + 1 == vector.size() ? "" : ", ");
    }
    std::cout << "]";
}

int main() {
    try {
        BallProjectionConfig cfg;
        cfg.validate();

        const Vector z{1.4, -0.8};
        std::cout << std::fixed << std::setprecision(8);
        std::cout << "Euclidean projection: ";
        printVector(euclideanProjectBall(z, cfg.thetaMax));
        std::cout << "\n";

        const Vector theta{1.05, 0.20};
        const Vector outwardY{0.8, 0.3};
        const Vector inwardY{-0.8, -0.3};
        std::cout << "f(theta): " << convexFunction(theta, cfg) << "\n";
        std::cout << "outward direction after projection: ";
        printVector(smoothProjection(theta, outwardY, cfg));
        std::cout << "\n";
        std::cout << "inward direction after projection: ";
        printVector(smoothProjection(theta, inwardY, cfg));
        std::cout << "\n";

        std::mt19937 generator(7U);
        std::normal_distribution<double> normal(0.0, 1.0);
        std::uniform_real_distribution<double> unit(0.0, 1.0);
        double worstLeftSide = -1.0e100;
        int violations = 0;

        for (int sample = 0; sample < 5000; ++sample) {
            Vector direction{normal(generator), normal(generator)};
            direction = scale(direction, 1.0 / norm(direction));
            const Vector sampledTheta = scale(direction, unit(generator) * cfg.outerRadius());

            Vector starDirection{normal(generator), normal(generator)};
            starDirection = scale(starDirection, 1.0 / norm(starDirection));
            const Vector thetaStar = scale(starDirection, unit(generator) * cfg.thetaMax);

            const Vector y{normal(generator), normal(generator)};
            const Vector projectedY = smoothProjection(sampledTheta, y, cfg);
            const double leftSide = dot(subtract(sampledTheta, thetaStar), subtract(projectedY, y));
            worstLeftSide = std::max(worstLeftSide, leftSide);
            if (leftSide > 1.0e-11) {
                ++violations;
            }
        }

        std::cout << "largest sampled Lyapunov left side: " << worstLeftSide << "\n";
        std::cout << "sampled inequality violations: " << violations << "\n";

        const double dt = 0.001;
        const double finalTime = 16.0;
        const int steps = static_cast<int>(std::lround(finalTime / dt));
        Vector simulatedTheta{0.15, -0.10};
        double maximumNorm = norm(simulatedTheta);

        for (int k = 0; k < steps; ++k) {
            const double time = k * dt;
            Vector y(2);
            if (time < 7.0) {
                y = {0.55, 0.30};
            } else if (time < 10.0) {
                y = {0.18 * std::cos(2.0 * time), 0.18 * std::sin(2.0 * time)};
            } else {
                y = {-0.65, -0.20};
            }
            const Vector thetaDot = smoothProjection(simulatedTheta, y, cfg);
            simulatedTheta = add(simulatedTheta, scale(thetaDot, dt));
            maximumNorm = std::max(maximumNorm, norm(simulatedTheta));
        }

        std::cout << "maximum simulated norm: " << maximumNorm << "\n";
        std::cout << "outer radius: " << cfg.outerRadius() << "\n";
        return violations == 0 ? 0 : 1;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << "\n";
        return 2;
    }
}

The Java implementation uses primitive arrays and explicit checks. EJML or Apache Commons Math can replace these helpers when the parameter dimension or matrix structure becomes larger.

Chapter9_Lesson2.java

// Chapter9_Lesson2.java
// Smooth projection-operator definitions and numerical property checks.
// Build and run: javac Chapter9_Lesson2.java && java Chapter9_Lesson2

import java.util.Arrays;
import java.util.Random;

public final class Chapter9_Lesson2 {
    private Chapter9_Lesson2() {}

    private static final class BallProjectionConfig {
        final double thetaMax;
        final double epsilon;

        BallProjectionConfig(double thetaMax, double epsilon) {
            if (thetaMax <= 0.0 || epsilon <= 0.0) {
                throw new IllegalArgumentException("thetaMax and epsilon must be positive");
            }
            this.thetaMax = thetaMax;
            this.epsilon = epsilon;
        }

        double outerRadius() {
            return thetaMax + epsilon;
        }

        double denominator() {
            return 2.0 * epsilon * thetaMax + epsilon * epsilon;
        }
    }

    private static double dot(double[] a, double[] b) {
        requireSameLength(a, b);
        double value = 0.0;
        for (int i = 0; i < a.length; i++) {
            value += a[i] * b[i];
        }
        return value;
    }

    private static double norm(double[] a) {
        return Math.sqrt(dot(a, a));
    }

    private static double[] add(double[] a, double[] b) {
        requireSameLength(a, b);
        double[] result = new double[a.length];
        for (int i = 0; i < a.length; i++) {
            result[i] = a[i] + b[i];
        }
        return result;
    }

    private static double[] subtract(double[] a, double[] b) {
        requireSameLength(a, b);
        double[] result = new double[a.length];
        for (int i = 0; i < a.length; i++) {
            result[i] = a[i] - b[i];
        }
        return result;
    }

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

    private static void requireSameLength(double[] a, double[] b) {
        if (a.length != b.length) {
            throw new IllegalArgumentException("vector dimensions must agree");
        }
    }

    private static double convexFunction(double[] theta, BallProjectionConfig cfg) {
        return (dot(theta, theta) - cfg.thetaMax * cfg.thetaMax) / cfg.denominator();
    }

    private static double[] convexGradient(double[] theta, BallProjectionConfig cfg) {
        return scale(theta, 2.0 / cfg.denominator());
    }

    private static double[] smoothProjection(
            double[] theta, double[] y, BallProjectionConfig cfg) {
        requireSameLength(theta, y);
        double fValue = convexFunction(theta, cfg);
        double[] gradient = convexGradient(theta, cfg);
        double gradientNormSquared = dot(gradient, gradient);
        double outwardMeasure = dot(y, gradient);

        if (fValue > 0.0 && outwardMeasure > 0.0) {
            if (gradientNormSquared <= 1.0e-14) {
                throw new ArithmeticException("active projection encountered a zero gradient");
            }
            double multiplier = outwardMeasure / gradientNormSquared * fValue;
            return subtract(y, scale(gradient, multiplier));
        }
        return Arrays.copyOf(y, y.length);
    }

    private static double[] euclideanProjectBall(double[] z, double radius) {
        double zNorm = norm(z);
        if (zNorm <= radius || zNorm == 0.0) {
            return Arrays.copyOf(z, z.length);
        }
        return scale(z, radius / zNorm);
    }

    public static void main(String[] args) {
        BallProjectionConfig cfg = new BallProjectionConfig(1.0, 0.2);

        double[] z = {1.4, -0.8};
        System.out.println("Euclidean projection: " + Arrays.toString(euclideanProjectBall(z, cfg.thetaMax)));

        double[] theta = {1.05, 0.20};
        double[] outwardY = {0.8, 0.3};
        double[] inwardY = {-0.8, -0.3};
        System.out.printf("f(theta): %.8f%n", convexFunction(theta, cfg));
        System.out.println("outward direction after projection: "
                + Arrays.toString(smoothProjection(theta, outwardY, cfg)));
        System.out.println("inward direction after projection: "
                + Arrays.toString(smoothProjection(theta, inwardY, cfg)));

        Random random = new Random(7L);
        double worstLeftSide = -Double.MAX_VALUE;
        int violations = 0;

        for (int sample = 0; sample < 5000; sample++) {
            double[] direction = {random.nextGaussian(), random.nextGaussian()};
            direction = scale(direction, 1.0 / norm(direction));
            double[] sampledTheta = scale(direction, random.nextDouble() * cfg.outerRadius());

            double[] starDirection = {random.nextGaussian(), random.nextGaussian()};
            starDirection = scale(starDirection, 1.0 / norm(starDirection));
            double[] thetaStar = scale(starDirection, random.nextDouble() * cfg.thetaMax);

            double[] y = {random.nextGaussian(), random.nextGaussian()};
            double[] projectedY = smoothProjection(sampledTheta, y, cfg);
            double leftSide = dot(subtract(sampledTheta, thetaStar), subtract(projectedY, y));
            worstLeftSide = Math.max(worstLeftSide, leftSide);
            if (leftSide > 1.0e-11) {
                violations++;
            }
        }

        System.out.printf("largest sampled Lyapunov left side: %.12e%n", worstLeftSide);
        System.out.println("sampled inequality violations: " + violations);

        double dt = 0.001;
        double finalTime = 16.0;
        int steps = (int) Math.round(finalTime / dt);
        double[] simulatedTheta = {0.15, -0.10};
        double maximumNorm = norm(simulatedTheta);

        for (int k = 0; k < steps; k++) {
            double time = k * dt;
            double[] y;
            if (time < 7.0) {
                y = new double[] {0.55, 0.30};
            } else if (time < 10.0) {
                y = new double[] {0.18 * Math.cos(2.0 * time), 0.18 * Math.sin(2.0 * time)};
            } else {
                y = new double[] {-0.65, -0.20};
            }
            double[] thetaDot = smoothProjection(simulatedTheta, y, cfg);
            simulatedTheta = add(simulatedTheta, scale(thetaDot, dt));
            maximumNorm = Math.max(maximumNorm, norm(simulatedTheta));
        }

        System.out.printf("maximum simulated norm: %.8f%n", maximumNorm);
        System.out.printf("outer radius: %.8f%n", cfg.outerRadius());

        if (violations != 0) {
            throw new IllegalStateException("Lyapunov inequality test failed");
        }
    }
}

The MATLAB script contains local functions for the convex boundary function, gradient, smooth direction projection, and static ball projection. For Simulink, place the smoothProjection logic in a MATLAB Function block and integrate its output with an Integrator block. The complete adaptive-loop interconnection is reserved for Lesson 3.

Chapter9_Lesson2.m

% Chapter9_Lesson2.m
% Projection-operator definitions and numerical property checks.
% No specialized toolbox is required.

clear; clc; close all;
rng(7);

thetaMax = 1.0;
epsilon = 0.2;
outerRadius = thetaMax + epsilon;

% Static nearest-point projection.
z = [1.4; -0.8];
fprintf('Euclidean projection: [% .8f, % .8f]\n', ...
    euclideanProjectBall(z, thetaMax));

% Direction-projection examples.
theta = [1.05; 0.20];
outwardY = [0.8; 0.3];
inwardY = [-0.8; -0.3];
fprintf('f(theta): %.8f\n', convexFunction(theta, thetaMax, epsilon));
disp('outward direction after projection:');
disp(smoothProjection(theta, outwardY, thetaMax, epsilon));
disp('inward direction after projection:');
disp(smoothProjection(theta, inwardY, thetaMax, epsilon));

% Numerical verification of the Lyapunov inequality.
worstLeftSide = -inf;
violations = 0;
for sample = 1:5000
    direction = randn(2, 1);
    direction = direction / norm(direction);
    sampledTheta = direction * (rand() * outerRadius);

    starDirection = randn(2, 1);
    starDirection = starDirection / norm(starDirection);
    thetaStar = starDirection * (rand() * thetaMax);

    y = randn(2, 1);
    projectedY = smoothProjection(sampledTheta, y, thetaMax, epsilon);
    leftSide = (sampledTheta - thetaStar).' * (projectedY - y);
    worstLeftSide = max(worstLeftSide, leftSide);
    if leftSide > 1.0e-11
        violations = violations + 1;
    end
end
fprintf('largest sampled Lyapunov left side: %.12e\n', worstLeftSide);
fprintf('sampled inequality violations: %d\n', violations);

% Integrate theta_dot = Proj(theta, y(t)).
dt = 0.001;
finalTime = 16.0;
time = 0:dt:finalTime;
thetaHistory = zeros(2, numel(time));
thetaHistory(:, 1) = [0.15; -0.10];

for k = 1:(numel(time) - 1)
    t = time(k);
    if t < 7.0
        y = [0.55; 0.30];
    elseif t < 10.0
        y = 0.18 * [cos(2.0 * t); sin(2.0 * t)];
    else
        y = [-0.65; -0.20];
    end
    thetaDot = smoothProjection(thetaHistory(:, k), y, thetaMax, epsilon);
    thetaHistory(:, k + 1) = thetaHistory(:, k) + dt * thetaDot;
end

parameterNorm = vecnorm(thetaHistory, 2, 1);
fprintf('maximum simulated norm: %.8f\n', max(parameterNorm));
fprintf('outer radius: %.8f\n', outerRadius);

angle = linspace(0, 2*pi, 500);
figure;
plot(thetaMax*cos(angle), thetaMax*sin(angle), '--', 'LineWidth', 1.2);
hold on;
plot(outerRadius*cos(angle), outerRadius*sin(angle), '-', 'LineWidth', 1.2);
plot(thetaHistory(1, :), thetaHistory(2, :), 'LineWidth', 1.2);
axis equal; grid on;
xlabel('\theta_1'); ylabel('\theta_2');
title('Smooth adaptive projection in parameter space');
legend('inner set boundary: f=0', 'outer invariant boundary: f=1', ...
    'parameter trajectory', 'Location', 'best');

figure;
plot(time, parameterNorm, 'LineWidth', 1.2);
hold on;
yline(thetaMax, '--', '\theta_{max}');
yline(outerRadius, '-', '\theta_{max}+\epsilon');
grid on;
xlabel('time'); ylabel('||\theta(t)||_2');
title('Forward-invariance numerical check');

% Simulink use:
% Place smoothProjection in a MATLAB Function block and feed it theta and y.
% Keep thetaMax and epsilon as block parameters. Integrate the returned
% projected direction with an Integrator block. Lesson 3 develops the full
% continuous-time adaptive-law interconnection.

function value = convexFunction(theta, thetaMax, epsilon)
    denominator = 2.0 * epsilon * thetaMax + epsilon^2;
    value = (theta.' * theta - thetaMax^2) / denominator;
end

function gradient = convexGradient(theta, thetaMax, epsilon)
    denominator = 2.0 * epsilon * thetaMax + epsilon^2;
    gradient = 2.0 * theta / denominator;
end

function projectedY = smoothProjection(theta, y, thetaMax, epsilon)
    if ~isequal(size(theta), size(y))
        error('theta and y must have the same dimensions.');
    end

    fValue = convexFunction(theta, thetaMax, epsilon);
    gradient = convexGradient(theta, thetaMax, epsilon);
    gradientNormSquared = gradient.' * gradient;
    outwardMeasure = y.' * gradient;

    if fValue > 0.0 && outwardMeasure > 0.0
        if gradientNormSquared <= 1.0e-14
            error('Active projection encountered a zero gradient.');
        end
        projectedY = y - gradient * ...
            (outwardMeasure / gradientNormSquared) * fValue;
    else
        projectedY = y;
    end
end

function projectedZ = euclideanProjectBall(z, radius)
    zNorm = norm(z);
    if zNorm <= radius || zNorm == 0.0
        projectedZ = z;
    else
        projectedZ = (radius / zNorm) * z;
    end
end

The notebook contains an executable input cell and a readable copy of the Wolfram Language source. It performs the same static projection, direction projection, random inequality test, numerical integration, and plotting.

Chapter9_Lesson2.nb

ClearAll["Global`*"];
thetaMax = 1.0; epsilon = 0.2; outerRadius = thetaMax + epsilon;
denominator = 2.0 epsilon thetaMax + epsilon^2;
convexFunction[theta_List] := (theta.theta - thetaMax^2)/denominator;
convexGradient[theta_List] := 2.0 theta/denominator;
smoothProjection[theta_List, y_List] := Module[
  {fValue, gradient, gradientNormSquared, outwardMeasure},
  If[Length[theta] != Length[y],
    Return[Failure["DimensionMismatch", <|"Message" -> "theta and y must have equal dimensions."|>]]
  ];
  fValue = convexFunction[theta]; gradient = convexGradient[theta];
  gradientNormSquared = gradient.gradient; outwardMeasure = y.gradient;
  If[fValue > 0.0 && outwardMeasure > 0.0,
    If[gradientNormSquared <= 1.0*^-14,
      Return[Failure["ZeroGradient", <|"Message" -> "Active projection encountered a zero gradient."|>]]
    ];
    y - gradient (outwardMeasure/gradientNormSquared) fValue, y
  ]
];
euclideanProjectBall[z_List, radius_?Positive] := Module[{zNorm = Norm[z]},
  If[zNorm <= radius || zNorm == 0.0, z, (radius/zNorm) z]
];
Print["Euclidean projection: ", euclideanProjectBall[{1.4, -0.8}, thetaMax]];
theta = {1.05, 0.20}; outwardY = {0.8, 0.3}; inwardY = {-0.8, -0.3};
Print["f(theta): ", N[convexFunction[theta], 10]];
Print["outward direction after projection: ", smoothProjection[theta, outwardY]];
Print["inward direction after projection: ", smoothProjection[theta, inwardY]];
SeedRandom[7];
samples = Table[
  Module[{direction, sampledTheta, starDirection, thetaStar, y, projectedY, leftSide},
    direction = Normalize[RandomVariate[NormalDistribution[0, 1], 2]];
    sampledTheta = direction RandomReal[{0, outerRadius}];
    starDirection = Normalize[RandomVariate[NormalDistribution[0, 1], 2]];
    thetaStar = starDirection RandomReal[{0, thetaMax}];
    y = RandomVariate[NormalDistribution[0, 1], 2];
    projectedY = smoothProjection[sampledTheta, y];
    leftSide = (sampledTheta - thetaStar).(projectedY - y);
    leftSide
  ], {5000}
];
Print["largest sampled Lyapunov left side: ", Max[samples]];
Print["sampled inequality violations: ", Count[samples, value_ /; value > 1.0*^-11]];
dt = 0.001; finalTime = 16.0; time = Range[0.0, finalTime, dt];
rawDirection[t_] := Which[
  t < 7.0, {0.55, 0.30},
  t < 10.0, 0.18 {Cos[2.0 t], Sin[2.0 t]},
  True, {-0.65, -0.20}
];
thetaHistory = FoldList[
  Function[{currentTheta, t}, currentTheta + dt smoothProjection[currentTheta, rawDirection[t]]],
  {0.15, -0.10}, Most[time]
];
parameterNorm = Norm /@ thetaHistory;
Print["maximum simulated norm: ", Max[parameterNorm]]; Print["outer radius: ", outerRadius];
trajectoryPlot = Show[
  ParametricPlot[thetaMax {Cos[a], Sin[a]}, {a, 0, 2 Pi}, PlotStyle -> Dashed],
  ParametricPlot[outerRadius {Cos[a], Sin[a]}, {a, 0, 2 Pi}],
  ListLinePlot[thetaHistory], AspectRatio -> 1,
  AxesLabel -> {"theta_1", "theta_2"},
  PlotLabel -> "Smooth adaptive projection in parameter space"
];
normPlot = ListLinePlot[
  Transpose[{time, parameterNorm}], Frame -> True,
  FrameLabel -> {"time", "||theta(t)||_2"},
  PlotLabel -> "Forward-invariance numerical check",
  Epilog -> {Dashed, InfiniteLine[{{0, thetaMax}, {1, thetaMax}}],
    InfiniteLine[{{0, outerRadius}, {1, outerRadius}}]}
];
Print[trajectoryPlot]; Print[normPlot];        

Problem 1 — Static projection onto a ball

Let \( \Omega=\{\boldsymbol{\theta}\in\mathbb{R}^{2}: \|\boldsymbol{\theta}\|_{2}\leq 1\} \) and \( \mathbf{z}=[1.4,-0.8]^{T} \). Compute \( P_{\Omega}(\mathbf{z}) \).

Solution:

\[ \|\mathbf{z}\|_{2}=\sqrt{1.4^{2}+(-0.8)^{2}} =\sqrt{2.6}\approx 1.61245>1. \]

Therefore the point is radially scaled to the boundary:

\[ P_{\Omega}(\mathbf{z})= \frac{\mathbf{z}}{\|\mathbf{z}\|_{2}} \approx \begin{bmatrix}0.86824\\-0.49614\end{bmatrix}. \]

Problem 2 — Evaluate a smooth projected direction

Use \( \theta_{\max}=1 \), \( \varepsilon=0.2 \), \( \widehat{\boldsymbol{\theta}}=[1.05,0.20]^{T} \), and \( \mathbf{y}=[0.8,0.3]^{T} \). Determine whether projection is active and calculate the projected direction.

Solution: The denominator of the boundary function is

\[ d=2\varepsilon\theta_{\max}+\varepsilon^{2} =2(0.2)(1)+(0.2)^{2}=0.44. \]

Thus

\[ f(\widehat{\boldsymbol{\theta}})= \frac{1.05^{2}+0.20^{2}-1}{0.44} \approx 0.323864>0, \]

\[ \nabla f= \frac{2}{0.44} \begin{bmatrix}1.05\\0.20\end{bmatrix} \approx \begin{bmatrix}4.77273\\0.90909\end{bmatrix}. \]

The outward test is positive:

\[ \mathbf{y}^{T}\nabla f \approx 0.8(4.77273)+0.3(0.90909) \approx 4.09091>0. \]

Therefore projection is active, and substitution gives

\[ \operatorname{Proj} (\widehat{\boldsymbol{\theta}},\mathbf{y},f) \approx \begin{bmatrix}0.53212\\0.24898\end{bmatrix}. \]

Only the outward normal component is reduced; a tangential component remains.

Problem 3 — Prove the Lyapunov inequality

Let \( f \) be convex and differentiable, \( \boldsymbol{\theta}^{*}\in\Omega_{0} \), and suppose projection is active. Prove

\[ (\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*})^{T} [\operatorname{Proj}(\widehat{\boldsymbol{\theta}},\mathbf{y},f)-\mathbf{y}] \leq 0. \]

Solution: Convexity implies

\[ f(\boldsymbol{\theta}^{*})\geq f(\widehat{\boldsymbol{\theta}})+ \nabla f(\widehat{\boldsymbol{\theta}})^{T} (\boldsymbol{\theta}^{*}-\widehat{\boldsymbol{\theta}}). \]

Hence

\[ (\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*})^{T} \nabla f(\widehat{\boldsymbol{\theta}}) \geq f(\widehat{\boldsymbol{\theta}})-f(\boldsymbol{\theta}^{*}) >0. \]

In the active branch,

\[ \operatorname{Proj}-\mathbf{y} =-\nabla f \frac{\nabla f^{T}\mathbf{y}}{\|\nabla f\|_{2}^{2}}f. \]

The scalar multiplying \( -\nabla f \) is positive, while the inner product of parameter error with \( \nabla f \) is positive. Their product is therefore nonpositive.

Problem 4 — Boundary tangency and invariance

Show that at \( f(\widehat{\boldsymbol{\theta}})=1 \), an outward raw direction produces a projected direction tangent to the boundary.

Solution: In the active branch,

\[ \nabla f^{T}\operatorname{Proj} =\nabla f^{T}\mathbf{y} -\nabla f^{T}\nabla f \frac{\nabla f^{T}\mathbf{y}}{\|\nabla f\|_{2}^{2}}f. \]

At \( f=1 \),

\[ \nabla f^{T}\operatorname{Proj} =\nabla f^{T}\mathbf{y}-\nabla f^{T}\mathbf{y}=0. \]

Therefore the projected velocity lies in the tangent hyperplane. Since an inward direction already gives \( \dot f\leq 0 \), no direction at the boundary can increase \( f \) beyond one.

Problem 5 — Derive the ellipsoidal gradient

For

\[ f(\boldsymbol{\theta})= \frac{(\boldsymbol{\theta}-\mathbf{c})^{T} W(\boldsymbol{\theta}-\mathbf{c})-r_{0}^{2}} {r_{1}^{2}-r_{0}^{2}}, \qquad W=W^{T}>0, \]

derive \( \nabla f \).

Solution: For symmetric \( W \),

\[ \nabla_{\boldsymbol{\theta}} (\boldsymbol{\theta}-\mathbf{c})^{T} W(\boldsymbol{\theta}-\mathbf{c}) =2W(\boldsymbol{\theta}-\mathbf{c}). \]

The constant denominator therefore gives

\[ \nabla f(\boldsymbol{\theta})= \frac{2W(\boldsymbol{\theta}-\mathbf{c})} {r_{1}^{2}-r_{0}^{2}}. \]

Problem 6 — Static clipping versus direction projection

A programmer integrates an unconstrained update and then applies theta = clip(theta, lower, upper). Explain why this is not a complete substitute for the projected direction used in a Lyapunov proof.

Solution:

Elementwise clipping computes a static nearest point in a box after a finite numerical step. The stability proof, however, usually differentiates a Lyapunov function along a continuous-time update law and relies on a specific inequality involving the projected velocity. Finite-step clipping can be a useful feasibility safeguard, but its effect on the discrete-time Lyapunov difference must be analyzed separately. The two operations coincide only in special limiting cases and should not be treated as mathematically identical.

Projection in adaptive control modifies an update direction, not merely an already computed parameter point. A convex boundary function defines an unmodified inner set, an invariant outer set, and a smooth transition layer. The operator removes only an outward normal component, preserves inward and tangential motion, does not increase the update norm, satisfies a key nonpositive parameter-error inequality, and makes the outer admissible set forward invariant. Spherical, ellipsoidal, box, and gain-weighted versions provide different geometries for physical parameter knowledge. Lesson 3 uses these definitions inside continuous-time adaptive laws and complete Lyapunov proofs.