Consider a continuous-time linear plant with \( n \) states, \( m \) control inputs, and \( q \) measured outputs:

\[ \dot{\mathbf{x} }(t)=\mathbf{A}(\boldsymbol{\theta})\mathbf{x}(t) +\mathbf{B}(\boldsymbol{\theta})\mathbf{u}(t),\qquad \mathbf{y}(t)=\mathbf{C}(\boldsymbol{\theta})\mathbf{x}(t), \]

where \( \mathbf{x}\in\mathbb{R}^{n} \), \( \mathbf{u}\in\mathbb{R}^{m} \), \( \mathbf{y}\in\mathbb{R}^{q} \), and \( \boldsymbol{\theta}\in\mathbb{R}^{p} \) is an unknown constant parameter vector. In a SISO model, an uncertain gain changes one scalar channel. In a MIMO model, an uncertain matrix can change gain, coupling, control direction, relative influence among channels, and even the subspace in which the uncertainty acts.

It is useful to separate a known nominal model from an uncertain perturbation:

\[ \mathbf{A}(\boldsymbol{\theta})=\mathbf{A}_{0} +\Delta\mathbf{A}(\boldsymbol{\theta}),\qquad \mathbf{B}(\boldsymbol{\theta})=\mathbf{B}_{0} +\Delta\mathbf{B}(\boldsymbol{\theta}). \]

The adaptive-control problem is structurally tractable when the unknown matrices are generated by a finite set of known basis matrices and unknown scalar coefficients. This is the MIMO analogue of the linear-in-the-parameters models introduced earlier in the course.

flowchart TD
  U["m input channels u"] --> P["MIMO plant with uncertain A, B, C"]
  P --> X["n coupled states x"]
  X --> Y["q measured outputs y"]
  X --> R["Matrix regressor Phi(x,u)"]
  U --> R
  R --> E["Parameter estimator theta_hat"]
  E --> M["Estimated uncertainty matrices"]
  M --> PRED["Adaptive predictor"]
  X --> PRED
  PRED --> ERR["Prediction error"]
  ERR --> E
        

An affine parameterization writes the uncertain matrices as

\[ \mathbf{A}(\boldsymbol{\theta}) =\mathbf{A}_{0}+\sum_{i=1}^{p_A}\theta_i\mathbf{A}_{i},\qquad \mathbf{B}(\boldsymbol{\theta}) =\mathbf{B}_{0}+\sum_{j=1}^{p_B}\theta_{p_A+j}\mathbf{B}_{j}, \]

where every basis matrix \( \mathbf{A}_{i}\in\mathbb{R}^{n\times n} \) and \( \mathbf{B}_{j}\in\mathbb{R}^{n\times m} \) is known. The total number of unknown coefficients is \( p=p_A+p_B \).

Substitution into the state equation gives

\[ \dot{\mathbf{x} } =\mathbf{A}_{0}\mathbf{x}+\mathbf{B}_{0}\mathbf{u} +\sum_{i=1}^{p_A}\theta_i\mathbf{A}_{i}\mathbf{x} +\sum_{j=1}^{p_B}\theta_{p_A+j}\mathbf{B}_{j}\mathbf{u}. \]

Define the nominal residual \( \mathbf{z} \) and matrix regressor \( \boldsymbol{\Phi} \) by

\[ \mathbf{z} :=\dot{\mathbf{x} }-\mathbf{A}_{0}\mathbf{x}-\mathbf{B}_{0}\mathbf{u}, \]

\[ \boldsymbol{\Phi}(\mathbf{x},\mathbf{u}) := \begin{bmatrix} \mathbf{A}_{1}\mathbf{x} & \cdots & \mathbf{A}_{p_A}\mathbf{x} & \mathbf{B}_{1}\mathbf{u} & \cdots & \mathbf{B}_{p_B}\mathbf{u} \end{bmatrix}\in\mathbb{R}^{n\times p}. \]

The uncertain dynamics then become the linear regression

\[ \boxed{\mathbf{z}(t)=\boldsymbol{\Phi}(t)\boldsymbol{\theta} }. \]

The regression output is now a vector, not a scalar. Each sample contributes up to \( n \) scalar equations, but those equations need not be independent. Coupling and trajectory geometry determine the actual information content.

2.1 Structural economy

Treating every entry of \( \Delta\mathbf{A} \) and \( \Delta\mathbf{B} \) as independent would require

\[ p_{\mathrm{full} }=n^{2}+nm \]

parameters. Physical structure often gives a much smaller \( p \). A common mass, inertia, aerodynamic derivative, stiffness, or actuator efficiency may affect several entries simultaneously. Using one coefficient for all entries generated by that physical parameter preserves consistency and reduces excitation requirements.

2.2 Shared parameters across matrices

The same physical coefficient may appear in both \( \mathbf{A} \) and \( \mathbf{B} \). The most general affine form is

\[ \mathbf{A}(\boldsymbol{\theta}) =\mathbf{A}_{0}+\sum_{k=1}^{p}\theta_k\mathbf{A}_{k},\qquad \mathbf{B}(\boldsymbol{\theta}) =\mathbf{B}_{0}+\sum_{k=1}^{p}\theta_k\mathbf{B}_{k}, \]

with regressor columns \( \mathbf{A}_{k}\mathbf{x}+\mathbf{B}_{k}\mathbf{u} \). This avoids duplicating a single physical unknown as two unrelated estimator states.

When an entire uncertain matrix is to be estimated, vectorization gives a systematic regression. For any \( \Delta\mathbf{A}\in\mathbb{R}^{n\times n} \),

\[ \Delta\mathbf{A}\mathbf{x} =(\mathbf{x}^{T}\otimes\mathbf{I}_{n}) \operatorname{vec}(\Delta\mathbf{A}). \]

Similarly,

\[ \Delta\mathbf{B}\mathbf{u} =(\mathbf{u}^{T}\otimes\mathbf{I}_{n}) \operatorname{vec}(\Delta\mathbf{B}). \]

Therefore, with

\[ \boldsymbol{\vartheta} := \begin{bmatrix} \operatorname{vec}(\Delta\mathbf{A})\\ \operatorname{vec}(\Delta\mathbf{B}) \end{bmatrix}, \]

the residual dynamics can be written as

\[ \mathbf{z} = \underbrace{ \begin{bmatrix} \mathbf{x}^{T}\otimes\mathbf{I}_{n} & \mathbf{u}^{T}\otimes\mathbf{I}_{n} \end{bmatrix} }_{\boldsymbol{\Omega}(\mathbf{x},\mathbf{u})} \boldsymbol{\vartheta}. \]

This construction is algebraically complete but often over-parameterized. A structured basis can be represented by a known selection matrix \( \mathbf{S} \) such that \( \boldsymbol{\vartheta}=\mathbf{S}\boldsymbol{\theta} \). Then

\[ \mathbf{z}=\boldsymbol{\Omega}\mathbf{S}\boldsymbol{\theta} =\boldsymbol{\Phi}\boldsymbol{\theta}. \]

This relationship explains why basis-matrix parameterization and vectorized parameterization are not competing models. The former is the physically structured reduction of the latter.

Proof of the vectorization identity

Let \( \Delta\mathbf{A}=[\mathbf{a}_{1}\ \cdots\ \mathbf{a}_{n}] \) and \( \mathbf{x}=[x_{1},\ldots,x_{n}]^{T} \). Then

\[ \Delta\mathbf{A}\mathbf{x} =\sum_{k=1}^{n}x_k\mathbf{a}_{k}. \]

Because column-wise vectorization gives \( \operatorname{vec}(\Delta\mathbf{A}) =[\mathbf{a}_{1}^{T},\ldots,\mathbf{a}_{n}^{T}]^{T} \), the block row \( \mathbf{x}^{T}\otimes\mathbf{I}_{n} \) forms exactly the weighted sum above. Hence \( \Delta\mathbf{A}\mathbf{x} =(\mathbf{x}^{T}\otimes\mathbf{I}_{n}) \operatorname{vec}(\Delta\mathbf{A}) \).

A central geometric question is whether the uncertain state dynamics act through directions available to the control input. Assume \( \mathbf{B}_{0} \) has full column rank. The orthogonal projector onto its column space is

\[ \mathbf{P}_{B} =\mathbf{B}_{0} (\mathbf{B}_{0}^{T}\mathbf{B}_{0})^{-1} \mathbf{B}_{0}^{T}. \]

Any state-matrix perturbation admits the decomposition

\[ \Delta\mathbf{A} =\underbrace{\mathbf{P}_{B}\Delta\mathbf{A} }_ {\Delta\mathbf{A}_{\mathrm{m} } } +\underbrace{(\mathbf{I}_{n}-\mathbf{P}_{B})\Delta\mathbf{A} }_ {\Delta\mathbf{A}_{\mathrm{u} } }. \]

The first term is matched because there exists a matrix \( \boldsymbol{\Lambda}_{A} \) satisfying

\[ \Delta\mathbf{A}_{\mathrm{m} } =\mathbf{B}_{0}\boldsymbol{\Lambda}_{A}. \]

The second term is unmatched because it lies outside \( \operatorname{im}(\mathbf{B}_{0}) \). Exact cancellation by an input of the form \( \mathbf{u}_{a}=-\boldsymbol{\Lambda}_{A}\mathbf{x} \) is possible only for the matched part.

\[ \Delta\mathbf{A}\ \text{is fully matched} \quad\Longleftrightarrow\quad (\mathbf{I}_{n}-\mathbf{P}_{B})\Delta\mathbf{A}=\mathbf{0}. \]

flowchart TD
  DA["State uncertainty DeltaA"] --> PROJ["Project onto input subspace col(B0)"]
  PROJ --> MATCH["Matched part PB DeltaA"]
  PROJ --> UNMATCH["Unmatched part (I-PB) DeltaA"]
  MATCH --> CAN["Can be represented as B0 LambdaA"]
  CAN --> ADAPT["Eligible for direct adaptive compensation"]
  UNMATCH --> CANNOT["Cannot be exactly canceled \nby available inputs"]
  CANNOT --> ROBUST["Requires stability margin, redesign, \nor robust augmentation"]
        

4.1 Why this distinction is more severe in MIMO systems

For a scalar input, the input direction is one-dimensional. For several inputs, the controller acts in the column space of \( \mathbf{B}_{0} \), whose orientation matters. A perturbation may be small in norm yet largely unmatched, or large in norm yet completely matched. Norm bounds alone do not reveal this geometry.

4.2 Input-matrix uncertainty

If \( \mathbf{B}=\mathbf{B}_{0}+\Delta\mathbf{B} \), the uncertainty changes the map from commanded inputs to state directions. A useful structured form is

\[ \mathbf{B}(\boldsymbol{\theta}) =\mathbf{B}_{0}\boldsymbol{\Lambda}(\boldsymbol{\theta}), \]

where \( \boldsymbol{\Lambda} \) is an unknown nonsingular input-effectiveness matrix. Diagonal \( \boldsymbol{\Lambda} \) models independent actuator effectiveness; a full matrix also models cross-actuation. Later MIMO MRAC designs require additional assumptions on the signs, principal minors, or factorization of this matrix.

Off-diagonal entries of \( \mathbf{A} \), \( \mathbf{B} \), and \( \mathbf{C} \) encode interaction among channels. Parameter uncertainty may therefore modify both individual loop gains and cross-coupling.

For a square plant whose first output derivative depends directly on the input, a local high-frequency gain matrix is often represented by

\[ \mathbf{K}_{p}=\mathbf{C}\mathbf{B}. \]

In SISO adaptive control, knowledge of the sign of the high-frequency gain is often enough. In MIMO control, there is no single sign for a matrix. Useful prior information may instead concern:

• nonsingularity of \( \mathbf{K}_{p} \);
• signs of leading principal minors;
• a known sign-definite symmetric part;
• an LDU or LDS factorization with known diagonal signs;
• known actuator direction with uncertain effectiveness.

These assumptions are not merely technical conveniences. They prevent the adaptive law from approaching singular controller matrices or reversing an unknown multivariable control direction. Full controller parameterization is postponed to the next lesson; here the important point is that uncertainty in a matrix direction needs more structural information than uncertainty in a scalar gain.

5.1 A norm bound does not imply a safe direction

Suppose \( \|\Delta\mathbf{B}\|_{2}\le\delta \). This bounds the perturbation magnitude but does not guarantee that \( \mathbf{B}_{0}+\Delta\mathbf{B} \) remains nonsingular. A sufficient condition is

\[ \|\mathbf{B}_{0}^{-1}\Delta\mathbf{B}\|_{2}<1, \]

because \( \mathbf{B} =\mathbf{B}_{0}(\mathbf{I}+\mathbf{B}_{0}^{-1}\Delta\mathbf{B}) \) and the Neumann series ensures invertibility of the second factor.

Given \( \mathbf{z}(t)=\boldsymbol{\Phi}(t)\boldsymbol{\theta} \), define the finite-window information matrix

\[ \mathbf{W}(t_{0},T) :=\int_{t_{0} }^{t_{0}+T} \boldsymbol{\Phi}^{T}(\tau)\boldsymbol{\Phi}(\tau)\,d\tau. \]

The parameter vector is uniquely identifiable from noiseless data on that window if

\[ \operatorname{rank}\mathbf{W}(t_{0},T)=p. \]

Persistent excitation strengthens this to a uniform condition: there exist constants \( T>0 \) and \( \alpha>0 \) such that

\[ \int_{t}^{t+T} \boldsymbol{\Phi}^{T}(\tau)\boldsymbol{\Phi}(\tau)\,d\tau \succeq \alpha\mathbf{I}_{p} \qquad\text{for every }t\ge 0. \]

6.1 Multichannel excitation is not the same as nonzero input

A MIMO input can be nonzero while still producing a rank-deficient regressor. For example, if \( u_{2}(t)=c\,u_{1}(t) \) and the two uncertain input-basis columns depend on those channels in the same state direction, their regressor columns may remain proportional. The estimator can then recover only a linear combination of parameters.

6.2 State and input trajectories jointly determine excitation

Since columns of \( \boldsymbol{\Phi} \) contain \( \mathbf{A}_{i}\mathbf{x} \) and \( \mathbf{B}_{j}\mathbf{u} \), informative input signals may still fail if the closed-loop state remains in an invariant subspace that annihilates some basis matrices. Thus PE is a property of the complete closed-loop trajectory, not of the command signal alone.

6.3 State convergence versus parameter convergence

Adaptive prediction or tracking error can converge to zero even when \( \boldsymbol{\theta} \) is not uniquely identified. If \( \boldsymbol{\Phi}(t)\tilde{\boldsymbol{\theta} }(t) \rightarrow\mathbf{0} \) but \( \boldsymbol{\Phi} \) lacks PE, a nonzero parameter error may lie in the unexcited nullspace. This distinction is especially important in MIMO systems because the nullspace may rotate with the state and input trajectories.

To estimate uncertainty without differentiating measured states, use the predictor

\[ \dot{\hat{\mathbf{x} } } =\mathbf{A}_{0}\mathbf{x}+\mathbf{B}_{0}\mathbf{u} +\boldsymbol{\Phi}(\mathbf{x},\mathbf{u}) \hat{\boldsymbol{\theta} } +\mathbf{L}(\mathbf{x}-\hat{\mathbf{x} }), \]

where \( -\mathbf{L} \) is Hurwitz. Define prediction and parameter errors

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

Their dynamics are

\[ \dot{\tilde{\mathbf{x} } } =-\mathbf{L}\tilde{\mathbf{x} } +\boldsymbol{\Phi}\tilde{\boldsymbol{\theta} }. \]

Choose \( \mathbf{P}=\mathbf{P}^{T}\succ0 \) as the solution of

\[ \mathbf{L}^{T}\mathbf{P}+\mathbf{P}\mathbf{L} =\mathbf{Q},\qquad \mathbf{Q}=\mathbf{Q}^{T}\succ0, \]

and use the matrix-gradient adaptive law

\[ \boxed{ \dot{\hat{\boldsymbol{\theta} } } =\boldsymbol{\Gamma}\boldsymbol{\Phi}^{T} \mathbf{P}\tilde{\mathbf{x} } },\qquad \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}\succ0. \]

7.1 Stability proof

Consider

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

Since the true parameters are constant, \( \dot{\tilde{\boldsymbol{\theta} } } =-\dot{\hat{\boldsymbol{\theta} } } \). Differentiation gives

\[ \begin{aligned} \dot V &= 2\tilde{\mathbf{x} }^{T}\mathbf{P} \left(-\mathbf{L}\tilde{\mathbf{x} } +\boldsymbol{\Phi}\tilde{\boldsymbol{\theta} }\right) +2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\tilde{\boldsymbol{\theta} } }\\ &= -\tilde{\mathbf{x} }^{T} (\mathbf{L}^{T}\mathbf{P}+\mathbf{P}\mathbf{L}) \tilde{\mathbf{x} } +2\tilde{\mathbf{x} }^{T}\mathbf{P} \boldsymbol{\Phi}\tilde{\boldsymbol{\theta} } -2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Phi}^{T}\mathbf{P}\tilde{\mathbf{x} }\\ &=-\tilde{\mathbf{x} }^{T}\mathbf{Q}\tilde{\mathbf{x} } \le 0. \end{aligned} \]

Hence \( \tilde{\mathbf{x} } \) and \( \tilde{\boldsymbol{\theta} } \) are bounded, and \( \tilde{\mathbf{x} }\in L_{2} \). Under bounded regressor signals and bounded derivatives, Barbalat's lemma yields \( \tilde{\mathbf{x} }(t)\rightarrow\mathbf{0} \). Parameter convergence additionally requires excitation conditions such as PE.

7.2 Interpretation of matrix dimensions

\[ \boldsymbol{\Phi}^{T}\mathbf{P}\tilde{\mathbf{x} } \in\mathbb{R}^{p}. \]

Thus all output components of the predictor error contribute to every parameter update according to the regressor geometry. This cross-channel information flow is a defining feature of multivariable adaptation.

Consider a two-state, two-input plant with four uncertain parameters:

\[ \mathbf{A}_{0}= \begin{bmatrix}-1.2&0.4\\-0.3&-0.9\end{bmatrix},\quad \mathbf{A}_{1}= \begin{bmatrix}0.6&0\\0&0\end{bmatrix},\quad \mathbf{A}_{2}= \begin{bmatrix}0&0.35\\-0.25&0\end{bmatrix}, \]

\[ \mathbf{B}_{0}= \begin{bmatrix}1&0.2\\0.1&0.9\end{bmatrix},\quad \mathbf{B}_{1}= \begin{bmatrix}0.4&0\\0&0\end{bmatrix},\quad \mathbf{B}_{2}= \begin{bmatrix}0&0\\0&0.3\end{bmatrix}. \]

Let \( \boldsymbol{\theta} =[0.5,-0.4,0.25,-0.3]^{T} \). Then

\[ \mathbf{A} =\mathbf{A}_{0}+0.5\mathbf{A}_{1}-0.4\mathbf{A}_{2} = \begin{bmatrix}-0.9&0.26\\-0.2&-0.9\end{bmatrix}, \]

\[ \mathbf{B} =\mathbf{B}_{0}+0.25\mathbf{B}_{1}-0.3\mathbf{B}_{2} = \begin{bmatrix}1.1&0.2\\0.1&0.81\end{bmatrix}. \]

For \( \mathbf{x}=[x_{1},x_{2}]^{T} \) and \( \mathbf{u}=[u_{1},u_{2}]^{T} \), the regressor is

\[ \boldsymbol{\Phi}= \begin{bmatrix} 0.6x_{1} & 0.35x_{2} & 0.4u_{1} & 0\\ 0 & -0.25x_{1} & 0 & 0.3u_{2} \end{bmatrix}. \]

This matrix reveals the information pathways directly. The first parameter is visible only through the first state equation; the fourth parameter is visible only through the second input channel. If \( u_{2} \) is never excited, the fourth parameter cannot be identified, although prediction error may still be reduced through the predictor feedback term.

The implementations below use \( \mathbf{L}=2\mathbf{I}_{2} \), \( \mathbf{P}=0.25\mathbf{I}_{2} \), and \( \boldsymbol{\Gamma}=10\mathbf{I}_{4} \). The two input channels contain different frequencies so that the accumulated information matrix is full rank.

Recommended libraries for larger projects are NumPy for matrix operations, SciPy for integration and Lyapunov equations, and python-control for state-space analysis. The downloadable example uses NumPy and an explicit RK4 integrator so that every adaptive equation remains visible.

Chapter23_Lesson1.py

"""Chapter23_Lesson1.py

MIMO affine-parametric uncertainty and a Lyapunov-motivated adaptive predictor.
The example uses a 2-state, 2-input plant with four unknown scalar parameters.

Dependencies:
    numpy
"""

from __future__ import annotations

import math
from typing import Tuple

import numpy as np


A0 = np.array([[-1.2, 0.4], [-0.3, -0.9]], dtype=float)
A1 = np.array([[0.6, 0.0], [0.0, 0.0]], dtype=float)
A2 = np.array([[0.0, 0.35], [-0.25, 0.0]], dtype=float)

B0 = np.array([[1.0, 0.2], [0.1, 0.9]], dtype=float)
B1 = np.array([[0.4, 0.0], [0.0, 0.0]], dtype=float)
B2 = np.array([[0.0, 0.0], [0.0, 0.3]], dtype=float)

THETA_TRUE = np.array([0.5, -0.4, 0.25, -0.3], dtype=float)

# Predictor injection and Lyapunov matrices:
# L = 2 I, P = 0.25 I satisfy L.T @ P + P @ L = I.
L = 2.0 * np.eye(2)
P = 0.25 * np.eye(2)
GAMMA = 10.0 * np.eye(4)


def input_signal(t: float) -> np.ndarray:
    """Two independent input channels used to excite the matrix regressor."""
    return np.array(
        [
            math.sin(0.7 * t) + 0.5 * math.sin(1.9 * t),
            math.cos(1.1 * t) + 0.3 * math.sin(2.3 * t),
        ],
        dtype=float,
    )


def true_matrices() -> Tuple[np.ndarray, np.ndarray]:
    """Construct A(theta) and B(theta) from the true parameter vector."""
    theta = THETA_TRUE
    a_true = A0 + theta[0] * A1 + theta[1] * A2
    b_true = B0 + theta[2] * B1 + theta[3] * B2
    return a_true, b_true


def regressor(x: np.ndarray, u: np.ndarray) -> np.ndarray:
    """Return Phi(x,u) in R^(2x4), where residual dynamics equal Phi theta."""
    return np.column_stack((A1 @ x, A2 @ x, B1 @ u, B2 @ u))


def dynamics(t: float, state: np.ndarray) -> np.ndarray:
    """Combined plant, adaptive predictor, and parameter-estimator dynamics."""
    x = state[0:2]
    x_hat = state[2:4]
    theta_hat = state[4:8]

    u = input_signal(t)
    phi = regressor(x, u)
    a_true, b_true = true_matrices()

    x_dot = a_true @ x + b_true @ u
    prediction_error = x - x_hat

    x_hat_dot = (
        A0 @ x
        + B0 @ u
        + phi @ theta_hat
        + L @ prediction_error
    )
    theta_hat_dot = GAMMA @ phi.T @ P @ prediction_error

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


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


def simulate(final_time: float = 40.0, dt: float = 0.002) -> None:
    """Simulate and print parameter, prediction-error, and excitation metrics."""
    if final_time <= 0.0 or dt <= 0.0:
        raise ValueError("final_time and dt must be positive")

    state = np.array([1.0, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    information_matrix = np.zeros((4, 4), dtype=float)

    steps = int(round(final_time / dt))
    for k in range(steps):
        t = k * dt
        state = rk4_step(t, state, dt)

        x = state[0:2]
        phi = regressor(x, input_signal(t))
        information_matrix += phi.T @ phi * dt

    x = state[0:2]
    x_hat = state[2:4]
    theta_hat = state[4:8]
    eigenvalues = np.linalg.eigvalsh(information_matrix)

    np.set_printoptions(precision=6, suppress=True)
    print("True theta:              ", THETA_TRUE)
    print("Estimated theta:         ", theta_hat)
    print("Final prediction error:  ", np.linalg.norm(x - x_hat))
    print("Information eigenvalues: ", eigenvalues)
    print("Information rank:        ", np.linalg.matrix_rank(information_matrix, tol=1e-8))


if __name__ == "__main__":
    simulate()

For production C++, Eigen is the standard choice for fixed and dynamic matrices, while Boost.Odeint provides numerical integration. This implementation uses fixed-size std::array operations and therefore compiles without third-party dependencies.

Chapter23_Lesson1.cpp

// Chapter23_Lesson1.cpp
// Standalone C++17 simulation of affine MIMO uncertainty and an adaptive predictor.
// No third-party library is required. Eigen or Boost.Odeint can replace the
// small fixed-size routines in production software.

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

using Vec2 = std::array<double, 2>;
using Vec4 = std::array<double, 4>;
using Vec8 = std::array<double, 8>;
using Mat2 = std::array<std::array<double, 2>, 2>;
using Mat24 = std::array<std::array<double, 4>, 2>;
using Mat4 = std::array<std::array<double, 4>, 4>;

const Mat2 A0{ { {-1.2, 0.4}, {-0.3, -0.9} } };
const Mat2 A1{ { {0.6, 0.0}, {0.0, 0.0} } };
const Mat2 A2{ { {0.0, 0.35}, {-0.25, 0.0} } };

const Mat2 B0{ { {1.0, 0.2}, {0.1, 0.9} } };
const Mat2 B1{ { {0.4, 0.0}, {0.0, 0.0} } };
const Mat2 B2{ { {0.0, 0.0}, {0.0, 0.3} } };

const Vec4 THETA_TRUE{ {0.5, -0.4, 0.25, -0.3} };
constexpr double L_GAIN = 2.0;
constexpr double P_GAIN = 0.25;
constexpr double GAMMA_GAIN = 10.0;

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

Mat2 add_scaled(const Mat2& base, const Mat2& m1, double s1,
                const Mat2& m2, double s2) {
    Mat2 result = base;
    for (std::size_t i = 0; i < 2; ++i) {
        for (std::size_t j = 0; j < 2; ++j) {
            result[i][j] += s1 * m1[i][j] + s2 * m2[i][j];
        }
    }
    return result;
}

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

Vec2 input_signal(double t) {
    return Vec2{ {
        std::sin(0.7 * t) + 0.5 * std::sin(1.9 * t),
        std::cos(1.1 * t) + 0.3 * std::sin(2.3 * t)
    } };
}

Mat24 regressor(const Vec2& x, const Vec2& u) {
    const Vec2 c1 = mat_vec(A1, x);
    const Vec2 c2 = mat_vec(A2, x);
    const Vec2 c3 = mat_vec(B1, u);
    const Vec2 c4 = mat_vec(B2, u);

    return Mat24{ {
        { {c1[0], c2[0], c3[0], c4[0]} },
        { {c1[1], c2[1], c3[1], c4[1]} }
    } };
}

Vec8 dynamics(double t, const Vec8& state) {
    const Vec2 x{ {state[0], state[1]} };
    const Vec2 x_hat{ {state[2], state[3]} };
    const Vec4 theta_hat{ {state[4], state[5], state[6], state[7]} };
    const Vec2 u = input_signal(t);

    const Mat2 a_true = add_scaled(A0, A1, THETA_TRUE[0], A2, THETA_TRUE[1]);
    const Mat2 b_true = add_scaled(B0, B1, THETA_TRUE[2], B2, THETA_TRUE[3]);
    const Mat24 phi = regressor(x, u);

    const Vec2 ax = mat_vec(a_true, x);
    const Vec2 bu = mat_vec(b_true, u);
    const Vec2 x_dot = add(ax, bu);

    const Vec2 nominal = add(mat_vec(A0, x), mat_vec(B0, u));
    const Vec2 error{ {x[0] - x_hat[0], x[1] - x_hat[1]} };

    Vec2 uncertain{ {0.0, 0.0} };
    for (std::size_t i = 0; i < 2; ++i) {
        for (std::size_t j = 0; j < 4; ++j) {
            uncertain[i] += phi[i][j] * theta_hat[j];
        }
    }

    Vec2 x_hat_dot{ {
        nominal[0] + uncertain[0] + L_GAIN * error[0],
        nominal[1] + uncertain[1] + L_GAIN * error[1]
    } };

    Vec4 theta_hat_dot{ {0.0, 0.0, 0.0, 0.0} };
    for (std::size_t j = 0; j < 4; ++j) {
        theta_hat_dot[j] = GAMMA_GAIN * P_GAIN *
            (phi[0][j] * error[0] + phi[1][j] * error[1]);
    }

    return Vec8{ {
        x_dot[0], x_dot[1],
        x_hat_dot[0], x_hat_dot[1],
        theta_hat_dot[0], theta_hat_dot[1],
        theta_hat_dot[2], theta_hat_dot[3]
    } };
}

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

Vec8 rk4_step(double t, const Vec8& state, double dt) {
    const Vec8 k1 = dynamics(t, state);
    const Vec8 k2 = dynamics(t + 0.5 * dt, add_scaled_state(state, k1, 0.5 * dt));
    const Vec8 k3 = dynamics(t + 0.5 * dt, add_scaled_state(state, k2, 0.5 * dt));
    const Vec8 k4 = dynamics(t + dt, add_scaled_state(state, k3, dt));

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

int matrix_rank(Mat4 a, double tolerance = 1e-8) {
    int rank = 0;
    for (int col = 0; col < 4; ++col) {
        int pivot = rank;
        for (int row = rank + 1; row < 4; ++row) {
            if (std::abs(a[row][col]) > std::abs(a[pivot][col])) {
                pivot = row;
            }
        }
        if (std::abs(a[pivot][col]) <= tolerance) {
            continue;
        }
        std::swap(a[pivot], a[rank]);
        const double pivot_value = a[rank][col];
        for (int j = col; j < 4; ++j) {
            a[rank][j] /= pivot_value;
        }
        for (int row = 0; row < 4; ++row) {
            if (row == rank) {
                continue;
            }
            const double factor = a[row][col];
            for (int j = col; j < 4; ++j) {
                a[row][j] -= factor * a[rank][j];
            }
        }
        ++rank;
    }
    return rank;
}

int main() {
    const double final_time = 40.0;
    const double dt = 0.002;
    if (final_time <= 0.0 || dt <= 0.0) {
        throw std::invalid_argument("final_time and dt must be positive");
    }

    Vec8 state{ {1.0, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0} };
    Mat4 information{};

    const int steps = static_cast<int>(std::llround(final_time / dt));
    for (int k = 0; k < steps; ++k) {
        const double t = k * dt;
        state = rk4_step(t, state, dt);

        const Vec2 x{ {state[0], state[1]} };
        const Mat24 phi = regressor(x, input_signal(t));
        for (std::size_t i = 0; i < 4; ++i) {
            for (std::size_t j = 0; j < 4; ++j) {
                information[i][j] +=
                    (phi[0][i] * phi[0][j] + phi[1][i] * phi[1][j]) * dt;
            }
        }
    }

    const double error0 = state[0] - state[2];
    const double error1 = state[1] - state[3];
    const double error_norm = std::sqrt(error0 * error0 + error1 * error1);

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "True theta:       ";
    for (double value : THETA_TRUE) {
        std::cout << std::setw(11) << value;
    }
    std::cout << "\nEstimated theta:  ";
    for (std::size_t i = 4; i < 8; ++i) {
        std::cout << std::setw(11) << state[i];
    }
    std::cout << "\nPrediction error: " << error_norm;
    std::cout << "\nInformation rank: " << matrix_rank(information) << "\n";

    return 0;
}

EJML and Apache Commons Math are suitable for general Java matrix and ODE work. The following code keeps the two-state example self-contained and implements the same RK4 equations used in Python and C++.

Chapter23_Lesson1.java

// Chapter23_Lesson1.java
// Standalone Java simulation of affine MIMO uncertainty and an adaptive predictor.
// EJML or Apache Commons Math can replace these fixed-size routines in larger projects.

import java.util.Arrays;

public final class Chapter23_Lesson1 {
    private Chapter23_Lesson1() {}

    private static final double[][] A0 = { {-1.2, 0.4}, {-0.3, -0.9} };
    private static final double[][] A1 = { {0.6, 0.0}, {0.0, 0.0} };
    private static final double[][] A2 = { {0.0, 0.35}, {-0.25, 0.0} };

    private static final double[][] B0 = { {1.0, 0.2}, {0.1, 0.9} };
    private static final double[][] B1 = { {0.4, 0.0}, {0.0, 0.0} };
    private static final double[][] B2 = { {0.0, 0.0}, {0.0, 0.3} };

    private static final double[] THETA_TRUE = {0.5, -0.4, 0.25, -0.3};
    private static final double L_GAIN = 2.0;
    private static final double P_GAIN = 0.25;
    private static final double GAMMA_GAIN = 10.0;

    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[][] addScaled(
            double[][] base,
            double[][] matrix1,
            double scale1,
            double[][] matrix2,
            double scale2) {
        double[][] result = new double[2][2];
        for (int i = 0; i < 2; i++) {
            for (int j = 0; j < 2; j++) {
                result[i][j] =
                    base[i][j] + scale1 * matrix1[i][j] + scale2 * matrix2[i][j];
            }
        }
        return result;
    }

    private static double[] inputSignal(double t) {
        return new double[] {
            Math.sin(0.7 * t) + 0.5 * Math.sin(1.9 * t),
            Math.cos(1.1 * t) + 0.3 * Math.sin(2.3 * t)
        };
    }

    private static double[][] regressor(double[] x, double[] u) {
        double[] c1 = matVec(A1, x);
        double[] c2 = matVec(A2, x);
        double[] c3 = matVec(B1, u);
        double[] c4 = matVec(B2, u);

        return new double[][] {
            {c1[0], c2[0], c3[0], c4[0]},
            {c1[1], c2[1], c3[1], c4[1]}
        };
    }

    private static double[] dynamics(double t, double[] state) {
        double[] x = {state[0], state[1]};
        double[] xHat = {state[2], state[3]};
        double[] thetaHat = {state[4], state[5], state[6], state[7]};
        double[] u = inputSignal(t);

        double[][] aTrue = addScaled(A0, A1, THETA_TRUE[0], A2, THETA_TRUE[1]);
        double[][] bTrue = addScaled(B0, B1, THETA_TRUE[2], B2, THETA_TRUE[3]);
        double[][] phi = regressor(x, u);

        double[] ax = matVec(aTrue, x);
        double[] bu = matVec(bTrue, u);
        double[] xDot = {ax[0] + bu[0], ax[1] + bu[1]};

        double[] a0x = matVec(A0, x);
        double[] b0u = matVec(B0, u);
        double[] error = {x[0] - xHat[0], x[1] - xHat[1]};

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

        double[] xHatDot = {
            a0x[0] + b0u[0] + uncertain[0] + L_GAIN * error[0],
            a0x[1] + b0u[1] + uncertain[1] + L_GAIN * error[1]
        };

        double[] thetaHatDot = new double[4];
        for (int j = 0; j < 4; j++) {
            thetaHatDot[j] = GAMMA_GAIN * P_GAIN
                * (phi[0][j] * error[0] + phi[1][j] * error[1]);
        }

        return new double[] {
            xDot[0], xDot[1],
            xHatDot[0], xHatDot[1],
            thetaHatDot[0], thetaHatDot[1], thetaHatDot[2], thetaHatDot[3]
        };
    }

    private static double[] addScaledState(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[] state, double dt) {
        double[] k1 = dynamics(t, state);
        double[] k2 = dynamics(t + 0.5 * dt, addScaledState(state, k1, 0.5 * dt));
        double[] k3 = dynamics(t + 0.5 * dt, addScaledState(state, k2, 0.5 * dt));
        double[] k4 = dynamics(t + dt, addScaledState(state, k3, dt));

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

    private static int matrixRank(double[][] source, double tolerance) {
        double[][] matrix = new double[source.length][source[0].length];
        for (int i = 0; i < source.length; i++) {
            matrix[i] = Arrays.copyOf(source[i], source[i].length);
        }

        int rank = 0;
        for (int col = 0; col < 4; col++) {
            int pivot = rank;
            for (int row = rank + 1; row < 4; row++) {
                if (Math.abs(matrix[row][col]) > Math.abs(matrix[pivot][col])) {
                    pivot = row;
                }
            }
            if (Math.abs(matrix[pivot][col]) <= tolerance) {
                continue;
            }

            double[] temporary = matrix[pivot];
            matrix[pivot] = matrix[rank];
            matrix[rank] = temporary;

            double pivotValue = matrix[rank][col];
            for (int j = col; j < 4; j++) {
                matrix[rank][j] /= pivotValue;
            }

            for (int row = 0; row < 4; row++) {
                if (row == rank) {
                    continue;
                }
                double factor = matrix[row][col];
                for (int j = col; j < 4; j++) {
                    matrix[row][j] -= factor * matrix[rank][j];
                }
            }
            rank++;
        }
        return rank;
    }

    public static void main(String[] args) {
        final double finalTime = 40.0;
        final double dt = 0.002;
        if (finalTime <= 0.0 || dt <= 0.0) {
            throw new IllegalArgumentException("finalTime and dt must be positive");
        }

        double[] state = {1.0, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
        double[][] information = new double[4][4];

        int steps = (int) Math.round(finalTime / dt);
        for (int k = 0; k < steps; k++) {
            double t = k * dt;
            state = rk4Step(t, state, dt);

            double[] x = {state[0], state[1]};
            double[][] phi = regressor(x, inputSignal(t));
            for (int i = 0; i < 4; i++) {
                for (int j = 0; j < 4; j++) {
                    information[i][j] +=
                        (phi[0][i] * phi[0][j] + phi[1][i] * phi[1][j]) * dt;
                }
            }
        }

        double error0 = state[0] - state[2];
        double error1 = state[1] - state[3];
        double errorNorm = Math.sqrt(error0 * error0 + error1 * error1);

        System.out.println("True theta:       " + Arrays.toString(THETA_TRUE));
        System.out.println(
            "Estimated theta:  "
                + Arrays.toString(Arrays.copyOfRange(state, 4, 8))
        );
        System.out.printf("Prediction error: %.6f%n", errorNorm);
        System.out.println("Information rank: " + matrixRank(information, 1e-8));
    }
}

MATLAB supplies ode45, matrix factorizations, and direct plotting. Control System Toolbox can be used for state-space analysis, while Simulink is useful for separating the plant, predictor, and adaptive-law subsystems. The script also creates a fixed-parameter Simulink realization of the uncertain two-input plant when Simulink is licensed.

Chapter23_Lesson1.m

% Chapter23_Lesson1.m
% MIMO affine-parametric uncertainty and a Lyapunov-motivated adaptive predictor.
% Requires base MATLAB. The optional model builder requires Simulink.

clear; clc;

A0 = [-1.2,  0.4;
      -0.3, -0.9];
A1 = [ 0.6,  0.0;
       0.0,  0.0];
A2 = [ 0.0,  0.35;
      -0.25, 0.0];

B0 = [1.0, 0.2;
      0.1, 0.9];
B1 = [0.4, 0.0;
      0.0, 0.0];
B2 = [0.0, 0.0;
      0.0, 0.3];

thetaTrue = [0.5; -0.4; 0.25; -0.3];
Atrue = A0 + thetaTrue(1) * A1 + thetaTrue(2) * A2;
Btrue = B0 + thetaTrue(3) * B1 + thetaTrue(4) * B2;

L = 2.0 * eye(2);
P = 0.25 * eye(2);     % L'P + PL = I
Gamma = 10.0 * eye(4);

z0 = [1.0; -0.5; zeros(2, 1); zeros(4, 1)];
timeSpan = [0.0, 40.0];

ode = @(t, z) combinedDynamics( ...
    t, z, A0, A1, A2, B0, B1, B2, Atrue, Btrue, L, P, Gamma);

options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[t, z] = ode45(ode, timeSpan, z0, options);

thetaHat = z(end, 5:8).';
predictionError = norm(z(end, 1:2).' - z(end, 3:4).');

information = zeros(4, 4);
for k = 1:(numel(t) - 1)
    x = z(k, 1:2).';
    u = inputSignal(t(k));
    Phi = makeRegressor(x, u, A1, A2, B1, B2);
    information = information + Phi.' * Phi * (t(k + 1) - t(k));
end

fprintf('True theta:\n');
disp(thetaTrue.');
fprintf('Estimated theta:\n');
disp(thetaHat.');
fprintf('Final prediction error: %.6e\n', predictionError);
fprintf('Information eigenvalues:\n');
disp(eig(information).');
fprintf('Information rank: %d\n', rank(information, 1e-8));

figure('Name', 'Chapter23 Lesson1 MIMO Adaptive Predictor');
subplot(2, 1, 1);
plot(t, z(:, 1:2), 'LineWidth', 1.1);
hold on;
plot(t, z(:, 3:4), '--', 'LineWidth', 1.1);
grid on;
xlabel('Time (s)');
ylabel('States');
legend('x_1', 'x_2', 'xHat_1', 'xHat_2', 'Location', 'best');
title('Plant and predictor states');

subplot(2, 1, 2);
plot(t, z(:, 5:8), 'LineWidth', 1.1);
grid on;
xlabel('Time (s)');
ylabel('Parameter estimates');
legend('thetaHat_1', 'thetaHat_2', 'thetaHat_3', 'thetaHat_4', ...
       'Location', 'best');
title('Online matrix-uncertainty estimates');

% Optional fixed-parameter Simulink realization of the uncertain MIMO plant.
% It creates two input channels, a vector State-Space block, and a workspace sink.
if license('test', 'Simulink')
    buildSimulinkPlant(Atrue, Btrue);
else
    fprintf('Simulink is not licensed; skipped model generation.\n');
end

function dz = combinedDynamics( ...
        t, z, A0, A1, A2, B0, B1, B2, Atrue, Btrue, L, P, Gamma)
    x = z(1:2);
    xHat = z(3:4);
    thetaHat = z(5:8);

    u = inputSignal(t);
    Phi = makeRegressor(x, u, A1, A2, B1, B2);

    xDot = Atrue * x + Btrue * u;
    error = x - xHat;
    xHatDot = A0 * x + B0 * u + Phi * thetaHat + L * error;
    thetaHatDot = Gamma * Phi.' * P * error;

    dz = [xDot; xHatDot; thetaHatDot];
end

function u = inputSignal(t)
    u = [sin(0.7 * t) + 0.5 * sin(1.9 * t);
         cos(1.1 * t) + 0.3 * sin(2.3 * t)];
end

function Phi = makeRegressor(x, u, A1, A2, B1, B2)
    Phi = [A1 * x, A2 * x, B1 * u, B2 * u];
end

function buildSimulinkPlant(Atrue, Btrue)
    modelName = 'Chapter23_Lesson1_Simulink';
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    if exist([modelName, '.slx'], 'file')
        delete([modelName, '.slx']);
    end

    new_system(modelName);
    open_system(modelName);

    add_block('simulink/Sources/Sine Wave', [modelName, '/Input1'], ...
        'Amplitude', '1', 'Frequency', '0.7', ...
        'Position', [40, 60, 90, 90]);
    add_block('simulink/Sources/Sine Wave', [modelName, '/Input2'], ...
        'Amplitude', '1', 'Frequency', '1.1', 'Phase', 'pi/2', ...
        'Position', [40, 130, 90, 160]);
    add_block('simulink/Signal Routing/Mux', [modelName, '/Mux'], ...
        'Inputs', '2', 'Position', [130, 75, 135, 150]);
    add_block('simulink/Continuous/State-Space', [modelName, '/UncertainPlant'], ...
        'A', 'Atrue', 'B', 'Btrue', 'C', 'eye(2)', 'D', 'zeros(2,2)', ...
        'X0', '[1;-0.5]', 'Position', [190, 85, 310, 145]);
    add_block('simulink/Sinks/To Workspace', [modelName, '/PlantStates'], ...
        'VariableName', 'simulinkStates', 'SaveFormat', 'Structure With Time', ...
        'Position', [370, 95, 475, 135]);

    add_line(modelName, 'Input1/1', 'Mux/1');
    add_line(modelName, 'Input2/1', 'Mux/2');
    add_line(modelName, 'Mux/1', 'UncertainPlant/1');
    add_line(modelName, 'UncertainPlant/1', 'PlantStates/1');

    set_param(modelName, 'StopTime', '20');
    save_system(modelName);
    fprintf('Created %s.slx\n', modelName);
end

Mathematica's NDSolveValue, symbolic matrix operations, NIntegrate, and MatrixRank provide a compact environment for analyzing the predictor and information matrix. The downloadable file is a notebook expression.

Chapter23_Lesson1.nb

Notebook[{
  Cell["Chapter23_Lesson1.nb", "Title"],
  Cell[
    "MIMO affine-parametric uncertainty and a Lyapunov-motivated adaptive predictor.",
    "Text"
  ],
  Cell[BoxData[
    RowBox[{
      RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}], ";"
    }]
  ], "Input"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"A0", "=", "{ {-1.2,0.4},{-0.3,-0.9} }"}], ";",
      RowBox[{"A1", "=", "{ {0.6,0.0},{0.0,0.0} }"}], ";",
      RowBox[{"A2", "=", "{ {0.0,0.35},{-0.25,0.0} }"}], ";",
      RowBox[{"B0", "=", "{ {1.0,0.2},{0.1,0.9} }"}], ";",
      RowBox[{"B1", "=", "{ {0.4,0.0},{0.0,0.0} }"}], ";",
      RowBox[{"B2", "=", "{ {0.0,0.0},{0.0,0.3} }"}], ";"
    }]
  ], "Input"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"thetaTrue", "=", "{0.5,-0.4,0.25,-0.3}"}], ";",
      RowBox[{"Atrue", "=", "A0+thetaTrue[[1]] A1+thetaTrue[[2]] A2"}], ";",
      RowBox[{"Btrue", "=", "B0+thetaTrue[[3]] B1+thetaTrue[[4]] B2"}], ";",
      RowBox[{"L", "=", "2 IdentityMatrix[2]"}], ";",
      RowBox[{"P", "=", "0.25 IdentityMatrix[2]"}], ";",
      RowBox[{"Gamma", "=", "10 IdentityMatrix[4]"}], ";"
    }]
  ], "Input"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"u", "[", "t_", "]"}], ":=",
      "{Sin[0.7 t]+0.5 Sin[1.9 t],Cos[1.1 t]+0.3 Sin[2.3 t]}"
    }]
  ], "Input"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"phi", "[", "x_", ",", "input_", "]"}], ":=",
      "Transpose[{A1.x,A2.x,B1.input,B2.input}]"
    }]
  ], "Input"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"solution", "=", "NDSolveValue["}],
      RowBox[{"{"}],
      RowBox[{
        RowBox[{"x'[t]", "==", "Atrue.x[t]+Btrue.u[t]"}], ",",
        RowBox[{
          "xh'[t]",
          "==",
          "A0.x[t]+B0.u[t]+phi[x[t],u[t]].th[t]+L.(x[t]-xh[t])"
        }], ",",
        RowBox[{
          "th'[t]",
          "==",
          "Gamma.Transpose[phi[x[t],u[t]]].P.(x[t]-xh[t])"
        }], ",",
        RowBox[{"x[0]", "==", "{1.0,-0.5}"}], ",",
        RowBox[{"xh[0]", "==", "{0.0,0.0}"}], ",",
        RowBox[{"th[0]", "==", "{0.0,0.0,0.0,0.0}"}]
      }],
      RowBox[{"}"}], ",",
      RowBox[{"{x,xh,th}"}], ",",
      RowBox[{"{t,0,40}"}],
      RowBox[{"]"}], ";"
    }]
  ], "Input"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"thetaEstimate", "=", "solution[[3]][40]"}], ";",
      RowBox[{"predictionError", "=", "Norm[solution[[1]][40]-solution[[2]][40]]"}], ";",
      RowBox[{"information", "=", "NIntegrate["}],
      RowBox[{
        "Transpose[phi[solution[[1]][s],u[s]]].phi[solution[[1]][s],u[s]]"
      }],
      RowBox[{",{s,0,40},Method->\"LocalAdaptive\""}],
      RowBox[{"]"}], ";",
      RowBox[{"Print[\"True theta: \",thetaTrue]"}], ";",
      RowBox[{"Print[\"Estimated theta: \",thetaEstimate]"}], ";",
      RowBox[{"Print[\"Prediction error: \",predictionError]"}], ";",
      RowBox[{"Print[\"Information eigenvalues: \",Eigenvalues[information]]"}], ";",
      RowBox[{"Print[\"Information rank: \",MatrixRank[information]]"}], ";"
    }]
  ], "Input"],
  Cell[BoxData[
    RowBox[{"Plot", "["}],
    RowBox[{
      "Evaluate[Join[solution[[1]][t],solution[[2]][t]]]"
    }],
    RowBox[{",{t,0,40},PlotLegends->{\"x1\",\"x2\",\"xHat1\",\"xHat2\"},AxesLabel->{\"t\",\"state\"}"}],
    RowBox[{"]"}]
  ], "Input"],
  Cell[BoxData[
    RowBox[{"Plot", "["}],
    RowBox[{"Evaluate[solution[[3]][t]]"}],
    RowBox[{",{t,0,40},PlotLegends->{\"thetaHat1\",\"thetaHat2\",\"thetaHat3\",\"thetaHat4\"},AxesLabel->{\"t\",\"estimate\"}"}],
    RowBox[{"]"}]
  ], "Input"]
}]

14.1 Do not estimate impossible independent entries

An over-parameterized matrix model may be algebraically valid but practically unidentifiable. Begin with physical parameters and derive basis matrices from them. Only use entry-wise uncertainty when no reliable structure is available.

14.2 Scaling must be matrix-aware

State and input channels with different engineering units can make columns of \( \boldsymbol{\Phi} \) differ by orders of magnitude. Normalize states, inputs, or regressor columns before selecting a common adaptation gain. Otherwise one parameter direction may adapt rapidly while another is numerically frozen.

14.3 Derivative-based regression amplifies noise

The direct equation \( \mathbf{z}=\dot{\mathbf{x} }-\mathbf{A}_{0}\mathbf{x} -\mathbf{B}_{0}\mathbf{u} \) requires state derivatives. Numerical differentiation is noise-sensitive. Predictors, stable filters, integral regressions, or observer-based formulations avoid direct differentiation.

14.4 Full-rank instantaneous regressors are unnecessary

Since \( \boldsymbol{\Phi}\in\mathbb{R}^{n\times p} \) and often \( p>n \), the instantaneous regressor cannot have column rank \( p \). Identification comes from temporal accumulation in \( \int\boldsymbol{\Phi}^{T}\boldsymbol{\Phi}dt \).

14.5 Projection and leakage remain relevant

The ideal proof assumes constant parameters and an exact model. In implementations, projection constrains physically meaningful matrix parameters, while leakage or robust modifications limit drift caused by noise and unmodeled dynamics. These modifications should preserve known matrix properties such as positive actuator effectiveness or nonsingularity.

14.6 Monitor geometry, not only parameter traces

Useful diagnostics include the smallest eigenvalue and condition number of the information matrix, the norm of the unmatched uncertainty estimate, singular values of the estimated input matrix, and the prediction-error norm.

Problem 1 (Constructing a Matrix Regressor): Consider

\[ \mathbf{A}=\mathbf{A}_{0}+\theta_{1} \begin{bmatrix}1&0\\0&0\end{bmatrix} +\theta_{2} \begin{bmatrix}0&1\\-1&0\end{bmatrix}, \qquad \mathbf{B}=\mathbf{B}_{0}+\theta_{3} \begin{bmatrix}0&0\\0&1\end{bmatrix}. \]

Derive \( \boldsymbol{\Phi}(\mathbf{x},\mathbf{u}) \).

Solution: The three columns are the action of each basis matrix on the appropriate signal:

\[ \boldsymbol{\Phi}= \begin{bmatrix} x_{1} & x_{2} & 0\\ 0 & -x_{1} & u_{2} \end{bmatrix}. \]

Therefore \( \dot{\mathbf{x} }-\mathbf{A}_{0}\mathbf{x} -\mathbf{B}_{0}\mathbf{u} =\boldsymbol{\Phi}\boldsymbol{\theta} \).

Problem 2 (Matched-Uncertainty Test): Let \( \mathbf{B}_{0}=[1,0]^{T} \) and

\[ \Delta\mathbf{A}= \begin{bmatrix}a&b\\c&d\end{bmatrix}. \]

Find the matched and unmatched parts.

Solution: Since

\[ \mathbf{P}_{B}= \begin{bmatrix}1&0\\0&0\end{bmatrix}, \]

we obtain

\[ \Delta\mathbf{A}_{\mathrm{m} } =\begin{bmatrix}a&b\\0&0\end{bmatrix},\qquad \Delta\mathbf{A}_{\mathrm{u} } =\begin{bmatrix}0&0\\c&d\end{bmatrix}. \]

Full matching holds only when \( c=d=0 \). The second row cannot be generated by the available input direction.

Problem 3 (Loss of Multichannel Identifiability): Suppose two unknown input gains produce regressor columns \( \boldsymbol{\phi}_{1}(t)=\mathbf{v}u_{1}(t) \) and \( \boldsymbol{\phi}_{2}(t)=\mathbf{v}u_{2}(t) \). Show that if \( u_{2}(t)=c\,u_{1}(t) \), the two parameters are not separately identifiable.

Solution: Under the stated relation,

\[ \boldsymbol{\phi}_{2}(t)=c\,\boldsymbol{\phi}_{1}(t). \]

Hence the two columns of \( \boldsymbol{\Phi} \) are linearly dependent for all time, and the corresponding two-by-two information block has rank one. The residual depends only on \( \theta_{1}+c\theta_{2} \), so infinitely many parameter pairs produce identical data.

Problem 4 (Adaptive-Predictor Lyapunov Derivative): For \( \dot{\tilde{\mathbf{x} } } =-\mathbf{L}\tilde{\mathbf{x} }+ \boldsymbol{\Phi}\tilde{\boldsymbol{\theta} } \), verify that the law \( \dot{\hat{\boldsymbol{\theta} } } =\boldsymbol{\Gamma}\boldsymbol{\Phi}^{T} \mathbf{P}\tilde{\mathbf{x} } \) cancels the cross term in \( \dot V \).

Solution: The state-error term produces

\[ 2\tilde{\mathbf{x} }^{T}\mathbf{P} \boldsymbol{\Phi}\tilde{\boldsymbol{\theta} }. \]

Since \( \dot{\tilde{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\Phi}^{T} \mathbf{P}\tilde{\mathbf{x} } \), the parameter term produces

\[ -2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Phi}^{T}\mathbf{P}\tilde{\mathbf{x} }. \]

The two scalars are transposes of one another and cancel exactly. Therefore \( \dot V=-\tilde{\mathbf{x} }^{T} \mathbf{Q}\tilde{\mathbf{x} } \).

Problem 5 (Kronecker Parameterization): For \( n=2 \), write \( \Delta\mathbf{A}\mathbf{x} \) as a regression in \( \operatorname{vec}(\Delta\mathbf{A}) \).

Solution: Let

\[ \Delta\mathbf{A}= \begin{bmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{bmatrix}, \qquad \operatorname{vec}(\Delta\mathbf{A}) =\begin{bmatrix}a_{11}&a_{21}&a_{12}&a_{22}\end{bmatrix}^{T}. \]

Then

\[ \mathbf{x}^{T}\otimes\mathbf{I}_{2} = \begin{bmatrix} x_{1}&0&x_{2}&0\\ 0&x_{1}&0&x_{2} \end{bmatrix}, \]

and multiplication yields

\[ \begin{bmatrix} a_{11}x_{1}+a_{12}x_{2}\\ a_{21}x_{1}+a_{22}x_{2} \end{bmatrix} =\Delta\mathbf{A}\mathbf{x}. \]

Problem 6 (Numerical Uncertainty Reconstruction): Using the matrices in Section 8 and \( \boldsymbol{\theta} =[0.5,-0.4,0.25,-0.3]^{T} \), calculate \( \mathbf{A} \) and \( \mathbf{B} \), and determine whether \( \mathbf{A} \) is Hurwitz.

Solution: Direct substitution gives

\[ \mathbf{A}= \begin{bmatrix}-0.9&0.26\\-0.2&-0.9\end{bmatrix},\qquad \mathbf{B}= \begin{bmatrix}1.1&0.2\\0.1&0.81\end{bmatrix}. \]

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

\[ \det(s\mathbf{I}-\mathbf{A}) =(s+0.9)^{2}+0.052. \]

Thus the eigenvalues are \( -0.9\pm j\sqrt{0.052} \), whose real parts are negative. The uncertain plant matrix is Hurwitz.

MIMO parametric uncertainty is naturally represented by affine basis matrices and a vector of unknown physical coefficients. This structure leads to a vector regression \( \mathbf{z}=\boldsymbol{\Phi}\boldsymbol{\theta} \) or, equivalently, a Kronecker-product regression for vectorized matrix entries. Unlike the SISO case, uncertainty must be classified by its direction relative to the input subspace. Matched components can be represented through the control matrix, whereas unmatched components cannot be exactly canceled by the available actuators. Parameter identifiability depends on the rank of the accumulated multichannel information matrix, and zero prediction error does not by itself imply parameter convergence. A Lyapunov-designed adaptive predictor provides bounded estimates and asymptotic prediction-error convergence under standard boundedness assumptions. Lesson 2 will use these uncertainty structures to parameterize multivariable feedback and feedforward controller matrices.