After completing this lesson, students should be able to:

• recognize a nonlinear system that is linear in its unknown parameters;
• construct the nonlinear regressor vector and parameter vector;
• state and verify the matched-uncertainty condition;
• derive the nonlinear MRAC tracking-error dynamics;
• design a gradient adaptive law by Lyapunov cross-term cancellation;
• prove boundedness and asymptotic tracking using Barbalat's lemma;
• explain why zero tracking error does not automatically imply exact parameter identification;
• implement the same controller in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The lesson uses material already established in earlier chapters: reference models, matching conditions, quadratic Lyapunov functions, gradient adaptation, normalization, projection, persistent excitation, and adaptive backstepping. It does not require new observer theory.

A system is linear in the parameters when the unknown constants multiply known functions of measured signals. The functions may be highly nonlinear in the state and input. A general control-affine model is

\[ \dot{\mathbf{x} } = \mathbf{f}_0(\mathbf{x}) + \mathbf{g}(\mathbf{x})u + \mathbf{Y}(\mathbf{x})\boldsymbol{\theta}, \qquad \mathbf{x}\in\mathbb{R}^n, \quad \boldsymbol{\theta}\in\mathbb{R}^p . \]

Here \( \mathbf{Y}(\mathbf{x}) \) is a known regressor matrix and \( \boldsymbol{\theta} \) is an unknown, constant parameter vector. The map \( \mathbf{Y}(\mathbf{x})\boldsymbol{\theta} \) is linear with respect to \( \boldsymbol{\theta} \), although it can be nonlinear with respect to \( \mathbf{x} \).

For a scalar uncertain term, the same structure is written as

\[ \Delta(\mathbf{x}) = \boldsymbol{\theta}^{T}\boldsymbol{\phi}(\mathbf{x}) = \sum_{i=1}^p \theta_i \phi_i(\mathbf{x}), \]

where \( \boldsymbol{\phi}(\mathbf{x}) \) is the nonlinear regressor vector. Examples include

\[ \Delta(\mathbf{x}) = \theta_1\sin x_1 + \theta_2 x_1x_2 + \theta_3 x_2|x_2|, \]

and the mechanical expression

\[ \tau = m\,\ddot q_r + b\,\dot q + g\sin q = \begin{bmatrix} \ddot q_r & \dot q & \sin q \end{bmatrix} \begin{bmatrix} m \\ b \\ g \end{bmatrix}. \]

The second formula is nonlinear in the configuration variable \(q\), but linear in the physical parameters \(m,b,g\). This distinction is the key reason Lyapunov gradient laws remain constructive for many nonlinear plants.

2.1 What does not satisfy the assumption?

Terms such as \( \sin(\theta x) \), \( e^{\theta x} \), or \( \theta_1\theta_2\phi(x) \) are not linear in the original unknown parameter vector. They may require reparameterization, overparameterization, bounds, convexity, immersion and invariance, neural approximators, or another nonlinear adaptive design.

The cleanest nonlinear MRAC extension occurs when the uncertain term acts through the same channel as the control. Consider

\[ \dot{\mathbf{x} } = \mathbf{A}\mathbf{x} + \mathbf{b} \left[ u + \boldsymbol{\theta}^{T}\boldsymbol{\phi}(\mathbf{x}) \right]. \]

The uncertainty is matched because both \(u\) and \(\boldsymbol{\theta}^{T}\boldsymbol{\phi}(\mathbf{x})\) multiply the same vector \(\mathbf{b}\). Select a stable reference model

\[ \dot{\mathbf{x} }_m = \mathbf{A}_m\mathbf{x}_m + \mathbf{b}_m r, \qquad \mathbf{A}_m \text{ Hurwitz}. \]

Assume constant ideal gains \(\mathbf{k}_x^\ast\) and \(k_r^\ast\) satisfy the linear matching equations

\[ \mathbf{A}+\mathbf{b}\mathbf{k}_x^{\ast T}=\mathbf{A}_m, \qquad \mathbf{b}k_r^\ast=\mathbf{b}_m. \]

The nonlinear uncertainty is then cancelled by an adaptive feedforward term:

\[ u = \mathbf{k}_x^{\ast T}\mathbf{x} + k_r^\ast r - \hat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}(\mathbf{x}). \]

flowchart TD
  R["Command r"] --> RM["Reference model"]
  RM --> XM["Desired state xm"]
  R --> C["MRAC control law"]
  X["Plant state x"] --> C
  X --> REG["Nonlinear regressor phi(x)"]
  REG --> C
  XM --> E["Tracking error e = x - xm"]
  X --> E
  E --> AD["Adaptive law"]
  REG --> AD
  AD --> TH["Parameter estimate thetaHat"]
  TH --> C
  C --> U["Control u"]
  U --> P["Nonlinear plant"]
  P --> X
        

Define the tracking error and parameter-estimation error as

\[ \mathbf{e} = \mathbf{x}-\mathbf{x}_m, \qquad \tilde{\boldsymbol{\theta} } = \hat{\boldsymbol{\theta} }-\boldsymbol{\theta}. \]

Substitute the controller into the plant:

\[ \begin{aligned} \dot{\mathbf{x} } &= \mathbf{A}\mathbf{x} + \mathbf{b} \left[ \mathbf{k}_x^{\ast T}\mathbf{x} + k_r^\ast r - \hat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} + \boldsymbol{\theta}^{T}\boldsymbol{\phi} \right] \\ &= \mathbf{A}_m\mathbf{x} + \mathbf{b}_m r - \mathbf{b}\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\phi}(\mathbf{x}). \end{aligned} \]

Subtracting the reference-model dynamics gives

\[ \boxed{ \dot{\mathbf{e} } = \mathbf{A}_m\mathbf{e} - \mathbf{b}\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\phi}(\mathbf{x}) } \]

This has the same structural form as linear state-feedback MRAC, except that the regressor now depends nonlinearly on the plant state.

Since \(\mathbf{A}_m\) is Hurwitz, for every symmetric positive-definite matrix \(\mathbf{Q}\) there exists a unique symmetric positive-definite matrix \(\mathbf{P}\) satisfying

\[ \mathbf{A}_m^T\mathbf{P} + \mathbf{P}\mathbf{A}_m = -\mathbf{Q}. \]

Choose the composite Lyapunov function

\[ V = \mathbf{e}^T\mathbf{P}\mathbf{e} + \tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta} }, \qquad \boldsymbol{\Gamma} = \boldsymbol{\Gamma}^T > \mathbf{0}. \]

Because the true parameter vector is constant, \(\dot{\tilde{\boldsymbol{\theta} } }= \dot{\hat{\boldsymbol{\theta} } }\). Differentiating \(V\) gives

\[ \begin{aligned} \dot V &= \dot{\mathbf{e} }^T\mathbf{P}\mathbf{e} + \mathbf{e}^T\mathbf{P}\dot{\mathbf{e} } + 2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta} } } \\ &= -\mathbf{e}^T\mathbf{Q}\mathbf{e} - 2\mathbf{e}^T\mathbf{P}\mathbf{b} \tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}(\mathbf{x}) + 2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta} } }. \end{aligned} \]

Select the adaptive law

\[ \boxed{ \dot{\hat{\boldsymbol{\theta} } } = \boldsymbol{\Gamma} \boldsymbol{\phi}(\mathbf{x}) \mathbf{e}^T\mathbf{P}\mathbf{b} } \]

so that

\[ 2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta} } } = 2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\phi}(\mathbf{x}) \mathbf{e}^T\mathbf{P}\mathbf{b}, \]

which cancels the mixed error-parameter term. Therefore,

\[ \boxed{ \dot V = -\mathbf{e}^T\mathbf{Q}\mathbf{e} \le 0 } \]

5.1 Consequences

Since \(V(t)\le V(0)\), both \(\mathbf{e}(t)\) and \(\tilde{\boldsymbol{\theta} }(t)\) are bounded. Moreover,

\[ \int_0^\infty \mathbf{e}^T(\tau)\mathbf{Q}\mathbf{e}(\tau)\,d\tau \le V(0), \]

so \(\mathbf{e}\in L_2\cap L_\infty\). If the command \(r\) and reference-model state are bounded, and if \(\boldsymbol{\phi}(\mathbf{x})\) remains bounded on bounded sets, then \(\dot{\mathbf{e} }\) is bounded. Hence \(\mathbf{e}\) is uniformly continuous. Barbalat's lemma yields

\[ \lim_{t→\infty}\mathbf{e}(t)=\mathbf{0}. \]

Thus the adaptive controller guarantees asymptotic model following under the stated boundedness assumptions.

The Lyapunov proof establishes \(\mathbf{e}(t)→\mathbf{0}\), but does not by itself establish \(\tilde{\boldsymbol{\theta} }(t)→\mathbf{0}\). Once the tracking error becomes small, the update signal \(\mathbf{e}^T\mathbf{P}\mathbf{b}\) also becomes small, so adaptation can stop at a parameter vector that produces correct behavior only along the experienced trajectory.

A standard sufficient excitation condition is that there exist constants \(T>0\) and \(\alpha>0\) such that

\[ \int_t^{t+T} \boldsymbol{\phi}(\mathbf{x}(\tau)) \boldsymbol{\phi}^T(\mathbf{x}(\tau)) \,d\tau \ge \alpha\mathbf{I} \qquad \text{for all } t\ge 0. \]

This is persistent excitation of the nonlinear regressor. Under suitable boundedness and identifiability conditions, PE can strengthen convergence from tracking alone to parameter convergence. Without PE, only the component of the parameter error visible through the trajectory is forced to vanish:

\[ \tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\phi}(\mathbf{x}(t)) → 0, \]

which need not imply \(\tilde{\boldsymbol{\theta} }→\mathbf{0}\).

6.1 Basis expansion and approximation error

An unknown nonlinear function may be approximated by known basis functions:

\[ \Delta(\mathbf{x}) = \boldsymbol{\theta}^{\ast T} \boldsymbol{\phi}(\mathbf{x}) + \varepsilon(\mathbf{x}). \]

If the residual \(\varepsilon(\mathbf{x})\) is nonzero, exact cancellation is impossible. Plain gradient MRAC may then exhibit parameter drift or only practical tracking. The robust modifications from Chapters 8 and 9—projection, normalization, dead zones, leakage, and sigma-modification—become relevant. Their use changes the strict asymptotic result into an ultimate-boundedness result unless additional structure is available.

Consider the nonlinear plant

\[ \begin{aligned} \dot x_1 &= x_2, \\ \dot x_2 &= \theta_1\sin x_1 + \theta_2 x_1x_2 + u. \end{aligned} \]

The uncertainty is linearly parameterized using

\[ \boldsymbol{\theta} = \begin{bmatrix}\theta_1\\\theta_2\end{bmatrix}, \qquad \boldsymbol{\phi}(\mathbf{x}) = \begin{bmatrix}\sin x_1\\x_1x_2\end{bmatrix}. \]

Select the reference model

\[ \begin{aligned} \dot x_{m1} &= x_{m2}, \\ \dot x_{m2} &= -4x_{m1}-3x_{m2}+4r. \end{aligned} \]

Thus

\[ \mathbf{A}_m = \begin{bmatrix}0&1\\-4&-3\end{bmatrix}, \qquad \mathbf{b} = \begin{bmatrix}0\\1\end{bmatrix}, \qquad \mathbf{b}_m = \begin{bmatrix}0\\4\end{bmatrix}. \]

The controller is

\[ u = -4x_1 - 3x_2 + 4r - \hat\theta_1\sin x_1 - \hat\theta_2x_1x_2. \]

With \(\mathbf{Q}=\mathbf{I}\), the Lyapunov equation gives

\[ \mathbf{P} = \begin{bmatrix} 29/24 & 1/8 \\ 1/8 & 5/24 \end{bmatrix}. \]

Because

\[ \mathbf{e}^T\mathbf{P}\mathbf{b} = \frac{1}{8}e_1 + \frac{5}{24}e_2, \]

the adaptive laws are

\[ \begin{aligned} \dot{\hat\theta}_1 &= \gamma_1\sin x_1 \left( \frac{1}{8}e_1+\frac{5}{24}e_2 \right), \\ \dot{\hat\theta}_2 &= \gamma_2x_1x_2 \left( \frac{1}{8}e_1+\frac{5}{24}e_2 \right). \end{aligned} \]

For the ideal continuous-time model, the closed-loop Lyapunov derivative is

\[ \dot V = -e_1^2-e_2^2. \]

The simulations below use \(\boldsymbol{\theta}=[1.2,-0.4]^T\), \(\boldsymbol{\Gamma}=\operatorname{diag}(8,8)\), and \(r(t)=0.8\sin(0.7t)+0.35\sin(1.9t)\).

The continuous-time proof assumes exact integration. Software must update the plant, reference model, and estimator consistently. A fourth-order Runge-Kutta integrator is used in the examples to reduce discretization error. In embedded control, the estimator is normally integrated at the same sample period as the controller or at a known submultiple.

flowchart TD
  S["Read x and command r"] --> M["Update reference model xm"]
  M --> E["Compute e = x - xm"]
  E --> G["Evaluate nonlinear regressor phi(x)"]
  G --> A["Compute thetaHatDot"]
  A --> I["Integrate parameter estimates"]
  I --> C["Compute adaptive control u"]
  C --> L["Apply limits and safety logic"]
  L --> P["Advance plant or send actuator command"]
  P --> N["Log x, xm, e, thetaHat, u"]
  N --> S
        

8.1 Discrete approximation

A forward-Euler estimator would use

\[ \hat{\boldsymbol{\theta} }_{k+1} = \hat{\boldsymbol{\theta} }_k + T_s \boldsymbol{\Gamma} \boldsymbol{\phi}(\mathbf{x}_k) \mathbf{e}_k^T\mathbf{P}\mathbf{b}. \]

Small sample time is necessary because the continuous-time Lyapunov proof does not directly prove stability of this discrete implementation. Saturation must also be treated explicitly: clipping the control without modifying the adaptive law breaks the exact cancellation used in the proof.

8.2 Recommended software checks

• verify that every state, regressor, estimate, and control value is finite;
• log the norm of the tracking error and the parameter estimates;
• use consistent units and scaling for each regressor component;
• compare decreasing sample periods to detect numerical artifacts;
• add projection or safety bounds before deploying on hardware;
• do not infer parameter correctness solely from good tracking.

This implementation uses NumPy for vector operations and a from-scratch RK4 integrator. It writes all trajectories to CSV.

Chapter18_Lesson2.py

"""
Chapter18_Lesson2.py
Linear-in-the-parameters nonlinear MRAC example.

Plant:
    x1_dot = x2
    x2_dot = theta1*sin(x1) + theta2*x1*x2 + u

Reference model:
    xm1_dot = xm2
    xm2_dot = -4*xm1 - 3*xm2 + 4*r

Controller:
    u = -4*x1 - 3*x2 + 4*r - theta_hat^T*phi(x)

Adaptive law:
    theta_hat_dot = Gamma*phi(x)*(e^T*P*b)
"""

from __future__ import annotations

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

import numpy as np


@dataclass(frozen=True)
class Config:
    dt: float = 0.001
    final_time: float = 20.0
    gamma1: float = 8.0
    gamma2: float = 8.0
    theta1: float = 1.2
    theta2: float = -0.4


def reference(t: float) -> float:
    """Bounded command with several frequencies to enrich the regressor."""
    return 0.8 * np.sin(0.7 * t) + 0.35 * np.sin(1.9 * t)


def derivative(t: float, z: np.ndarray, cfg: Config, p: np.ndarray) -> np.ndarray:
    """Combined plant, reference-model, and estimator dynamics."""
    x = z[0:2]
    xm = z[2:4]
    theta_hat = z[4:6]

    r = reference(t)
    phi = np.array([np.sin(x[0]), x[0] * x[1]], dtype=float)
    e = x - xm
    b = np.array([0.0, 1.0])
    sigma = float(e @ p @ b)

    u = -4.0 * x[0] - 3.0 * x[1] + 4.0 * r - float(theta_hat @ phi)

    x_dot = np.array(
        [
            x[1],
            cfg.theta1 * np.sin(x[0]) + cfg.theta2 * x[0] * x[1] + u,
        ],
        dtype=float,
    )
    xm_dot = np.array([xm[1], -4.0 * xm[0] - 3.0 * xm[1] + 4.0 * r])
    theta_hat_dot = np.array(
        [cfg.gamma1 * phi[0] * sigma, cfg.gamma2 * phi[1] * sigma]
    )

    return np.concatenate((x_dot, xm_dot, theta_hat_dot))


def rk4_step(
    f: Callable[[float, np.ndarray, Config, np.ndarray], np.ndarray],
    t: float,
    z: np.ndarray,
    h: float,
    cfg: Config,
    p: np.ndarray,
) -> np.ndarray:
    k1 = f(t, z, cfg, p)
    k2 = f(t + 0.5 * h, z + 0.5 * h * k1, cfg, p)
    k3 = f(t + 0.5 * h, z + 0.5 * h * k2, cfg, p)
    k4 = f(t + h, z + h * k3, cfg, p)
    return z + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def main() -> None:
    cfg = Config()

    # For A_m = [[0, 1], [-4, -3]], Q = I, solve A_m^T P + P A_m = -I.
    p = np.array([[1.2083333333333333, 0.125], [0.125, 0.2083333333333333]])

    steps = int(round(cfg.final_time / cfg.dt))
    z = np.array([0.8, -0.3, 0.0, 0.0, 0.0, 0.0], dtype=float)

    rows: list[list[float]] = []
    for k in range(steps + 1):
        t = k * cfg.dt
        x = z[0:2]
        xm = z[2:4]
        theta_hat = z[4:6]
        phi = np.array([np.sin(x[0]), x[0] * x[1]])
        u = -4.0 * x[0] - 3.0 * x[1] + 4.0 * reference(t) - float(theta_hat @ phi)
        e = x - xm
        rows.append(
            [
                t,
                x[0],
                x[1],
                xm[0],
                xm[1],
                e[0],
                e[1],
                theta_hat[0],
                theta_hat[1],
                u,
            ]
        )
        if k < steps:
            z = rk4_step(derivative, t, z, cfg.dt, cfg, p)
            if not np.all(np.isfinite(z)):
                raise FloatingPointError("Simulation became non-finite.")

    output = Path("Chapter18_Lesson2_results_python.csv")
    with output.open("w", newline="", encoding="utf-8") as file:
        writer = csv.writer(file)
        writer.writerow(
            [
                "t",
                "x1",
                "x2",
                "xm1",
                "xm2",
                "e1",
                "e2",
                "theta_hat1",
                "theta_hat2",
                "u",
            ]
        )
        writer.writerows(rows)

    final = rows[-1]
    print(f"Saved {len(rows)} samples to {output}")
    print(f"Final tracking-error norm: {np.hypot(final[5], final[6]):.6e}")
    print(f"Final parameter estimates: [{final[7]:.6f}, {final[8]:.6f}]")


if __name__ == "__main__":
    main()

The C++ version uses only the standard library, fixed-size arrays, explicit finite-value checks, and CSV output.

Chapter18_Lesson2.cpp

// Chapter18_Lesson2.cpp
// Linear-in-the-parameters nonlinear MRAC simulation using RK4.

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

using State = std::array<double, 6>;

struct Config {
    double dt = 0.001;
    double final_time = 20.0;
    double gamma1 = 8.0;
    double gamma2 = 8.0;
    double theta1 = 1.2;
    double theta2 = -0.4;
};

double reference_signal(double t) {
    return 0.8 * std::sin(0.7 * t) + 0.35 * std::sin(1.9 * t);
}

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

State derivative(double t, const State& z, const Config& cfg) {
    constexpr double p12 = 0.125;
    constexpr double p22 = 0.2083333333333333;

    const double x1 = z[0];
    const double x2 = z[1];
    const double xm1 = z[2];
    const double xm2 = z[3];
    const double theta_hat1 = z[4];
    const double theta_hat2 = z[5];

    const double r = reference_signal(t);
    const double phi1 = std::sin(x1);
    const double phi2 = x1 * x2;
    const double e1 = x1 - xm1;
    const double e2 = x2 - xm2;
    const double sigma = p12 * e1 + p22 * e2;

    const double u =
        -4.0 * x1 - 3.0 * x2 + 4.0 * r
        - theta_hat1 * phi1 - theta_hat2 * phi2;

    State dz{};
    dz[0] = x2;
    dz[1] = cfg.theta1 * phi1 + cfg.theta2 * phi2 + u;
    dz[2] = xm2;
    dz[3] = -4.0 * xm1 - 3.0 * xm2 + 4.0 * r;
    dz[4] = cfg.gamma1 * phi1 * sigma;
    dz[5] = cfg.gamma2 * phi2 * sigma;
    return dz;
}

State rk4_step(double t, const State& z, double h, const Config& cfg) {
    const State k1 = derivative(t, z, cfg);
    const State k2 = derivative(t + 0.5 * h, add_scaled(z, k1, 0.5 * h), cfg);
    const State k3 = derivative(t + 0.5 * h, add_scaled(z, k2, 0.5 * h), cfg);
    const State k4 = derivative(t + h, add_scaled(z, k3, h), cfg);

    State next{};
    for (std::size_t i = 0; i < z.size(); ++i) {
        next[i] = z[i] + (h / 6.0) *
            (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        if (!std::isfinite(next[i])) {
            throw std::runtime_error("Simulation became non-finite.");
        }
    }
    return next;
}

int main() {
    try {
        const Config cfg;
        const int steps = static_cast<int>(std::llround(cfg.final_time / cfg.dt));
        State z{0.8, -0.3, 0.0, 0.0, 0.0, 0.0};

        std::ofstream file("Chapter18_Lesson2_results_cpp.csv");
        if (!file) {
            throw std::runtime_error("Cannot open output CSV file.");
        }
        file << "t,x1,x2,xm1,xm2,e1,e2,theta_hat1,theta_hat2,u\n";
        file << std::setprecision(15);

        for (int k = 0; k <= steps; ++k) {
            const double t = k * cfg.dt;
            const double phi1 = std::sin(z[0]);
            const double phi2 = z[0] * z[1];
            const double u =
                -4.0 * z[0] - 3.0 * z[1] + 4.0 * reference_signal(t)
                - z[4] * phi1 - z[5] * phi2;
            const double e1 = z[0] - z[2];
            const double e2 = z[1] - z[3];

            file << t << ',' << z[0] << ',' << z[1] << ','
                 << z[2] << ',' << z[3] << ',' << e1 << ',' << e2 << ','
                 << z[4] << ',' << z[5] << ',' << u << '\n';

            if (k < steps) {
                z = rk4_step(t, z, cfg.dt, cfg);
            }
        }

        const double e1 = z[0] - z[2];
        const double e2 = z[1] - z[3];
        std::cout << "Final tracking-error norm: "
                  << std::hypot(e1, e2) << '\n';
        std::cout << "Final parameter estimates: ["
                  << z[4] << ", " << z[5] << "]\n";
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

The Java implementation mirrors the same equations and integration order, using only the Java standard library.

Chapter18_Lesson2.java

// Chapter18_Lesson2.java
// Linear-in-the-parameters nonlinear MRAC simulation using RK4.

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

public final class Chapter18_Lesson2 {
    private static final double DT = 0.001;
    private static final double FINAL_TIME = 20.0;
    private static final double GAMMA1 = 8.0;
    private static final double GAMMA2 = 8.0;
    private static final double THETA1 = 1.2;
    private static final double THETA2 = -0.4;

    private Chapter18_Lesson2() {
    }

    private static double referenceSignal(double t) {
        return 0.8 * Math.sin(0.7 * t) + 0.35 * Math.sin(1.9 * t);
    }

    private static double[] derivative(double t, double[] z) {
        final double p12 = 0.125;
        final double p22 = 0.2083333333333333;

        double x1 = z[0];
        double x2 = z[1];
        double xm1 = z[2];
        double xm2 = z[3];
        double thetaHat1 = z[4];
        double thetaHat2 = z[5];

        double r = referenceSignal(t);
        double phi1 = Math.sin(x1);
        double phi2 = x1 * x2;
        double e1 = x1 - xm1;
        double e2 = x2 - xm2;
        double sigma = p12 * e1 + p22 * e2;

        double u = -4.0 * x1 - 3.0 * x2 + 4.0 * r
                - thetaHat1 * phi1 - thetaHat2 * phi2;

        return new double[] {
            x2,
            THETA1 * phi1 + THETA2 * phi2 + u,
            xm2,
            -4.0 * xm1 - 3.0 * xm2 + 4.0 * r,
            GAMMA1 * phi1 * sigma,
            GAMMA2 * phi2 * sigma
        };
    }

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

    private static double[] rk4Step(double t, double[] z, double h) {
        double[] k1 = derivative(t, z);
        double[] k2 = derivative(t + 0.5 * h, addScaled(z, k1, 0.5 * h));
        double[] k3 = derivative(t + 0.5 * h, addScaled(z, k2, 0.5 * h));
        double[] k4 = derivative(t + h, addScaled(z, k3, h));

        double[] next = new double[z.length];
        for (int i = 0; i < z.length; i++) {
            next[i] = z[i] + (h / 6.0)
                    * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
            if (!Double.isFinite(next[i])) {
                throw new IllegalStateException("Simulation became non-finite.");
            }
        }
        return next;
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);
        int steps = (int) Math.round(FINAL_TIME / DT);
        double[] z = {0.8, -0.3, 0.0, 0.0, 0.0, 0.0};
        Path output = Path.of("Chapter18_Lesson2_results_java.csv");

        try (BufferedWriter writer = Files.newBufferedWriter(
                output, StandardCharsets.UTF_8)) {
            writer.write("t,x1,x2,xm1,xm2,e1,e2,theta_hat1,theta_hat2,u");
            writer.newLine();

            for (int k = 0; k <= steps; k++) {
                double t = k * DT;
                double phi1 = Math.sin(z[0]);
                double phi2 = z[0] * z[1];
                double u = -4.0 * z[0] - 3.0 * z[1]
                        + 4.0 * referenceSignal(t)
                        - z[4] * phi1 - z[5] * phi2;
                double e1 = z[0] - z[2];
                double e2 = z[1] - z[3];

                writer.write(String.format(
                        Locale.ROOT,
                        "%.15g,%.15g,%.15g,%.15g,%.15g,%.15g,%.15g,%.15g,%.15g,%.15g",
                        t, z[0], z[1], z[2], z[3], e1, e2, z[4], z[5], u));
                writer.newLine();

                if (k < steps) {
                    z = rk4Step(t, z, DT);
                }
            }
        } catch (IOException ex) {
            System.err.println("I/O error: " + ex.getMessage());
            System.exit(1);
        }

        double e1 = z[0] - z[2];
        double e2 = z[1] - z[3];
        System.out.printf(Locale.ROOT, "Final tracking-error norm: %.6e%n",
                Math.hypot(e1, e2));
        System.out.printf(Locale.ROOT, "Final parameter estimates: [%.6f, %.6f]%n",
                z[4], z[5]);
    }
}

The MATLAB script integrates the augmented closed-loop system, writes a table to CSV, and plots tracking and parameter estimates.

Chapter18_Lesson2.m

% Chapter18_Lesson2.m
% Linear-in-the-parameters nonlinear MRAC simulation.
clear; clc; close all;

dt = 1e-3;
tf = 20;
t = (0:dt:tf).';
z = zeros(numel(t), 6);
z(1, :) = [0.8, -0.3, 0, 0, 0, 0];

cfg.gamma = diag([8, 8]);
cfg.theta = [1.2; -0.4];
cfg.P = [1.2083333333333333, 0.125; ...
         0.125, 0.2083333333333333];

for k = 1:numel(t)-1
    zk = z(k, :).';
    h = dt;
    k1 = dynamics(t(k), zk, cfg);
    k2 = dynamics(t(k) + h/2, zk + h*k1/2, cfg);
    k3 = dynamics(t(k) + h/2, zk + h*k2/2, cfg);
    k4 = dynamics(t(k) + h, zk + h*k3, cfg);
    z(k+1, :) = (zk + h*(k1 + 2*k2 + 2*k3 + k4)/6).';

    if any(~isfinite(z(k+1, :)))
        error('Simulation became non-finite.');
    end
end

x = z(:, 1:2);
xm = z(:, 3:4);
e = x - xm;
thetaHat = z(:, 5:6);
u = zeros(numel(t), 1);

for k = 1:numel(t)
    r = referenceSignal(t(k));
    phi = [sin(x(k,1)); x(k,1)*x(k,2)];
    u(k) = -4*x(k,1) - 3*x(k,2) + 4*r - thetaHat(k,:)*phi;
end

result = table(t, x(:,1), x(:,2), xm(:,1), xm(:,2), e(:,1), e(:,2), ...
    thetaHat(:,1), thetaHat(:,2), u, ...
    'VariableNames', {'t','x1','x2','xm1','xm2','e1','e2', ...
    'theta_hat1','theta_hat2','u'});
writetable(result, 'Chapter18_Lesson2_results_matlab.csv');

fprintf('Final tracking-error norm: %.6e\n', norm(e(end,:)));
fprintf('Final parameter estimates: [%.6f, %.6f]\n', ...
    thetaHat(end,1), thetaHat(end,2));

figure;
plot(t, x(:,1), t, xm(:,1), '--', 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Position state');
legend('x_1', 'x_{m1}', 'Location', 'best');
title('Nonlinear MRAC tracking');

figure;
plot(t, thetaHat(:,1), t, thetaHat(:,2), 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Parameter estimate');
legend('\hat{\theta}_1', '\hat{\theta}_2', 'Location', 'best');
title('Adaptive parameters');

function dz = dynamics(t, z, cfg)
    x = z(1:2);
    xm = z(3:4);
    thetaHat = z(5:6);

    r = referenceSignal(t);
    phi = [sin(x(1)); x(1)*x(2)];
    e = x - xm;
    b = [0; 1];
    sigma = e.'*cfg.P*b;

    u = -4*x(1) - 3*x(2) + 4*r - thetaHat.'*phi;

    xDot = [x(2); cfg.theta.'*phi + u];
    xmDot = [xm(2); -4*xm(1) - 3*xm(2) + 4*r];
    thetaHatDot = cfg.gamma*phi*sigma;
    dz = [xDot; xmDot; thetaHatDot];
end

function r = referenceSignal(t)
    r = 0.8*sin(0.7*t) + 0.35*sin(1.9*t);
end

% Simulink mapping:
% 1. Use two Integrator blocks for x and two for xm.
% 2. Form phi = [sin(x1); x1*x2] with Trigonometric Function and Product blocks.
% 3. Implement u = -4*x1 - 3*x2 + 4*r - thetaHat'*phi.
% 4. Integrate thetaHatDot = Gamma*phi*(e'*P*b) with two Integrator blocks.
% 5. Use Saturation or Projection only when the design explicitly includes it.

12.1 Simulink construction

A direct Simulink realization uses four state Integrator blocks for \(x_1,x_2,x_{m1},x_{m2}\) and two estimator Integrator blocks for \(\hat\theta_1,\hat\theta_2\). Build \(\sin x_1\) with a Trigonometric Function block and \(x_1x_2\) with a Product block. A MATLAB Function block may compute the controller and update signal. Scope or To Workspace blocks should record the tracking error, control input, and estimates.

For hardware-oriented models, place actuator saturation after the controller and add an adaptation-enabling signal or projection operator. Do not silently insert saturation while retaining the unmodified continuous-time proof.

The Wolfram Language version uses the same explicit RK4 recursion and exports a CSV file for direct comparison with the other implementations.

Chapter18_Lesson2.nb

(* Chapter18_Lesson2.nb *)
Notebook[{
  Cell["Chapter 18, Lesson 2: Linear-in-the-Parameters Nonlinear MRAC", "Title"],
  Cell[BoxData[
"ClearAll[\"Global`*\"];
dt = 0.001; tf = 20.0; gamma = DiagonalMatrix[{8.0, 8.0}];
thetaTrue = {1.2, -0.4}; b = {0.0, 1.0};
p = { {1.2083333333333333, 0.125}, {0.125, 0.2083333333333333} };
referenceSignal[t_?NumericQ] := 0.8 Sin[0.7 t] + 0.35 Sin[1.9 t];
rhs[t_?NumericQ, z_List] := Module[
  {x, xm, thetaHat, r, phi, e, sigma, u, xDot, xmDot, thetaHatDot},
  x = z[[1 ;; 2]]; xm = z[[3 ;; 4]]; thetaHat = z[[5 ;; 6]];
  r = referenceSignal[t]; phi = {Sin[x[[1]]], x[[1]] x[[2]]}; e = x - xm;
  sigma = e . p . b;
  u = -4.0 x[[1]] - 3.0 x[[2]] + 4.0 r - thetaHat . phi;
  xDot = {x[[2]], thetaTrue . phi + u};
  xmDot = {xm[[2]], -4.0 xm[[1]] - 3.0 xm[[2]] + 4.0 r};
  thetaHatDot = gamma . phi sigma;
  Join[xDot, xmDot, thetaHatDot]
];
rk4Step[t_?NumericQ, z_List, h_?NumericQ] := Module[{k1, k2, k3, k4},
  k1 = rhs[t, z]; k2 = rhs[t + h/2.0, z + h k1/2.0];
  k3 = rhs[t + h/2.0, z + h k2/2.0]; k4 = rhs[t + h, z + h k3];
  z + h (k1 + 2.0 k2 + 2.0 k3 + k4)/6.0
];
times = N@Range[0.0, tf, dt];
states = Rest@FoldList[rk4Step[#2 - dt, #1, dt] &,
  {0.8, -0.3, 0.0, 0.0, 0.0, 0.0}, Rest[times]];
states = Prepend[states, {0.8, -0.3, 0.0, 0.0, 0.0, 0.0}];
rows = MapThread[Function[{t, z}, Module[{x, xm, thetaHat, phi, u, e},
  x = z[[1 ;; 2]]; xm = z[[3 ;; 4]]; thetaHat = z[[5 ;; 6]];
  phi = {Sin[x[[1]]], x[[1]] x[[2]]};
  u = -4.0 x[[1]] - 3.0 x[[2]] + 4.0 referenceSignal[t] - thetaHat . phi;
  e = x - xm; Join[{t}, x, xm, e, thetaHat, {u}]
]], {times, states}];

Export[\"Chapter18_Lesson2_results_mathematica.csv\", Prepend[rows,
  {\"t\", \"x1\", \"x2\", \"xm1\", \"xm2\", \"e1\", \"e2\",
   \"theta_hat1\", \"theta_hat2\", \"u\"}]];

final = Last[rows];
Print[\"Final tracking-error norm: \", NumberForm[Norm[final[[6 ;; 7]]], {8, 6}]];
Print[\"Final parameter estimates: \", NumberForm[final[[8 ;; 9]], {8, 6}]];

trackingPlot = ListLinePlot[
  {Transpose[{times, states[[All, 1]]}], Transpose[{times, states[[All, 3]]}]},
  PlotLegends -> {\"x1\", \"xm1\"}, AxesLabel -> {\"t\", \"state\"},
  PlotLabel -> \"Nonlinear MRAC tracking\", ImageSize -> Large];

parameterPlot = ListLinePlot[
  {Transpose[{times, states[[All, 5]]}], Transpose[{times, states[[All, 6]]}]},
  PlotLegends -> {\"thetaHat1\", \"thetaHat2\"}, AxesLabel -> {\"t\", \"estimate\"},
  PlotLabel -> \"Adaptive parameters\", ImageSize -> Large];

Column[{trackingPlot, parameterPlot}]"], "Input"]
},
WindowTitle -> "Chapter18_Lesson2",
StyleDefinitions -> "Default.nb"
]        

Problem 1 — Identify a linear parameterization. Consider

\[ \dot x_2 = -a\sin x_1 - b x_2|x_2| + u. \]

Construct \(\boldsymbol{\theta}\) and \(\boldsymbol{\phi}(\mathbf{x})\) so that the uncertain term is \(\boldsymbol{\theta}^T\boldsymbol{\phi}\).

Solution:

\[ \boldsymbol{\theta} = \begin{bmatrix}a\\b\end{bmatrix}, \qquad \boldsymbol{\phi}(\mathbf{x}) = \begin{bmatrix}-\sin x_1\\-x_2|x_2|\end{bmatrix}. \]

Then \(\boldsymbol{\theta}^T\boldsymbol{\phi} =-a\sin x_1-bx_2|x_2|\). The plant is nonlinear in the state but affine in the unknown constants.

Problem 2 — Derive the error model. For

\[ \dot{\mathbf{x} } = \mathbf{A}\mathbf{x} + \mathbf{b} \left[u+\boldsymbol{\theta}^T\boldsymbol{\phi}(\mathbf{x})\right], \]

use \(u=\mathbf{k}_x^{\ast T}\mathbf{x}+k_r^\ast r- \hat{\boldsymbol{\theta} }^T\boldsymbol{\phi}(\mathbf{x})\) and derive the tracking-error dynamics.

Solution:

Apply the matching equations \(\mathbf{A}+\mathbf{b}\mathbf{k}_x^{\ast T} =\mathbf{A}_m\) and \(\mathbf{b}k_r^\ast=\mathbf{b}_m\). Then

\[ \dot{\mathbf{x} } = \mathbf{A}_m\mathbf{x} + \mathbf{b}_m r - \mathbf{b} (\hat{\boldsymbol{\theta} }-\boldsymbol{\theta})^T \boldsymbol{\phi}(\mathbf{x}). \]

Subtracting \(\dot{\mathbf{x} }_m=\mathbf{A}_m\mathbf{x}_m+\mathbf{b}_m r\) yields

\[ \dot{\mathbf{e} } = \mathbf{A}_m\mathbf{e} - \mathbf{b}\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\phi}(\mathbf{x}). \]

Problem 3 — Complete the Lyapunov cancellation. Starting with

\[ V = \mathbf{e}^T\mathbf{P}\mathbf{e} + \tilde{\boldsymbol{\theta} }^T \boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta} }, \]

derive an adaptive law that makes \(\dot V=-\mathbf{e}^T\mathbf{Q}\mathbf{e}\).

Solution:

\[ \dot V = -\mathbf{e}^T\mathbf{Q}\mathbf{e} - 2\mathbf{e}^T\mathbf{P}\mathbf{b} \tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} + 2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta} } }. \]

Choose

\[ \dot{\hat{\boldsymbol{\theta} } } = \boldsymbol{\Gamma}\boldsymbol{\phi} \mathbf{e}^T\mathbf{P}\mathbf{b}. \]

The last two terms are equal and opposite scalars, leaving \(\dot V=-\mathbf{e}^T\mathbf{Q}\mathbf{e}\le0\).

Problem 4 — Solve the Lyapunov equation. For \(\mathbf{A}_m=\begin{bmatrix}0&1\\-4&-3\end{bmatrix}\) and \(\mathbf{Q}=\mathbf{I}\), find the symmetric matrix \(\mathbf{P}\).

Solution:

Let

\[ \mathbf{P} = \begin{bmatrix}p_{11}&p_{12}\\p_{12}&p_{22}\end{bmatrix}. \]

Equating coefficients in \(\mathbf{A}_m^T\mathbf{P}+\mathbf{P}\mathbf{A}_m=-\mathbf{I}\) gives

\[ -8p_{12}=-1, \qquad p_{11}-3p_{12}-4p_{22}=0, \qquad 2p_{12}-6p_{22}=-1. \]

Therefore,

\[ p_{12}=\frac{1}{8}, \qquad p_{22}=\frac{5}{24}, \qquad p_{11}=\frac{29}{24}, \]

and

\[ \mathbf{P} = \begin{bmatrix} 29/24&1/8\\ 1/8&5/24 \end{bmatrix} > \mathbf{0}. \]

Problem 5 — Explain nonconvergent parameters with perfect tracking. Suppose the reference command is zero and the plant reaches the origin. Can the two parameters in the worked example necessarily be identified?

Solution:

No. At the origin, \(\boldsymbol{\phi}(\mathbf{0})=[0,0]^T\). Therefore the update law becomes zero for every parameter estimate. The tracking error can converge to zero while the estimates retain arbitrary constant errors. The regressor is not persistently exciting, so the data contain no information that distinguishes the true parameter vector from other vectors.

Problem 6 — Analyze unmatched uncertainty. Consider

\[ \dot x_1=x_2+\theta\sin x_1, \qquad \dot x_2=u. \]

Explain why the controller from this lesson cannot directly cancel the uncertainty.

Solution:

The control acts in the second state equation, while the uncertainty acts in the first. Their vector directions are different, so the uncertainty is unmatched. A term in \(u\) cannot algebraically cancel a term in \(\dot x_1\). One needs a recursive nonlinear design such as adaptive backstepping, a suitable coordinate transformation, or another structure-specific controller.

Problem 7 — Programming exercise. Modify the simulation to compare \(\boldsymbol{\Gamma}=\operatorname{diag}(2,2)\), \(\operatorname{diag}(8,8)\), and \(\operatorname{diag}(30,30)\). Record peak tracking error, final tracking error, and peak control magnitude.

Solution guidance:

Run the same initial conditions and command for each gain. A larger adaptation gain usually changes the transient more rapidly, but it does not guarantee a uniformly better response and may amplify measurement noise or numerical error. The comparison should be based on all three metrics rather than parameter speed alone. Save the modified files as Chapter18_Lesson2_Ex1.py, Chapter18_Lesson2_Ex1.cpp, Chapter18_Lesson2_Ex1.java, Chapter18_Lesson2_Ex1.m, or Chapter18_Lesson2_Ex1.nb.

A nonlinear plant can retain a tractable adaptive structure when its unknown constants enter through a known linear parameterization. Under matched uncertainty, the nonlinear regressor can be inserted into the familiar MRAC error model. A composite state-parameter Lyapunov function then produces a gradient law that cancels the mixed term and guarantees boundedness and asymptotic tracking. Exact parameter convergence is a stronger property and generally requires persistent excitation. The result is powerful but structural: unmatched uncertainty, nonlinear parameterization, approximation residuals, actuator saturation, and discrete implementation require additional analysis or modifications.