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

• formulate a coupled 2×2 plant and a stable 2×2 reference model;
• write a direct MIMO MRAC law with matrix-valued adaptive gains;
• derive the exact matching conditions for feedback and command gains;
• construct a trace-based Lyapunov function for matrix parameter errors;
• derive adaptive laws that cancel all state-parameter cross terms;
• explain why asymptotic tracking does not automatically imply convergence of every gain estimate;
• implement and numerically validate the controller in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The example uses full-state feedback and a known nonsingular input matrix. These assumptions isolate the essential multivariable adaptation mechanism. Unknown high-frequency gain matrices, output feedback, relative-degree structure, and interactor matrices require additional machinery and are not silently assumed away.

Consider the continuous-time plant \( \mathbf{x}\in\mathbb{R}^{2} \), \( \mathbf{u}\in\mathbb{R}^{2} \), and \( \mathbf{y}=\mathbf{x} \):

\[ \dot{\mathbf{x} }(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}\mathbf{u}(t), \qquad \mathbf{B}=\mathbf{I}_{2}. \]

The off-diagonal elements of \( \mathbf{A} \) create interaction: input-driven motion in one state can affect both output channels. The desired closed-loop behavior is generated by

\[ \dot{\mathbf{x} }_{m}(t) = \mathbf{A}_{m}\mathbf{x}_{m}(t) + \mathbf{B}_{m}\mathbf{r}(t), \]

where \( \mathbf{A}_{m} \) is Hurwitz and \( \mathbf{r}\in\mathbb{R}^{2} \) is the two-channel command. Define the tracking error

\[ \mathbf{e}(t)=\mathbf{x}(t)-\mathbf{x}_{m}(t). \]

Because the plant has two inputs and two measured states, the adaptive controller can use a full matrix of feedback gains rather than two independent scalar loops.

Choose the direct adaptive controller

\[ \mathbf{u}(t) = \widehat{\mathbf{K} }_{x}(t)\mathbf{x}(t) + \widehat{\mathbf{K} }_{r}(t)\mathbf{r}(t), \]

with \( \widehat{\mathbf{K} }_{x}, \widehat{\mathbf{K} }_{r}\in\mathbb{R}^{2\times 2} \). The feedback matrix contains four adaptable state-feedback coefficients, while the command matrix contains four adaptable feedforward coefficients.

Suppose constant ideal matrices \( \mathbf{K}_{x}^{*} \) and \( \mathbf{K}_{r}^{*} \) exist such that

\[ \mathbf{A} + \mathbf{B}\mathbf{K}_{x}^{*} = \mathbf{A}_{m}, \qquad \mathbf{B}\mathbf{K}_{r}^{*} = \mathbf{B}_{m}. \]

These are the plant-model matching conditions. Since \( \mathbf{B}=\mathbf{I}_{2} \), the ideal gains are

\[ \mathbf{K}_{x}^{*} = \mathbf{A}_{m}-\mathbf{A}, \qquad \mathbf{K}_{r}^{*} = \mathbf{B}_{m}. \]

Define matrix estimation errors

\[ \widetilde{\mathbf{K} }_{x} = \widehat{\mathbf{K} }_{x}-\mathbf{K}_{x}^{*}, \qquad \widetilde{\mathbf{K} }_{r} = \widehat{\mathbf{K} }_{r}-\mathbf{K}_{r}^{*}. \]

Substituting the adaptive controller into the plant and subtracting the reference-model dynamics gives

\[ \dot{\mathbf{e} } = \mathbf{A}_{m}\mathbf{e} + \widetilde{\mathbf{K} }_{x}\mathbf{x} + \widetilde{\mathbf{K} }_{r}\mathbf{r}. \]

The final two terms are multivariable parameter-error forcing terms. The adaptive laws will be selected to cancel their contribution to the derivative of a Lyapunov function.

flowchart TD
  R["Command vector r(t)"] --> RM["Reference model: xm_dot = Am xm + Bm r"]
  R --> C["Adaptive controller: u = Kx_hat x + Kr_hat r"]
  C --> P["Coupled 2x2 plant: x_dot = A x + u"]
  P --> X["Measured state x"]
  RM --> E["Tracking error e = x - xm"]
  X --> E
  E --> U["Matrix update laws"]
  X --> U
  R --> U
  U --> C
        

Unlike decentralized adaptation, every row of each gain matrix is updated from the common error vector. Therefore, an error in one output can alter coefficients multiplying both state channels and both command channels.

Select any symmetric matrix \( \mathbf{Q}=\mathbf{Q}^{T}\succ 0 \). Since \( \mathbf{A}_{m} \) is Hurwitz, the Lyapunov equation

\[ \mathbf{A}_{m}^{T}\mathbf{P} + \mathbf{P}\mathbf{A}_{m} = -\mathbf{Q} \]

has a unique symmetric solution \( \mathbf{P}=\mathbf{P}^{T}\succ 0 \). Let \( \boldsymbol{\Gamma}_{x}\succ 0 \) and \( \boldsymbol{\Gamma}_{r}\succ 0 \) be adaptation-gain matrices. Consider

\[ \begin{aligned} V &= \mathbf{e}^{T}\mathbf{P}\mathbf{e} + \operatorname{tr}\!\left( \widetilde{\mathbf{K} }_{x}^{T} \boldsymbol{\Gamma}_{x}^{-1} \widetilde{\mathbf{K} }_{x} \right)\\ &\quad+ \operatorname{tr}\!\left( \widetilde{\mathbf{K} }_{r}^{T} \boldsymbol{\Gamma}_{r}^{-1} \widetilde{\mathbf{K} }_{r} \right). \end{aligned} \]

The trace terms are squared Frobenius norms under weighted inner products. Differentiating and using the error dynamics yields

\[ \begin{aligned} \dot{V} &= -\mathbf{e}^{T}\mathbf{Q}\mathbf{e} + 2\mathbf{e}^{T}\mathbf{P} \widetilde{\mathbf{K} }_{x}\mathbf{x} + 2\mathbf{e}^{T}\mathbf{P} \widetilde{\mathbf{K} }_{r}\mathbf{r}\\ &\quad+ 2\operatorname{tr}\!\left( \widetilde{\mathbf{K} }_{x}^{T} \boldsymbol{\Gamma}_{x}^{-1} \dot{\widehat{\mathbf{K} } }_{x} \right) + 2\operatorname{tr}\!\left( \widetilde{\mathbf{K} }_{r}^{T} \boldsymbol{\Gamma}_{r}^{-1} \dot{\widehat{\mathbf{K} } }_{r} \right). \end{aligned} \]

The scalar-trace identity

\[ \mathbf{e}^{T}\mathbf{P} \widetilde{\mathbf{K} }_{x}\mathbf{x} = \operatorname{tr}\!\left( \widetilde{\mathbf{K} }_{x}^{T} \mathbf{P}\mathbf{e}\mathbf{x}^{T} \right) \]

and the analogous identity for the command term suggest the adaptive laws

\[ \boxed{ \dot{\widehat{\mathbf{K} } }_{x} = -\boldsymbol{\Gamma}_{x} \mathbf{P}\mathbf{e}\mathbf{x}^{T} }, \qquad \boxed{ \dot{\widehat{\mathbf{K} } }_{r} = -\boldsymbol{\Gamma}_{r} \mathbf{P}\mathbf{e}\mathbf{r}^{T} }. \]

Substitution gives the exact cancellation

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

Consequently, \( V(t)\leq V(0) \), so \( \mathbf{e} \), \( \widetilde{\mathbf{K} }_{x} \), and \( \widetilde{\mathbf{K} }_{r} \) are bounded. Moreover,

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

so \( \mathbf{e}\in\mathcal{L}_{2} \). For bounded commands, the closed-loop equations imply bounded \( \dot{\mathbf{e} } \). Barbalat's lemma then gives

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

This proof establishes asymptotic tracking. It does not by itself prove that every element of the two gain matrices converges to its ideal value.

Use the coupled, open-loop-unstable plant matrix

\[ \mathbf{A} = \begin{bmatrix} 0.2 & 1.0\\ -1.5 & 0.3 \end{bmatrix}, \qquad \mathbf{B}=\mathbf{I}_{2}. \]

The desired model is

\[ \mathbf{A}_{m} = \begin{bmatrix} -1.5 & 0.4\\ -0.3 & -1.2 \end{bmatrix}, \qquad \mathbf{B}_{m} = \begin{bmatrix} 1.0 & 0.2\\ 0.1 & 1.0 \end{bmatrix}. \]

The characteristic polynomial of \( \mathbf{A}_{m} \) is

\[ \det(s\mathbf{I}_{2}-\mathbf{A}_{m}) = s^{2}+2.7s+1.92, \]

whose roots have negative real parts. Therefore, \( \mathbf{A}_{m} \) is Hurwitz. The ideal gains are

\[ \mathbf{K}_{x}^{*} = \begin{bmatrix} -1.7 & -0.6\\ 1.2 & -1.5 \end{bmatrix}, \qquad \mathbf{K}_{r}^{*} = \begin{bmatrix} 1.0 & 0.2\\ 0.1 & 1.0 \end{bmatrix}. \]

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

\[ \mathbf{P} = \begin{bmatrix} 0.33275463 & 0.00289352\\ 0.00289352 & 0.41763117 \end{bmatrix}. \]

Use

\[ \boldsymbol{\Gamma}_{x}=4\mathbf{I}_{2}, \qquad \boldsymbol{\Gamma}_{r}=2\mathbf{I}_{2}, \]

and the bounded command

\[ \mathbf{r}(t) = \begin{bmatrix} 1+0.35\sin(0.7t)\\ -0.6+0.25\cos(0.5t) \end{bmatrix}. \]

The initial conditions are \( \mathbf{x}(0)=[0.8,-0.5]^{T} \), \( \mathbf{x}_{m}(0)=\mathbf{0} \), and both estimated gain matrices initially zero.

flowchart TD
  S["Initialize x, xm, Kx_hat, Kr_hat"] --> R["Evaluate two-channel command r(t)"]
  R --> E["Compute e = x - xm"]
  E --> U["Compute u = Kx_hat x + Kr_hat r"]
  U --> D["Evaluate plant and model derivatives"]
  E --> A["Evaluate matrix gain derivatives"]
  D --> I["Integrate all 12 states with RK4 or adaptive ODE solver"]
  A --> I
  I --> L["Log states, errors, and gain entries"]
  L --> C{"Final time reached?"}
  C -->|no| R
  C -->|yes| O["Plot tracking and report final error"]
        

The augmented numerical state has twelve entries: two plant states, two model states, four entries of \( \widehat{\mathbf{K} }_{x} \), and four entries of \( \widehat{\mathbf{K} }_{r} \). The fixed-step C++, Java, and Mathematica implementations use classical fourth-order Runge-Kutta integration. Python and MATLAB use adaptive ODE solvers with a restricted maximum step.

Chapter23_Lesson5.py

# Chapter23_Lesson5.py
# 2x2 full-state MIMO model-reference adaptive controller (MRAC).
# Dependencies: numpy, scipy, matplotlib
# Install: python -m pip install numpy scipy matplotlib

from __future__ import annotations

import csv
from dataclasses import dataclass
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import solve_continuous_lyapunov


@dataclass(frozen=True)
class Design:
    A: np.ndarray
    A_m: np.ndarray
    B_m: np.ndarray
    Q: np.ndarray
    P: np.ndarray
    Gamma_x: np.ndarray
    Gamma_r: np.ndarray


def build_design() -> Design:
    # Actual coupled plant: x_dot = A x + u.  B = I_2.
    A = np.array([[0.2, 1.0], [-1.5, 0.3]], dtype=float)

    # Desired stable reference model: x_m_dot = A_m x_m + B_m r.
    A_m = np.array([[-1.5, 0.4], [-0.3, -1.2]], dtype=float)
    B_m = np.array([[1.0, 0.2], [0.1, 1.0]], dtype=float)

    Q = np.eye(2)
    # SciPy solves A X + X A^T = Q.  Passing A_m^T and -Q gives
    # A_m^T P + P A_m = -Q.
    P = solve_continuous_lyapunov(A_m.T, -Q)

    Gamma_x = 4.0 * np.eye(2)
    Gamma_r = 2.0 * np.eye(2)
    return Design(A, A_m, B_m, Q, P, Gamma_x, Gamma_r)


def command(t: float) -> np.ndarray:
    """Bounded two-channel command with distinct frequencies."""
    return np.array(
        [1.0 + 0.35 * np.sin(0.7 * t), -0.6 + 0.25 * np.cos(0.5 * t)],
        dtype=float,
    )


def unpack(z: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    x = z[0:2]
    x_m = z[2:4]
    k_x = z[4:8].reshape(2, 2)
    k_r = z[8:12].reshape(2, 2)
    return x, x_m, k_x, k_r


def closed_loop_rhs(t: float, z: np.ndarray, d: Design) -> np.ndarray:
    x, x_m, k_x, k_r = unpack(z)
    r = command(t)
    e = x - x_m

    # Adaptive MIMO control: u = K_x_hat x + K_r_hat r.
    u = k_x @ x + k_r @ r

    x_dot = d.A @ x + u
    x_m_dot = d.A_m @ x_m + d.B_m @ r

    # Matrix gradient laws that cancel the Lyapunov cross terms.
    k_x_dot = -d.Gamma_x @ d.P @ np.outer(e, x)
    k_r_dot = -d.Gamma_r @ d.P @ np.outer(e, r)

    return np.concatenate((x_dot, x_m_dot, k_x_dot.ravel(), k_r_dot.ravel()))


def simulate(t_final: float = 30.0) -> tuple[Design, np.ndarray, np.ndarray]:
    d = build_design()

    z0 = np.zeros(12)
    z0[0:2] = np.array([0.8, -0.5])  # plant initial state
    z0[2:4] = np.array([0.0, 0.0])   # model initial state
    # K_x_hat(0) = 0 and K_r_hat(0) = 0.

    t_eval = np.linspace(0.0, t_final, 6001)
    sol = solve_ivp(
        fun=lambda t, z: closed_loop_rhs(t, z, d),
        t_span=(0.0, t_final),
        y0=z0,
        t_eval=t_eval,
        method="RK45",
        rtol=1.0e-8,
        atol=1.0e-10,
        max_step=0.01,
    )
    if not sol.success:
        raise RuntimeError(sol.message)

    return d, sol.t, sol.y.T


def save_results(d: Design, t: np.ndarray, z: np.ndarray) -> None:
    x = z[:, 0:2]
    x_m = z[:, 2:4]
    e = x - x_m

    csv_path = Path("Chapter23_Lesson5_results_python.csv")
    with csv_path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(
            [
                "t",
                "x1",
                "x2",
                "xm1",
                "xm2",
                "e1",
                "e2",
                "error_norm",
                "Kx11",
                "Kx12",
                "Kx21",
                "Kx22",
                "Kr11",
                "Kr12",
                "Kr21",
                "Kr22",
            ]
        )
        for ti, zi, ei in zip(t, z, e):
            writer.writerow(
                [
                    ti,
                    zi[0],
                    zi[1],
                    zi[2],
                    zi[3],
                    ei[0],
                    ei[1],
                    np.linalg.norm(ei),
                    *zi[4:12],
                ]
            )

    fig = plt.figure(figsize=(10, 7))
    ax1 = fig.add_subplot(2, 1, 1)
    ax1.plot(t, x[:, 0], label="x1")
    ax1.plot(t, x_m[:, 0], "--", label="xm1")
    ax1.plot(t, x[:, 1], label="x2")
    ax1.plot(t, x_m[:, 1], "--", label="xm2")
    ax1.set_ylabel("States")
    ax1.grid(True)
    ax1.legend(ncol=2)

    ax2 = fig.add_subplot(2, 1, 2)
    ax2.semilogy(t, np.maximum(np.linalg.norm(e, axis=1), 1.0e-12))
    ax2.set_xlabel("Time [s]")
    ax2.set_ylabel("||e||2")
    ax2.grid(True)

    fig.tight_layout()
    fig.savefig("Chapter23_Lesson5_python.png", dpi=180)
    plt.close(fig)

    k_x_star = d.A_m - d.A
    k_r_star = d.B_m
    print("P =\n", d.P)
    print("Ideal K_x* =\n", k_x_star)
    print("Ideal K_r* =\n", k_r_star)
    print("Final K_x_hat =\n", z[-1, 4:8].reshape(2, 2))
    print("Final K_r_hat =\n", z[-1, 8:12].reshape(2, 2))
    print("Final tracking error norm =", np.linalg.norm(e[-1]))
    print("Saved:", csv_path)
    print("Saved: Chapter23_Lesson5_python.png")


if __name__ == "__main__":
    design, time, trajectory = simulate()
    save_results(design, time, trajectory)

The Python version uses scipy.linalg.solve_continuous_lyapunov for \( \mathbf{P} \), scipy.integrate.solve_ivp for numerical integration, NumPy for matrix operations, and Matplotlib for plots.

Chapter23_Lesson5.cpp

// Chapter23_Lesson5.cpp
// 2x2 full-state MIMO MRAC using only the C++ standard library.
// Build: g++ -std=c++17 -O2 Chapter23_Lesson5.cpp -o Chapter23_Lesson5
// Run:   ./Chapter23_Lesson5

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

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

static Vec2 mat_vec(const Mat2& M, const Vec2& v) {
    return {
        M[0][0] * v[0] + M[0][1] * v[1],
        M[1][0] * v[0] + M[1][1] * v[1]
    };
}

static Vec2 add(const Vec2& a, const Vec2& b) {
    return {a[0] + b[0], a[1] + b[1]};
}

static Vec2 sub(const Vec2& a, const Vec2& b) {
    return {a[0] - b[0], a[1] - b[1]};
}

static Vec2 command(double t) {
    return {
        1.0 + 0.35 * std::sin(0.7 * t),
        -0.6 + 0.25 * std::cos(0.5 * t)
    };
}

static Vec12 rhs(double t, const Vec12& z) {
    const Mat2 A{ { {0.2, 1.0}, {-1.5, 0.3} } };
    const Mat2 Am{ { {-1.5, 0.4}, {-0.3, -1.2} } };
    const Mat2 Bm{ { {1.0, 0.2}, {0.1, 1.0} } };
    const Mat2 P{ { {0.3327546296, 0.0028935185},
                  {0.0028935185, 0.4176311728} } };
    const double gamma_x = 4.0;
    const double gamma_r = 2.0;

    const Vec2 x{z[0], z[1]};
    const Vec2 xm{z[2], z[3]};
    const Mat2 Kx{ { {z[4], z[5]}, {z[6], z[7]} } };
    const Mat2 Kr{ { {z[8], z[9]}, {z[10], z[11]} } };

    const Vec2 r = command(t);
    const Vec2 e = sub(x, xm);
    const Vec2 u = add(mat_vec(Kx, x), mat_vec(Kr, r));
    const Vec2 xdot = add(mat_vec(A, x), u);
    const Vec2 xmdot = add(mat_vec(Am, xm), mat_vec(Bm, r));
    const Vec2 Pe = mat_vec(P, e);

    Vec12 dz{};
    dz[0] = xdot[0];
    dz[1] = xdot[1];
    dz[2] = xmdot[0];
    dz[3] = xmdot[1];

    // Kx_dot = -gamma_x (P e) x^T, row-major.
    dz[4] = -gamma_x * Pe[0] * x[0];
    dz[5] = -gamma_x * Pe[0] * x[1];
    dz[6] = -gamma_x * Pe[1] * x[0];
    dz[7] = -gamma_x * Pe[1] * x[1];

    // Kr_dot = -gamma_r (P e) r^T, row-major.
    dz[8] = -gamma_r * Pe[0] * r[0];
    dz[9] = -gamma_r * Pe[0] * r[1];
    dz[10] = -gamma_r * Pe[1] * r[0];
    dz[11] = -gamma_r * Pe[1] * r[1];
    return dz;
}

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

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

    Vec12 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]);
        if (!std::isfinite(next[i])) {
            throw std::runtime_error("Non-finite state detected.");
        }
    }
    return next;
}

int main() {
    constexpr double dt = 0.002;
    constexpr double t_final = 30.0;
    const int steps = static_cast<int>(std::round(t_final / dt));

    Vec12 z{};
    z[0] = 0.8;
    z[1] = -0.5;

    std::ofstream csv("Chapter23_Lesson5_results_cpp.csv");
    if (!csv) {
        std::cerr << "Cannot create output CSV.\n";
        return 1;
    }

    csv << "t,x1,x2,xm1,xm2,e1,e2,error_norm,"
           "Kx11,Kx12,Kx21,Kx22,Kr11,Kr12,Kr21,Kr22\n";
    csv << std::setprecision(12);

    for (int k = 0; k <= steps; ++k) {
        const double t = k * dt;
        const double e1 = z[0] - z[2];
        const double e2 = z[1] - z[3];
        csv << t
            << ',' << z[0] << ',' << z[1]
            << ',' << z[2] << ',' << z[3]
            << ',' << e1 << ',' << e2
            << ',' << std::sqrt(e1 * e1 + e2 * e2);
        for (std::size_t i = 4; i < z.size(); ++i) {
            csv << ',' << z[i];
        }
        csv << '\n';

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

    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 Kx_hat = [[" << z[4] << ", " << z[5] << "], ["
              << z[6] << ", " << z[7] << "]]\n";
    std::cout << "Final Kr_hat = [[" << z[8] << ", " << z[9] << "], ["
              << z[10] << ", " << z[11] << "]]\n";
    std::cout << "Saved: Chapter23_Lesson5_results_cpp.csv\n";
    return 0;
}

This implementation uses fixed-size arrays and a from-scratch RK4 integrator, so no external numerical library is required.

Chapter23_Lesson5.java

// Chapter23_Lesson5.java
// 2x2 full-state MIMO MRAC using only the Java standard library.
// Build: javac Chapter23_Lesson5.java
// Run:   java Chapter23_Lesson5

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

public final class Chapter23_Lesson5 {
    private static final double[][] A = {
        {0.2, 1.0},
        {-1.5, 0.3}
    };

    private static final double[][] AM = {
        {-1.5, 0.4},
        {-0.3, -1.2}
    };

    private static final double[][] BM = {
        {1.0, 0.2},
        {0.1, 1.0}
    };

    private static final double[][] P = {
        {0.3327546296, 0.0028935185},
        {0.0028935185, 0.4176311728}
    };

    private static final double GAMMA_X = 4.0;
    private static final double GAMMA_R = 2.0;

    private Chapter23_Lesson5() {
    }

    private static double[] matVec(double[][] matrix, double[] vector) {
        return new double[] {
            matrix[0][0] * vector[0] + matrix[0][1] * vector[1],
            matrix[1][0] * vector[0] + matrix[1][1] * vector[1]
        };
    }

    private static double[] add(double[] a, double[] b) {
        return new double[] {a[0] + b[0], a[1] + b[1]};
    }

    private static double[] command(double t) {
        return new double[] {
            1.0 + 0.35 * Math.sin(0.7 * t),
            -0.6 + 0.25 * Math.cos(0.5 * t)
        };
    }

    private static double[] rhs(double t, double[] z) {
        double[] x = {z[0], z[1]};
        double[] xm = {z[2], z[3]};
        double[][] kx = {
            {z[4], z[5]},
            {z[6], z[7]}
        };
        double[][] kr = {
            {z[8], z[9]},
            {z[10], z[11]}
        };

        double[] r = command(t);
        double[] e = {x[0] - xm[0], x[1] - xm[1]};
        double[] u = add(matVec(kx, x), matVec(kr, r));
        double[] xDot = add(matVec(A, x), u);
        double[] xmDot = add(matVec(AM, xm), matVec(BM, r));
        double[] pe = matVec(P, e);

        double[] dz = new double[12];
        dz[0] = xDot[0];
        dz[1] = xDot[1];
        dz[2] = xmDot[0];
        dz[3] = xmDot[1];

        dz[4] = -GAMMA_X * pe[0] * x[0];
        dz[5] = -GAMMA_X * pe[0] * x[1];
        dz[6] = -GAMMA_X * pe[1] * x[0];
        dz[7] = -GAMMA_X * pe[1] * x[1];

        dz[8] = -GAMMA_R * pe[0] * r[0];
        dz[9] = -GAMMA_R * pe[0] * r[1];
        dz[10] = -GAMMA_R * pe[1] * r[0];
        dz[11] = -GAMMA_R * pe[1] * r[1];
        return dz;
    }

    private static double[] scaledAdd(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 dt) {
        double[] k1 = rhs(t, z);
        double[] k2 = rhs(t + 0.5 * dt, scaledAdd(z, k1, 0.5 * dt));
        double[] k3 = rhs(t + 0.5 * dt, scaledAdd(z, k2, 0.5 * dt));
        double[] k4 = rhs(t + dt, scaledAdd(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]);
            if (!Double.isFinite(next[i])) {
                throw new IllegalStateException("Non-finite state detected.");
            }
        }
        return next;
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);

        final double dt = 0.002;
        final double tFinal = 30.0;
        final int steps = (int) Math.round(tFinal / dt);

        double[] z = new double[12];
        z[0] = 0.8;
        z[1] = -0.5;

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter23_Lesson5_results_java.csv"))) {
            writer.write(
                "t,x1,x2,xm1,xm2,e1,e2,error_norm,"
                + "Kx11,Kx12,Kx21,Kx22,Kr11,Kr12,Kr21,Kr22\n"
            );

            for (int k = 0; k <= steps; ++k) {
                double t = k * dt;
                double e1 = z[0] - z[2];
                double e2 = z[1] - z[3];

                writer.write(String.format(
                    Locale.US,
                    "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g",
                    t,
                    z[0],
                    z[1],
                    z[2],
                    z[3],
                    e1,
                    e2,
                    Math.hypot(e1, e2)
                ));
                for (int i = 4; i < z.length; ++i) {
                    writer.write(String.format(Locale.US, ",%.12g", z[i]));
                }
                writer.newLine();

                if (k < steps) {
                    z = rk4Step(t, z, dt);
                }
            }
        }

        double finalError = Math.hypot(z[0] - z[2], z[1] - z[3]);
        System.out.printf(Locale.US, "Final tracking error norm = %.8g%n", finalError);
        System.out.printf(
            Locale.US,
            "Final Kx_hat = [[%.6f, %.6f], [%.6f, %.6f]]%n",
            z[4], z[5], z[6], z[7]
        );
        System.out.printf(
            Locale.US,
            "Final Kr_hat = [[%.6f, %.6f], [%.6f, %.6f]]%n",
            z[8], z[9], z[10], z[11]
        );
        System.out.println("Saved: Chapter23_Lesson5_results_java.csv");
    }
}

The Java implementation mirrors the C++ state ordering and produces a CSV file with the same state, error, and gain columns.

Chapter23_Lesson5.m

% Chapter23_Lesson5.m
% 2x2 full-state MIMO model-reference adaptive controller (MRAC).
% Requires MATLAB. The lyap call requires Control System Toolbox.
% Run from the MATLAB command window:
%   Chapter23_Lesson5

clear; clc; close all;

A = [0.2, 1.0; -1.5, 0.3];
Am = [-1.5, 0.4; -0.3, -1.2];
Bm = [1.0, 0.2; 0.1, 1.0];

Q = eye(2);
P = lyap(Am', Q);  % Am' P + P Am = -Q.

GammaX = 4.0 * eye(2);
GammaR = 2.0 * eye(2);

z0 = zeros(12, 1);
z0(1:2) = [0.8; -0.5];

options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10, 'MaxStep', 0.01);
[t, z] = ode45(@(t, z) rhs(t, z, A, Am, Bm, P, GammaX, GammaR), ...
               [0, 30], z0, options);

x = z(:, 1:2);
xm = z(:, 3:4);
e = x - xm;
errorNorm = vecnorm(e, 2, 2);

figure('Name', 'Chapter 23 Lesson 5: MIMO MRAC');
subplot(2, 1, 1);
plot(t, x(:, 1), 'LineWidth', 1.1); hold on;
plot(t, xm(:, 1), '--', 'LineWidth', 1.1);
plot(t, x(:, 2), 'LineWidth', 1.1);
plot(t, xm(:, 2), '--', 'LineWidth', 1.1);
grid on;
ylabel('States');
legend('x_1', 'x_{m1}', 'x_2', 'x_{m2}', 'Location', 'best');

subplot(2, 1, 2);
semilogy(t, max(errorNorm, 1e-12), 'LineWidth', 1.1);
grid on;
xlabel('Time [s]');
ylabel('||e||_2');

exportgraphics(gcf, 'Chapter23_Lesson5_matlab.png', 'Resolution', 180);

results = array2table([t, z, e, errorNorm], 'VariableNames', ...
    {'t', 'x1', 'x2', 'xm1', 'xm2', ...
     'Kx11', 'Kx21', 'Kx12', 'Kx22', ...
     'Kr11', 'Kr21', 'Kr12', 'Kr22', ...
     'e1', 'e2', 'errorNorm'});
writetable(results, 'Chapter23_Lesson5_results_matlab.csv');

KxFinal = reshape(z(end, 5:8), 2, 2);
KrFinal = reshape(z(end, 9:12), 2, 2);
KxStar = Am - A;
KrStar = Bm;

disp('P ='); disp(P);
disp('Ideal Kx* ='); disp(KxStar);
disp('Ideal Kr* ='); disp(KrStar);
disp('Final Kx_hat ='); disp(KxFinal);
disp('Final Kr_hat ='); disp(KrFinal);
fprintf('Final tracking error norm = %.8g\n', errorNorm(end));

% Simulink mapping:
% 1. Use two State-Space blocks for (A,I,C=I,D=0) and (Am,Bm,C=I,D=0).
% 2. Compute e = x - xm.
% 3. Form -GammaX*P*e*x' and -GammaR*P*e*r' with Matrix Multiply blocks.
% 4. Integrate the four entries of each gain matrix.
% 5. Form u = KxHat*x + KrHat*r and feed it to the plant.

function dz = rhs(t, z, A, Am, Bm, P, GammaX, GammaR)
    x = z(1:2);
    xm = z(3:4);
    Kx = reshape(z(5:8), 2, 2);
    Kr = reshape(z(9:12), 2, 2);

    r = [1.0 + 0.35*sin(0.7*t); ...
        -0.6 + 0.25*cos(0.5*t)];
    e = x - xm;

    u = Kx*x + Kr*r;
    xDot = A*x + u;
    xmDot = Am*xm + Bm*r;

    KxDot = -GammaX*P*e*x';
    KrDot = -GammaR*P*e*r';

    dz = [xDot; xmDot; KxDot(:); KrDot(:)];
end

In Simulink, use separate State-Space blocks for the plant and reference model, Sum blocks for the tracking error, Matrix Multiply blocks for the outer products, and Integrator blocks for the eight adaptive parameters. Preserve matrix dimensions explicitly: the error and command are 2×1 vectors, while each adaptive gain is 2×2.

Chapter23_Lesson5.nb

Notebook[{
Cell["Chapter 23, Lesson 5: 2x2 Adaptive MIMO Controller", "Title"],
Cell["Chapter23_Lesson5.nb", "Subtitle"],
Cell[BoxData[
"ClearAll[\"Global`*\"];\n\
A = { {0.2, 1.0}, {-1.5, 0.3} };\n\
Am = { {-1.5, 0.4}, {-0.3, -1.2} };\n\
Bm = { {1.0, 0.2}, {0.1, 1.0} };\n\
P = { {0.3327546296, 0.0028935185}, {0.0028935185, 0.4176311728} };\n\
gammaX = 4.0;\n\
gammaR = 2.0;\n\
dt = 0.002;\n\
tFinal = 30.0;\n\
nSteps = Round[tFinal/dt];\n\
command[t_] := {1.0 + 0.35 Sin[0.7 t], -0.6 + 0.25 Cos[0.5 t]};\n\
rhs[t_, z_List] := Module[{x, xm, kx, kr, r, e, u, xDot, xmDot, pe, kxDot, krDot},\n\
  x = z[[1 ;; 2]];\n\
  xm = z[[3 ;; 4]];\n\
  kx = Partition[z[[5 ;; 8]], 2];\n\
  kr = Partition[z[[9 ;; 12]], 2];\n\
  r = command[t];\n\
  e = x - xm;\n\
  u = kx.x + kr.r;\n\
  xDot = A.x + u;\n\
  xmDot = Am.xm + Bm.r;\n\
  pe = P.e;\n\
  kxDot = -gammaX Outer[Times, pe, x];\n\
  krDot = -gammaR Outer[Times, pe, r];\n\
  Join[xDot, xmDot, Flatten[kxDot], Flatten[krDot]]\n\
];\n\
rk4Step[t_, z_List] := Module[{k1, k2, k3, k4},\n\
  k1 = rhs[t, z];\n\
  k2 = rhs[t + dt/2, z + (dt/2) k1];\n\
  k3 = rhs[t + dt/2, z + (dt/2) k2];\n\
  k4 = rhs[t + dt, z + dt k3];\n\
  z + (dt/6) (k1 + 2 k2 + 2 k3 + k4)\n\
];\n\
times = N@Range[0, tFinal, dt];\n\
states = ConstantArray[0.0, {Length[times], 12}];\n\
states[[1, 1 ;; 2]] = {0.8, -0.5};\n\
Do[states[[k + 1]] = rk4Step[times[[k]], states[[k]]], {k, 1, nSteps}];\n\
errors = states[[All, 1 ;; 2]] - states[[All, 3 ;; 4]];\n\
errorNorm = Norm /@ errors;\n\
statePlot = ListLinePlot[\n\
  {\n\
    Transpose[{times, states[[All, 1]]}],\n\
    Transpose[{times, states[[All, 3]]}],\n\
    Transpose[{times, states[[All, 2]]}],\n\
    Transpose[{times, states[[All, 4]]]}\n\
  },\n\
  PlotLegends -> {\"x1\", \"xm1\", \"x2\", \"xm2\"},\n\
  Frame -> True,\n\
  FrameLabel -> {\"Time [s]\", \"States\"},\n\
  PlotRange -> All,\n\
  ImageSize -> Large\n\
];\n\
errorPlot = ListLogPlot[\n\
  Transpose[{times, Clip[errorNorm, {10^-12, Infinity}]}],\n\
  Frame -> True,\n\
  FrameLabel -> {\"Time [s]\", \"||e||2\"},\n\
  PlotRange -> All,\n\
  ImageSize -> Large\n\
];\n\
Print[statePlot];\n\
Print[errorPlot];\n\
Export[\"Chapter23_Lesson5_mathematica.png\", Column[{statePlot, errorPlot}]];\n\
header = {\"t\", \"x1\", \"x2\", \"xm1\", \"xm2\", \"Kx11\", \"Kx12\", \"Kx21\", \"Kx22\", \"Kr11\", \"Kr12\", \"Kr21\", \"Kr22\", \"e1\", \"e2\", \"errorNorm\"};\n\
rows = MapThread[Join[{ #1 }, #2, #3, { #4 }] &, {times, states, errors, errorNorm}];\n\
Export[\"Chapter23_Lesson5_results_mathematica.csv\", Prepend[rows, header]];\n\
Print[\"Ideal Kx* = \", MatrixForm[Am - A]];\n\
Print[\"Ideal Kr* = \", MatrixForm[Bm]];\n\
Print[\"Final KxHat = \", MatrixForm[Partition[states[[-1, 5 ;; 8]], 2]]];\n\
Print[\"Final KrHat = \", MatrixForm[Partition[states[[-1, 9 ;; 12]], 2]]];\n\
Print[\"Final tracking error norm = \", Last[errorNorm]];"
], "Input"]
},
WindowSize -> {1200, 850},
StyleDefinitions -> "Default.nb"
]

The notebook implements the same fixed-step RK4 method, exports tracking data, and plots both state trajectories and the logarithmic error norm.

All five implementations represent the same closed-loop differential equations. With the stated initial conditions and gains, the tested Python, C++, and Java implementations agree to numerical precision and produce a tracking-error norm of approximately \( 4.47\times 10^{-2} \) at \( t=30\,\text{s} \).

The error need not decrease monotonically. The Lyapunov result states that the total energy-like function is nonincreasing, while individual state-error components may exchange energy through the coupled dynamics.

The estimated matrices at a finite time generally differ from \( \mathbf{K}_{x}^{*} \) and \( \mathbf{K}_{r}^{*} \). This does not contradict the tracking theorem. Parameter convergence requires sufficient excitation of the combined regressor

\[ \boldsymbol{\omega}(t) = \begin{bmatrix} \mathbf{x}(t)\\ \mathbf{r}(t) \end{bmatrix}. \]

A standard persistent-excitation condition is the existence of constants \( T>0 \) and \( \alpha>0 \) such that

\[ \int_{t}^{t+T} \boldsymbol{\omega}(\tau) \boldsymbol{\omega}^{T}(\tau)\,d\tau \succeq \alpha\mathbf{I} \quad \text{for all }t\geq 0. \]

If the command and resulting state trajectory fail to excite every independent parameter direction, multiple gain matrices can generate nearly identical input signals along the observed trajectory.

The adaptive controller is not merely two scalar MRAC loops. For example,

\[ u_{1} = \widehat{k}_{x,11}x_{1} + \widehat{k}_{x,12}x_{2} + \widehat{k}_{r,11}r_{1} + \widehat{k}_{r,12}r_{2}. \]

Thus the first actuator can compensate the influence of the second state and second command. The off-diagonal coefficients perform adaptive interaction compensation. Imposing diagonal gain matrices would reduce the parameter count, but exact model matching could then fail whenever the plant or desired model requires cross-channel action.

A structural restriction is acceptable only when it is consistent with the matching equations. In general, the question is not whether a decentralized controller is simpler, but whether the chosen restricted parameterization contains an ideal controller.

The nominal proof assumes exact matching, full state measurement, bounded commands, exact continuous-time adaptation, and no actuator constraints. Practical implementations should consider:

• Projection: constrain gain matrices to a known admissible set and prevent numerical drift.
• Normalization: divide updates by a term such as \( 1+\|\mathbf{x}\|^{2}+\|\mathbf{r}\|^{2} \) to limit update magnitude under large regressors.
• Leakage or sigma modification: add a small restoring term to improve robustness in the presence of noise and unmodeled dynamics.
• Actuator saturation logic: stop or modify adaptation when the commanded input cannot be realized.
• Sampled-data verification: check that the integration and control sampling rates are fast relative to both the reference model and adaptation dynamics.

These modifications alter the exact derivative \( \dot{V}=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e} \). Their stability arguments normally establish uniform ultimate boundedness or convergence to a residual set rather than the same ideal asymptotic result.

Problem 1 — Verify the Lyapunov Matrix: For the numerical \( \mathbf{A}_{m} \) and \( \mathbf{Q}=\mathbf{I}_{2} \), verify that the stated matrix \( \mathbf{P} \) satisfies the Lyapunov equation.

Solution: Direct multiplication gives

\[ \mathbf{A}_{m}^{T}\mathbf{P} + \mathbf{P}\mathbf{A}_{m} \approx \begin{bmatrix} -1 & 0\\ 0 & -1 \end{bmatrix} = -\mathbf{I}_{2}, \]

up to rounding in the displayed entries of \( \mathbf{P} \). Its leading principal minors are positive, so \( \mathbf{P}\succ 0 \).

Problem 2 — Derive One Gain-Entry Update: Expand the differential equation for \( \widehat{k}_{x,12} \) when \( \boldsymbol{\Gamma}_{x}=\gamma_{x}\mathbf{I}_{2} \).

Solution: Since

\[ \dot{\widehat{\mathbf{K} } }_{x} = -\gamma_{x}(\mathbf{P}\mathbf{e})\mathbf{x}^{T}, \]

the first-row, second-column element is

\[ \dot{\widehat{k} }_{x,12} = -\gamma_{x} \left(p_{11}e_{1}+p_{12}e_{2}\right)x_{2}. \]

Hence the update depends on a weighted combination of both output errors and on the second state regressor.

Problem 3 — Restricted Diagonal Feedback: Suppose the controller is forced to use diagonal \( \mathbf{K}_{x} \). Can it exactly satisfy the numerical matching equation?

Solution: No. Exact matching requires

\[ \mathbf{K}_{x}^{*} = \begin{bmatrix} -1.7 & -0.6\\ 1.2 & -1.5 \end{bmatrix}, \]

whose off-diagonal entries are nonzero. A diagonal parameterization excludes the ideal controller, leaving an irreducible model-matching error.

Problem 4 — Why Tracking Can Converge Without Gain Convergence: Explain how \( \mathbf{e}(t)\rightarrow\mathbf{0} \) can occur while the estimated gains do not converge to their ideal matrices.

Solution: The parameter errors enter through

\[ \widetilde{\mathbf{K} }_{x}\mathbf{x} + \widetilde{\mathbf{K} }_{r}\mathbf{r}. \]

Along a non-persistently-exciting trajectory, nonzero parameter-error matrices can lie in directions that annihilate the observed regressor. The input generated along that trajectory can therefore be correct even though the individual matrix entries are not uniquely identified.

Problem 5 — General Known Input Matrix: Let the plant be \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+ \mathbf{B}\mathbf{u} \), where \( \mathbf{B} \) is known and nonsingular. State the matching equations and modify the adaptation laws.

Solution: The matching equations remain

\[ \mathbf{A}+\mathbf{B}\mathbf{K}_{x}^{*} = \mathbf{A}_{m}, \qquad \mathbf{B}\mathbf{K}_{r}^{*} = \mathbf{B}_{m}. \]

The error equation contains

\[ \mathbf{B}\widetilde{\mathbf{K} }_{x}\mathbf{x} + \mathbf{B}\widetilde{\mathbf{K} }_{r}\mathbf{r}. \]

A compatible pair of update laws is therefore

\[ \dot{\widehat{\mathbf{K} } }_{x} = -\boldsymbol{\Gamma}_{x} \mathbf{B}^{T}\mathbf{P}\mathbf{e}\mathbf{x}^{T}, \qquad \dot{\widehat{\mathbf{K} } }_{r} = -\boldsymbol{\Gamma}_{r} \mathbf{B}^{T}\mathbf{P}\mathbf{e}\mathbf{r}^{T}. \]

The factor \( \mathbf{B}^{T} \) appears because the Lyapunov cross term is \( \mathbf{e}^{T}\mathbf{P}\mathbf{B} \widetilde{\mathbf{K} }\boldsymbol{\phi} \).

A 2×2 adaptive MIMO controller requires matrix-valued feedback and command gains. Under exact plant-model matching and a known input direction, the error dynamics are linear in the matrix parameter errors. A Lyapunov function combining the tracking energy with weighted Frobenius norms leads directly to outer-product adaptive laws. These laws guarantee bounded closed-loop signals and asymptotic state tracking for bounded commands. The worked example also demonstrates a central adaptive-control distinction: successful tracking is weaker than identification of all ideal controller parameters.