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

• define a closed convex admissible parameter set and its Euclidean projection;
• distinguish projection of a point from projection of a continuous-time update direction;
• construct scalar, box, ball, and smooth-boundary projection formulas;
• prove that projection preserves feasibility and does not add an adverse parameter-error term to a Lyapunov derivative;
• implement numerically safe projected updates in continuous- and discrete-time computations.

In the previous lesson, a generic parameter estimate \( \hat{\boldsymbol{\theta} } \in \mathbb{R}^{p} \) was moved in a descent direction. Here the unconstrained update is denoted by \( \mathbf{v}(t) \):

\[ \dot{\hat{\boldsymbol{\theta} } }=\mathbf{v}(t). \]

The vector \( \mathbf{v} \) may be a negative gradient, a Lyapunov-designed direction, or any other candidate update. Projection modifies only those components that would drive the estimate out of the admissible set.

Physical knowledge often supplies bounds before adaptation begins. A mass is positive, an aerodynamic coefficient may lie in a certified interval, and a normalized efficiency factor may be known to remain between two limits. If this information is ignored, a transient gradient direction can generate implausible estimates even when the state response remains bounded.

Projection serves three distinct purposes:

1. Feasibility: the estimate remains physically meaningful.
2. Analytical boundedness: the parameter estimate belongs to a prescribed compact set.
3. Numerical protection: finite precision and large temporary update directions cannot create arbitrarily large estimates.

Projection is not the same as actuator saturation. Actuator saturation constrains a plant input, while parameter projection constrains an internal estimate. The two operations have different stability implications.

Let \( \Omega \subset \mathbb{R}^{p} \) denote the set of allowable estimates. The standard projection theory assumes that \( \Omega \) is nonempty, closed, and convex. Convexity means

\[ \boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2}\in\Omega,\quad \lambda\in[0,1] \quad\Longrightarrow\quad \lambda\boldsymbol{\theta}_{1}+(1-\lambda)\boldsymbol{\theta}_{2}\in\Omega. \]

Common choices are:

\[ \begin{aligned} \text{interval: }&\quad \Omega=[\ell,u],\\ \text{box: }&\quad \Omega=\{\boldsymbol{\theta}:\boldsymbol{\ell}\le \boldsymbol{\theta}\le\mathbf{u}\},\\ \text{ball: }&\quad \Omega=\{\boldsymbol{\theta}:\|\boldsymbol{\theta}-\mathbf{c}\|_{2}\le R\},\\ \text{ellipsoid: }&\quad \Omega=\{\boldsymbol{\theta}: (\boldsymbol{\theta}-\mathbf{c})^{\top}\mathbf{Q} (\boldsymbol{\theta}-\mathbf{c})\le 1\},\quad \mathbf{Q}=\mathbf{Q}^{\top}>0. \end{aligned} \]

Compactness is obtained when the set is also bounded. For Lyapunov analysis, one normally assumes that the unknown true parameter satisfies \( \boldsymbol{\theta}^{*}\in\Omega \). If this assumption is false, feasibility of the estimate can still be enforced, but the usual projection inequality relative to the true parameter is no longer available.

For any \( \mathbf{x}\in\mathbb{R}^{p} \), the Euclidean projection onto \( \Omega \) is the closest feasible point:

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

Closedness gives existence, and convexity together with strict convexity of the squared norm gives uniqueness. Let \( \mathbf{p}=\Pi_{\Omega}(\mathbf{x}) \). For any \( \mathbf{z}\in\Omega \), the segment \( \mathbf{p}+s(\mathbf{z}-\mathbf{p}) \) is feasible for \( s\in[0,1] \). Define

\[ h(s)=\frac{1}{2}\|\mathbf{p}+s(\mathbf{z}-\mathbf{p})-\mathbf{x}\|_{2}^{2}. \]

Because \( s=0 \) is a minimum along this segment,

\[ h'(0^{+})=(\mathbf{p}-\mathbf{x})^{\top}(\mathbf{z}-\mathbf{p})\ge 0, \]

or equivalently the fundamental variational inequality

\[ (\mathbf{x}-\mathbf{p})^{\top}(\mathbf{z}-\mathbf{p})\le 0, \qquad \forall\mathbf{z}\in\Omega. \]

Taking \( \mathbf{z}=\boldsymbol{\theta}^{*}\in\Omega \) and expanding norms gives

\[ \|\mathbf{p}-\boldsymbol{\theta}^{*}\|_{2}^{2} \le \|\mathbf{x}-\boldsymbol{\theta}^{*}\|_{2}^{2} -\|\mathbf{x}-\mathbf{p}\|_{2}^{2} \le \|\mathbf{x}-\boldsymbol{\theta}^{*}\|_{2}^{2}. \]

Thus a point projection cannot increase the distance to any feasible true parameter. This fact is central in discrete-time projected adaptation.

In continuous time, projecting the point after it has already left the set is conceptually weaker than preventing an outward velocity. At a feasible point \( \hat{\boldsymbol{\theta} }\in\Omega \), the tangent cone contains all instantaneous feasible directions:

\[ T_{\Omega}(\hat{\boldsymbol{\theta} }) =\operatorname{cl}\left\{\alpha(\mathbf{z}-\hat{\boldsymbol{\theta} }): \alpha\ge 0,\;\mathbf{z}\in\Omega\right\}. \]

The corresponding normal cone is

\[ N_{\Omega}(\hat{\boldsymbol{\theta} }) =\left\{\mathbf{n}:\mathbf{n}^{\top} (\mathbf{z}-\hat{\boldsymbol{\theta} })\le 0, \;\forall\mathbf{z}\in\Omega\right\}. \]

The projected update direction is the closest vector in the tangent cone:

\[ \operatorname{Proj}_{\Omega} (\hat{\boldsymbol{\theta} },\mathbf{v}) =\underset{\mathbf{w}\in T_{\Omega}(\hat{\boldsymbol{\theta} })} {\operatorname{arg\,min} }\; \frac{1}{2}\|\mathbf{w}-\mathbf{v}\|_{2}^{2}. \]

flowchart TD
  A["Current estimate inside admissible set"] --> B["Compute raw update direction v"]
  B --> C{"Would v leave the set?"}
  C -->|"No"| D["Use the raw direction"]
  C -->|"Yes"| E["Remove the outward normal component"]
  D --> F["Feasible parameter motion"]
  E --> F
        

Geometrically, the raw direction decomposes into a feasible tangential part and an outward normal part. Projection keeps the former and removes the latter. In the interior, the tangent cone is all of \( \mathbb{R}^{p} \), so the projection leaves the update unchanged.

Suppose the admissible set is described by a differentiable convex function:

\[ \Omega=\{\boldsymbol{\theta}:f(\boldsymbol{\theta})\le 0\}. \]

At a regular boundary point where \( f(\hat{\boldsymbol{\theta} })=0 \) and \( \nabla f(\hat{\boldsymbol{\theta} })\ne\mathbf{0} \), the gradient is an outward normal. Therefore

\[ \operatorname{Proj}_{\Omega}(\hat{\boldsymbol{\theta} },\mathbf{v}) =\begin{cases} \mathbf{v}, & f(\hat{\boldsymbol{\theta} })<0,\\ \mathbf{v}, & f(\hat{\boldsymbol{\theta} })=0\text{ and } \nabla f^{\top}\mathbf{v}\le 0,\\ \mathbf{v}-\dfrac{\nabla f\,\nabla f^{\top} } {\|\nabla f\|_{2}^{2} }\mathbf{v}, & f(\hat{\boldsymbol{\theta} })=0\text{ and } \nabla f^{\top}\mathbf{v}>0. \end{cases} \]

The last line subtracts exactly the outward normal component. It does not reverse the update and does not alter tangential motion.

For a scalar interval \( \Omega=[\ell,u] \), point projection is

\[ \Pi_{[\ell,u]}(x)=\min\{u,\max\{\ell,x\}\}. \]

The continuous-time direction projection is

\[ \operatorname{Proj}_{[\ell,u]}(\hat\theta,v)= \begin{cases} 0, & \hat\theta=\ell\text{ and }v<0,\\ 0, & \hat\theta=u\text{ and }v>0,\\ v, & \text{otherwise}. \end{cases} \]

For a box \( \Omega=\prod_{i=1}^{p}[\ell_i,u_i] \), both operations act componentwise:

\[ [\Pi_{\Omega}(\mathbf{x})]_i =\min\{u_i,\max\{\ell_i,x_i\}\}, \]

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

In floating-point software, exact boundary tests should be replaced by a small tolerance. A component is regarded as being on the lower boundary when \( \hat\theta_i\le\ell_i+\varepsilon \), and similarly for the upper boundary.

For the ball \( \Omega=\{\boldsymbol{\theta}:\|\boldsymbol{\theta}-\mathbf{c}\|_2\le R\} \), the point projection is

\[ \Pi_{\Omega}(\mathbf{x})= \begin{cases} \mathbf{x}, & \|\mathbf{x}-\mathbf{c}\|_2\le R,\\ \mathbf{c}+R\dfrac{\mathbf{x}-\mathbf{c} } {\|\mathbf{x}-\mathbf{c}\|_2}, & \|\mathbf{x}-\mathbf{c}\|_2>R. \end{cases} \]

At a boundary point, define \( \mathbf{r}=\hat{\boldsymbol{\theta} }-\mathbf{c} \). If \( \mathbf{r}^{\top}\mathbf{v}>0 \), the direction is outward and

\[ \operatorname{Proj}_{\Omega}(\hat{\boldsymbol{\theta} },\mathbf{v}) =\mathbf{v}-\frac{\mathbf{r}^{\top}\mathbf{v} } {\mathbf{r}^{\top}\mathbf{r} }\mathbf{r}. \]

For an ellipsoid with \( g(\boldsymbol{\theta})=(\boldsymbol{\theta}-\mathbf{c})^{\top} \mathbf{Q}(\boldsymbol{\theta}-\mathbf{c})-1 \), the outward normal is

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

The smooth-boundary formula in Section 6 then gives the projected direction. Unlike box or ball point projection, Euclidean point projection onto a general ellipsoid is not usually obtained by simple componentwise clipping; it may require solving a scalar multiplier equation or a convex quadratic program.

Consider the projected differential equation

\[ \dot{\hat{\boldsymbol{\theta} } } =\operatorname{Proj}_{\Omega} (\hat{\boldsymbol{\theta} },\mathbf{v}(t)), \qquad \hat{\boldsymbol{\theta} }(0)\in\Omega. \]

At every boundary point, the right-hand side belongs to the tangent cone. Hence the vector field never points strictly outside the set. Under the usual existence assumptions for the projected differential equation, this gives

\[ \hat{\boldsymbol{\theta} }(0)\in\Omega \quad\Longrightarrow\quad \hat{\boldsymbol{\theta} }(t)\in\Omega, \qquad \forall t\ge 0. \]

This statement is called positive invariance. An important implementation condition is visible here: the initial estimate must be feasible. If it is not, first apply the point projection \( \hat{\boldsymbol{\theta} }(0)= \Pi_{\Omega}(\hat{\boldsymbol{\theta} }_{ ext{raw} }(0)) \).

Define the parameter error using the convention \( \tilde{\boldsymbol{\theta} }= \hat{\boldsymbol{\theta} }-\boldsymbol{\theta}^{*} \), where \( \boldsymbol{\theta}^{*}\in\Omega \). Let

\[ \mathbf{w}=\operatorname{Proj}_{\Omega} (\hat{\boldsymbol{\theta} },\mathbf{v}), \qquad \mathbf{c}=\mathbf{w}-\mathbf{v}. \]

The correction \( \mathbf{c} \) is zero in the interior. At an active boundary, it is the negative of an outward normal component. Therefore, because the true parameter is feasible,

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

Proof. At an active boundary write

\[ \mathbf{w}=\mathbf{v}-\mathbf{n}, \qquad \mathbf{n}\in N_{\Omega}(\hat{\boldsymbol{\theta} }). \]

By the definition of the normal cone and feasibility of the true parameter,

\[ \mathbf{n}^{\top} (\boldsymbol{\theta}^{*}-\hat{\boldsymbol{\theta} })\le 0 \quad\Longrightarrow\quad (\hat{\boldsymbol{\theta} }-\boldsymbol{\theta}^{*})^{\top} \mathbf{n}\ge 0. \]

Since \( \mathbf{c}=-\mathbf{n} \),

\[ \tilde{\boldsymbol{\theta} }^{\top}\mathbf{c} =-\tilde{\boldsymbol{\theta} }^{\top}\mathbf{n}\le 0. \]

Now choose the parameter-error Lyapunov term

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

With a projection defined in the corresponding \( \boldsymbol{\Gamma}^{-1} \) metric and an unconstrained candidate update \( \boldsymbol{\Gamma}\boldsymbol{\tau} \),

\[ \dot V_{\theta} =\tilde{\boldsymbol{\theta} }^{\top}\boldsymbol{\Gamma}^{-1} \operatorname{Proj}^{\boldsymbol{\Gamma} }_{\Omega} (\hat{\boldsymbol{\theta} },\boldsymbol{\Gamma}\boldsymbol{\tau}) \le \tilde{\boldsymbol{\theta} }^{\top}\boldsymbol{\tau}. \]

Thus projection does not introduce a positive extra term compared with the unconstrained update. This is the precise reason projection can be inserted into many Lyapunov-based adaptive laws without destroying the sign estimate already obtained for the derivative.

For a differentiable objective \( J(\boldsymbol{\theta}) \), a constrained gradient flow is

\[ \dot{\hat{\boldsymbol{\theta} } } =\operatorname{Proj}_{\Omega} (\hat{\boldsymbol{\theta} },-\nabla J(\hat{\boldsymbol{\theta} })). \]

Let \( \mathbf{w} \) denote the projected direction. The Moreau decomposition of the raw direction into tangent and normal components gives orthogonality between those components. Consequently,

\[ \dot J=\nabla J^{\top}\mathbf{w}=-\|\mathbf{w}\|_{2}^{2}\le 0. \]

A projected equilibrium need not satisfy \( \nabla J=\mathbf{0} \). At a constrained optimum, the negative gradient may point outside the set and be cancelled by a normal. The equilibrium condition is

\[ \mathbf{0}\in\nabla J(\hat{\boldsymbol{\theta} }) +N_{\Omega}(\hat{\boldsymbol{\theta} }). \]

This is the first-order optimality condition for minimization over a convex set.

A direct Euler discretization of a projected continuous-time direction is

\[ \hat{\boldsymbol{\theta} }_{k+1}^{\text{candidate} } =\hat{\boldsymbol{\theta} }_{k} +h\operatorname{Proj}_{\Omega} (\hat{\boldsymbol{\theta} }_{k},\mathbf{v}_{k}). \]

Even if the direction is tangent at the current point, a finite step can overshoot a curved boundary or cross a box boundary when the state is close to, but not exactly on, that boundary. A robust sampled implementation uses

\[ \hat{\boldsymbol{\theta} }_{k+1} =\Pi_{\Omega}\!\left( \hat{\boldsymbol{\theta} }_{k} +h\operatorname{Proj}_{\Omega} (\hat{\boldsymbol{\theta} }_{k},\mathbf{v}_{k})\right). \]

flowchart TD
  A["Read current parameter estimate"] --> B["Compute raw update"]
  B --> C["Project update onto feasible tangent directions"]
  C --> D["Take one numerical integration step"]
  D --> E["Project candidate point back into admissible set"]
  E --> F["Store constrained estimate"]
  F --> A
        

Recommended safeguards are:

• project the initial estimate before starting;
• validate that every lower bound is not greater than its upper bound;
• use a scale-aware boundary tolerance;
• avoid excessively large integration steps;
• log how often each constraint is active, because permanent boundary activity may indicate incorrect prior bounds.

The downloadable programs use closed-form box and ball formulas from scratch so that the geometric operation is visible. General polyhedral or nonlinear convex sets can instead be handled by a quadratic-programming solver.

The program uses NumPy, compares unconstrained and constrained gradient flows, and numerically verifies the non-positive Lyapunov correction term.

Chapter3_Lesson4.py

"""Chapter 3, Lesson 4: Projection operators and parameter constraints.

The example compares an unconstrained gradient flow with a projected flow.
It also implements Euclidean point projection and tangent-direction projection
for boxes and Euclidean balls.
"""

from __future__ import annotations

import numpy as np


def project_point_box(x: np.ndarray, lower: np.ndarray, upper: np.ndarray) -> np.ndarray:
    """Euclidean projection of x onto an axis-aligned box."""
    x = np.asarray(x, dtype=float)
    lower = np.asarray(lower, dtype=float)
    upper = np.asarray(upper, dtype=float)
    if x.shape != lower.shape or x.shape != upper.shape:
        raise ValueError("x, lower, and upper must have the same shape")
    if np.any(lower > upper):
        raise ValueError("every lower bound must not exceed its upper bound")
    return np.minimum(np.maximum(x, lower), upper)


def project_direction_box(
    theta: np.ndarray,
    direction: np.ndarray,
    lower: np.ndarray,
    upper: np.ndarray,
    tolerance: float = 1.0e-12,
) -> np.ndarray:
    """Project a continuous-time update direction onto the box tangent cone."""
    theta = np.asarray(theta, dtype=float)
    direction = np.asarray(direction, dtype=float)
    lower = np.asarray(lower, dtype=float)
    upper = np.asarray(upper, dtype=float)
    if not (theta.shape == direction.shape == lower.shape == upper.shape):
        raise ValueError("all vectors must have the same shape")

    projected = direction.copy()
    outward_at_lower = (theta <= lower + tolerance) & (direction < 0.0)
    outward_at_upper = (theta >= upper - tolerance) & (direction > 0.0)
    projected[outward_at_lower | outward_at_upper] = 0.0
    return projected


def project_point_ball(x: np.ndarray, centre: np.ndarray, radius: float) -> np.ndarray:
    """Euclidean projection of x onto a closed Euclidean ball."""
    x = np.asarray(x, dtype=float)
    centre = np.asarray(centre, dtype=float)
    if x.shape != centre.shape:
        raise ValueError("x and centre must have the same shape")
    if radius <= 0.0:
        raise ValueError("radius must be positive")

    displacement = x - centre
    norm = np.linalg.norm(displacement)
    if norm <= radius:
        return x.copy()
    return centre + radius * displacement / norm


def project_direction_ball(
    theta: np.ndarray,
    direction: np.ndarray,
    centre: np.ndarray,
    radius: float,
    tolerance: float = 1.0e-12,
) -> np.ndarray:
    """Project a direction onto the tangent cone of a Euclidean ball."""
    theta = np.asarray(theta, dtype=float)
    direction = np.asarray(direction, dtype=float)
    centre = np.asarray(centre, dtype=float)
    if not (theta.shape == direction.shape == centre.shape):
        raise ValueError("all vectors must have the same shape")
    if radius <= 0.0:
        raise ValueError("radius must be positive")

    radial = theta - centre
    radial_norm = np.linalg.norm(radial)
    outward_measure = float(radial @ direction)
    on_boundary = radial_norm >= radius - tolerance
    if on_boundary and outward_measure > 0.0:
        return direction - (outward_measure / (radial @ radial)) * radial
    return direction.copy()


def verify_lyapunov_projection_inequality() -> float:
    """Return the projection correction term; it should be non-positive."""
    lower = np.array([-1.0, 0.25])
    upper = np.array([1.5, 2.0])
    theta_hat = np.array([1.5, 0.25])
    theta_star = np.array([1.0, 0.5])
    raw_direction = np.array([0.8, -0.6])
    projected = project_direction_box(theta_hat, raw_direction, lower, upper)
    correction = projected - raw_direction
    return float((theta_hat - theta_star) @ correction)


def simulate_projected_gradient() -> tuple[np.ndarray, np.ndarray]:
    """Compare unconstrained and projected gradient flows."""
    lower = np.array([-1.0, 0.25])
    upper = np.array([1.5, 2.0])
    theta_initial = np.array([0.0, 1.0])

    # The unconstrained minimiser lies outside the admissible box.
    target = np.array([2.2, -0.4])
    gain = np.diag([1.2, 0.8])
    step = 0.01
    iterations = 800

    theta_free = theta_initial.copy()
    theta_projected = theta_initial.copy()

    for _ in range(iterations):
        free_direction = -gain @ (theta_free - target)
        theta_free = theta_free + step * free_direction

        raw_direction = -gain @ (theta_projected - target)
        feasible_direction = project_direction_box(
            theta_projected, raw_direction, lower, upper
        )
        # The final point projection protects a sampled implementation against
        # finite-step overshoot while preserving the continuous-time idea.
        theta_projected = project_point_box(
            theta_projected + step * feasible_direction, lower, upper
        )

    if np.any(theta_projected < lower - 1.0e-10) or np.any(
        theta_projected > upper + 1.0e-10
    ):
        raise RuntimeError("projected estimate escaped the admissible set")

    return theta_free, theta_projected


def main() -> None:
    free, constrained = simulate_projected_gradient()
    correction_term = verify_lyapunov_projection_inequality()

    print("Unconstrained final estimate:", free)
    print("Projected final estimate:  ", constrained)
    print("Lyapunov correction term:  ", correction_term)
    print("Inequality satisfied:       ", correction_term <= 1.0e-12)

    point = np.array([3.0, 4.0])
    ball_projection = project_point_ball(point, np.zeros(2), radius=2.0)
    print("Projection of [3, 4] onto radius-2 ball:", ball_projection)


if __name__ == "__main__":
    main()
      

This fixed-size implementation uses only the C++ standard library and is suitable as a starting point for embedded or real-time code. Eigen can replace the small helper functions when parameter dimension is larger.

Chapter3_Lesson4.cpp

// Chapter 3, Lesson 4: Projection operators and parameter constraints.
// Standard-library implementation for a two-parameter example.

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

using Vec2 = std::array<double, 2>;

Vec2 projectPointBox(const Vec2& x, const Vec2& lower, const Vec2& upper) {
    Vec2 result{};
    for (std::size_t i = 0; i < x.size(); ++i) {
        if (lower[i] > upper[i]) {
            throw std::invalid_argument("lower bound exceeds upper bound");
        }
        result[i] = std::min(std::max(x[i], lower[i]), upper[i]);
    }
    return result;
}

Vec2 projectDirectionBox(
    const Vec2& theta,
    const Vec2& direction,
    const Vec2& lower,
    const Vec2& upper,
    double tolerance = 1.0e-12
) {
    Vec2 result = direction;
    for (std::size_t i = 0; i < theta.size(); ++i) {
        const bool outwardAtLower =
            theta[i] <= lower[i] + tolerance && direction[i] < 0.0;
        const bool outwardAtUpper =
            theta[i] >= upper[i] - tolerance && direction[i] > 0.0;
        if (outwardAtLower || outwardAtUpper) {
            result[i] = 0.0;
        }
    }
    return result;
}

double dot(const Vec2& a, const Vec2& b) {
    return a[0] * b[0] + a[1] * b[1];
}

double norm(const Vec2& x) {
    return std::sqrt(dot(x, x));
}

Vec2 projectPointBall(const Vec2& x, const Vec2& centre, double radius) {
    if (radius <= 0.0) {
        throw std::invalid_argument("radius must be positive");
    }
    const Vec2 displacement{x[0] - centre[0], x[1] - centre[1]};
    const double displacementNorm = norm(displacement);
    if (displacementNorm <= radius) {
        return x;
    }
    const double scale = radius / displacementNorm;
    return Vec2{
        centre[0] + scale * displacement[0],
        centre[1] + scale * displacement[1]
    };
}

Vec2 projectDirectionBall(
    const Vec2& theta,
    const Vec2& direction,
    const Vec2& centre,
    double radius,
    double tolerance = 1.0e-12
) {
    if (radius <= 0.0) {
        throw std::invalid_argument("radius must be positive");
    }
    const Vec2 radial{theta[0] - centre[0], theta[1] - centre[1]};
    const double radialNorm = norm(radial);
    const double outwardMeasure = dot(radial, direction);
    if (radialNorm >= radius - tolerance && outwardMeasure > 0.0) {
        const double denominator = dot(radial, radial);
        const double scale = outwardMeasure / denominator;
        return Vec2{
            direction[0] - scale * radial[0],
            direction[1] - scale * radial[1]
        };
    }
    return direction;
}

double verifyLyapunovProjectionInequality() {
    const Vec2 lower{-1.0, 0.25};
    const Vec2 upper{1.5, 2.0};
    const Vec2 thetaHat{1.5, 0.25};
    const Vec2 thetaStar{1.0, 0.5};
    const Vec2 rawDirection{0.8, -0.6};
    const Vec2 projected =
        projectDirectionBox(thetaHat, rawDirection, lower, upper);
    const Vec2 error{thetaHat[0] - thetaStar[0], thetaHat[1] - thetaStar[1]};
    const Vec2 correction{
        projected[0] - rawDirection[0],
        projected[1] - rawDirection[1]
    };
    return dot(error, correction);
}

std::pair<Vec2, Vec2> simulateProjectedGradient() {
    const Vec2 lower{-1.0, 0.25};
    const Vec2 upper{1.5, 2.0};
    const Vec2 target{2.2, -0.4};
    const Vec2 gain{1.2, 0.8};
    const double step = 0.01;
    const int iterations = 800;

    Vec2 thetaFree{0.0, 1.0};
    Vec2 thetaProjected = thetaFree;

    for (int k = 0; k < iterations; ++k) {
        Vec2 freeDirection{
            -gain[0] * (thetaFree[0] - target[0]),
            -gain[1] * (thetaFree[1] - target[1])
        };
        thetaFree[0] += step * freeDirection[0];
        thetaFree[1] += step * freeDirection[1];

        const Vec2 rawDirection{
            -gain[0] * (thetaProjected[0] - target[0]),
            -gain[1] * (thetaProjected[1] - target[1])
        };
        const Vec2 feasibleDirection = projectDirectionBox(
            thetaProjected, rawDirection, lower, upper
        );
        const Vec2 candidate{
            thetaProjected[0] + step * feasibleDirection[0],
            thetaProjected[1] + step * feasibleDirection[1]
        };
        thetaProjected = projectPointBox(candidate, lower, upper);
    }

    return {thetaFree, thetaProjected};
}

void printVector(const char* label, const Vec2& x) {
    std::cout << label << "[" << x[0] << ", " << x[1] << "]\n";
}

int main() {
    try {
        const auto [freeEstimate, projectedEstimate] = simulateProjectedGradient();
        const double correctionTerm = verifyLyapunovProjectionInequality();

        std::cout << std::fixed << std::setprecision(8);
        printVector("Unconstrained final estimate: ", freeEstimate);
        printVector("Projected final estimate:   ", projectedEstimate);
        std::cout << "Lyapunov correction term:   " << correctionTerm << '\n';
        std::cout << "Inequality satisfied:        "
                  << std::boolalpha << (correctionTerm <= 1.0e-12) << '\n';

        const Vec2 point{3.0, 4.0};
        const Vec2 ballProjection = projectPointBall(point, Vec2{0.0, 0.0}, 2.0);
        printVector("Projection onto radius-2 ball: ", ballProjection);
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
    return 0;
}
      

The Java version uses standard arrays. For larger models, EJML can provide matrix and decomposition operations while the projection logic remains unchanged.

Chapter3_Lesson4.java

import java.util.Arrays;

/**
 * Chapter 3, Lesson 4: Projection operators and parameter constraints.
 * Standard-library implementation for a two-parameter example.
 */
public final class Chapter3_Lesson4 {
    private Chapter3_Lesson4() {
    }

    static double[] projectPointBox(double[] x, double[] lower, double[] upper) {
        requireSameLength(x, lower, upper);
        double[] result = new double[x.length];
        for (int i = 0; i < x.length; i++) {
            if (lower[i] > upper[i]) {
                throw new IllegalArgumentException("lower bound exceeds upper bound");
            }
            result[i] = Math.min(Math.max(x[i], lower[i]), upper[i]);
        }
        return result;
    }

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

    static double[] projectPointBall(double[] x, double[] centre, double radius) {
        requireSameLength(x, centre);
        if (radius <= 0.0) {
            throw new IllegalArgumentException("radius must be positive");
        }
        double[] displacement = subtract(x, centre);
        double displacementNorm = norm(displacement);
        if (displacementNorm <= radius) {
            return x.clone();
        }
        double scale = radius / displacementNorm;
        return new double[] {
            centre[0] + scale * displacement[0],
            centre[1] + scale * displacement[1]
        };
    }

    static double[] projectDirectionBall(
            double[] theta,
            double[] direction,
            double[] centre,
            double radius,
            double tolerance) {
        requireSameLength(theta, direction, centre);
        if (radius <= 0.0) {
            throw new IllegalArgumentException("radius must be positive");
        }
        double[] radial = subtract(theta, centre);
        double radialNorm = norm(radial);
        double outwardMeasure = dot(radial, direction);
        if (radialNorm >= radius - tolerance && outwardMeasure > 0.0) {
            double scale = outwardMeasure / dot(radial, radial);
            return new double[] {
                direction[0] - scale * radial[0],
                direction[1] - scale * radial[1]
            };
        }
        return direction.clone();
    }

    static double verifyLyapunovProjectionInequality() {
        double[] lower = {-1.0, 0.25};
        double[] upper = {1.5, 2.0};
        double[] thetaHat = {1.5, 0.25};
        double[] thetaStar = {1.0, 0.5};
        double[] rawDirection = {0.8, -0.6};
        double[] projected = projectDirectionBox(
                thetaHat, rawDirection, lower, upper, 1.0e-12);
        double[] error = subtract(thetaHat, thetaStar);
        double[] correction = subtract(projected, rawDirection);
        return dot(error, correction);
    }

    static double[][] simulateProjectedGradient() {
        double[] lower = {-1.0, 0.25};
        double[] upper = {1.5, 2.0};
        double[] target = {2.2, -0.4};
        double[] gain = {1.2, 0.8};
        double step = 0.01;
        int iterations = 800;

        double[] thetaFree = {0.0, 1.0};
        double[] thetaProjected = thetaFree.clone();

        for (int k = 0; k < iterations; k++) {
            double[] freeDirection = {
                -gain[0] * (thetaFree[0] - target[0]),
                -gain[1] * (thetaFree[1] - target[1])
            };
            thetaFree[0] += step * freeDirection[0];
            thetaFree[1] += step * freeDirection[1];

            double[] rawDirection = {
                -gain[0] * (thetaProjected[0] - target[0]),
                -gain[1] * (thetaProjected[1] - target[1])
            };
            double[] feasibleDirection = projectDirectionBox(
                    thetaProjected, rawDirection, lower, upper, 1.0e-12);
            double[] candidate = {
                thetaProjected[0] + step * feasibleDirection[0],
                thetaProjected[1] + step * feasibleDirection[1]
            };
            thetaProjected = projectPointBox(candidate, lower, upper);
        }
        return new double[][] {thetaFree, thetaProjected};
    }

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

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

    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;
    }

    static void requireSameLength(double[] first, double[]... remaining) {
        for (double[] vector : remaining) {
            if (vector.length != first.length) {
                throw new IllegalArgumentException("all vectors must have the same length");
            }
        }
    }

    public static void main(String[] args) {
        double[][] estimates = simulateProjectedGradient();
        double correctionTerm = verifyLyapunovProjectionInequality();

        System.out.println("Unconstrained final estimate: " + Arrays.toString(estimates[0]));
        System.out.println("Projected final estimate:   " + Arrays.toString(estimates[1]));
        System.out.println("Lyapunov correction term:   " + correctionTerm);
        System.out.println("Inequality satisfied:        " + (correctionTerm <= 1.0e-12));

        double[] ballProjection = projectPointBall(
                new double[] {3.0, 4.0}, new double[] {0.0, 0.0}, 2.0);
        System.out.println("Projection onto radius-2 ball: "
                + Arrays.toString(ballProjection));
    }
}
      

The local MATLAB functions implement point and direction projection. In Simulink, place projectDirectionBox inside a MATLAB Function block before an Integrator. For a discrete design, apply projectPointBox to the candidate estimate before a Unit Delay or Discrete-Time Integrator.

Chapter3_Lesson4.m

%% Chapter 3, Lesson 4: Projection Operators and Parameter Constraints
% This script compares an unconstrained gradient flow with a projected flow.
% The local functions are also suitable for use inside a Simulink MATLAB
% Function block when the vector dimensions are fixed.

clear; clc;

lower = [-1.0; 0.25];
upper = [ 1.5; 2.00];
thetaInitial = [0.0; 1.0];
target = [2.2; -0.4];       % Unconstrained minimiser is outside the box.
gain = diag([1.2, 0.8]);
step = 0.01;
iterations = 800;

thetaFree = thetaInitial;
thetaProjected = thetaInitial;

for k = 1:iterations
    freeDirection = -gain * (thetaFree - target);
    thetaFree = thetaFree + step * freeDirection;

    rawDirection = -gain * (thetaProjected - target);
    feasibleDirection = projectDirectionBox( ...
        thetaProjected, rawDirection, lower, upper, 1.0e-12);

    % A sampled implementation can overshoot during one finite step, so a
    % point projection is applied after the tangent-direction projection.
    candidate = thetaProjected + step * feasibleDirection;
    thetaProjected = projectPointBox(candidate, lower, upper);
end

fprintf('Unconstrained final estimate: [% .8f, % .8f]\n', ...
    thetaFree(1), thetaFree(2));
fprintf('Projected final estimate:   [% .8f, % .8f]\n', ...
    thetaProjected(1), thetaProjected(2));

% Verify the Lyapunov-compatible correction inequality.
thetaHat = [1.5; 0.25];
thetaStar = [1.0; 0.5];
rawDirection = [0.8; -0.6];
projectedDirection = projectDirectionBox( ...
    thetaHat, rawDirection, lower, upper, 1.0e-12);
correctionTerm = (thetaHat - thetaStar).' * ...
    (projectedDirection - rawDirection);

fprintf('Lyapunov correction term:   % .8f\n', correctionTerm);
fprintf('Inequality satisfied:        %d\n', correctionTerm <= 1.0e-12);

ballProjection = projectPointBall([3.0; 4.0], [0.0; 0.0], 2.0);
fprintf('Projection onto radius-2 ball: [% .8f, % .8f]\n', ...
    ballProjection(1), ballProjection(2));

%% Simulink use
% In a MATLAB Function block, call:
%
%   thetaDot = projectDirectionBox(thetaHat, rawUpdate, lower, upper, 1e-9);
%
% Integrate thetaDot with an Integrator block. For a discrete model, compute
% the candidate estimate and pass it through projectPointBox before storing
% it in a Unit Delay or Discrete-Time Integrator block.

function projected = projectPointBox(x, lower, upper)
%PROJECTPOINTBOX Euclidean projection onto an axis-aligned box.
    if any(size(x) ~= size(lower)) || any(size(x) ~= size(upper))
        error('x, lower, and upper must have the same size.');
    end
    if any(lower > upper)
        error('Every lower bound must not exceed its upper bound.');
    end
    projected = min(max(x, lower), upper);
end

function projected = projectDirectionBox(theta, direction, lower, upper, tolerance)
%PROJECTDIRECTIONBOX Project a velocity onto the tangent cone of a box.
    if nargin < 5
        tolerance = 1.0e-12;
    end
    projected = direction;
    outwardAtLower = (theta <= lower + tolerance) & (direction < 0.0);
    outwardAtUpper = (theta >= upper - tolerance) & (direction > 0.0);
    projected(outwardAtLower | outwardAtUpper) = 0.0;
end

function projected = projectPointBall(x, centre, radius)
%PROJECTPOINTBALL Euclidean projection onto a closed Euclidean ball.
    if radius <= 0.0
        error('Radius must be positive.');
    end
    displacement = x - centre;
    displacementNorm = norm(displacement, 2);
    if displacementNorm <= radius
        projected = x;
    else
        projected = centre + radius * displacement / displacementNorm;
    end
end

function projected = projectDirectionBall(theta, direction, centre, radius, tolerance)
%PROJECTDIRECTIONBALL Project a velocity onto a ball tangent cone.
    if nargin < 5
        tolerance = 1.0e-12;
    end
    if radius <= 0.0
        error('Radius must be positive.');
    end
    radial = theta - centre;
    outwardMeasure = radial.' * direction;
    if norm(radial, 2) >= radius - tolerance && outwardMeasure > 0.0
        projected = direction - ...
            (outwardMeasure / (radial.' * radial)) * radial;
    else
        projected = direction;
    end
end
      

The notebook stores the complete implementation in an input cell. Built-in region functions can later replace the closed-form formulas for more general convex sets.

Chapter3_Lesson4.nb

Notebook[{
 Cell["Chapter 3, Lesson 4: Projection Operators and Parameter Constraints", "Title"],
 Cell["The input cell implements Euclidean point projection and tangent-direction projection for boxes and Euclidean balls, then compares unconstrained and projected gradient flows.", "Text"],
 Cell[BoxData[
"ClearAll[projectPointBox, projectDirectionBox, projectPointBall,
  projectDirectionBall, simulateProjectedGradient];
projectPointBox[x_List, lower_List, upper_List] := Module[{},
  If[Length[x] != Length[lower] || Length[x] != Length[upper],
    Print[\"All vectors must have the same length.\"]; Abort[]];
  If[AnyTrue[MapThread[Greater, {lower, upper}], TrueQ],
    Print[\"Every lower bound must not exceed its upper bound.\"]; Abort[]];
  MapThread[Clip[ #1, { #2, #3}] &, {N[x], N[lower], N[upper]}]
];
projectDirectionBox[theta_List, direction_List, lower_List, upper_List,
  tolerance_: 10^-12] := Module[{projected = N[direction]},
  Do[If[(theta[[i]] <= lower[[i]] + tolerance && direction[[i]] < 0) ||
        (theta[[i]] >= upper[[i]] - tolerance && direction[[i]] > 0),
      projected[[i]] = 0.0], {i, Length[theta]}];
  projected
];
projectPointBall[x_List, centre_List, radius_?NumericQ] := Module[
  {displacement, displacementNorm},
  If[radius <= 0, Print[\"Radius must be positive.\"]; Abort[]];
  displacement = N[x - centre]; displacementNorm = Norm[displacement];
  If[displacementNorm <= radius, N[x],
    N[centre] + radius displacement/displacementNorm]
];
projectDirectionBall[theta_List, direction_List, centre_List,
  radius_?NumericQ, tolerance_: 10^-12] := Module[{radial, outwardMeasure},
  If[radius <= 0, Print[\"Radius must be positive.\"]; Abort[]];
  radial = N[theta - centre]; outwardMeasure = radial . N[direction];
  If[Norm[radial] >= radius - tolerance && outwardMeasure > 0,
    N[direction] - (outwardMeasure/(radial . radial)) radial, N[direction]]
];
simulateProjectedGradient[] := Module[
  {lower = {-1.0, 0.25}, upper = {1.5, 2.0}, target = {2.2, -0.4},
   gain = DiagonalMatrix[{1.2, 0.8}], step = 0.01, iterations = 800,
   thetaFree = {0.0, 1.0}, thetaProjected = {0.0, 1.0}, freeDirection,
   rawDirection, feasibleDirection, candidate},
  Do[freeDirection = -gain . (thetaFree - target);
    thetaFree = thetaFree + step freeDirection;
    rawDirection = -gain . (thetaProjected - target);
    feasibleDirection = projectDirectionBox[thetaProjected, rawDirection, lower, upper];
    candidate = thetaProjected + step feasibleDirection;
    thetaProjected = projectPointBox[candidate, lower, upper], {iterations}];
  {thetaFree, thetaProjected}
];
{thetaFree, thetaProjected} = simulateProjectedGradient[];
lower = {-1.0, 0.25}; upper = {1.5, 2.0};
thetaHat = {1.5, 0.25}; thetaStar = {1.0, 0.5}; rawDirection = {0.8, -0.6};
projectedDirection = projectDirectionBox[thetaHat, rawDirection, lower, upper];
correctionTerm = (thetaHat - thetaStar) . (projectedDirection - rawDirection);
Print[\"Unconstrained final estimate: \", thetaFree];
Print[\"Projected final estimate:   \", thetaProjected];
Print[\"Lyapunov correction term:   \", correctionTerm];
Print[\"Inequality satisfied:        \", correctionTerm <= 10^-12];
Print[\"Projection onto radius-2 ball: \",
  projectPointBall[{3.0, 4.0}, {0.0, 0.0}, 2.0]];"
 ], "Input"]
}, WindowTitle -> "Chapter3_Lesson4"]
      

The programs minimize the quadratic objective

\[ J(\boldsymbol{\theta})= rac{1}{2} (\boldsymbol{\theta}-\mathbf{q})^{\top}\mathbf{H} (\boldsymbol{\theta}-\mathbf{q}), \]

with

\[ \mathbf{q}=\begin{bmatrix}2.2\\-0.4\end{bmatrix},\qquad \mathbf{H}=\begin{bmatrix}1.2&0\\0&0.8\end{bmatrix},\qquad \Omega=[-1,1.5]\times[0.25,2]. \]

The unconstrained minimizer is \( \mathbf{q} \), which is infeasible in both coordinates. Because the Hessian is diagonal and the set is a box, the constrained optimum is obtained by clipping each coordinate:

\[ \boldsymbol{\theta}_{\Omega}^{*} =\begin{bmatrix}1.5\\0.25\end{bmatrix}. \]

The simulation converges to this boundary point. The verification example also produces a correction term of \( -0.55 \), confirming that the projection modification contributes a non-positive term to the Euclidean parameter-error Lyapunov derivative.

Problem 1 (Scalar Direction Projection): Let \( \Omega=[-2,3] \). Evaluate \( \operatorname{Proj}_{\Omega}(\hat\theta,v) \) for (a) \( \hat\theta=1,v=4 \), (b) \( \hat\theta=3,v=4 \), and (c) \( \hat\theta=-2,v=-1 \).

Solution:

(a) The estimate is in the interior, so every scalar direction is tangent and the result is \( 4 \). (b) The estimate is at the upper bound and the direction is outward, so the result is \( 0 \). (c) The estimate is at the lower bound and the direction is outward, so the result is again \( 0 \).

Problem 2 (Box Point Projection): Project \( \mathbf{x}=[-3,\;0.4,\;7]^{\top} \) onto

\[ \Omega=[-1,2]\times[0,1]\times[4,6]. \]

Solution: Clip each coordinate independently:

\[ \Pi_{\Omega}(\mathbf{x})= \begin{bmatrix}-1\\0.4\\6\end{bmatrix}. \]

The first coordinate is raised to its lower bound, the second is unchanged, and the third is reduced to its upper bound.

Problem 3 (Ball Direction Projection): Let \( \Omega=\{\boldsymbol{\theta}:\|\boldsymbol{\theta}\|_2\le 2\} \), \( \hat{\boldsymbol{\theta} }=[2,0]^{\top} \), and \( \mathbf{v}=[3,4]^{\top} \). Find the projected direction.

Solution: The radial vector is \( \mathbf{r}=[2,0]^{\top} \), and

\[ \mathbf{r}^{\top}\mathbf{v}=6>0, \]

so the direction has an outward component. Hence

\[ \begin{aligned} \operatorname{Proj}_{\Omega}(\hat{\boldsymbol{\theta} },\mathbf{v}) &=\mathbf{v}-\frac{\mathbf{r}^{\top}\mathbf{v} } {\mathbf{r}^{\top}\mathbf{r} }\mathbf{r}\\ &=\begin{bmatrix}3\\4\end{bmatrix} -\frac{6}{4}\begin{bmatrix}2\\0\end{bmatrix} =\begin{bmatrix}0\\4\end{bmatrix}. \end{aligned} \]

The remaining velocity is tangent to the circle.

Problem 4 (Distance Reduction by Point Projection): Let \( \mathbf{p}=\Pi_{\Omega}(\mathbf{x}) \) and \( \boldsymbol{\theta}^{*}\in\Omega \). Prove

\[ \|\mathbf{p}-\boldsymbol{\theta}^{*}\|_2^2 \le\|\mathbf{x}-\boldsymbol{\theta}^{*}\|_2^2 -\|\mathbf{x}-\mathbf{p}\|_2^2. \]

Solution: Expand the right-hand distance:

\[ \begin{aligned} \|\mathbf{x}-\boldsymbol{\theta}^{*}\|_2^2 &=\|\mathbf{x}-\mathbf{p}+\mathbf{p}-\boldsymbol{\theta}^{*}\|_2^2\\ &=\|\mathbf{x}-\mathbf{p}\|_2^2 +\|\mathbf{p}-\boldsymbol{\theta}^{*}\|_2^2 +2(\mathbf{x}-\mathbf{p})^{\top} (\mathbf{p}-\boldsymbol{\theta}^{*}). \end{aligned} \]

The projection variational inequality with \( \mathbf{z}=\boldsymbol{\theta}^{*} \) gives

\[ (\mathbf{x}-\mathbf{p})^{\top} (\boldsymbol{\theta}^{*}-\mathbf{p})\le 0, \]

so \( (\mathbf{x}-\mathbf{p})^{\top} (\mathbf{p}-\boldsymbol{\theta}^{*})\ge 0 \). Substitution and rearrangement produce the required inequality.

Problem 5 (Lyapunov Correction): For a scalar upper bound \( \hat\theta\le u \), suppose \( \hat\theta=u \), \( \theta^{*}<u \), and \( v>0 \). Show directly that the projection correction is non-positive in the parameter-error product.

Solution: At the upper boundary with an outward direction,

\[ \operatorname{Proj}_{\Omega}(u,v)=0, \qquad c=0-v=-v<0. \]

Also \( \tilde\theta=u-\theta^{*}>0 \). Therefore

\[ \tilde\theta c=(u-\theta^{*})(-v)<0. \]

Projection strictly improves the parameter-error contribution whenever the raw direction pushes outward and the true parameter lies strictly inside.

Problem 6 (Finite-Step Overshoot): Consider the interval \( [0,1] \), current estimate \( \hat\theta_k=0.99 \), update \( v_k=0.5 \), and step \( h=0.1 \). The direction projection does not activate because the estimate is not exactly on the boundary. Compute the Euler candidate and give a safe corrected update.

Solution:

\[ \hat\theta_{k+1}^{\text{candidate} } =0.99+0.1(0.5)=1.04. \]

The candidate is infeasible. Post-step point projection gives

\[ \hat\theta_{k+1}=\Pi_{[0,1]}(1.04)=1. \]

This example explains why sampled implementations often use both direction projection and point projection.

Projection incorporates reliable prior parameter information into an update law. Euclidean point projection maps an arbitrary candidate to the nearest feasible point, while continuous-time direction projection removes only the outward normal component of a raw update. Closed convex sets provide unique point projections, tangent and normal cones provide the correct boundary geometry, and the projection correction satisfies a non-positive parameter-error inequality. Consequently, feasibility can be enforced without worsening the standard Lyapunov derivative estimate. Box and ball constraints admit simple formulas; general convex sets can be handled by numerical optimization. In a sampled implementation, point projection after integration protects against finite-step overshoot.