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

• distinguish state nonlinearity from parameter nonlinearity;
• test whether an uncertain model is linear in the unknown parameters;
• derive the exact Lyapunov cross-term cancellation available in linear-in-the-parameters models;
• derive the Taylor remainder that obstructs the same argument for nonlinear parameterizations;
• analyze local identifiability using sensitivity matrices and information Gramians;
• recognize sign ambiguity, periodic aliases, singular sensitivities, and nonconvex prediction costs;
• evaluate whether reparameterization preserves physical identifiability; and
• interpret simulations without confusing good tracking with correct parameter recovery.

The lesson assumes the Lyapunov, gradient, projection, persistent-excitation, self-tuning, and adaptive-backstepping concepts developed in earlier chapters. It does not yet present a universal nonlinear adaptive controller; instead, it establishes the structural obstacles that any valid design must overcome.

Consider an uncertain nonlinear plant

\[ \dot{\mathbf{x}} = \mathbf{f}_0(\mathbf{x},\mathbf{u}) + \mathbf{F}(\mathbf{x},\mathbf{u},\boldsymbol{\theta}), \qquad \mathbf{y}=\mathbf{h}(\mathbf{x}), \qquad \boldsymbol{\theta}\in\Theta\subset\mathbb{R}^{p}. \]

The plant may be highly nonlinear in the measurable signals while remaining linear in the unknown parameters. The defining condition is the existence of a known regressor matrix independent of the unknown parameter:

\[ \mathbf{F}(\mathbf{x},\mathbf{u},\boldsymbol{\theta}) = \mathbf{Y}(\mathbf{x},\mathbf{u})\boldsymbol{\theta}. \]

By contrast, a model is nonlinearly parameterized when no such known finite-dimensional regressor representation is available without introducing additional unknowns, constraints, or approximations:

\[ \mathbf{F}(\mathbf{x},\mathbf{u},\boldsymbol{\theta}) \neq \mathbf{Y}(\mathbf{x},\mathbf{u})\boldsymbol{\theta} \quad \text{for every admissible known } \mathbf{Y}. \]

flowchart TD
  A["Uncertain nonlinear plant"] --> B["Ask how unknown parameters enter"]
  B --> C["Known regressor times \nparameter vector"]
  B --> D["Powers, products, sin, exp, \nratios, unknown coordinates"]
  C --> E["Linear in parameters"]
  D --> F["Nonlinear in parameters"]
  E --> G["Exact Lyapunov cross-term \ncancellation may be available"]
  F --> H["Taylor remainder, aliases, \nnonconvexity, local guarantees"]
        

Therefore, the phrase nonlinear adaptive control is ambiguous unless one separately specifies nonlinearity in the state and nonlinearity in the parameter.

A useful diagnostic is differentiation with respect to the parameter. For a linear-in-parameter model, the sensitivity is independent of the parameter:

\[ \frac{\partial}{\partial\boldsymbol{\theta}} \left[\mathbf{Y}(\mathbf{x},\mathbf{u})\boldsymbol{\theta}\right] = \mathbf{Y}(\mathbf{x},\mathbf{u}). \]

Parameter-dependent sensitivity is a strong indication of nonlinear parameterization, although algebraic reparameterization may sometimes expose a different linear parameter vector.

Consider a scalar tracking-error model obtained after nominal feedback design:

\[ \dot e = -k e + \mathbf{Y}(z)\boldsymbol{\theta} -\mathbf{Y}(z)\hat{\boldsymbol{\theta}}, \qquad k>0, \]

where \(z\) collects known closed-loop signals. Define

\[ \tilde{\boldsymbol{\theta}}= \hat{\boldsymbol{\theta}}-\boldsymbol{\theta}, \qquad \dot e=-ke-\mathbf{Y}(z)\tilde{\boldsymbol{\theta}}. \]

Use the standard composite Lyapunov function

\[ V(e,\tilde{\boldsymbol{\theta}})=\frac{1}{2}e^2+ \frac{1}{2}\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta}}, \qquad \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{\mathsf T}>0. \]

Because the true parameter is constant,

\[ \dot V=-ke^2-e\mathbf{Y}(z)\tilde{\boldsymbol{\theta}}+ \tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta}}}. \]

Selecting

\[ \dot{\hat{\boldsymbol{\theta}}}= \boldsymbol{\Gamma}\mathbf{Y}^{\mathsf T}(z)e \]

gives the exact identity

\[ \dot V=-ke^2\le 0. \]

The essential property is not merely smoothness. It is the exact factorization of the uncertainty mismatch into a known regressor multiplied by the parameter error.

Now let the error model be

\[ \dot e=-ke+F(z,\boldsymbol{\theta})- F(z,\hat{\boldsymbol{\theta}}). \]

Define the sensitivity row vector evaluated at the estimate:

\[ \mathbf{J}(z,\hat{\boldsymbol{\theta}})= \left.\frac{\partial F(z,\boldsymbol{\vartheta})} {\partial\boldsymbol{\vartheta}}\right|_{ \boldsymbol{\vartheta}=\hat{\boldsymbol{\theta}}}. \]

A second-order Taylor expansion around \(\hat{\boldsymbol{\theta}}\) gives

\[ F(z,\boldsymbol{\theta})-F(z,\hat{\boldsymbol{\theta}}) =-\mathbf{J}(z,\hat{\boldsymbol{\theta}}) \tilde{\boldsymbol{\theta}}+R_2(z,\tilde{\boldsymbol{\theta}}). \]

If the parameter Hessian satisfies \(\|\nabla_{\theta}^{2}F\|\le L_F\) on the segment connecting \(\boldsymbol{\theta}\) and \(\hat{\boldsymbol{\theta}}\), then

\[ |R_2(z,\tilde{\boldsymbol{\theta}})|\le \frac{L_F}{2}\|\tilde{\boldsymbol{\theta}}\|^2. \]

The natural Jacobian update

\[ \dot{\hat{\boldsymbol{\theta}}}= \boldsymbol{\Gamma}\mathbf{J}^{\mathsf T} (z,\hat{\boldsymbol{\theta}})e \]

cancels only the first-order term. The Lyapunov derivative becomes

\[ \dot V=-ke^2+eR_2(z,\tilde{\boldsymbol{\theta}}) \]

and therefore only satisfies the bound

\[ \dot V\le -ke^2+\frac{L_F}{2}|e| \|\tilde{\boldsymbol{\theta}}\|^2. \]

The last term has no fixed sign. This is the central mathematical reason why replacing a known regressor by a local Jacobian does not automatically produce a globally stable adaptive controller.

Taylor expansion is local, but the fundamental theorem of calculus yields an exact identity:

\[ F(z,\boldsymbol{\theta})-F(z,\hat{\boldsymbol{\theta}}) =-\left[\int_{0}^{1}\mathbf{J} \left(z,\hat{\boldsymbol{\theta}}-s \tilde{\boldsymbol{\theta}}\right)ds\right] \tilde{\boldsymbol{\theta}}. \]

Define the averaged Jacobian

\[ \bar{\mathbf{J}}(z,\boldsymbol{\theta}, \hat{\boldsymbol{\theta}})=\int_0^1 \mathbf{J}\left(z,\hat{\boldsymbol{\theta}}-s \tilde{\boldsymbol{\theta}}\right)ds. \]

The mismatch is linear in the parameter error with regressor \(\bar{\mathbf{J}}\), but that regressor depends on the unknown true parameter. Thus, the algebraic factorization exists mathematically while remaining unusable for direct cancellation unless additional monotonicity, bounding, filtering, or structural information is available.

For a realized known signal trajectory \(z(t)\), two parameters are indistinguishable on \([0,T]\) if

\[ F(z(t),\boldsymbol{\theta}_1)=F(z(t),\boldsymbol{\theta}_2) \quad \text{for all }t\in[0,T]. \]

Global identifiability requires the implication

\[ F(z(t),\boldsymbol{\theta}_1)=F(z(t),\boldsymbol{\theta}_2) \ \forall t \quad \Longrightarrow \quad \boldsymbol{\theta}_1=\boldsymbol{\theta}_2. \]

Three important counterexamples are:

1. Even dependence: \(F(x,\theta)=\theta^2x\) gives \(F(x,\theta)=F(x,-\theta)\). The sign of \(\theta\) is not identifiable.
2. Product ambiguity: \(F(x,\theta_1,\theta_2)= \theta_1\theta_2\phi(x)\). Infinitely many parameter pairs produce the same product.
3. Periodic aliasing: \(F(x,\theta)=\sin(\theta x)\). Over a limited set of state values, separated parameter values may generate almost indistinguishable outputs.

Consequently, no adaptive law can recover information that the measured trajectory does not contain. Controller design cannot overcome structural non-identifiability.

Let the prediction error be \(\boldsymbol{\varepsilon}(t,\hat{\boldsymbol{\theta}})= \mathbf{F}(z(t),\hat{\boldsymbol{\theta}})- \mathbf{F}(z(t),\boldsymbol{\theta})\). A local sensitivity matrix is

\[ \mathbf{S}(t,\hat{\boldsymbol{\theta}})= \frac{\partial\mathbf{F}(z(t),\hat{\boldsymbol{\theta}})} {\partial\hat{\boldsymbol{\theta}}}. \]

The finite-horizon local information Gramian is

\[ \mathbf{G}_T(\hat{\boldsymbol{\theta}})= \int_0^T \mathbf{S}^{\mathsf T}(t,\hat{\boldsymbol{\theta}}) \mathbf{S}(t,\hat{\boldsymbol{\theta}})dt. \]

Positive definiteness of \(\mathbf{G}_T\) indicates local distinguishability of infinitesimal parameter perturbations. It does not prove global injectivity. For \(F(x,\theta)=\sin(\theta x)\),

\[ S(t,\hat\theta)=x(t)\cos(\hat\theta x(t)), \qquad G_T(\hat\theta)=\int_0^T x^2(t)\cos^2(\hat\theta x(t))dt. \]

Information is weak when \(x(t)\) remains near zero or when \(\cos(\hat\theta x(t))\) remains near zero. Thus excitation must be assessed through the nonlinear sensitivity, not only through the raw state amplitude.

Define a batch prediction cost

\[ J_T(\hat{\boldsymbol{\theta}})=\frac{1}{2} \int_0^T \|\boldsymbol{\varepsilon} (t,\hat{\boldsymbol{\theta}})\|^2dt. \]

Its gradient is

\[ \nabla J_T=\int_0^T \mathbf{S}^{\mathsf T} \boldsymbol{\varepsilon}\,dt. \]

For a scalar prediction error, the Hessian is

\[ \nabla^2J_T=\int_0^T \left[\mathbf{S}^{\mathsf T}\mathbf{S}+ \varepsilon\nabla_{\theta}^{2}F\right]dt. \]

The first term is positive semidefinite, but the second can be indefinite. Therefore, the cost may contain local minima, maxima, saddle points, or flat directions. A zero gradient can mean correct identification, an alias, a local extremum, or zero sensitivity.

A mapping \(\mathbf{M}:\Theta\mapsto\mathbb{R}^p\) is strongly monotone on \(\Theta\) if there is \(m>0\) such that

\[ (\boldsymbol{\theta}_1-\boldsymbol{\theta}_2)^{\mathsf T} [\mathbf{M}(\boldsymbol{\theta}_1)- \mathbf{M}(\boldsymbol{\theta}_2)]\ge m\|\boldsymbol{\theta}_1-\boldsymbol{\theta}_2\|^2. \]

If the measurable estimation error were exactly \(\mathbf{M}(\hat{\boldsymbol{\theta}})- \mathbf{M}(\boldsymbol{\theta})\), then the update

\[ \dot{\hat{\boldsymbol{\theta}}}=-\boldsymbol{\Gamma} [\mathbf{M}(\hat{\boldsymbol{\theta}})- \mathbf{M}(\boldsymbol{\theta})] \]

would yield

\[ \dot V_{\theta}=-\tilde{\boldsymbol{\theta}}^{\mathsf T} [\mathbf{M}(\hat{\boldsymbol{\theta}})- \mathbf{M}(\boldsymbol{\theta})] \le -m\|\tilde{\boldsymbol{\theta}}\|^2. \]

The difficulty is constructing such a measurable monotone error from closed-loop data. Convexity, concavity, triangular structure, filtered transformations, and invariant-manifold methods are examples of extra assumptions used in the research literature. They are not consequences of smoothness alone.

Some nonlinear parameterizations can be made linear by defining new parameters. For example,

\[ F(x,\theta)=\theta^2\phi(x), \qquad \alpha=\theta^2, \qquad F(x,\alpha)=\alpha\phi(x). \]

This is useful for control if only \(\alpha\) matters, but it cannot recover the sign of the physical parameter. A more complicated example is

\[ F(x,\boldsymbol{\theta})=\theta_1\phi_1(x)+ \theta_1\theta_2\phi_2(x)+\theta_2^2\phi_3(x). \]

Define

\[ \alpha_1=\theta_1, \qquad \alpha_2=\theta_1\theta_2, \qquad \alpha_3=\theta_2^2. \]

The model becomes linear in \(\boldsymbol{\alpha}\), but admissible augmented parameters satisfy the nonlinear constraint

\[ \alpha_2^2=\alpha_1^2\alpha_3, \qquad \alpha_3\ge 0. \]

An unconstrained adaptive law for \(\boldsymbol{\alpha}\) can leave this physical manifold. Overparameterization may simplify stability analysis while increasing dimension, destroying uniqueness, worsening excitation requirements, and producing estimates that cannot be mapped back to physical parameters.

flowchart TD
  A["Write uncertainty F of signals and theta"] --> B["Can it be expressed as known Y times theta?"]
  B -->|yes| C["Use linear parameter adaptive tools"]
  B -->|no| D["Check injectivity on admissible parameter set"]
  D --> E["Compute sensitivity and locate singular regions"]
  E --> F["Check convexity, monotonicity, or triangular structure"]
  F --> G["Test reparameterization and physical constraints"]
  G --> H["Specify local, semiglobal, or global objective"]
  H --> I["Add projection, bounds, and robustness safeguards"]
  I --> J["Validate multiple initial estimates and excitation levels"]
        

This assessment prevents a common error: writing a gradient update first and only later discovering that the parameter is not identifiable, the sensitivity vanishes, or the Lyapunov remainder has an uncontrolled sign.

Consider the scalar plant

\[ \dot x=-ax+u+\sin(\theta x), \qquad a>0, \]

with reference \(r(t)\), error \(e=x-r\), and certainty-equivalent control

\[ u=\dot r+ax-\sin(\hat\theta x)-ke. \]

The exact error dynamics are

\[ \dot e=-ke+\sin(\theta x)-\sin(\hat\theta x). \]

The local sensitivity is

\[ \psi(x,\hat\theta)=\frac{\partial}{\partial\hat\theta} \sin(\hat\theta x)=x\cos(\hat\theta x). \]

The pedagogical update used in the code is

\[ \dot{\hat\theta}=\operatorname{Proj}_{[-4,4]} \{\gamma e\psi(x,\hat\theta)\}. \]

This update mimics the linear-parameter cancellation locally. It is intentionally not presented as a global theorem. Running several initial estimates reveals that tracking may remain acceptable while the estimate converges to a biased value or becomes trapped at a projection boundary.

Chapter18_Lesson1.py

# Chapter18_Lesson1.py
# Local Jacobian adaptation for a nonlinearly parameterized scalar plant.
# Dependencies: numpy, matplotlib

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
from typing import Iterable

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    a: float = 1.1
    k: float = 3.0
    gamma: float = 2.5
    theta_true: float = 1.4
    theta_min: float = -4.0
    theta_max: float = 4.0
    dt: float = 0.002
    t_final: float = 30.0


def reference(t: float) -> tuple[float, float]:
    r = 0.8 * np.sin(0.7 * t) + 0.3 * np.sin(1.3 * t)
    r_dot = 0.56 * np.cos(0.7 * t) + 0.39 * np.cos(1.3 * t)
    return float(r), float(r_dot)


def project_rate(theta_hat: float, raw_rate: float, cfg: Config) -> float:
    if theta_hat <= cfg.theta_min and raw_rate < 0.0:
        return 0.0
    if theta_hat >= cfg.theta_max and raw_rate > 0.0:
        return 0.0
    return raw_rate


def dynamics(t: float, state: np.ndarray, cfg: Config) -> np.ndarray:
    x, theta_hat = state
    r, r_dot = reference(t)
    error = x - r

    q_true = np.sin(cfg.theta_true * x)
    q_hat = np.sin(theta_hat * x)
    sensitivity = x * np.cos(theta_hat * x)

    control = r_dot + cfg.a * x - q_hat - cfg.k * error
    x_dot = -cfg.a * x + control + q_true

    raw_theta_dot = cfg.gamma * error * sensitivity
    theta_dot = project_rate(theta_hat, raw_theta_dot, cfg)
    return np.array([x_dot, theta_dot], dtype=float)


def rk4_step(t: float, state: np.ndarray, dt: float, cfg: Config) -> np.ndarray:
    k1 = dynamics(t, state, cfg)
    k2 = dynamics(t + 0.5 * dt, state + 0.5 * dt * k1, cfg)
    k3 = dynamics(t + 0.5 * dt, state + 0.5 * dt * k2, cfg)
    k4 = dynamics(t + dt, state + dt * k3, cfg)
    next_state = state + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
    next_state[1] = np.clip(next_state[1], cfg.theta_min, cfg.theta_max)
    return next_state


def simulate(theta_hat_0: float, cfg: Config) -> dict[str, np.ndarray | float]:
    time = np.arange(0.0, cfg.t_final + cfg.dt, cfg.dt)
    state = np.array([0.0, theta_hat_0], dtype=float)

    x = np.empty_like(time)
    r = np.empty_like(time)
    error = np.empty_like(time)
    theta_hat = np.empty_like(time)
    control = np.empty_like(time)

    for i, t in enumerate(time):
        r_i, r_dot_i = reference(float(t))
        x[i] = state[0]
        r[i] = r_i
        error[i] = state[0] - r_i
        theta_hat[i] = state[1]
        control[i] = (
            r_dot_i
            + cfg.a * state[0]
            - np.sin(state[1] * state[0])
            - cfg.k * error[i]
        )
        if i + 1 < time.size:
            state = rk4_step(float(t), state, cfg.dt, cfg)

    rms_error = float(np.sqrt(np.mean(error**2)))
    return {
        "time": time,
        "x": x,
        "r": r,
        "error": error,
        "theta_hat": theta_hat,
        "control": control,
        "rms_error": rms_error,
        "theta_final": float(theta_hat[-1]),
    }


def save_csv(result: dict[str, np.ndarray | float], name: str) -> None:
    data = np.column_stack(
        [
            result["time"],
            result["x"],
            result["r"],
            result["error"],
            result["theta_hat"],
            result["control"],
        ]
    )
    np.savetxt(
        name,
        data,
        delimiter=",",
        header="time,x,reference,error,theta_hat,control",
        comments="",
    )


def main(initial_estimates: Iterable[float] = (-2.5, 0.0, 2.5)) -> None:
    cfg = Config()
    results: list[tuple[float, dict[str, np.ndarray | float]]] = []

    for theta_0 in initial_estimates:
        result = simulate(theta_0, cfg)
        results.append((theta_0, result))
        tag = str(theta_0).replace("-", "m").replace(".", "p")
        save_csv(result, f"Chapter18_Lesson1_theta0_{tag}.csv")
        print(
            f"theta_hat(0)={theta_0:5.2f} | "
            f"theta_hat(T)={result['theta_final']:8.4f} | "
            f"RMS(e)={result['rms_error']:.6f}"
        )

    best = min(results, key=lambda item: float(item[1]["rms_error"]))
    _, best_result = best

    plt.figure(figsize=(9, 4.8))
    plt.plot(best_result["time"], best_result["r"], label="reference r")
    plt.plot(best_result["time"], best_result["x"], label="state x")
    plt.xlabel("Time [s]")
    plt.ylabel("State")
    plt.title("Tracking for the best tested initialization")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter18_Lesson1_tracking.png", dpi=180)

    plt.figure(figsize=(9, 4.8))
    for theta_0, result in results:
        plt.plot(result["time"], result["theta_hat"], label=f"theta_hat(0)={theta_0}")
    plt.axhline(cfg.theta_true, linestyle="--", label="true theta")
    plt.xlabel("Time [s]")
    plt.ylabel("Parameter estimate")
    plt.title("Initialization sensitivity of local Jacobian adaptation")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter18_Lesson1_parameters.png", dpi=180)
    plt.show()


if __name__ == "__main__":
    main()

The script uses fourth-order Runge-Kutta integration, writes one CSV file per initialization, and generates tracking and parameter-estimate plots. Required libraries are numpy and matplotlib.

Chapter18_Lesson1.cpp

// Chapter18_Lesson1.cpp
// Compile: g++ -std=c++17 -O2 Chapter18_Lesson1.cpp -o Chapter18_Lesson1

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

struct Config {
    double a = 1.1;
    double k = 3.0;
    double gamma = 2.5;
    double thetaTrue = 1.4;
    double thetaMin = -4.0;
    double thetaMax = 4.0;
    double dt = 0.002;
    double tFinal = 30.0;
};

struct State {
    double x;
    double thetaHat;
};

struct Sample {
    double time;
    double x;
    double reference;
    double error;
    double thetaHat;
    double control;
};

std::pair<double, double> referenceSignal(double t) {
    const double r = 0.8 * std::sin(0.7 * t) + 0.3 * std::sin(1.3 * t);
    const double rDot = 0.56 * std::cos(0.7 * t) + 0.39 * std::cos(1.3 * t);
    return {r, rDot};
}

double projectRate(double thetaHat, double rawRate, const Config& cfg) {
    if (thetaHat <= cfg.thetaMin && rawRate < 0.0) return 0.0;
    if (thetaHat >= cfg.thetaMax && rawRate > 0.0) return 0.0;
    return rawRate;
}

State dynamics(double t, const State& state, const Config& cfg) {
    const auto [r, rDot] = referenceSignal(t);
    const double error = state.x - r;
    const double qTrue = std::sin(cfg.thetaTrue * state.x);
    const double qHat = std::sin(state.thetaHat * state.x);
    const double sensitivity = state.x * std::cos(state.thetaHat * state.x);

    const double control = rDot + cfg.a * state.x - qHat - cfg.k * error;
    const double xDot = -cfg.a * state.x + control + qTrue;
    const double rawThetaDot = cfg.gamma * error * sensitivity;
    const double thetaDot = projectRate(state.thetaHat, rawThetaDot, cfg);
    return {xDot, thetaDot};
}

State addScaled(const State& a, const State& b, double scale) {
    return {a.x + scale * b.x, a.thetaHat + scale * b.thetaHat};
}

State rk4Step(double t, const State& state, const Config& cfg) {
    const State k1 = dynamics(t, state, cfg);
    const State k2 = dynamics(t + 0.5 * cfg.dt, addScaled(state, k1, 0.5 * cfg.dt), cfg);
    const State k3 = dynamics(t + 0.5 * cfg.dt, addScaled(state, k2, 0.5 * cfg.dt), cfg);
    const State k4 = dynamics(t + cfg.dt, addScaled(state, k3, cfg.dt), cfg);

    State next{
        state.x + (cfg.dt / 6.0) * (k1.x + 2.0 * k2.x + 2.0 * k3.x + k4.x),
        state.thetaHat + (cfg.dt / 6.0) *
            (k1.thetaHat + 2.0 * k2.thetaHat + 2.0 * k3.thetaHat + k4.thetaHat)
    };
    next.thetaHat = std::clamp(next.thetaHat, cfg.thetaMin, cfg.thetaMax);
    return next;
}

std::string tagFor(double value) {
    std::ostringstream out;
    out << std::fixed << std::setprecision(1) << value;
    std::string tag = out.str();
    std::replace(tag.begin(), tag.end(), '-', 'm');
    std::replace(tag.begin(), tag.end(), '.', 'p');
    return tag;
}

std::vector<Sample> simulate(double thetaHat0, const Config& cfg) {
    const int steps = static_cast<int>(std::llround(cfg.tFinal / cfg.dt));
    std::vector<Sample> samples;
    samples.reserve(static_cast<std::size_t>(steps + 1));
    State state{0.0, thetaHat0};

    for (int i = 0; i <= steps; ++i) {
        const double t = i * cfg.dt;
        const auto [r, rDot] = referenceSignal(t);
        const double error = state.x - r;
        const double control = rDot + cfg.a * state.x
            - std::sin(state.thetaHat * state.x) - cfg.k * error;
        samples.push_back({t, state.x, r, error, state.thetaHat, control});
        if (i < steps) state = rk4Step(t, state, cfg);
    }
    return samples;
}

void writeCsv(const std::vector<Sample>& samples, const std::string& filename) {
    std::ofstream file(filename);
    if (!file) throw std::runtime_error("Cannot open " + filename);
    file << "time,x,reference,error,theta_hat,control\n";
    file << std::setprecision(12);
    for (const auto& s : samples) {
        file << s.time << ',' << s.x << ',' << s.reference << ',' << s.error << ','
             << s.thetaHat << ',' << s.control << '\n';
    }
}

int main() {
    try {
        const Config cfg;
        const std::array<double, 3> initialEstimates{-2.5, 0.0, 2.5};

        for (double theta0 : initialEstimates) {
            const auto samples = simulate(theta0, cfg);
            double sumSquares = 0.0;
            for (const auto& s : samples) sumSquares += s.error * s.error;
            const double rms = std::sqrt(sumSquares / static_cast<double>(samples.size()));
            const double thetaFinal = samples.back().thetaHat;

            const std::string filename =
                "Chapter18_Lesson1_theta0_" + tagFor(theta0) + ".csv";
            writeCsv(samples, filename);
            std::cout << std::fixed << std::setprecision(6)
                      << "thetaHat(0)=" << std::setw(8) << theta0
                      << " | thetaHat(T)=" << std::setw(10) << thetaFinal
                      << " | RMS(e)=" << rms << '\n';
        }
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
    return 0;
}

The C++17 implementation uses only the standard library and writes CSV output for plotting in any external tool.

Chapter18_Lesson1.java

// Chapter18_Lesson1.java
// Compile: javac Chapter18_Lesson1.java
// Run:     java Chapter18_Lesson1

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;

public final class Chapter18_Lesson1 {
    private static final class Config {
        final double a = 1.1;
        final double k = 3.0;
        final double gamma = 2.5;
        final double thetaTrue = 1.4;
        final double thetaMin = -4.0;
        final double thetaMax = 4.0;
        final double dt = 0.002;
        final double tFinal = 30.0;
    }

    private record State(double x, double thetaHat) {}

    private record Sample(
        double time,
        double x,
        double reference,
        double error,
        double thetaHat,
        double control
    ) {}

    private static double[] referenceSignal(double t) {
        double r = 0.8 * Math.sin(0.7 * t) + 0.3 * Math.sin(1.3 * t);
        double rDot = 0.56 * Math.cos(0.7 * t) + 0.39 * Math.cos(1.3 * t);
        return new double[] {r, rDot};
    }

    private static double projectRate(double thetaHat, double rawRate, Config cfg) {
        if (thetaHat <= cfg.thetaMin && rawRate < 0.0) return 0.0;
        if (thetaHat >= cfg.thetaMax && rawRate > 0.0) return 0.0;
        return rawRate;
    }

    private static State dynamics(double t, State state, Config cfg) {
        double[] ref = referenceSignal(t);
        double r = ref[0];
        double rDot = ref[1];
        double error = state.x() - r;
        double qTrue = Math.sin(cfg.thetaTrue * state.x());
        double qHat = Math.sin(state.thetaHat() * state.x());
        double sensitivity = state.x() * Math.cos(state.thetaHat() * state.x());

        double control = rDot + cfg.a * state.x() - qHat - cfg.k * error;
        double xDot = -cfg.a * state.x() + control + qTrue;
        double rawThetaDot = cfg.gamma * error * sensitivity;
        double thetaDot = projectRate(state.thetaHat(), rawThetaDot, cfg);
        return new State(xDot, thetaDot);
    }

    private static State addScaled(State a, State b, double scale) {
        return new State(a.x() + scale * b.x(), a.thetaHat() + scale * b.thetaHat());
    }

    private static State rk4Step(double t, State state, Config cfg) {
        State k1 = dynamics(t, state, cfg);
        State k2 = dynamics(t + 0.5 * cfg.dt, addScaled(state, k1, 0.5 * cfg.dt), cfg);
        State k3 = dynamics(t + 0.5 * cfg.dt, addScaled(state, k2, 0.5 * cfg.dt), cfg);
        State k4 = dynamics(t + cfg.dt, addScaled(state, k3, cfg.dt), cfg);

        double xNext = state.x() + (cfg.dt / 6.0)
            * (k1.x() + 2.0 * k2.x() + 2.0 * k3.x() + k4.x());
        double thetaNext = state.thetaHat() + (cfg.dt / 6.0)
            * (k1.thetaHat() + 2.0 * k2.thetaHat() + 2.0 * k3.thetaHat() + k4.thetaHat());
        thetaNext = Math.max(cfg.thetaMin, Math.min(cfg.thetaMax, thetaNext));
        return new State(xNext, thetaNext);
    }

    private static List<Sample> simulate(double thetaHat0, Config cfg) {
        int steps = (int) Math.round(cfg.tFinal / cfg.dt);
        List<Sample> samples = new ArrayList<>(steps + 1);
        State state = new State(0.0, thetaHat0);

        for (int i = 0; i <= steps; i++) {
            double t = i * cfg.dt;
            double[] ref = referenceSignal(t);
            double error = state.x() - ref[0];
            double control = ref[1] + cfg.a * state.x()
                - Math.sin(state.thetaHat() * state.x()) - cfg.k * error;
            samples.add(new Sample(t, state.x(), ref[0], error, state.thetaHat(), control));
            if (i < steps) state = rk4Step(t, state, cfg);
        }
        return samples;
    }

    private static String tagFor(double value) {
        return String.format(Locale.US, "%.1f", value)
            .replace("-", "m")
            .replace(".", "p");
    }

    private static void writeCsv(List<Sample> samples, String filename) throws IOException {
        try (BufferedWriter writer = Files.newBufferedWriter(Path.of(filename))) {
            writer.write("time,x,reference,error,theta_hat,control\n");
            for (Sample s : samples) {
                writer.write(String.format(
                    Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                    s.time(), s.x(), s.reference(), s.error(), s.thetaHat(), s.control()
                ));
            }
        }
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        Config cfg = new Config();
        double[] initialEstimates = {-2.5, 0.0, 2.5};

        try {
            for (double theta0 : initialEstimates) {
                List<Sample> samples = simulate(theta0, cfg);
                double sumSquares = 0.0;
                for (Sample s : samples) sumSquares += s.error() * s.error();
                double rms = Math.sqrt(sumSquares / samples.size());
                double thetaFinal = samples.get(samples.size() - 1).thetaHat();

                String filename = "Chapter18_Lesson1_theta0_" + tagFor(theta0) + ".csv";
                writeCsv(samples, filename);
                System.out.printf(
                    "thetaHat(0)=%8.3f | thetaHat(T)=%10.6f | RMS(e)=%.6f%n",
                    theta0, thetaFinal, rms
                );
            }
        } catch (IOException ex) {
            System.err.println("File error: " + ex.getMessage());
            System.exit(1);
        }
    }
}

The Java implementation uses records for immutable numerical samples and has no external dependencies.

Chapter18_Lesson1.m

% Chapter18_Lesson1.m
% Local Jacobian adaptation for q(x,theta) = sin(theta*x).
% This example is intentionally local: different initial estimates can
% converge differently or become trapped at the projection boundary.

clear; clc; close all;

cfg.a = 1.1;
cfg.k = 3.0;
cfg.gamma = 2.5;
cfg.thetaTrue = 1.4;
cfg.thetaMin = -4.0;
cfg.thetaMax = 4.0;
cfg.dt = 0.002;
cfg.tFinal = 30.0;
initialEstimates = [-2.5, 0.0, 2.5];

results = cell(numel(initialEstimates), 1);
for j = 1:numel(initialEstimates)
    results{j} = simulateCase(initialEstimates(j), cfg);
    tag = strrep(strrep(sprintf('%.1f', initialEstimates(j)), '-', 'm'), '.', 'p');
    fileName = ['Chapter18_Lesson1_theta0_' tag '.csv'];
    writetable(results{j}.table, fileName);
    fprintf('thetaHat(0)=%8.3f | thetaHat(T)=%10.6f | RMS(e)=%.6f\n', ...
        initialEstimates(j), results{j}.thetaFinal, results{j}.rmsError);
end

[~, bestIndex] = min(cellfun(@(s) s.rmsError, results));
best = results{bestIndex};

figure;
plot(best.table.time, best.table.reference, '--', 'LineWidth', 1.2); hold on;
plot(best.table.time, best.table.x, 'LineWidth', 1.2);
grid on; xlabel('Time [s]'); ylabel('State');
legend('reference r', 'state x', 'Location', 'best');
title('Tracking for the best tested initialization');

figure; hold on;
for j = 1:numel(results)
    plot(results{j}.table.time, results{j}.table.theta_hat, 'LineWidth', 1.2);
end
yline(cfg.thetaTrue, '--', 'true theta', 'LineWidth', 1.2);
grid on; xlabel('Time [s]'); ylabel('Parameter estimate');
legend(compose('thetaHat(0)=%.1f', initialEstimates), 'Location', 'best');
title('Initialization sensitivity of local Jacobian adaptation');

function result = simulateCase(thetaHat0, cfg)
    time = (0:cfg.dt:cfg.tFinal).';
    n = numel(time);
    state = [0.0; thetaHat0];

    x = zeros(n, 1);
    reference = zeros(n, 1);
    error = zeros(n, 1);
    theta_hat = zeros(n, 1);
    control = zeros(n, 1);

    for i = 1:n
        [r, rDot] = referenceSignal(time(i));
        x(i) = state(1);
        reference(i) = r;
        error(i) = state(1) - r;
        theta_hat(i) = state(2);
        control(i) = rDot + cfg.a * state(1) ...
            - sin(state(2) * state(1)) - cfg.k * error(i);

        if i < n
            state = rk4Step(time(i), state, cfg);
        end
    end

    result.table = table(time, x, reference, error, theta_hat, control);
    result.rmsError = sqrt(mean(error.^2));
    result.thetaFinal = theta_hat(end);
end

function next = rk4Step(t, state, cfg)
    k1 = dynamics(t, state, cfg);
    k2 = dynamics(t + 0.5 * cfg.dt, state + 0.5 * cfg.dt * k1, cfg);
    k3 = dynamics(t + 0.5 * cfg.dt, state + 0.5 * cfg.dt * k2, cfg);
    k4 = dynamics(t + cfg.dt, state + cfg.dt * k3, cfg);
    next = state + (cfg.dt / 6.0) * (k1 + 2*k2 + 2*k3 + k4);
    next(2) = min(max(next(2), cfg.thetaMin), cfg.thetaMax);
end

function dz = dynamics(t, state, cfg)
    x = state(1);
    thetaHat = state(2);
    [r, rDot] = referenceSignal(t);
    error = x - r;

    qTrue = sin(cfg.thetaTrue * x);
    qHat = sin(thetaHat * x);
    sensitivity = x * cos(thetaHat * x);

    control = rDot + cfg.a*x - qHat - cfg.k*error;
    xDot = -cfg.a*x + control + qTrue;

    rawThetaDot = cfg.gamma * error * sensitivity;
    thetaDot = projectRate(thetaHat, rawThetaDot, cfg);
    dz = [xDot; thetaDot];
end

function rate = projectRate(thetaHat, rawRate, cfg)
    if (thetaHat <= cfg.thetaMin && rawRate < 0) || ...
       (thetaHat >= cfg.thetaMax && rawRate > 0)
        rate = 0.0;
    else
        rate = rawRate;
    end
end

function [r, rDot] = referenceSignal(t)
    r = 0.8*sin(0.7*t) + 0.3*sin(1.3*t);
    rDot = 0.56*cos(0.7*t) + 0.39*cos(1.3*t);
end

The next script creates a Simulink block diagram using standard continuous-time, mathematical, signal-routing, and sink blocks.

Chapter18_Lesson1_Simulink.m

% Chapter18_Lesson1_Simulink.m
% Programmatically builds a continuous-time Simulink model for the lesson.
% Requires Simulink. Run this script, then press Run in the generated model.

clear; clc;
model = 'Chapter18_Lesson1_Model';
if bdIsLoaded(model)
    close_system(model, 0);
end
if exist([model '.slx'], 'file')
    delete([model '.slx']);
end

new_system(model);
open_system(model);
set_param(model, 'StopTime', '30', 'Solver', 'ode45');

% Reference r and derivative rDot.
add_block('simulink/Sources/Sine Wave', [model '/r1'], ...
    'Amplitude', '0.8', 'Frequency', '0.7', 'Position', [40 40 90 70]);
add_block('simulink/Sources/Sine Wave', [model '/r2'], ...
    'Amplitude', '0.3', 'Frequency', '1.3', 'Position', [40 90 90 120]);
add_block('simulink/Math Operations/Sum', [model '/r'], ...
    'Inputs', '++', 'Position', [130 50 160 110]);

add_block('simulink/Sources/Sine Wave', [model '/rd1'], ...
    'Amplitude', '0.56', 'Frequency', '0.7', 'Phase', 'pi/2', ...
    'Position', [40 160 90 190]);
add_block('simulink/Sources/Sine Wave', [model '/rd2'], ...
    'Amplitude', '0.39', 'Frequency', '1.3', 'Phase', 'pi/2', ...
    'Position', [40 210 90 240]);
add_block('simulink/Math Operations/Sum', [model '/rDot'], ...
    'Inputs', '++', 'Position', [130 170 160 230]);

% Plant state and tracking error.
add_block('simulink/Continuous/Integrator', [model '/x'], ...
    'InitialCondition', '0', 'Position', [680 90 710 120]);
add_block('simulink/Math Operations/Sum', [model '/error'], ...
    'Inputs', '+-', 'Position', [230 45 260 105]);

% True and estimated nonlinear terms.
add_block('simulink/Math Operations/Fcn', [model '/qTrue'], ...
    'Expr', 'sin(1.4*u)', 'Position', [450 210 540 240]);
add_block('simulink/Signal Routing/Mux', [model '/Mux_qHat'], ...
    'Inputs', '2', 'Position', [330 250 335 305]);
add_block('simulink/Math Operations/Fcn', [model '/qHat'], ...
    'Expr', 'sin(u(1)*u(2))', 'Position', [370 255 470 285]);

% Controller: u = rDot + a*x - qHat - k*e.
add_block('simulink/Math Operations/Gain', [model '/aGain'], ...
    'Gain', '1.1', 'Position', [300 125 350 155]);
add_block('simulink/Math Operations/Gain', [model '/kGain'], ...
    'Gain', '3.0', 'Position', [300 65 350 95]);
add_block('simulink/Math Operations/Sum', [model '/controller'], ...
    'Inputs', '++--', 'Position', [500 70 530 160]);

% Plant derivative: xDot = -a*x + u + qTrue.
add_block('simulink/Math Operations/Gain', [model '/minusAGain'], ...
    'Gain', '-1.1', 'Position', [480 175 530 205]);
add_block('simulink/Math Operations/Sum', [model '/plantSum'], ...
    'Inputs', '+++', 'Position', [600 85 630 175]);

% Sensitivity psi = x*cos(thetaHat*x), adaptation, and projection by saturation.
add_block('simulink/Signal Routing/Mux', [model '/Mux_sensitivity'], ...
    'Inputs', '2', 'Position', [330 340 335 395]);
add_block('simulink/Math Operations/Fcn', [model '/sensitivity'], ...
    'Expr', 'u(1)*cos(u(2)*u(1))', 'Position', [370 345 500 375]);
add_block('simulink/Math Operations/Product', [model '/eTimesPsi'], ...
    'Position', [530 325 560 380]);
add_block('simulink/Math Operations/Gain', [model '/gamma'], ...
    'Gain', '2.5', 'Position', [590 335 640 365]);
add_block('simulink/Continuous/Integrator Limited', [model '/thetaHat'], ...
    'InitialCondition', '0', 'UpperSaturationLimit', '4', ...
    'LowerSaturationLimit', '-4', 'Position', [680 330 720 365]);

% Scopes.
add_block('simulink/Signal Routing/Mux', [model '/Mux_tracking'], ...
    'Inputs', '2', 'Position', [760 45 765 105]);
add_block('simulink/Sinks/Scope', [model '/Tracking Scope'], ...
    'Position', [805 45 855 105]);
add_block('simulink/Sinks/Scope', [model '/Parameter Scope'], ...
    'Position', [805 325 855 375]);

% Connections.
add_line(model, 'r1/1', 'r/1');
add_line(model, 'r2/1', 'r/2');
add_line(model, 'rd1/1', 'rDot/1');
add_line(model, 'rd2/1', 'rDot/2');
add_line(model, 'x/1', 'error/1');
add_line(model, 'r/1', 'error/2');
add_line(model, 'error/1', 'kGain/1');
add_line(model, 'x/1', 'aGain/1');
add_line(model, 'rDot/1', 'controller/1');
add_line(model, 'aGain/1', 'controller/2');
add_line(model, 'qHat/1', 'controller/3');
add_line(model, 'kGain/1', 'controller/4');
add_line(model, 'x/1', 'minusAGain/1');
add_line(model, 'controller/1', 'plantSum/1');
add_line(model, 'minusAGain/1', 'plantSum/2');
add_line(model, 'qTrue/1', 'plantSum/3');
add_line(model, 'plantSum/1', 'x/1');
add_line(model, 'x/1', 'qTrue/1');
add_line(model, 'thetaHat/1', 'Mux_qHat/1');
add_line(model, 'x/1', 'Mux_qHat/2');
add_line(model, 'Mux_qHat/1', 'qHat/1');
add_line(model, 'x/1', 'Mux_sensitivity/1');
add_line(model, 'thetaHat/1', 'Mux_sensitivity/2');
add_line(model, 'Mux_sensitivity/1', 'sensitivity/1');
add_line(model, 'error/1', 'eTimesPsi/1');
add_line(model, 'sensitivity/1', 'eTimesPsi/2');
add_line(model, 'eTimesPsi/1', 'gamma/1');
add_line(model, 'gamma/1', 'thetaHat/1');
add_line(model, 'r/1', 'Mux_tracking/1');
add_line(model, 'x/1', 'Mux_tracking/2');
add_line(model, 'Mux_tracking/1', 'Tracking Scope/1');
add_line(model, 'thetaHat/1', 'Parameter Scope/1');

save_system(model);
Simulink.BlockDiagram.arrangeSystem(model);
open_system(model);
disp(['Created ' model '.slx']);

The Simulink limited integrator implements the parameter interval constraint. Its hard saturation is a practical approximation to the directional projection used in the numerical scripts, so boundary behavior can differ slightly.

Chapter18_Lesson1.nb

Notebook[{
 Cell["Chapter 18, Lesson 1: Nonlinear Parameterization and Challenges", "Title"],
 Cell["Local Jacobian adaptation for q(x, theta) = Sin[theta x]. Different initial estimates demonstrate the local nature of the update.", "Text"],
 Cell[BoxData[
  RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"a", "=", "1.1"}], ";",
   RowBox[{"k", "=", "3.0"}], ";",
   RowBox[{"gamma", "=", "2.5"}], ";",
   RowBox[{"thetaTrue", "=", "1.4"}], ";",
   RowBox[{"thetaMin", "=", RowBox[{"-", "4.0"}]}], ";",
   RowBox[{"thetaMax", "=", "4.0"}], ";",
   RowBox[{"tFinal", "=", "30.0"}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"r", "[", "t_", "]"}], ":=",
   RowBox[{RowBox[{"0.8", " ", RowBox[{"Sin", "[", RowBox[{"0.7", " t"}], "]"}]}], "+",
           RowBox[{"0.3", " ", RowBox[{"Sin", "[", RowBox[{"1.3", " t"}], "]"}]}]}], ";",
   RowBox[{"rd", "[", "t_", "]"}], ":=",
   RowBox[{RowBox[{"0.56", " ", RowBox[{"Cos", "[", RowBox[{"0.7", " t"}], "]"}]}], "+",
           RowBox[{"0.39", " ", RowBox[{"Cos", "[", RowBox[{"1.3", " t"}], "]"}]}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"projectedRate", "[", RowBox[{"th_", ",", "raw_"}], "]"}], ":=",
   RowBox[{"Which", "[",
    RowBox[{
     RowBox[{"th", "<=", "thetaMin", "&&", "raw", "<", "0"}], ",", "0",
     ",", RowBox[{"th", ">=", "thetaMax", "&&", "raw", ">", "0"}], ",", "0",
     ",", "True", ",", "raw"}], "]"}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"simulate", "[", "theta0_", "]"}], ":=",
   RowBox[{"Module", "[",
    RowBox[{
     RowBox[{"{", RowBox[{"sol", ",", "x", ",", "th", ",", "e", ",", "psi"}], "}"}], ",",
     RowBox[{
      RowBox[{"sol", "=",
       RowBox[{"NDSolveValue", "[",
        RowBox[{
         RowBox[{"{",
          RowBox[{
           RowBox[{RowBox[{"x", "'"}], "[", "t", "]"}], "==",
            RowBox[{RowBox[{"-", "k"}], " ", RowBox[{"(", RowBox[{RowBox[{"x", "[", "t", "]"}], "-", RowBox[{"r", "[", "t", "]"}]}], ")"}], "+",
             RowBox[{"Sin", "[", RowBox[{"thetaTrue", " ", RowBox[{"x", "[", "t", "]"}]}], "]"}], "-",
             RowBox[{"Sin", "[", RowBox[{RowBox[{"th", "[", "t", "]"}], " ", RowBox[{"x", "[", "t", "]"}]}], "]"}]}], ",",
           RowBox[{RowBox[{"th", "'"}], "[", "t", "]"}], "==",
            RowBox[{"projectedRate", "[",
             RowBox[{RowBox[{"th", "[", "t", "]"}], ",",
              RowBox[{"gamma", " ", RowBox[{"(", RowBox[{RowBox[{"x", "[", "t", "]"}], "-", RowBox[{"r", "[", "t", "]"}]}], ")"}], " ",
               RowBox[{"x", "[", "t", "]"}], " ",
               RowBox[{"Cos", "[", RowBox[{RowBox[{"th", "[", "t", "]"}], " ", RowBox[{"x", "[", "t", "]"}]}], "]"}]}]}], "]"}], ",",
           RowBox[{RowBox[{"x", "[", "0", "]"}], "==", "0"}], ",",
           RowBox[{RowBox[{"th", "[", "0", "]"}], "==", "theta0"}]
          }], "}"}], ",", RowBox[{"{", RowBox[{"x", ",", "th"}], "}"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "tFinal"}], "}"}],
         ",", RowBox[{"MaxStepFraction", "->", RowBox[{"1", "/", "500"}]}]}], "]"}]}], ";",
      RowBox[{"<|", RowBox[{
       RowBox[{"\"theta0\"", "->", "theta0"}], ",",
       RowBox[{"\"x\"", "->", RowBox[{"sol", "[[", "1", "]]"}]}], ",",
       RowBox[{"\"theta\"", "->", RowBox[{"sol", "[[", "2", "]]"}]}]
      }], "|>"}]
     }]
    }], "]"}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"initials", "=", RowBox[{"{", RowBox[{RowBox[{"-", "2.5"}], ",", "0.0", ",", "2.5"}], "}"}]}], ";",
   RowBox[{"runs", "=", RowBox[{"simulate", "/@", "initials"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{"Plot", "[",
   RowBox[{
    RowBox[{"Evaluate", "[", RowBox[{RowBox[{"Through", "[", RowBox[{RowBox[{"runs", "[[", RowBox[{"All", ",", "\"theta\""}], "]]"}], "[", "t", "]"}], "]"}], "~", "Join", "~", RowBox[{"{", "thetaTrue", "}"}]}], "]"}], ",",
    RowBox[{"{", RowBox[{"t", ",", "0", ",", "tFinal"}], "}"}], ",",
    RowBox[{"PlotLegends", "->", RowBox[{"{", RowBox[{"\"theta0=-2.5\"", ",", "\"theta0=0\"", ",", "\"theta0=2.5\"", ",", "\"true theta\""}], "}"}]}], ",",
    RowBox[{"PlotRange", "->", "All"}], ",", RowBox[{"AxesLabel", "->", RowBox[{"{", RowBox[{"\"t\"", ",", "\"thetaHat\""}], "}"}]}]
   }], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"Table", "[",
   RowBox[{
    RowBox[{"With", "[",
     RowBox[{
      RowBox[{"{", RowBox[{"run", "=", RowBox[{"runs", "[[", "i", "]]"}]}], "}"}], ",",
      RowBox[{"Export", "[",
       RowBox[{
        RowBox[{"\"Chapter18_Lesson1_theta0_\"", "<>", RowBox[{"StringReplace", "[", RowBox[{RowBox[{"ToString", "[", RowBox[{"initials", "[[", "i", "]]"}], "]"}], ",", RowBox[{"{", RowBox[{RowBox[{"\"-\"", "->", "\"m\""}], ",", RowBox[{"\".\"", "->", "\"p\""}]}], "}"}]}], "]"}], "<>", "\".csv\""}], ",",
        RowBox[{"Prepend", "[",
         RowBox[{
          RowBox[{"Table", "[",
           RowBox[{
            RowBox[{"{", RowBox[{"tt", ",", RowBox[{RowBox[{"run", "[", "\"x\"", "]"}], "[", "tt", "]"}], ",", RowBox[{"r", "[", "tt", "]"}], ",", RowBox[{RowBox[{"run", "[", "\"x\"", "]"}], "[", "tt", "]"}], "-", RowBox[{"r", "[", "tt", "]"}], ",", RowBox[{RowBox[{"run", "[", "\"theta\"", "]"}], "[", "tt", "]"}]}], "}"}], ",",
            RowBox[{"{", RowBox[{"tt", ",", "0", ",", "tFinal", ",", "0.01"}], "}"}]}], "]"}], ",",
          RowBox[{"{", RowBox[{"\"time\"", ",", "\"x\"", ",", "\"reference\"", ",", "\"error\"", ",", "\"theta_hat\""}], "}"}]}], "]"}], ",", "\"CSV\""}], "]"}]
     }], "]"}], ",", RowBox[{"{", RowBox[{"i", ",", RowBox[{"Length", "[", "runs", "]"}]}], "}"}]}], "]"}]], "Input"]
}, WindowSize->{1100, 800}, StyleDefinitions->"Default.nb"]

Save the displayed notebook expression as Chapter18_Lesson1.nb and open it in Wolfram Mathematica. The downloadable package already includes this file.

The simulations should be interpreted through four separate questions:

1. Boundedness: do \(x\), \(u\), and \(\hat\theta\) remain bounded?
2. Tracking: is the RMS tracking error small?
3. Identification: does \(\hat\theta\) approach the true parameter?
4. Robustness of outcome: are the conclusions insensitive to initialization?

In a nonlinear parameterization, these properties need not coincide. A controller can track well because feedback suppresses the error while the estimate remains inaccurate. Conversely, aggressive adaptation can disturb tracking when sensitivity changes sign or becomes large. The projection interval prevents unbounded drift but cannot create missing information or remove false stationary points.

Problem 1: Classifying Parameterizations

Classify each uncertainty with respect to the physical parameter vector:

\[ F_1=\theta_1x^3+\theta_2\sin x, \quad F_2=\sin(\theta x), \quad F_3=\theta^2x, \quad F_4=\frac{\theta_1x}{1+\theta_2x^2}. \]

Solution:

\(F_1=[x^3\ \ \sin x]\boldsymbol{\theta}\), so it is linear in the parameters despite being nonlinear in the state. \(F_2\) is nonlinearly parameterized because its sensitivity \(x\cos(\theta x)\) depends on the unknown parameter. \(F_3\) is nonlinear in physical \(\theta\), although the reparameterization \(\alpha=\theta^2\) makes it linear in \(\alpha\). \(F_4\) is linear in \(\theta_1\) only when \(\theta_2\) is known; with both unknown it is nonlinearly parameterized.

Problem 2: Exact Cancellation for a Linear-in-Parameter Model

For

\[ \dot e=-ke-\boldsymbol{\phi}^{\mathsf T}(t) \tilde{\boldsymbol{\theta}}, \]

prove that \(\dot{\hat{\boldsymbol{\theta}}}= \boldsymbol{\Gamma}\boldsymbol{\phi}e\) gives \(\dot V=-ke^2\).

Solution:

\[ V=\frac12e^2+\frac12\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta}}. \]

\[ \dot V=e(-ke-\boldsymbol{\phi}^{\mathsf T} \tilde{\boldsymbol{\theta}})+ \tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1} (\boldsymbol{\Gamma}\boldsymbol{\phi}e). \]

Since the last two scalar cross terms are negatives of one another, they cancel, yielding \(\dot V=-ke^2\le0\).

Problem 3: Taylor Remainder for a Nonlinear Parameterization

Let \(F(\theta)=\sin(\theta x)\) with fixed measured \(x\). Expand \(F(\theta)-F(\hat\theta)\) around \(\hat\theta\) and bound the remainder.

Solution:

With \(\tilde\theta=\hat\theta-\theta\),

\[ \sin(\theta x)-\sin(\hat\theta x)= -x\cos(\hat\theta x)\tilde\theta+R_2. \]

The second derivative with respect to \(\theta\) is \(-x^2\sin(\theta x)\), whose magnitude is bounded by \(x^2\). Hence

\[ |R_2|\le\frac{x^2}{2}|\tilde\theta|^2. \]

The local Jacobian update cancels the term linear in \(\tilde\theta\), but the quadratic remainder persists in \(\dot V\).

Problem 4: Structural Non-Identifiability

Show that \(F(x,\theta)=\theta^2x\) cannot globally identify the sign of \(\theta\), regardless of the excitation in \(x(t)\).

Solution:

\[ F(x(t),\theta)=\theta^2x(t)=(-\theta)^2x(t)=F(x(t),-\theta) \]

for every trajectory and every time. Thus, the input-output map is invariant under \(\theta\mapsto-\theta\). No excitation can distinguish two parameters that produce exactly the same model. Only \(\alpha=\theta^2\) is identifiable.

Problem 5: Sensitivity Gramian

For \(F(x,\theta)=\sin(\theta x)\), derive the scalar information Gramian and state two conditions under which it is small.

Solution:

\[ S(t,\hat\theta)=x(t)\cos(\hat\theta x(t)), \qquad G_T=\int_0^T x^2(t)\cos^2(\hat\theta x(t))dt. \]

The Gramian is small when the state amplitude is small, or when the trajectory remains near sensitivity-null sets satisfying \(\hat\theta x(t)\approx\pi/2+k\pi\). Large state amplitude alone does not guarantee useful parameter information.

Problem 6: Local Lyapunov Region

Suppose \(|R_2|\le(L_F/2)|\tilde\theta|^2\). Find a sufficient relation among \(|e|\), \(|\tilde\theta|\), \(k\), and \(L_F\) that ensures \(\dot V\le0\).

Solution:

\[ \dot V\le-k e^2+\frac{L_F}{2}|e||\tilde\theta|^2. \]

For \(e\neq0\), a sufficient condition is

\[ \frac{L_F}{2}|\tilde\theta|^2\le k|e|. \]

This condition is state dependent and is not invariant by itself. It illustrates why the naive Jacobian law usually provides, at best, a local argument that requires additional bounds or a more sophisticated design.

Problem 7: Overparameterization Constraint

For \(\alpha_1=\theta_1\), \(\alpha_2=\theta_1\theta_2\), and \(\alpha_3=\theta_2^2\), derive the constraint among the augmented parameters.

Solution:

\[ \alpha_2^2=(\theta_1\theta_2)^2= \theta_1^2\theta_2^2=\alpha_1^2\alpha_3, \qquad \alpha_3\ge0. \]

Therefore, an arbitrary estimate in \(\mathbb{R}^3\) need not correspond to any physical pair \((\theta_1,\theta_2)\). A constrained estimator or a direct nonlinear parameter estimator is needed if physical recovery matters.

Linear-in-the-parameters models permit exact separation of uncertainty mismatch into a known regressor and a parameter error, enabling standard Lyapunov cancellation. Nonlinear parameterizations replace this identity with an unknown averaged Jacobian or a local Jacobian plus a higher-order remainder. The resulting challenges include sign-indefinite Lyapunov terms, local sensitivity loss, global aliases, nonconvex prediction costs, initialization dependence, and physically inconsistent overparameterization. Successful nonlinear adaptive control therefore requires structural assumptions beyond smoothness, together with explicit statements about identifiability and the domain of the stability guarantee.