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

• formulate a full-state MIMO MRAC problem with matched parametric uncertainty;
• test exact plant-model matching by range-space and left-annihilator conditions;
• derive a matrix tracking-error model using one stacked regressor;
• construct a trace-form Lyapunov function and cancel all matrix cross terms;
• separate tracking convergence from adaptive-parameter convergence;
• explain why an unknown scalar control sign becomes a matrix-factorization problem in MIMO systems;
• recognize when unmatched uncertainty or insufficient actuation invalidates classical MRAC.

The lesson assumes the scalar and SISO MRAC material of Chapters 4–10, robust modifications of Chapters 8–9, and the MIMO uncertainty and controller parameterizations introduced in Lessons 1 and 2 of this chapter.

Consider the uncertain plant

\[ \dot{\mathbf{x} } = \mathbf{A}\mathbf{x} + \mathbf{B}\boldsymbol{\Lambda} \left(\mathbf{u}+\boldsymbol{\Theta}^{*T}\boldsymbol{\phi}(\mathbf{x})\right), \qquad \mathbf{y}=\mathbf{C}\mathbf{x}. \]

Here \( \mathbf{x}\in\mathbb{R}^{n} \), \( \mathbf{u}\in\mathbb{R}^{m} \), \( \boldsymbol{\phi}\in\mathbb{R}^{p} \), and \( \boldsymbol{\Theta}^{*}\in\mathbb{R}^{p\times m} \). The matrix \( \boldsymbol{\Lambda}\in\mathbb{R}^{m\times m} \) represents uncertain input effectiveness and interchannel scaling. For the basic theorem below, it is assumed constant, nonsingular, symmetric, and positive definite. This is the matrix counterpart of knowing the control direction in a SISO plant.

The desired reference model is

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

The core assumptions are:

1. Full-state availability: \( \mathbf{x} \) is measured.
2. Known input distribution: \( \mathbf{B} \) is known and has full column rank.
3. Matched uncertainty: the unknown term enters through \( \mathbf{B}\boldsymbol{\Lambda} \).
4. Known basis: \( \boldsymbol{\phi}(\mathbf{x}) \) is known and locally Lipschitz.
5. Exact matching: constant ideal controller matrices exist.

flowchart TD
  R["Command vector r"] --> RM["Stable reference model"]
  RM --> XM["Reference state xm"]
  X["Measured plant state x"] --> E["Tracking error e = x - xm"]
  XM --> E
  X --> REG["Stacked regressor omega"]
  R --> REG
  REG --> CTRL["Adaptive matrix controller"]
  E --> LAW["Matrix update law"]
  LAW --> CTRL
  CTRL --> U["Control vector u"]
  U --> PLANT["Coupled MIMO plant"]
  PLANT --> X
        

Choose the ideal control law

\[ \mathbf{u}^{*}=\mathbf{K}_x^{*T}\mathbf{x} +\mathbf{K}_r^{*T}\mathbf{r} -\boldsymbol{\Theta}^{*T}\boldsymbol{\phi}(\mathbf{x}). \]

Substitution into the plant gives the desired model only if

\[ \mathbf{A}+\mathbf{B}\boldsymbol{\Lambda}\mathbf{K}_x^{*T}=\mathbf{A}_m, \qquad \mathbf{B}\boldsymbol{\Lambda}\mathbf{K}_r^{*T}=\mathbf{B}_m. \]

Because \( \boldsymbol{\Lambda} \) is nonsingular, these equations require the columns of \( \mathbf{A}_m-\mathbf{A} \) and \( \mathbf{B}_m \) to lie in the range of \( \mathbf{B} \):

\[ \operatorname{im}(\mathbf{A}_m-\mathbf{A})\subseteq\operatorname{im}(\mathbf{B}), \qquad \operatorname{im}(\mathbf{B}_m)\subseteq\operatorname{im}(\mathbf{B}). \]

Let \( \mathbf{B}_{\perp} \) be any full-column-rank basis for the left nullspace of \( \mathbf{B} \), so that \( \mathbf{B}_{\perp}^{T}\mathbf{B}=\mathbf{0} \). Necessary and sufficient matching tests are

\[ \mathbf{B}_{\perp}^{T}(\mathbf{A}_m-\mathbf{A})=\mathbf{0}, \qquad \mathbf{B}_{\perp}^{T}\mathbf{B}_m=\mathbf{0}. \]

This result is decisive for underactuated plants. When \( m<n \), an arbitrary stable \( \mathbf{A}_m \) cannot generally be assigned. Adaptation can estimate unknown coefficients inside the actuated subspace, but it cannot create control authority in an unactuated state direction.

If \( \mathbf{B} \) is square and nonsingular, the ideal gains are

\[ \mathbf{K}_x^{*T}=\boldsymbol{\Lambda}^{-1}\mathbf{B}^{-1} (\mathbf{A}_m-\mathbf{A}), \qquad \mathbf{K}_r^{*T}=\boldsymbol{\Lambda}^{-1}\mathbf{B}^{-1}\mathbf{B}_m. \]

The adaptive control law is

\[ \mathbf{u}=\widehat{\mathbf{K} }_x^{T}\mathbf{x} +\widehat{\mathbf{K} }_r^{T}\mathbf{r} -\widehat{\boldsymbol{\Theta} }^{T}\boldsymbol{\phi}(\mathbf{x}). \]

Define the stacked regressor and controller matrix

\[ \boldsymbol{\omega}=\operatorname{col} \left(\mathbf{x},\mathbf{r},\boldsymbol{\phi}(\mathbf{x})\right) \in\mathbb{R}^{s},\qquad s=n+q+p, \]

\[ \widehat{\mathbf{W} }= \begin{bmatrix} \widehat{\mathbf{K} }_x\\ \widehat{\mathbf{K} }_r\\ -\widehat{\boldsymbol{\Theta} } \end{bmatrix}\in\mathbb{R}^{s\times m}, \qquad \mathbf{u}=\widehat{\mathbf{W} }^{T}\boldsymbol{\omega}. \]

The compact form is not merely notation. It exposes the central MIMO adaptation object: every regressor component has one coefficient for each input channel. Thus an \( s \)-component regressor produces \( sm \) adaptive coefficients. For large systems, this quadratic growth motivates structured, sparse, decentralized, or reduced-order parameterizations.

With

\[ \widetilde{\mathbf{W} }=\widehat{\mathbf{W} }-\mathbf{W}^{*}, \qquad \mathbf{e}=\mathbf{x}-\mathbf{x}_m, \]

the exact matching equations yield the error model

\[ \dot{\mathbf{e} }=\mathbf{A}_m\mathbf{e} +\mathbf{B}\boldsymbol{\Lambda}\widetilde{\mathbf{W} }^{T} \boldsymbol{\omega}. \]

Select arbitrary \( \mathbf{Q}=\mathbf{Q}^{T}\succ\mathbf{0} \) and solve

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

Let \( \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}\succ\mathbf{0} \) be the adaptation-gain matrix. Use the update law

\[ \dot{\widehat{\mathbf{W} } } =-\boldsymbol{\Gamma}\boldsymbol{\omega}\mathbf{e}^{T} \mathbf{P}\mathbf{B}. \]

Theorem. Under the assumptions of Section 2, for bounded reference input, all closed-loop signals are bounded and the tracking error satisfies \( \lim_{t\to\infty}\mathbf{e}(t)=\mathbf{0} \).

Proof. Consider

\[ V=\mathbf{e}^{T}\mathbf{P}\mathbf{e} +\operatorname{tr}\!\left( \widetilde{\mathbf{W} }^{T}\boldsymbol{\Gamma}^{-1} \widetilde{\mathbf{W} }\boldsymbol{\Lambda}\right). \]

The parameter term is positive because it can be written as

\[ \operatorname{tr}\!\left( \widetilde{\mathbf{W} }^{T}\boldsymbol{\Gamma}^{-1} \widetilde{\mathbf{W} }\boldsymbol{\Lambda}\right) =\left\|\boldsymbol{\Gamma}^{-1/2} \widetilde{\mathbf{W} }\boldsymbol{\Lambda}^{1/2}\right\|_{F}^{2}. \]

Differentiating and using the error equation gives

\[ \begin{aligned} \dot V={}&-\mathbf{e}^{T}\mathbf{Q}\mathbf{e} +2\mathbf{e}^{T}\mathbf{P}\mathbf{B}\boldsymbol{\Lambda} \widetilde{\mathbf{W} }^{T}\boldsymbol{\omega}\\ &+2\operatorname{tr}\!\left( \widetilde{\mathbf{W} }^{T}\boldsymbol{\Gamma}^{-1} \dot{\widehat{\mathbf{W} } }\boldsymbol{\Lambda}\right). \end{aligned} \]

The scalar cross term can be expressed through the trace identity as

\[ \mathbf{e}^{T}\mathbf{P}\mathbf{B}\boldsymbol{\Lambda} \widetilde{\mathbf{W} }^{T}\boldsymbol{\omega} =\operatorname{tr}\!\left( \widetilde{\mathbf{W} }^{T}\boldsymbol{\omega}\mathbf{e}^{T} \mathbf{P}\mathbf{B}\boldsymbol{\Lambda}\right). \]

Substituting the adaptive law cancels the two matrix cross terms exactly:

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

Therefore \( \mathbf{e} \) and \( \widetilde{\mathbf{W} } \) are bounded and \( \mathbf{e}\in L_2 \). Since \( \mathbf{A}_m \) is Hurwitz and \( \mathbf{r} \) is bounded, \( \mathbf{x}_m \) is bounded. Hence \( \mathbf{x}=\mathbf{e}+\mathbf{x}_m \), the regressor, control, and error derivative are bounded. Barbalat's lemma then gives \( \mathbf{e}(t)\to\mathbf{0} \). \(\square\)

6.1 Matrix control direction

In SISO MRAC, the known sign of the high-frequency gain determines the sign of the update law. In MIMO systems, the analogous object is a nonsingular matrix. Its off-diagonal entries mix input channels, so there is generally no single sign that can be inserted into the gradient law.

6.2 Cross-coupled adaptive coefficients

The coefficient associated with one state can affect every actuator. A diagonal adaptive gain does not imply a decentralized controller because \( \widehat{\mathbf{W} } \) itself is generally dense. Ignoring this structure can leave uncancelled cross terms in the Lyapunov derivative.

6.3 Rank and right-invertibility

Exact tracking requires sufficient control authority and compatible relative degree. State-feedback matching is expressed by the range conditions in Section 3. Output-feedback MRAC additionally requires suitable minimum-phase, interactor, and high-frequency-gain properties.

6.4 Overparameterization and nonuniqueness

If the same input-output behavior can be produced by several controller matrices, the ideal parameter is not unique. Tracking may converge while the parameter matrix approaches a set rather than one physical value. Redundant basis functions worsen this effect.

6.5 Unmatched uncertainty

Suppose the actual plant contains

\[ \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} +\mathbf{B}\boldsymbol{\Lambda}(\mathbf{u}+\boldsymbol{\Theta}^{*T}\boldsymbol{\phi}) +\mathbf{D}\boldsymbol{\Delta}(\mathbf{x}), \]

where \( \operatorname{im}(\mathbf{D}) \) is not contained in \( \operatorname{im}(\mathbf{B}) \). The additional term cannot be represented by \( \widetilde{\mathbf{W} }^{T}\boldsymbol{\omega} \), so the classical cancellation proof fails. Robust modifications may establish ultimate boundedness, but exact asymptotic model following is no longer automatic.

The symmetric positive-definite assumption on \( \boldsymbol{\Lambda} \) makes the trace Lyapunov term positive and allows direct cancellation. A general nonsymmetric high-frequency gain matrix requires additional structure. One influential route uses an LDS factorization

\[ \mathbf{K}_p=\mathbf{L}\mathbf{D}\mathbf{S}, \]

where \( \mathbf{L} \) is unit triangular, \( \mathbf{D} \) is diagonal with known signs, and \( \mathbf{S}=\mathbf{S}^{T}\succ\mathbf{0} \). Under nonzero leading-principal-minor conditions, this decomposition converts the unknown matrix direction into triangularly ordered scalar-like directions and an SPD factor suitable for Lyapunov analysis.

The resulting adaptive controller is more complicated than the basic law of Section 5: triangular coupling terms must be parameterized and updated in an order consistent with \( \mathbf{L} \). The essential lesson is structural: for a general MIMO plant, stability depends on a usable factorization or sign pattern of the high-frequency gain matrix, not only on its nonsingularity.

The Lyapunov proof establishes

\[ \mathbf{e}(t)\to\mathbf{0}, \qquad \widehat{\mathbf{W} }(t)\in L_{\infty}, \]

but it does not by itself establish

\[ \widetilde{\mathbf{W} }(t)\to\mathbf{0}. \]

A standard matrix persistent-excitation condition is

\[ \int_{t}^{t+T}\boldsymbol{\omega}(\tau) \boldsymbol{\omega}^{T}(\tau)\,d\tau\succeq\alpha\mathbf{I}_{s}, \qquad T>0,\;\alpha>0. \]

Even when this condition holds, parameter convergence additionally requires an identifiable parameterization: different columns or basis functions must not generate indistinguishable closed-loop effects. In practical multivariable systems, commands designed only for good tracking often fail to excite every coupled parameter direction. Consequently, a simulation should not claim identification merely because the tracking error becomes small.

The implementations use

\[ \mathbf{A}=\begin{bmatrix}0&1\\-1&-0.4\end{bmatrix},\quad \mathbf{B}=\mathbf{I}_2,\quad \boldsymbol{\Lambda}=\operatorname{diag}(1.2,0.8), \]

\[ \mathbf{A}_m=\begin{bmatrix}-1&0.3\\-0.2&-1.4\end{bmatrix},\quad \mathbf{B}_m=\mathbf{I}_2,\quad \boldsymbol{\Theta}^{*}=\begin{bmatrix}0.25&-0.10\\0.15&0.20\end{bmatrix}. \]

The nonlinear basis and multisine reference are

\[ \boldsymbol{\phi}(\mathbf{x})= \begin{bmatrix}\sin x_1\\x_2/(1+x_2^2)\end{bmatrix}, \]

\[ \mathbf{r}(t)=\begin{bmatrix} 0.8\sin(0.6t)+0.35\sin(1.7t)\\ 0.7\cos(0.4t)+0.25\sin(1.3t) \end{bmatrix}. \]

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

\[ \mathbf{P}=\begin{bmatrix} 0.49372146&0.03139269\\0.03139269&0.36386986 \end{bmatrix}. \]

All parameter estimates start at zero and \( \boldsymbol{\Gamma}=4\mathbf{I}_6 \). The example is intentionally coupled, but it satisfies the SPD input-gain assumption. Its purpose is to verify the stability mechanism, not to prove that the selected finite-time command is persistently exciting for every controller coefficient.

Chapter23_Lesson3.py

"""Chapter23_Lesson3.py

Full-state MIMO MRAC simulation for a coupled 2x2 plant with matched
nonlinear uncertainty and an unknown positive-definite input-gain matrix.

Dependencies:
    pip install numpy matplotlib
"""

from __future__ import annotations

import numpy as np
import matplotlib.pyplot as plt

A = np.array([[0.0, 1.0], [-1.0, -0.4]], dtype=float)
B = np.eye(2)
LAMBDA = np.diag([1.2, 0.8])
A_M = np.array([[-1.0, 0.3], [-0.2, -1.4]], dtype=float)
B_M = np.eye(2)
THETA_STAR = np.array([[0.25, -0.10], [0.15, 0.20]], dtype=float)
P = np.array([[0.49372146, 0.03139269],
              [0.03139269, 0.36386986]], dtype=float)
GAMMA = 4.0 * np.eye(6)


def reference(t: float) -> np.ndarray:
    """Bounded multisine command with both input channels active."""
    return np.array([
        0.8 * np.sin(0.6 * t) + 0.35 * np.sin(1.7 * t),
        0.7 * np.cos(0.4 * t) + 0.25 * np.sin(1.3 * t),
    ])


def basis(x: np.ndarray) -> np.ndarray:
    """Known basis multiplying the unknown matched coefficient matrix."""
    return np.array([np.sin(x[0]), x[1] / (1.0 + x[1] ** 2)])


def rhs(t: float, z: np.ndarray) -> np.ndarray:
    """Augmented plant, reference-model, and adaptive-parameter dynamics."""
    x = z[0:2]
    x_m = z[2:4]
    w_hat = z[4:].reshape(6, 2)

    r = reference(t)
    phi = basis(x)
    omega = np.concatenate((x, r, phi))
    u = w_hat.T @ omega

    x_dot = A @ x + B @ LAMBDA @ (u + THETA_STAR.T @ phi)
    x_m_dot = A_M @ x_m + B_M @ r

    e = x - x_m
    w_hat_dot = -GAMMA @ np.outer(omega, e.T @ P @ B)
    return np.concatenate((x_dot, x_m_dot, w_hat_dot.reshape(-1)))


def rk4_step(t: float, z: np.ndarray, dt: float) -> np.ndarray:
    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):
    steps = int(round(t_final / dt))
    z = np.zeros(16)
    z[0:2] = np.array([0.6, -0.4])

    time = np.linspace(0.0, t_final, steps + 1)
    x_hist = np.zeros((steps + 1, 2))
    xm_hist = np.zeros((steps + 1, 2))
    u_hist = np.zeros((steps + 1, 2))
    w_norm_hist = np.zeros(steps + 1)

    for k, t in enumerate(time):
        x = z[0:2]
        x_m = z[2:4]
        w_hat = z[4:].reshape(6, 2)
        omega = np.concatenate((x, reference(t), basis(x)))

        x_hist[k] = x
        xm_hist[k] = x_m
        u_hist[k] = w_hat.T @ omega
        w_norm_hist[k] = np.linalg.norm(w_hat)

        if k < steps:
            z = rk4_step(t, z, dt)

    return time, x_hist, xm_hist, u_hist, w_norm_hist


def main() -> None:
    time, x, x_m, u, w_norm = simulate()
    error = x - x_m

    print(f"Final tracking-error norm: {np.linalg.norm(error[-1]):.6f}")
    print(f"Maximum tracking-error norm: {np.max(np.linalg.norm(error, axis=1)):.6f}")
    print(f"Final adaptive-matrix Frobenius norm: {w_norm[-1]:.6f}")

    fig, axes = plt.subplots(3, 1, figsize=(10, 10), sharex=True)
    axes[0].plot(time, x[:, 0], label="x1")
    axes[0].plot(time, x_m[:, 0], "--", label="xm1")
    axes[0].plot(time, x[:, 1], label="x2")
    axes[0].plot(time, x_m[:, 1], "--", label="xm2")
    axes[0].set_ylabel("states")
    axes[0].grid(True)
    axes[0].legend(ncol=2)

    axes[1].plot(time, error[:, 0], label="e1")
    axes[1].plot(time, error[:, 1], label="e2")
    axes[1].plot(time, np.linalg.norm(error, axis=1), label="||e||")
    axes[1].set_ylabel("tracking error")
    axes[1].grid(True)
    axes[1].legend()

    axes[2].plot(time, u[:, 0], label="u1")
    axes[2].plot(time, u[:, 1], label="u2")
    axes[2].set_xlabel("time [s]")
    axes[2].set_ylabel("control")
    axes[2].grid(True)
    axes[2].legend()

    fig.tight_layout()
    fig.savefig("Chapter23_Lesson3_results.png", dpi=180)
    plt.show()


if __name__ == "__main__":
    main()

Chapter23_Lesson3.cpp

// Chapter23_Lesson3.cpp
// Full-state 2x2 MIMO MRAC with matched uncertainty and RK4 integration.

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

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

constexpr double DT = 0.002;
constexpr double T_FINAL = 20.0;
constexpr double GAMMA = 4.0;

constexpr double A[2][2] = { {0.0, 1.0}, {-1.0, -0.4} };
constexpr double AM[2][2] = { {-1.0, 0.3}, {-0.2, -1.4} };
constexpr double LAMBDA_DIAG[2] = {1.2, 0.8};
constexpr double THETA_STAR[2][2] = { {0.25, -0.10}, {0.15, 0.20} };
constexpr double P[2][2] = { {0.49372146, 0.03139269},
                            {0.03139269, 0.36386986} };

std::array<double, 2> reference(double t) {
    return {0.8 * std::sin(0.6 * t) + 0.35 * std::sin(1.7 * t),
            0.7 * std::cos(0.4 * t) + 0.25 * std::sin(1.3 * t)};
}

std::array<double, 2> basis(const std::array<double, 2>& x) {
    return {std::sin(x[0]), x[1] / (1.0 + x[1] * x[1])};
}

State rhs(double t, const State& z) {
    State dz{};
    const std::array<double, 2> x = {z[0], z[1]};
    const std::array<double, 2> xm = {z[2], z[3]};
    const auto r = reference(t);
    const auto phi = basis(x);
    const std::array<double, 6> omega = {x[0], x[1], r[0], r[1], phi[0], phi[1]};

    std::array<double, 2> u{};
    for (int j = 0; j < 2; ++j) {
        for (int i = 0; i < 6; ++i) {
            u[j] += z[4 + 2 * i + j] * omega[i];
        }
    }

    std::array<double, 2> uncertainty{};
    for (int j = 0; j < 2; ++j) {
        for (int i = 0; i < 2; ++i) {
            uncertainty[j] += THETA_STAR[i][j] * phi[i];
        }
    }

    for (int i = 0; i < 2; ++i) {
        dz[i] = A[i][0] * x[0] + A[i][1] * x[1]
                + LAMBDA_DIAG[i] * (u[i] + uncertainty[i]);
        dz[2 + i] = AM[i][0] * xm[0] + AM[i][1] * xm[1] + r[i];
    }

    const std::array<double, 2> e = {x[0] - xm[0], x[1] - xm[1]};
    std::array<double, 2> eTP{};
    for (int j = 0; j < 2; ++j) {
        eTP[j] = e[0] * P[0][j] + e[1] * P[1][j];
    }

    for (int i = 0; i < 6; ++i) {
        for (int j = 0; j < 2; ++j) {
            dz[4 + 2 * i + j] = -GAMMA * omega[i] * eTP[j];
        }
    }
    return dz;
}

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 rk4_step(double t, const State& z) {
    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 < z.size(); ++i) {
        next[i] = z[i] + (DT / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
    }
    return next;
}

int main() {
    State z{};
    z[0] = 0.6;
    z[1] = -0.4;

    std::ofstream csv("Chapter23_Lesson3_cpp.csv");
    if (!csv) {
        throw std::runtime_error("Cannot open output CSV file.");
    }
    csv << "t,x1,x2,xm1,xm2,e1,e2,u1,u2\n";

    const int steps = static_cast<int>(std::round(T_FINAL / DT));
    double max_error = 0.0;

    for (int k = 0; k <= steps; ++k) {
        const double t = k * DT;
        const auto r = reference(t);
        const std::array<double, 2> x = {z[0], z[1]};
        const auto phi = basis(x);
        const std::array<double, 6> omega = {x[0], x[1], r[0], r[1], phi[0], phi[1]};
        std::array<double, 2> u{};
        for (int j = 0; j < 2; ++j) {
            for (int i = 0; i < 6; ++i) {
                u[j] += z[4 + 2 * i + j] * omega[i];
            }
        }

        const double e1 = z[0] - z[2];
        const double e2 = z[1] - z[3];
        max_error = std::max(max_error, std::hypot(e1, e2));

        if (k % 10 == 0) {
            csv << std::setprecision(10) << t << ',' << z[0] << ',' << z[1] << ','
                << z[2] << ',' << z[3] << ',' << e1 << ',' << e2 << ','
                << u[0] << ',' << u[1] << '\n';
        }
        if (k < steps) {
            z = rk4_step(t, z);
        }
    }

    const double final_error = std::hypot(z[0] - z[2], z[1] - z[3]);
    std::cout << "Final tracking-error norm: " << final_error << '\n';
    std::cout << "Maximum tracking-error norm: " << max_error << '\n';
    std::cout << "Wrote Chapter23_Lesson3_cpp.csv\n";
    return 0;
}

Chapter23_Lesson3.java

// Chapter23_Lesson3.java
// Full-state 2x2 MIMO MRAC with matched uncertainty and RK4 integration.

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Locale;

public final class Chapter23_Lesson3 {
    private static final double DT = 0.002;
    private static final double T_FINAL = 20.0;
    private static final double GAMMA = 4.0;

    private static final double[][] A = { {0.0, 1.0}, {-1.0, -0.4} };
    private static final double[][] AM = { {-1.0, 0.3}, {-0.2, -1.4} };
    private static final double[] LAMBDA_DIAG = {1.2, 0.8};
    private static final double[][] THETA_STAR = { {0.25, -0.10}, {0.15, 0.20} };
    private static final double[][] P = {
        {0.49372146, 0.03139269},
        {0.03139269, 0.36386986}
    };

    private Chapter23_Lesson3() {}

    private static double[] reference(double t) {
        return new double[] {
            0.8 * Math.sin(0.6 * t) + 0.35 * Math.sin(1.7 * t),
            0.7 * Math.cos(0.4 * t) + 0.25 * Math.sin(1.3 * t)
        };
    }

    private static double[] basis(double[] x) {
        return new double[] {Math.sin(x[0]), x[1] / (1.0 + x[1] * x[1])};
    }

    private static double[] rhs(double t, double[] z) {
        double[] dz = new double[16];
        double[] x = {z[0], z[1]};
        double[] xm = {z[2], z[3]};
        double[] r = reference(t);
        double[] phi = basis(x);
        double[] omega = {x[0], x[1], r[0], r[1], phi[0], phi[1]};

        double[] u = new double[2];
        for (int j = 0; j < 2; j++) {
            for (int i = 0; i < 6; i++) {
                u[j] += z[4 + 2 * i + j] * omega[i];
            }
        }

        double[] uncertainty = new double[2];
        for (int j = 0; j < 2; j++) {
            for (int i = 0; i < 2; i++) {
                uncertainty[j] += THETA_STAR[i][j] * phi[i];
            }
        }

        for (int i = 0; i < 2; i++) {
            dz[i] = A[i][0] * x[0] + A[i][1] * x[1]
                    + LAMBDA_DIAG[i] * (u[i] + uncertainty[i]);
            dz[2 + i] = AM[i][0] * xm[0] + AM[i][1] * xm[1] + r[i];
        }

        double[] e = {x[0] - xm[0], x[1] - xm[1]};
        double[] eTP = new double[2];
        for (int j = 0; j < 2; j++) {
            eTP[j] = e[0] * P[0][j] + e[1] * P[1][j];
        }

        for (int i = 0; i < 6; i++) {
            for (int j = 0; j < 2; j++) {
                dz[4 + 2 * i + j] = -GAMMA * omega[i] * eTP[j];
            }
        }
        return dz;
    }

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

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        double[] z = new double[16];
        z[0] = 0.6;
        z[1] = -0.4;
        int steps = (int) Math.round(T_FINAL / DT);
        double maxError = 0.0;

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter23_Lesson3_java.csv"))) {
            writer.write("t,x1,x2,xm1,xm2,e1,e2,u1,u2\n");

            for (int k = 0; k <= steps; k++) {
                double t = k * DT;
                double[] r = reference(t);
                double[] x = {z[0], z[1]};
                double[] phi = basis(x);
                double[] omega = {x[0], x[1], r[0], r[1], phi[0], phi[1]};
                double[] u = new double[2];
                for (int j = 0; j < 2; j++) {
                    for (int i = 0; i < 6; i++) {
                        u[j] += z[4 + 2 * i + j] * omega[i];
                    }
                }

                double e1 = z[0] - z[2];
                double e2 = z[1] - z[3];
                maxError = Math.max(maxError, Math.hypot(e1, e2));

                if (k % 10 == 0) {
                    writer.write(String.format(Locale.US,
                        "%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f%n",
                        t, z[0], z[1], z[2], z[3], e1, e2, u[0], u[1]));
                }
                if (k < steps) {
                    z = rk4Step(t, z);
                }
            }
        }

        double finalError = Math.hypot(z[0] - z[2], z[1] - z[3]);
        System.out.printf(Locale.US, "Final tracking-error norm: %.6f%n", finalError);
        System.out.printf(Locale.US, "Maximum tracking-error norm: %.6f%n", maxError);
        System.out.println("Wrote Chapter23_Lesson3_java.csv");
    }
}

Chapter23_Lesson3.m

function Chapter23_Lesson3
% Chapter23_Lesson3.m
% Full-state 2x2 MIMO MRAC with matched uncertainty and RK4 integration.

A = [0, 1; -1, -0.4];
B = eye(2);
Lambda = diag([1.2, 0.8]);
Am = [-1, 0.3; -0.2, -1.4];
Bm = eye(2);
ThetaStar = [0.25, -0.10; 0.15, 0.20];
P = [0.49372146, 0.03139269; 0.03139269, 0.36386986];
gamma = 4.0;

dt = 0.002;
tFinal = 20.0;
time = 0:dt:tFinal;
z = zeros(16, 1);
z(1:2) = [0.6; -0.4];

xHistory = zeros(2, numel(time));
xmHistory = zeros(2, numel(time));
uHistory = zeros(2, numel(time));
wNormHistory = zeros(1, numel(time));

for k = 1:numel(time)
    t = time(k);
    x = z(1:2);
    xm = z(3:4);
    What = reshape(z(5:end), 6, 2);
    r = referenceSignal(t);
    phi = basisVector(x);
    omega = [x; r; phi];

    xHistory(:, k) = x;
    xmHistory(:, k) = xm;
    uHistory(:, k) = What' * omega;
    wNormHistory(k) = norm(What, 'fro');

    if k < numel(time)
        k1 = dynamics(t, z);
        k2 = dynamics(t + dt/2, z + dt*k1/2);
        k3 = dynamics(t + dt/2, z + dt*k2/2);
        k4 = dynamics(t + dt, z + dt*k3);
        z = z + dt*(k1 + 2*k2 + 2*k3 + k4)/6;
    end
end

errorHistory = xHistory - xmHistory;
fprintf('Final tracking-error norm: %.6f\n', norm(errorHistory(:, end)));
fprintf('Maximum tracking-error norm: %.6f\n', max(vecnorm(errorHistory)));
fprintf('Final adaptive-matrix Frobenius norm: %.6f\n', wNormHistory(end));

figure('Name', 'Chapter 23 Lesson 3: MIMO MRAC');
subplot(3,1,1);
plot(time, xHistory(1,:), time, xmHistory(1,:), '--', ...
     time, xHistory(2,:), time, xmHistory(2,:), '--');
grid on; ylabel('states'); legend('x_1','x_{m1}','x_2','x_{m2}');

subplot(3,1,2);
plot(time, errorHistory(1,:), time, errorHistory(2,:), ...
     time, vecnorm(errorHistory));
grid on; ylabel('tracking error'); legend('e_1','e_2','||e||');

subplot(3,1,3);
plot(time, uHistory(1,:), time, uHistory(2,:));
grid on; xlabel('time [s]'); ylabel('control'); legend('u_1','u_2');

saveas(gcf, 'Chapter23_Lesson3_matlab.png');

    function dz = dynamics(t, zLocal)
        xLocal = zLocal(1:2);
        xmLocal = zLocal(3:4);
        WhatLocal = reshape(zLocal(5:end), 6, 2);
        rLocal = referenceSignal(t);
        phiLocal = basisVector(xLocal);
        omegaLocal = [xLocal; rLocal; phiLocal];
        uLocal = WhatLocal' * omegaLocal;

        xDot = A*xLocal + B*Lambda*(uLocal + ThetaStar'*phiLocal);
        xmDot = Am*xmLocal + Bm*rLocal;
        eLocal = xLocal - xmLocal;
        WhatDot = -gamma * omegaLocal * (eLocal' * P * B);
        dz = [xDot; xmDot; WhatDot(:)];
    end

    function r = referenceSignal(t)
        r = [0.8*sin(0.6*t) + 0.35*sin(1.7*t); ...
             0.7*cos(0.4*t) + 0.25*sin(1.3*t)];
    end

    function phi = basisVector(x)
        phi = [sin(x(1)); x(2)/(1 + x(2)^2)];
    end
end

Simulink realization

A direct Simulink implementation uses four subsystems:

1. Reference Model: a State-Space block with \( \mathbf{A}_m,\mathbf{B}_m,\mathbf{I},\mathbf{0} \).
2. Plant: a MATLAB Function block computing \( \mathbf{A}\mathbf{x}+\mathbf{B}\boldsymbol{\Lambda}(\mathbf{u}+\boldsymbol{\Theta}^{*T}\boldsymbol{\phi}) \), followed by an Integrator.
3. Adaptive Controller: form \( \boldsymbol{\omega} \), reshape the adaptive state into \( \widehat{\mathbf{W} } \), and compute \( \mathbf{u}=\widehat{\mathbf{W} }^T\boldsymbol{\omega} \).
4. Update Law: compute \( -\boldsymbol{\Gamma}\boldsymbol{\omega}\mathbf{e}^{T}\mathbf{P}\mathbf{B} \), vectorize it, and integrate its \( sm \) entries.

Use a continuous solver such as ode45 for an initial study. For a real-time fixed-step implementation, reduce the adaptation gain or add normalization and projection before decreasing the sample period.

Chapter23_Lesson3.nb

Notebook[{
 Cell["Chapter 23, Lesson 3: MRAC for MIMO Plants", "Title"],
 Cell["Full-state 2x2 MIMO MRAC with matched uncertainty and RK4 integration.", "Text"],
 Cell[BoxData[
"ClearAll[\"Global`*\"];
A = { {0., 1.}, {-1., -0.4} };
B = IdentityMatrix[2];
lambda = DiagonalMatrix[{1.2, 0.8}];
Am = { {-1., 0.3}, {-0.2, -1.4} };
Bm = IdentityMatrix[2];
thetaStar = { {0.25, -0.10}, {0.15, 0.20} };
pMat = { {0.49372146, 0.03139269}, {0.03139269, 0.36386986} };
gamma = 4.0 IdentityMatrix[6];

reference[t_] := {0.8 Sin[0.6 t] + 0.35 Sin[1.7 t],
  0.7 Cos[0.4 t] + 0.25 Sin[1.3 t]};
basis[x_] := {Sin[x[[1]]], x[[2]]/(1 + x[[2]]^2)};

rhs[t_, z_] := Module[{x, xm, wHat, r, phi, omega, u, xDot, xmDot,
   e, wDot},
  x = z[[1 ;; 2]];
  xm = z[[3 ;; 4]];
  wHat = ArrayReshape[z[[5 ;;]], {6, 2}];
  r = reference[t];
  phi = basis[x];
  omega = Join[x, r, phi];
  u = Transpose[wHat].omega;
  xDot = A.x + B.lambda.(u + Transpose[thetaStar].phi);
  xmDot = Am.xm + Bm.r;
  e = x - xm;
  wDot = -gamma.Outer[Times, omega, e.pMat.B];
  Join[xDot, xmDot, Flatten[wDot]]
];

rk4Step[t_, z_, dt_] := Module[{k1, k2, k3, k4},
  k1 = rhs[t, z];
  k2 = rhs[t + dt/2, z + dt k1/2];
  k3 = rhs[t + dt/2, z + dt k2/2];
  k4 = rhs[t + dt, z + dt k3];
  z + dt (k1 + 2 k2 + 2 k3 + k4)/6
];

dt = 0.002; tFinal = 20.;
times = Range[0., tFinal, dt];
z0 = Join[{0.6, -0.4}, {0., 0.}, ConstantArray[0., 12]];
states = Rest@FoldList[rk4Step[#2 - dt, #1, dt] &, z0, Rest[times]];
states = Prepend[states, z0];

xHistory = states[[All, 1 ;; 2]];
xmHistory = states[[All, 3 ;; 4]];
errorHistory = xHistory - xmHistory;
uHistory = MapThread[
  Function[{t, z}, Module[{x, wHat, omega},
    x = z[[1 ;; 2]];
    wHat = ArrayReshape[z[[5 ;;]], {6, 2}];
    omega = Join[x, reference[t], basis[x]];
    Transpose[wHat].omega]], {times, states}];

Print[\"Final tracking-error norm: \", Norm[Last[errorHistory]]];
Print[\"Maximum tracking-error norm: \", Max[Norm /@ errorHistory]];

GraphicsGrid[{ {
  ListLinePlot[{Transpose[{times, xHistory[[All, 1]]}],
    Transpose[{times, xmHistory[[All, 1]]}],
    Transpose[{times, xHistory[[All, 2]]}],
    Transpose[{times, xmHistory[[All, 2]]}]},
   PlotLegends -> {\"x1\", \"xm1\", \"x2\", \"xm2\"},
   GridLines -> Automatic, ImageSize -> Large]}, {
  ListLinePlot[{Transpose[{times, errorHistory[[All, 1]]}],
    Transpose[{times, errorHistory[[All, 2]]}],
    Transpose[{times, Norm /@ errorHistory}]},
   PlotLegends -> {\"e1\", \"e2\", \"norm e\"},
   GridLines -> Automatic, ImageSize -> Large]}, {
  ListLinePlot[{Transpose[{times, uHistory[[All, 1]]}],
    Transpose[{times, uHistory[[All, 2]]}]},
   PlotLegends -> {\"u1\", \"u2\"}, GridLines -> Automatic,
   ImageSize -> Large]} }]
"], "Input"]
}, WindowTitle -> "Chapter23_Lesson3"]

Problem 1 — Matching test for an underactuated plant. Let

\[ \mathbf{B}=\begin{bmatrix}1\\0\end{bmatrix},\qquad \mathbf{A}_m-\mathbf{A}=\begin{bmatrix}a&b\\c&d\end{bmatrix}. \]

Find the exact state-feedback matching condition.

Solution. A basis for the left nullspace is \( \mathbf{B}_{\perp}=[0\;1]^T \). Therefore

\[ \mathbf{B}_{\perp}^{T}(\mathbf{A}_m-\mathbf{A}) =\begin{bmatrix}c&d\end{bmatrix}. \]

Exact matching requires \( c=d=0 \). The controller can change only the first row of the state matrix because its action lies in \( \operatorname{im}(\mathbf{B}) \).

Problem 2 — Derive the trace cancellation. Starting from the Lyapunov function in Section 5, prove that the adaptive law cancels the cross term.

Solution. Since \( \mathbf{W}^{*} \) is constant, \( \dot{\widetilde{\mathbf{W} } }=\dot{\widehat{\mathbf{W} } } \). Then

\[ \frac{d}{dt}\operatorname{tr}\!\left( \widetilde{\mathbf{W} }^{T}\boldsymbol{\Gamma}^{-1} \widetilde{\mathbf{W} }\boldsymbol{\Lambda}\right) =2\operatorname{tr}\!\left( \widetilde{\mathbf{W} }^{T}\boldsymbol{\Gamma}^{-1} \dot{\widehat{\mathbf{W} } }\boldsymbol{\Lambda}\right). \]

Substitution of the update law gives

\[ -2\operatorname{tr}\!\left( \widetilde{\mathbf{W} }^{T}\boldsymbol{\omega}\mathbf{e}^{T} \mathbf{P}\mathbf{B}\boldsymbol{\Lambda}\right), \]

which is the negative of the error-energy cross term after cyclic permutation under the trace. Hence only \( -\mathbf{e}^{T}\mathbf{Q}\mathbf{e} \) remains.

Problem 3 — Why positive eigenvalues are insufficient. Suppose a real nonsymmetric input-gain matrix has positive eigenvalues. Explain why the Lyapunov function of Section 5 is not automatically positive definite.

Solution. Positivity of eigenvalues does not imply symmetry or positive definiteness of the quadratic form. The trace term depends on

\[ \operatorname{tr}\!\left( \widetilde{\mathbf{W} }^{T}\boldsymbol{\Gamma}^{-1} \widetilde{\mathbf{W} }\boldsymbol{\Lambda}\right). \]

For nonsymmetric \( \boldsymbol{\Lambda} \), this expression cannot generally be represented as a squared Frobenius norm and may fail to be positive. A structural factorization such as LDS, or another known matrix direction property, is required.

Problem 4 — Parameter convergence. Assume tracking error converges to zero but the command is constant after a transient. Can all entries of \( \widehat{\mathbf{W} } \) be expected to converge to their ideal values?

Solution. Generally no. A constant command and converged state may produce a regressor confined to a low-dimensional subspace, so

\[ \int_{t}^{t+T}\boldsymbol{\omega}(\tau) \boldsymbol{\omega}^{T}(\tau)\,d\tau \]

is rank deficient. The adaptive parameters can stop at any values that produce zero prediction in the excited directions. Tracking convergence is therefore compatible with parameter error remaining nonzero.

Problem 5 — Unmatched disturbance bound. Add a bounded unmatched term \( \mathbf{d}(t) \) to the error equation:

\[ \dot{\mathbf{e} }=\mathbf{A}_m\mathbf{e} +\mathbf{B}\boldsymbol{\Lambda}\widetilde{\mathbf{W} }^{T}\boldsymbol{\omega} +\mathbf{d}(t),\qquad \|\mathbf{d}(t)\|\le\bar d. \]

Using the same adaptive law, derive an ultimate-bound argument.

Solution. The derivative becomes

\[ \dot V=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e} +2\mathbf{e}^{T}\mathbf{P}\mathbf{d}. \]

Using norm inequalities,

\[ \dot V\le-\lambda_{\min}(\mathbf{Q})\|\mathbf{e}\|^2 +2\|\mathbf{P}\|\bar d\|\mathbf{e}\|. \]

Thus \( \dot V<0 \) whenever

\[ \|\mathbf{e}\|> \frac{2\|\mathbf{P}\|\bar d}{\lambda_{\min}(\mathbf{Q})}. \]

The argument establishes ultimate boundedness, not asymptotic tracking. This is the precise stability consequence of an uncompensated unmatched term.

MIMO MRAC requires structural feasibility before adaptation begins. Exact model matching constrains the reference model to the plant's actuated subspace. Under matched uncertainty and an SPD input-effectiveness matrix, a stacked controller matrix and a trace-form Lyapunov function produce a clean stability proof: all signals remain bounded and the tracking error converges to zero. These conclusions do not imply parameter convergence without persistent excitation and identifiability. General nonsymmetric high-frequency gain matrices require additional sign or factorization information, while unmatched uncertainty changes asymptotic tracking into a bounded-error problem.