After completing this lesson, students should be able to:

• define the high-frequency gain of a SISO plant and identify its counterpart in controllable canonical form;
• distinguish an unknown gain magnitude with known sign from a truly unknown control direction;
• derive the MRAC error model and sign-corrected parameter update law;
• complete the Lyapunov proof without using the unknown gain magnitude in the implemented controller;
• explain mathematically why using the wrong gain sign invalidates the standard proof; and
• reproduce the same adaptive closed loop in several programming environments.

The wording “unknown high-frequency gain sign” in the lesson title must be interpreted together with the parenthetical qualification sign-definite case. In this lesson, the numerical gain is unknown, but all admissible gains have the same known sign. The case in which even that sign is unavailable requires a different mechanism and is reserved for Chapter 16.

flowchart TD
  A["Unknown plant input gain"] --> B["Is the admissible gain set sign-definite?"]
  B -->|"yes"| C["Store known sign s_b = +1 or -1"]
  C --> D["Use sign-corrected Lyapunov update law"]
  D --> E["Prove boundedness and tracking"]
  B -->|"no"| F["Standard canonical-form MRAC law is not justified"]
  F --> G["Defer unknown-control-direction methods \nto Chapter 16"]
        

Consider a proper SISO transfer function \( G_p(s)=N_p(s)/D_p(s) \) with relative degree \( \rho \). Its high-frequency gain is the leading input-output coefficient

\[ k_p = \lim_{s\to\infty}s^\rho G_p(s). \]

The sign of \( k_p \) is the control direction: at sufficiently high frequency, it determines whether a positive input initially drives the measured output in the positive or negative direction. In the full-state controllable-canonical representation used in this chapter,

\[ \dot{\mathbf{x}} = A_p\mathbf{x} + b_p\mathbf{b}_0 u, \qquad \mathbf{b}_0 = \begin{bmatrix} 0 & \cdots & 0 & 1 \end{bmatrix}^{T}, \]

the scalar \( b_p \) plays the same structural role. Its magnitude is unknown. The sign-definite assumption states that the uncertainty set is entirely positive or entirely negative:

\[ b_p\in\mathcal{B}_{+} = [\underline b,\overline b] \quad\text{or}\quad b_p\in\mathcal{B}_{-} = [-\overline b,-\underline b], \qquad 0<\underline b\le\overline b. \]

Therefore the controller knows the constant metadata \( s_b=\operatorname{sgn}(b_p)\in\{-1,+1\} \), but it does not know \( |b_p| \). Equivalently,

\[ b_p=s_b|b_p|, \qquad s_b\ \text{known}, \qquad |b_p|\ \text{unknown}. \]

This is a structural assumption, not an online estimate. It normally comes from actuator wiring, force/torque orientation, valve convention, motor polarity, or a validated local model. A controller should not infer it casually from noisy transient data.

Let the uncertain plant and desired reference model be

\[ \begin{aligned} \dot{\mathbf{x}} &= A_p\mathbf{x}+b_p\mathbf{b}_0u,\\ \dot{\mathbf{x}}_m &= A_m\mathbf{x}_m+b_m\mathbf{b}_0r, \end{aligned} \]

where \( A_m \) is Hurwitz, \( r(t) \) is bounded and piecewise continuous, and the states are available. The canonical-form matching equations require constant ideal controller parameters \( \boldsymbol{\theta}_x^\ast \) and \( \theta_r^\ast \) satisfying

\[ A_p+b_p\mathbf{b}_0{\boldsymbol{\theta}_x^\ast}^{T}=A_m, \qquad b_p\theta_r^\ast=b_m. \]

Define the combined ideal parameter, estimate, and regressor:

\[ \boldsymbol{\theta}^{\ast} = \begin{bmatrix} \boldsymbol{\theta}_x^\ast\\ \theta_r^\ast \end{bmatrix}, \qquad \hat{\boldsymbol{\theta}} = \begin{bmatrix} \hat{\boldsymbol{\theta}}_x\\ \hat{\theta}_r \end{bmatrix}, \qquad \boldsymbol{\omega} = \begin{bmatrix} \mathbf{x}\\ r \end{bmatrix}. \]

The direct adaptive control law is linear in its estimates:

\[ u = \hat{\boldsymbol{\theta}}^{T}\boldsymbol{\omega} = \hat{\boldsymbol{\theta}}_x^{T}\mathbf{x} + \hat{\theta}_r r. \]

With tracking error \( \mathbf{e}=\mathbf{x}-\mathbf{x}_m \) and parameter error \( \tilde{\boldsymbol{\theta}} =\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{\ast} \) , substitution of the matching equations gives

\[ \dot{\mathbf{e}} = A_m\mathbf{e} + b_p\mathbf{b}_0 \tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega}. \]

The unknown gain multiplies the entire parameter-error channel. This is exactly why its sign enters the adaptive law.

Choose any symmetric matrix \( Q=Q^{T}>0 \). Since \( A_m \) is Hurwitz, the Lyapunov equation

\[ A_m^{T}P+PA_m=-Q \]

has a unique symmetric solution \( P=P^{T}>0 \). Let the adaptation-gain matrix satisfy \( \Gamma=\Gamma^{T}>0 \). Define the scalar filtered error

\[ \sigma_e = \mathbf{e}^{T}P\mathbf{b}_0. \]

The sign-corrected update law is

\[ \boxed{ \dot{\hat{\boldsymbol{\theta}}} = -s_b\Gamma\boldsymbol{\omega}\sigma_e }. \]

The unknown magnitude \( |b_p| \) does not appear in this differential equation. Only the known sign \( s_b \) is implemented.

flowchart TD
  R["Command r"] --> RM["Reference model"]
  R --> REG["Regressor omega = [x; r]"]
  P["Plant state x"] --> REG
  P --> ERR["Tracking error e = x - xm"]
  RM --> ERR
  ERR --> SIG["sigma_e = e^T P b0"]
  REG --> LAW["theta_hat_dot = -s_b Gamma omega sigma_e"]
  SIG --> LAW
  LAW --> INT["Parameter integrator"]
  INT --> CTRL["u = theta_hat^T omega"]
  REG --> CTRL
  CTRL --> P
        

Theorem. Suppose that:

1. \( A_m \) is Hurwitz;
2. the constant matching parameters exist;
3. \( b_p\neq0 \) is constant and \( s_b=\operatorname{sgn}(b_p) \) is known;
4. \( r(t) \) is bounded and piecewise continuous; and
5. the plant and model states used by the controller are available.

Then the adaptive law in Section 4 makes all closed-loop signals bounded and gives \( \lim_{t\to\infty}\mathbf{e}(t)=\mathbf{0} \).

Proof. Use the Lyapunov candidate

\[ V = \frac{1}{2}\mathbf{e}^{T}P\mathbf{e} + \frac{|b_p|}{2} \tilde{\boldsymbol{\theta}}^{T} \Gamma^{-1} \tilde{\boldsymbol{\theta}}. \]

Although \( |b_p| \) is unknown, it is a fixed positive number and may be used in a proof. It is not needed by the controller. Since the ideal parameters are constant, \( \dot{\tilde{\boldsymbol{\theta}}} =\dot{\hat{\boldsymbol{\theta}}} \) . Differentiation gives

\[ \begin{aligned} \dot V &= \frac{1}{2} \mathbf{e}^{T} (A_m^{T}P+PA_m) \mathbf{e}\\ &\quad+ b_p \mathbf{e}^{T}P\mathbf{b}_0 \tilde{\boldsymbol{\theta}}^{T} \boldsymbol{\omega}\\ &\quad+ |b_p| \tilde{\boldsymbol{\theta}}^{T} \Gamma^{-1} \dot{\hat{\boldsymbol{\theta}}}. \end{aligned} \]

Substitute the Lyapunov equation and adaptive law:

\[ \begin{aligned} \dot V &= -\frac{1}{2}\mathbf{e}^{T}Q\mathbf{e} + b_p\sigma_e \tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega}\\ &\quad- |b_p|s_b \tilde{\boldsymbol{\theta}}^{T} \boldsymbol{\omega}\sigma_e. \end{aligned} \]

Since \( |b_p|s_b=b_p \), the two mixed terms cancel:

\[ \boxed{ \dot V = -\frac{1}{2}\mathbf{e}^{T}Q\mathbf{e} \le0 }. \]

Consequently, \( V(t)\le V(0) \), so \( \mathbf{e} \) and \( \tilde{\boldsymbol{\theta}} \) are bounded. A bounded command passed through the stable reference model yields bounded \( \mathbf{x}_m \); therefore \( \mathbf{x}=\mathbf{e}+\mathbf{x}_m \), \( \boldsymbol{\omega} \), and \( u=\hat{\boldsymbol{\theta}}^{T}\boldsymbol{\omega} \) are bounded. Integration of \( \dot V \) gives

\[ \int_{0}^{\infty} \mathbf{e}^{T}Q\mathbf{e}\,dt \le 2V(0), \]

hence \( \mathbf{e}\in L_2 \). The error dynamics then imply bounded \( \dot{\mathbf{e}} \). By Barbalat’s lemma, which was introduced in Chapter 3,

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

This theorem guarantees tracking, not exact recovery of every ideal parameter. The additional signal condition needed for parameter convergence is developed later in Chapter 10.

Suppose an implemented sign \( \hat s_b\in\{-1,+1\} \) is used:

\[ \dot{\hat{\boldsymbol{\theta}}} = -\hat s_b\Gamma\boldsymbol{\omega}\sigma_e. \]

Repeating the derivative calculation gives

\[ \dot V = -\frac{1}{2}\mathbf{e}^{T}Q\mathbf{e} + \left( b_p-|b_p|\hat s_b \right) \sigma_e \tilde{\boldsymbol{\theta}}^{T} \boldsymbol{\omega}. \]

If \( \hat s_b=s_b \), the coefficient of the mixed term is zero. If the sign is reversed, \( \hat s_b=-s_b \), then

\[ b_p-|b_p|\hat s_b = b_p+|b_p|s_b = 2b_p, \]

and therefore

\[ \dot V = -\frac{1}{2}\mathbf{e}^{T}Q\mathbf{e} + 2b_p\sigma_e \tilde{\boldsymbol{\theta}}^{T} \boldsymbol{\omega}. \]

The remaining term is sign-indefinite and can inject energy into the combined tracking-parameter dynamics. Thus the classical proof does not merely become conservative; its central cancellation fails. A sign that may change with operating condition also violates the constant sign-definite assumption.

The implementations use the following plant and reference model:

\[ A_p = \begin{bmatrix} 0 & 1\\ -1 & -1.2 \end{bmatrix}, \qquad b_p=2, \qquad \mathbf{b}_0 = \begin{bmatrix} 0\\1 \end{bmatrix}, \]

\[ A_m = \begin{bmatrix} 0 & 1\\ -4 & -4 \end{bmatrix}, \qquad b_m=4. \]

The controller does not use the magnitude \( b_p=2 \) in its adaptive update. The simulation plant necessarily uses it to generate the true state trajectory. Solving the matching equations gives

\[ {\boldsymbol{\theta}_x^\ast}^{T} = \begin{bmatrix} -1.5 & -1.4 \end{bmatrix}, \qquad \theta_r^\ast=2. \]

For \( Q=I_2 \), the Lyapunov solution is

\[ P = \begin{bmatrix} 9/8 & 1/8\\ 1/8 & 5/32 \end{bmatrix}, \qquad P\mathbf{b}_0 = \begin{bmatrix} 1/8\\ 5/32 \end{bmatrix}. \]

Hence

\[ \sigma_e = \frac{1}{8}e_1+\frac{5}{32}e_2. \]

With \( \Gamma=\operatorname{diag}(8,8,5) \) and \( s_b=+1 \), the update equations are

\[ \begin{aligned} \dot{\hat\theta}_{x1} &=-8x_1\sigma_e,\\ \dot{\hat\theta}_{x2} &=-8x_2\sigma_e,\\ \dot{\hat\theta}_{r} &=-5r\sigma_e. \end{aligned} \]

The common initial condition is \( \mathbf{x}(0)=[0.7,-0.4]^{T} \), \( \mathbf{x}_m(0)=\mathbf{0} \), \( \hat{\boldsymbol{\theta}}(0)=\mathbf{0} \) , and the command is a unit step.

All implementations integrate one augmented state containing the plant, reference model, and parameter estimates:

\[ \mathbf{z} = \begin{bmatrix} \mathbf{x}^{T} & \mathbf{x}_m^{T} & \hat{\boldsymbol{\theta}}^{T} \end{bmatrix}^{T}. \]

A single right-hand-side function is preferable because the controller and plant are dynamically coupled. Python, C++, and Java use a fourth-order Runge–Kutta method with \( \Delta t=0.002\ \mathrm{s} \). MATLAB uses ode45. The Simulink script builds the same equations from continuous State-Space, Integrator, Gain, Product, and Dot Product blocks.

The programs export a common CSV schema: t, x1, x2, xm1, xm2, e1, e2, theta_x1, theta_x2, theta_r, u. This enables cross-language regression testing. A suitable acceptance test is a small final tracking-error norm together with finite states, control, and parameter estimates.

The estimates are not expected to equal the matching parameters under every bounded command. The theorem requires tracking convergence only.

Required libraries: numpy for vector algebra and matplotlib for plots. The CSV module is part of the Python standard library.

Chapter7_Lesson4.py

"""Chapter 7, Lesson 4: MRAC with a known sign-definite input gain.

The plant input-gain magnitude is unknown to the controller. Its constant
sign is known and is the only gain information used by the adaptive law.
"""

from __future__ import annotations

import csv
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


A_P = np.array([[0.0, 1.0], [-1.0, -1.2]])
A_M = np.array([[0.0, 1.0], [-4.0, -4.0]])
B_0 = np.array([0.0, 1.0])
B_P = 2.0                 # Used only by the simulated plant.
B_M = 4.0
SIGN_B = 1.0              # Known sign of the admissible input-gain set.
GAMMA = np.array([8.0, 8.0, 5.0])
P = np.array([[9.0 / 8.0, 1.0 / 8.0],
              [1.0 / 8.0, 5.0 / 32.0]])
PB0 = P @ B_0


def reference(t: float) -> float:
    """Unit-step command."""
    return 1.0 if t >= 0.0 else 0.0


def rhs(t: float, z: np.ndarray) -> np.ndarray:
    """Augmented plant, model, and adaptive-parameter dynamics."""
    x = z[0:2]
    xm = z[2:4]
    theta_hat = z[4:7]

    r = reference(t)
    omega = np.array([x[0], x[1], r])
    u = float(theta_hat @ omega)

    x_dot = A_P @ x + B_P * B_0 * u
    xm_dot = A_M @ xm + B_M * B_0 * r

    e = x - xm
    sigma = float(e @ PB0)
    theta_dot = -SIGN_B * GAMMA * omega * sigma

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


def rk4_step(t: float, z: np.ndarray, dt: float) -> np.ndarray:
    """One classical fourth-order Runge-Kutta step."""
    k1 = rhs(t, z)
    k2 = rhs(t + 0.5 * dt, z + 0.5 * dt * k1)
    k3 = rhs(t + 0.5 * dt, z + 0.5 * dt * k2)
    k4 = rhs(t + dt, z + dt * k3)
    return z + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(t_final: float = 20.0, dt: float = 0.002) -> tuple[np.ndarray, np.ndarray]:
    """Run the MRAC simulation and return time and state histories."""
    if SIGN_B not in (-1.0, 1.0):
        raise ValueError("SIGN_B must be +1 or -1.")
    if dt <= 0.0 or t_final <= 0.0:
        raise ValueError("dt and t_final must be positive.")

    steps = int(round(t_final / dt))
    time = np.linspace(0.0, steps * dt, steps + 1)
    history = np.zeros((steps + 1, 7))
    history[0] = np.array([0.7, -0.4, 0.0, 0.0, 0.0, 0.0, 0.0])

    for k in range(steps):
        history[k + 1] = rk4_step(time[k], history[k], dt)
        if not np.all(np.isfinite(history[k + 1])):
            raise FloatingPointError("Simulation became non-finite.")

    return time, history


def write_csv(path: Path, time: np.ndarray, history: np.ndarray) -> None:
    """Write the common cross-language output format."""
    with path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(
            ["t", "x1", "x2", "xm1", "xm2", "e1", "e2",
             "theta_x1", "theta_x2", "theta_r", "u"]
        )
        for t, z in zip(time, history):
            x = z[0:2]
            xm = z[2:4]
            theta_hat = z[4:7]
            r = reference(float(t))
            omega = np.array([x[0], x[1], r])
            u = float(theta_hat @ omega)
            e = x - xm
            writer.writerow(
                [t, x[0], x[1], xm[0], xm[1], e[0], e[1],
                 theta_hat[0], theta_hat[1], theta_hat[2], u]
            )


def main() -> None:
    time, history = simulate()
    output = Path("Chapter7_Lesson4_python.csv")
    write_csv(output, time, history)

    e_final = history[-1, 0:2] - history[-1, 2:4]
    print(f"Final tracking-error norm: {np.linalg.norm(e_final):.6e}")
    print(f"Final parameter estimate: {history[-1, 4:7]}")
    print("CSV:", output.resolve())

    plt.figure()
    plt.plot(time, history[:, 0], label="x1")
    plt.plot(time, history[:, 2], "--", label="xm1")
    plt.xlabel("Time (s)")
    plt.ylabel("Position-like state")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter7_Lesson4_tracking.png", dpi=180)

    plt.figure()
    plt.plot(time, history[:, 4], label="theta_x1")
    plt.plot(time, history[:, 5], label="theta_x2")
    plt.plot(time, history[:, 6], label="theta_r")
    plt.xlabel("Time (s)")
    plt.ylabel("Adaptive parameters")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter7_Lesson4_parameters.png", dpi=180)
    plt.show()


if __name__ == "__main__":
    main()

This version uses only the C++17 standard library. Compile with g++ -std=c++17 -O2 Chapter7_Lesson4.cpp -o Chapter7_Lesson4.

Chapter7_Lesson4.cpp

// Chapter 7, Lesson 4: MRAC with a known sign-definite input gain.
// Build: g++ -std=c++17 -O2 Chapter7_Lesson4.cpp -o Chapter7_Lesson4

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

namespace {
using State = std::array<double, 7>;

constexpr double BP = 2.0;      // Used only by the simulated plant.
constexpr double BM = 4.0;
constexpr double SIGN_B = 1.0;  // Known sign of the admissible gain set.
constexpr std::array<double, 3> GAMMA{8.0, 8.0, 5.0};
constexpr std::array<double, 2> PB0{1.0 / 8.0, 5.0 / 32.0};

double reference(double t) {
    return t >= 0.0 ? 1.0 : 0.0;
}

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

State rhs(double t, const State& z) {
    const double x1 = z[0];
    const double x2 = z[1];
    const double xm1 = z[2];
    const double xm2 = z[3];
    const double th1 = z[4];
    const double th2 = z[5];
    const double thr = z[6];

    const double r = reference(t);
    const std::array<double, 3> omega{x1, x2, r};
    const double u = th1 * x1 + th2 * x2 + thr * r;

    const double x1_dot = x2;
    const double x2_dot = -x1 - 1.2 * x2 + BP * u;
    const double xm1_dot = xm2;
    const double xm2_dot = -4.0 * xm1 - 4.0 * xm2 + BM * r;

    const double e1 = x1 - xm1;
    const double e2 = x2 - xm2;
    const double sigma = e1 * PB0[0] + e2 * PB0[1];

    State dz{};
    dz[0] = x1_dot;
    dz[1] = x2_dot;
    dz[2] = xm1_dot;
    dz[3] = xm2_dot;
    for (std::size_t i = 0; i < 3; ++i) {
        dz[4 + i] = -SIGN_B * GAMMA[i] * omega[i] * sigma;
    }
    return dz;
}

State rk4_step(double t, const State& z, double dt) {
    const State k1 = rhs(t, z);
    const State k2 = rhs(t + 0.5 * dt, add_scaled(z, k1, 0.5 * dt));
    const State k3 = rhs(t + 0.5 * dt, add_scaled(z, k2, 0.5 * dt));
    const State k4 = rhs(t + dt, add_scaled(z, k3, dt));

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

void write_row(std::ofstream& out, double t, const State& z) {
    const double r = reference(t);
    const double e1 = z[0] - z[2];
    const double e2 = z[1] - z[3];
    const double u = z[4] * z[0] + z[5] * z[1] + z[6] * r;
    out << t << ',' << z[0] << ',' << z[1] << ','
        << z[2] << ',' << z[3] << ',' << e1 << ',' << e2 << ','
        << z[4] << ',' << z[5] << ',' << z[6] << ',' << u << '\n';
}
}  // namespace

int main() {
    try {
        if (SIGN_B != 1.0 && SIGN_B != -1.0) {
            throw std::invalid_argument("SIGN_B must be +1 or -1.");
        }

        constexpr double dt = 0.002;
        constexpr double t_final = 20.0;
        const int steps = static_cast<int>(std::lround(t_final / dt));

        State z{0.7, -0.4, 0.0, 0.0, 0.0, 0.0, 0.0};
        std::ofstream out("Chapter7_Lesson4_cpp.csv");
        if (!out) {
            throw std::runtime_error("Could not open CSV output.");
        }
        out << std::setprecision(15);
        out << "t,x1,x2,xm1,xm2,e1,e2,theta_x1,theta_x2,theta_r,u\n";
        write_row(out, 0.0, z);

        for (int k = 0; k < steps; ++k) {
            const double t = k * dt;
            z = rk4_step(t, z, dt);
            for (double value : z) {
                if (!std::isfinite(value)) {
                    throw std::runtime_error("Simulation became non-finite.");
                }
            }
            write_row(out, (k + 1) * dt, z);
        }

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

This implementation uses the Java standard library and writes the same CSV fields as the other programs.

Chapter7_Lesson4.java

// Chapter 7, Lesson 4: MRAC with a known sign-definite input gain.
// Build: javac Chapter7_Lesson4.java
// Run:   java Chapter7_Lesson4

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 Chapter7_Lesson4 {
    private static final double BP = 2.0;      // Simulated plant only.
    private static final double BM = 4.0;
    private static final double SIGN_B = 1.0;  // Known admissible-set sign.
    private static final double[] GAMMA = {8.0, 8.0, 5.0};
    private static final double[] PB0 = {1.0 / 8.0, 5.0 / 32.0};

    private Chapter7_Lesson4() {
    }

    private static double reference(double t) {
        return t >= 0.0 ? 1.0 : 0.0;
    }

    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[] rhs(double t, double[] z) {
        double x1 = z[0];
        double x2 = z[1];
        double xm1 = z[2];
        double xm2 = z[3];
        double th1 = z[4];
        double th2 = z[5];
        double thr = z[6];

        double r = reference(t);
        double[] omega = {x1, x2, r};
        double u = th1 * x1 + th2 * x2 + thr * r;

        double e1 = x1 - xm1;
        double e2 = x2 - xm2;
        double sigma = e1 * PB0[0] + e2 * PB0[1];

        double[] dz = new double[7];
        dz[0] = x2;
        dz[1] = -x1 - 1.2 * x2 + BP * u;
        dz[2] = xm2;
        dz[3] = -4.0 * xm1 - 4.0 * xm2 + BM * r;
        for (int i = 0; i < 3; i++) {
            dz[4 + i] = -SIGN_B * GAMMA[i] * omega[i] * sigma;
        }
        return dz;
    }

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

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

    private static void writeRow(BufferedWriter writer, double t, double[] z)
            throws IOException {
        double r = reference(t);
        double e1 = z[0] - z[2];
        double e2 = z[1] - z[3];
        double u = z[4] * z[0] + z[5] * z[1] + z[6] * r;
        writer.write(String.format(
            Locale.US,
            "%.15g,%.15g,%.15g,%.15g,%.15g,%.15g,%.15g,"
                + "%.15g,%.15g,%.15g,%.15g%n",
            t, z[0], z[1], z[2], z[3], e1, e2,
            z[4], z[5], z[6], u
        ));
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        if (SIGN_B != 1.0 && SIGN_B != -1.0) {
            throw new IllegalArgumentException("SIGN_B must be +1 or -1.");
        }

        final double dt = 0.002;
        final double tFinal = 20.0;
        final int steps = (int) Math.round(tFinal / dt);
        double[] z = {0.7, -0.4, 0.0, 0.0, 0.0, 0.0, 0.0};

        Path output = Path.of("Chapter7_Lesson4_java.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(
                output, StandardCharsets.UTF_8)) {
            writer.write(
                "t,x1,x2,xm1,xm2,e1,e2,theta_x1,theta_x2,theta_r,u\n"
            );
            writeRow(writer, 0.0, z);

            for (int k = 0; k < steps; k++) {
                double t = k * dt;
                z = rk4Step(t, z, dt);
                for (double value : z) {
                    if (!Double.isFinite(value)) {
                        throw new ArithmeticException(
                            "Simulation became non-finite."
                        );
                    }
                }
                writeRow(writer, (k + 1) * dt, z);
            }
        } catch (IOException error) {
            System.err.println("I/O error: " + error.getMessage());
            System.exit(1);
        }

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

The script uses ode45, creates tracking and parameter plots, and exports a table. No adaptive-control toolbox block is required because the update law is implemented explicitly.

Chapter7_Lesson4.m

%% Chapter 7, Lesson 4: MRAC with a known sign-definite input gain
clear; clc; close all;

Ap = [0 1; -1 -1.2];
Am = [0 1; -4 -4];
b0 = [0; 1];
bp = 2.0;                 % Used only by the simulated plant
bm = 4.0;
sign_b = 1.0;             % Known sign of the admissible input-gain set
Gamma = diag([8 8 5]);
P = [9/8 1/8; 1/8 5/32]; % Solves Am'*P + P*Am = -I
Pb0 = P*b0;

if ~ismember(sign_b, [-1 1])
    error('sign_b must be +1 or -1.');
end

z0 = [0.7; -0.4; 0; 0; 0; 0; 0];
tspan = [0 20];
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[t, z] = ode45(@closedLoop, tspan, z0, options, ...
    Ap, Am, b0, bp, bm, sign_b, Gamma, Pb0);

x = z(:, 1:2);
xm = z(:, 3:4);
thetaHat = z(:, 5:7);
e = x - xm;
r = ones(size(t));
u = sum(thetaHat .* [x r], 2);

result = table(t, x(:,1), x(:,2), xm(:,1), xm(:,2), ...
    e(:,1), e(:,2), thetaHat(:,1), thetaHat(:,2), ...
    thetaHat(:,3), u, ...
    'VariableNames', {'t','x1','x2','xm1','xm2','e1','e2', ...
    'theta_x1','theta_x2','theta_r','u'});
writetable(result, 'Chapter7_Lesson4_matlab.csv');

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

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

figure;
plot(t, thetaHat, 'LineWidth', 1.2);
grid on; xlabel('Time (s)'); ylabel('Adaptive parameters');
legend('\theta_{x1}', '\theta_{x2}', '\theta_r', 'Location', 'best');
title('Parameter-estimate histories');

function dz = closedLoop(t, z, Ap, Am, b0, bp, bm, sign_b, Gamma, Pb0)
    x = z(1:2);
    xm = z(3:4);
    thetaHat = z(5:7);

    r = double(t >= 0);
    omega = [x; r];
    u = thetaHat.' * omega;

    xDot = Ap*x + bp*b0*u;
    xmDot = Am*xm + bm*b0*r;

    e = x - xm;
    sigma = e.' * Pb0;
    thetaDot = -sign_b * Gamma * omega * sigma;

    dz = [xDot; xmDot; thetaDot];
end

The following MATLAB script creates and simulates Chapter7_Lesson4_Simulink_Model.slx. It requires Simulink. The gain magnitude appears in the simulated Plant block, while only sign_b enters the adaptive law.

Chapter7_Lesson4_Simulink.m

%% Chapter 7, Lesson 4: programmatic Simulink model
% Requires Simulink. The script creates a continuous-time block diagram
% whose adaptive law uses only the known sign of the plant input gain.
clear; clc;

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

Ap = [0 1; -1 -1.2];
Am = [0 1; -4 -4];
b0 = [0; 1];
bp = 2.0;                 % Simulated plant only
bm = 4.0;
sign_b = 1.0;             % Known sign-definite assumption
Gamma = diag([8 8 5]);
P = [9/8 1/8; 1/8 5/32];
Pb0 = P*b0;

new_system(model);
open_system(model);

add_block('simulink/Sources/Step', [model '/Reference r'], ...
    'Time', '0', 'Before', '0', 'After', '1', ...
    'Position', [30 95 60 125]);

add_block('simulink/Continuous/State-Space', ...
    [model '/Reference Model'], ...
    'A', 'Am', 'B', 'bm*b0', 'C', 'eye(2)', ...
    'D', 'zeros(2,1)', 'X0', '[0;0]', ...
    'Position', [120 45 270 105]);

add_block('simulink/Continuous/State-Space', [model '/Plant'], ...
    'A', 'Ap', 'B', 'bp*b0', 'C', 'eye(2)', ...
    'D', 'zeros(2,1)', 'X0', '[0.7;-0.4]', ...
    'Position', [600 155 750 215]);

add_block('simulink/Math Operations/Sum', [model '/Tracking Error'], ...
    'Inputs', '+-', 'IconShape', 'rectangular', ...
    'Position', [805 80 850 135]);

add_block('simulink/Signal Routing/Mux', [model '/Regressor omega'], ...
    'Inputs', '2', 'Position', [285 150 290 215]);

add_block('simulink/Math Operations/Dot Product', ...
    [model '/Control u = theta dot omega'], ...
    'Position', [465 155 545 205]);

add_block('simulink/Continuous/Integrator', ...
    [model '/Adaptive Parameters'], ...
    'InitialCondition', '[0;0;0]', ...
    'Position', [330 280 380 330]);

add_block('simulink/Math Operations/Gain', [model '/Sigma'], ...
    'Gain', 'transpose(Pb0)', 'Multiplication', 'Matrix(K*u)', ...
    'Position', [905 85 980 130]);

add_block('simulink/Math Operations/Product', ...
    [model '/omega times sigma'], ...
    'Inputs', '**', 'Multiplication', 'Element-wise(.*)', ...
    'Position', [300 385 355 435]);

add_block('simulink/Math Operations/Gain', ...
    [model '/Minus sign Gamma'], ...
    'Gain', '-sign_b*Gamma', 'Multiplication', 'Matrix(K*u)', ...
    'Position', [420 380 540 440]);

add_block('simulink/Sinks/To Workspace', [model '/x workspace'], ...
    'VariableName', 'x_sim', 'SaveFormat', 'Structure With Time', ...
    'Position', [805 180 895 210]);

add_block('simulink/Sinks/To Workspace', [model '/xm workspace'], ...
    'VariableName', 'xm_sim', 'SaveFormat', 'Structure With Time', ...
    'Position', [310 40 410 70]);

add_block('simulink/Sinks/To Workspace', [model '/theta workspace'], ...
    'VariableName', 'theta_sim', 'SaveFormat', 'Structure With Time', ...
    'Position', [420 285 520 315]);

add_block('simulink/Sinks/To Workspace', [model '/error workspace'], ...
    'VariableName', 'error_sim', 'SaveFormat', 'Structure With Time', ...
    'Position', [905 25 1005 55]);

add_line(model, 'Reference r/1', 'Reference Model/1');
add_line(model, 'Reference r/1', 'Regressor omega/2');
add_line(model, 'Reference Model/1', 'Tracking Error/2');
add_line(model, 'Reference Model/1', 'xm workspace/1');
add_line(model, 'Plant/1', 'Tracking Error/1');
add_line(model, 'Plant/1', 'Regressor omega/1');
add_line(model, 'Plant/1', 'x workspace/1');
add_line(model, 'Regressor omega/1', 'Control u = theta dot omega/1');
add_line(model, 'Adaptive Parameters/1', ...
    'Control u = theta dot omega/2');
add_line(model, 'Control u = theta dot omega/1', 'Plant/1');
add_line(model, 'Tracking Error/1', 'Sigma/1');
add_line(model, 'Tracking Error/1', 'error workspace/1');
add_line(model, 'Regressor omega/1', 'omega times sigma/1');
add_line(model, 'Sigma/1', 'omega times sigma/2');
add_line(model, 'omega times sigma/1', 'Minus sign Gamma/1');
add_line(model, 'Minus sign Gamma/1', 'Adaptive Parameters/1');
add_line(model, 'Adaptive Parameters/1', 'theta workspace/1');

set_param(model, 'Solver', 'ode45', 'StopTime', '20');
save_system(model);
simOut = sim(model);

eFinal = simOut.error_sim.signals.values(end, :);
fprintf('Final tracking-error norm: %.6e\n', norm(eFinal));
fprintf('Created model: %s.slx\n', model);

This notebook uses NDSolveValue to solve the coupled nonlinear closed-loop equations and produces tracking and parameter plots.

Chapter7_Lesson4.nb

Notebook[{
  Cell["Chapter 7, Lesson 4: Sign-Definite High-Frequency Gain", "Title"],
  Cell["MRAC simulation with unknown gain magnitude and known constant sign.", "Text"],
  Cell[BoxData["ClearAll[\"Global`*\"];\n\
Ap = {{0., 1.}, {-1., -1.2}};\n\
Am = {{0., 1.}, {-4., -4.}};\n\
b0 = {0., 1.}; bp = 2.; bm = 4.; signB = 1.;\n\
gamma = {8., 8., 5.};\n\
p = {{9./8., 1./8.}, {1./8., 5./32.}};\n\
pb0 = p.b0;\n\
tFinal = 20.;\n\
sol = NDSolveValue[{\n\
   x1'[t] == x2[t],\n\
   x2'[t] == -x1[t] - 1.2 x2[t] + bp u[t],\n\
   xm1'[t] == xm2[t],\n\
   xm2'[t] == -4 xm1[t] - 4 xm2[t] + bm,\n\
   th1'[t] == -signB gamma[[1]] x1[t] sigma[t],\n\
   th2'[t] == -signB gamma[[2]] x2[t] sigma[t],\n\
   thr'[t] == -signB gamma[[3]] sigma[t],\n\
   u[t] == th1[t] x1[t] + th2[t] x2[t] + thr[t],\n\
   sigma[t] == pb0.{x1[t] - xm1[t], x2[t] - xm2[t]},\n\
   x1[0] == 0.7, x2[0] == -0.4,\n\
   xm1[0] == 0., xm2[0] == 0.,\n\
   th1[0] == 0., th2[0] == 0., thr[0] == 0.\n\
   }, {x1, x2, xm1, xm2, th1, th2, thr}, {t, 0, tFinal},\n\
   Method -> {\"EquationSimplification\" -> \"Residual\"}];\n\
{x1f, x2f, xm1f, xm2f, th1f, th2f, thrf} = sol;\n\
eFinal = {x1f[tFinal] - xm1f[tFinal], x2f[tFinal] - xm2f[tFinal]};\n\
Print[\"Final tracking-error norm: \", ScientificForm[Norm[eFinal], 6]];\n\
Print[\"Final parameter estimate: \", {th1f[tFinal], th2f[tFinal], thrf[tFinal]}];\n\
trackingPlot = Plot[{x1f[t], xm1f[t]}, {t, 0, tFinal},\n\
   PlotLegends -> {\"x1\", \"xm1\"}, AxesLabel -> {\"t\", \"state\"},\n\
   PlotRange -> All, GridLines -> Automatic];\n\
parameterPlot = Plot[{th1f[t], th2f[t], thrf[t]}, {t, 0, tFinal},\n\
   PlotLegends -> {\"theta_x1\", \"theta_x2\", \"theta_r\"},\n\
   AxesLabel -> {\"t\", \"estimate\"}, PlotRange -> All,\n\
   GridLines -> Automatic];\n\
GraphicsRow[{trackingPlot, parameterPlot}, ImageSize -> Large]"], "Input"]
},
WindowTitle -> "Chapter7_Lesson4",
StyleDefinitions -> "Default.nb"
]

The sign-definite assumption should be treated as a verified interface contract. Before deployment:

1. document the physical direction associated with positive controller output;
2. verify sensor and actuator coordinate conventions independently;
3. test both the model equation and the software constant SIGN_B or sign_b;
4. enforce numerical monitoring for non-finite states and excessive control; and
5. reject operation if maintenance or rewiring can reverse actuator polarity without updating controller configuration.

The pure MRAC law developed here intentionally excludes projection, normalization, leakage, dead zones, and saturation compensation. Those additions are introduced in subsequent chapters after the nominal cancellation mechanism is understood.

Problem 1 — Matching Parameters. For

\[ A_p= \begin{bmatrix} 0 & 1\\ -a_1 & -a_2 \end{bmatrix}, \quad A_m= \begin{bmatrix} 0 & 1\\ -a_{m1} & -a_{m2} \end{bmatrix}, \quad b_p\neq0, \]

derive \( \boldsymbol{\theta}_x^\ast \) and \( \theta_r^\ast \).

Solution. The matching equation is

\[ A_p+b_p \begin{bmatrix} 0\\1 \end{bmatrix} \begin{bmatrix} \theta_{x1}^{\ast} & \theta_{x2}^{\ast} \end{bmatrix} = A_m. \]

Only the last row changes. Equating its two entries gives

\[ -a_1+b_p\theta_{x1}^{\ast}=-a_{m1}, \qquad -a_2+b_p\theta_{x2}^{\ast}=-a_{m2}. \]

Therefore

\[ \theta_{x1}^{\ast} = \frac{a_1-a_{m1}}{b_p}, \qquad \theta_{x2}^{\ast} = \frac{a_2-a_{m2}}{b_p}, \qquad \theta_r^\ast=\frac{b_m}{b_p}. \]

Problem 2 — Cancellation for a Negative Gain. Assume \( b_p<0 \). Show that the same Lyapunov derivative is obtained when \( s_b=-1 \).

Solution. Since

\[ b_p=-|b_p|, \qquad s_b=-1, \]

the adaptive law becomes

\[ \dot{\hat{\boldsymbol{\theta}}} = +\Gamma\boldsymbol{\omega}\sigma_e. \]

The parameter part of \( \dot V \) is then

\[ |b_p| \tilde{\boldsymbol{\theta}}^{T} \Gamma^{-1} \dot{\hat{\boldsymbol{\theta}}} = |b_p| \tilde{\boldsymbol{\theta}}^{T} \boldsymbol{\omega}\sigma_e. \]

The error-dynamics mixed term is \( b_p\sigma_e\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega} = -|b_p|\sigma_e \tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega} \) . The terms cancel, leaving \( \dot V=-\mathbf{e}^{T}Q\mathbf{e}/2 \).

Problem 3 — Wrong Sign. Let \( b_p>0 \), but implement \( \hat s_b=-1 \). Derive the uncancelled term.

Solution. Substituting the wrong sign gives

\[ \dot V = -\frac{1}{2}\mathbf{e}^{T}Q\mathbf{e} + \left(b_p+|b_p|\right) \sigma_e \tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega}. \]

Since \( b_p=|b_p| \),

\[ \dot V = -\frac{1}{2}\mathbf{e}^{T}Q\mathbf{e} + 2b_p\sigma_e \tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega}. \]

The second term has no fixed sign, so negative semidefiniteness cannot be established.

Problem 4 — Verify the Lyapunov Matrix. Verify that

\[ A_m= \begin{bmatrix} 0 & 1\\ -4 & -4 \end{bmatrix}, \qquad P= \begin{bmatrix} 9/8 & 1/8\\ 1/8 & 5/32 \end{bmatrix} \]

satisfy \( A_m^{T}P+PA_m=-I_2 \).

Solution. Direct multiplication gives

\[ A_m^{T}P = \begin{bmatrix} -1/2 & -5/8\\ 5/8 & -1/2 \end{bmatrix}, \qquad PA_m = \begin{bmatrix} -1/2 & 5/8\\ -5/8 & -1/2 \end{bmatrix}. \]

Adding the matrices yields

\[ A_m^{T}P+PA_m = \begin{bmatrix} -1 & 0\\ 0 & -1 \end{bmatrix} = -I_2. \]

Problem 5 — Simulation Study. Run one implementation with the correct sign. Record the final tracking-error norm and final parameter estimates. Explain why small tracking error does not require the estimates to equal \( [-1.5,-1.4,2]^{T} \).

Solution. With the supplied settings, the numerical solution produces a very small final tracking error while the final estimate generally differs from the ideal matching vector. The Lyapunov proof establishes bounded parameters and asymptotic tracking; it does not establish unique parameter identification for an arbitrary bounded command. Multiple estimate values can generate the required behavior along the realized trajectory. The stronger condition that removes this ambiguity is studied in Chapter 10.

In canonical-form SISO MRAC, the unknown input gain multiplies the parameter-error channel. When the gain belongs to a sign-definite set, its known sign is inserted into the update law: \( \dot{\hat{\boldsymbol{\theta}}} =-s_b\Gamma\boldsymbol{\omega} \mathbf{e}^{T}P\mathbf{b}_0 \) . A Lyapunov function weighted by the unknown but constant magnitude \( |b_p| \) proves boundedness and asymptotic tracking without requiring that magnitude in software. The correct sign is essential: reversing it leaves a sign-indefinite mixed term and destroys the standard cancellation argument. Truly unknown control direction is a separate adaptive-control problem.