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

• distinguish loop interaction from ordinary model uncertainty;
• compute a steady-state interaction measure and a static decoupler;
• derive a normalized matrix-gradient estimator for an unknown coupling matrix;
• construct a regularized adaptive inverse and explain why regularization is required;
• derive nominal and perturbed tracking-error dynamics for adaptive decoupling;
• interpret off-diagonal adaptive gains as interaction-compensation channels; and
• identify practical failure modes caused by poor excitation, saturation, noise, and near singularity.

The lesson assumes familiarity with linear state-space systems, transfer matrices, eigenvalues, singular values, Lyapunov equations, and the MIMO MRAC structure introduced in the preceding lesson. Students are not assumed to know advanced robust multivariable synthesis.

Consider a square linear plant with \(m\) inputs and \(m\) measured outputs:

\[ \mathbf{y}(s)=\mathbf{G}(s)\mathbf{u}(s), \qquad \mathbf{G}(s)=\begin{bmatrix} G_{11}(s)&\cdots&G_{1m}(s)\\ \vdots&\ddots&\vdots\\ G_{m1}(s)&\cdots&G_{mm}(s) \end{bmatrix}. \]

The diagonal terms describe intended input-output paths under a chosen pairing. The off-diagonal terms describe cross paths. A command applied to \(u_j\) therefore changes not only \(y_j\), but generally several outputs. Interaction is especially important when independently designed SISO loops compete with each other through these cross paths.

A decoupled closed-loop objective is commonly represented by a diagonal reference transfer matrix:

\[ \mathbf{y}(s)\approx \mathbf{M}(s)\mathbf{r}(s), \qquad \mathbf{M}(s)=\operatorname{diag}\{M_1(s),\ldots,M_m(s)\}. \]

Diagonal behavior does not require the physical plant itself to become diagonal. It requires the combined plant-controller map to suppress or compensate cross channels over the frequency range relevant to tracking and disturbance rejection.

Let the nonsingular steady-state gain matrix be \(\mathbf{K}=\mathbf{G}(0)\). The relative gain array is

\[ \boldsymbol{\Lambda}=\mathbf{K}\circ\mathbf{K}^{-T}, \]

where \(\circ\) denotes elementwise multiplication. Its element \(\lambda_{ij}\) compares the open-loop gain from \(u_j\) to \(y_i\) with the effective gain when all other outputs are regulated. For a nonsingular square matrix, each row and each column of \(\boldsymbol{\Lambda}\) sums to one.

A positive value near one suggests a relatively independent pairing. A large magnitude indicates strong interaction, while a negative value warns that closing the remaining loops may reverse the apparent sign of the selected path. The relative gain array is a steady-state diagnostic; it does not by itself capture dynamic phase, delays, nonminimum-phase zeros, or frequency-dependent interaction.

\[ \eta_{\mathrm{RGA} }=\left\|\boldsymbol{\Lambda}-\mathbf{I}\right\|_F \]

is one possible scalar summary for a proposed diagonal pairing. It should be interpreted together with the condition number of \(\mathbf{K}\) and the plant dynamics.

A static precompensator selects \(\mathbf{u}=\mathbf{D}_0\mathbf{v}\). If \(\mathbf{K}\) is known and nonsingular, one may choose

\[ \mathbf{D}_0=\mathbf{K}^{-1}\mathbf{K}_d, \]

where \(\mathbf{K}_d\) is diagonal. Then

\[ \mathbf{G}(0)\mathbf{D}_0=\mathbf{K}_d. \]

This removes steady-state interaction but generally not transient interaction. Exact dynamic decoupling would seek a transfer matrix \(\mathbf{D}(s)\) satisfying

\[ \mathbf{G}(s)\mathbf{D}(s)=\mathbf{Q}_d(s), \]

with diagonal \(\mathbf{Q}_d(s)\). Exact inversion can be noncausal, improper, unstable, or excessively high gain when the plant has delays, right-half-plane zeros, uncertain relative degree, or small singular values. Adaptive decoupling therefore usually targets approximate diagonalization over an operating region rather than unrestricted exact inversion.

A useful conceptual architecture separates a virtual diagonal controller from an adaptive interaction model. The virtual controller produces the behavior desired for independent channels. The adaptive block maps the virtual command into physical actuator commands.

flowchart TD
  R["Reference vector r"] --> C["Diagonal virtual controller"]
  Y["Measured output vector y"] --> C
  C --> V["Virtual command v"]
  V --> D["Adaptive decoupler Dhat"]
  D --> U["Physical input vector u"]
  U --> P["Interacting MIMO plant"]
  P --> Y
  U --> E["Online coupling estimator"]
  Y --> E
  E --> D
  S["Projection, regularization, saturation"] --> D
        

The architecture contains three distinct feedback mechanisms:

1. ordinary output feedback drives the tracking error toward zero;
2. parameter adaptation reduces prediction error in the coupling model; and
3. regularization and projection prevent a numerically dangerous inverse.

To isolate the main idea, consider a first-order square MIMO model

\[ \dot{\mathbf{x} }=-\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u}+\mathbf{d}, \qquad \mathbf{y}=\mathbf{x}, \]

where \(\mathbf{A}\) is known and stable, while the nonsingular input matrix \(\mathbf{B}\) contains uncertain direct and cross-channel gains. Define the regression output

\[ \mathbf{z}=\dot{\mathbf{x} }+\mathbf{A}\mathbf{x} =\mathbf{B}\mathbf{u}+\mathbf{d}. \]

The prediction and prediction error are

\[ \hat{\mathbf{z} }=\hat{\mathbf{B} }\mathbf{u}, \qquad \boldsymbol{\varepsilon}_B=\mathbf{z}-\hat{\mathbf{B} }\mathbf{u}. \]

For the instantaneous loss \(J=\tfrac{1}{2}\boldsymbol{\varepsilon}_B^T \boldsymbol{\varepsilon}_B\), matrix differentiation gives

\[ \frac{\partial J}{\partial\hat{\mathbf{B} } } =-\boldsymbol{\varepsilon}_B\mathbf{u}^T. \]

A normalized gradient estimator is therefore

\[ \dot{\hat{\mathbf{B} } } =\gamma\frac{\boldsymbol{\varepsilon}_B\mathbf{u}^T} {\epsilon+\mathbf{u}^T\mathbf{u} }, \qquad \gamma > 0,\quad \epsilon > 0. \]

Normalization limits the update magnitude when the actuator vector becomes large. Projection can additionally constrain \(\hat{\mathbf{B} }\) to a physically meaningful compact set. In a real implementation, direct differentiation of noisy states is usually replaced by stable filtering, an observer, or an integral regression.

Directly using \(\hat{\mathbf{B} }^{-1}\) is unsafe when the estimate is inaccurate or nearly singular. A Tikhonov-regularized right inverse is

\[ \mathbf{D}_{\rho}(\hat{\mathbf{B} }) =\hat{\mathbf{B} }^T \left(\hat{\mathbf{B} }\hat{\mathbf{B} }^T+\rho\mathbf{I}\right)^{-1}, \qquad \rho > 0. \]

If the singular value decomposition is \(\hat{\mathbf{B} }=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^T\), then

\[ \mathbf{D}_{\rho} =\mathbf{V}\operatorname{diag}\left\{ \frac{\sigma_i}{\sigma_i^2+\rho}\right\}\mathbf{U}^T. \]

Every inverse singular gain is bounded. For \(f(\sigma)=\sigma/(\sigma^2+\rho)\), differentiation gives

\[ f'(\sigma)=\frac{\rho-\sigma^2}{(\sigma^2+\rho)^2}. \]

The maximum occurs at \(\sigma=\sqrt{\rho}\), yielding

\[ \left\|\mathbf{D}_{\rho}\right\|_2 \leq \frac{1}{2\sqrt{\rho} }. \]

Thus \(\rho\) trades decoupling accuracy for actuator-gain limitation and numerical robustness. Small regularization improves inversion when the estimate is well conditioned; larger regularization is safer near singularity but leaves more residual interaction.

Let the virtual command be

\[ \mathbf{v}=\mathbf{A}\mathbf{r} +\mathbf{K}(\mathbf{r}-\mathbf{x}), \qquad \mathbf{K}=\mathbf{K}^T > 0, \]

and apply

\[ \mathbf{u}=\mathbf{D}_{\rho}(\hat{\mathbf{B} })\mathbf{v}. \]

For a constant reference and ideal inversion \(\mathbf{B}\mathbf{D}_{\rho}=\mathbf{I}\), the tracking error \(\mathbf{e}=\mathbf{x}-\mathbf{r}\) obeys

\[ \dot{\mathbf{e} }=-(\mathbf{A}+\mathbf{K})\mathbf{e}. \]

If \(\mathbf{A}+\mathbf{K}\) is positive definite, the channels have independent exponentially stable first-order error dynamics. With imperfect inversion, define

\[ \boldsymbol{\Delta}_D =\mathbf{B}\mathbf{D}_{\rho}(\hat{\mathbf{B} })-\mathbf{I}. \]

The error equation becomes

\[ \dot{\mathbf{e} } =-(\mathbf{A}+\mathbf{K})\mathbf{e} +\boldsymbol{\Delta}_D\mathbf{v}+\mathbf{d}. \]

Using \(V=\tfrac{1}{2}\mathbf{e}^T\mathbf{e}\) and \(\alpha=\lambda_{\min}(\mathbf{A}+\mathbf{K})\),

\[ \dot V\leq -\alpha\|\mathbf{e}\|^2 +\|\mathbf{e}\|\left( \|\boldsymbol{\Delta}_D\|\|\mathbf{v}\|+\|\mathbf{d}\|\right). \]

If \(\|\mathbf{v}\|\leq \bar v\) and \(\|\mathbf{d}\|\leq \bar d\), then \(\dot V<0\) whenever

\[ \|\mathbf{e}\| > \frac{\|\boldsymbol{\Delta}_D\|\bar v+\bar d}{\alpha}. \]

This establishes an ultimate-bound interpretation: better coupling estimation and better-conditioned inversion shrink the residual tracking neighborhood, while stronger nominal feedback enlarges \(\alpha\). The result is not a complete proof for the coupled estimator-controller dynamics, but it isolates the central robustness mechanism.

Instead of estimating a plant matrix and then inverting it, a direct MIMO adaptive controller updates the controller matrices themselves. Consider

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

with control law

\[ \mathbf{u}=\hat{\mathbf{K} }_x\mathbf{x} +\hat{\mathbf{K} }_r\mathbf{r}. \]

Ideal gains satisfy the matching equations

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

Choosing diagonal stable \(\mathbf{A}_m\) specifies the desired decoupled reference dynamics. The off-diagonal entries of \(\hat{\mathbf{K} }_x\) and \(\hat{\mathbf{K} }_r\) are then adaptive interaction compensators, not unwanted terms.

For pedagogical purposes, suppose \(\mathbf{B}\) is known and \(\mathbf{P}>0\) solves

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

Let \(\tilde{\mathbf{K} }_x=\hat{\mathbf{K} }_x-\mathbf{K}_x^*\) and similarly for \(\tilde{\mathbf{K} }_r\). Use

\[ V=\mathbf{e}^T\mathbf{P}\mathbf{e} +\frac{1}{\gamma_x}\|\tilde{\mathbf{K} }_x\|_F^2 +\frac{1}{\gamma_r}\|\tilde{\mathbf{K} }_r\|_F^2. \]

The matrix update laws

\[ \dot{\hat{\mathbf{K} } }_x =-\gamma_x\mathbf{B}^T\mathbf{P}\mathbf{e}\mathbf{x}^T, \qquad \dot{\hat{\mathbf{K} } }_r =-\gamma_r\mathbf{B}^T\mathbf{P}\mathbf{e}\mathbf{r}^T \]

cancel the parameter-error cross terms and give

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

In a genuine uncertain MIMO plant, \(\mathbf{B}\) is not generally available for these update laws. A known control-direction structure, sign-definite high-frequency gain, triangular factorization, or another structural assumption is required. This is why MIMO adaptive control needs stronger matching and sign information than a collection of unrelated SISO adaptive loops.

Tracking can be satisfactory even when the entire coupling matrix is not identified. Parameter convergence requires sufficient excitation of all input directions. For the matrix regression \(\mathbf{z}=\mathbf{B}\mathbf{u}\), a finite-window excitation condition is

\[ \int_t^{t+T}\mathbf{u}(\tau)\mathbf{u}^T(\tau)\,d\tau \geq \alpha_u\mathbf{I}, \qquad T>0,\quad \alpha_u>0. \]

If commands repeatedly excite only one actuator combination, different columns of \(\mathbf{B}\) cannot be separated. A small bounded probing signal can improve identifiability, but it introduces tracking ripple and may be unacceptable in safety-critical operation. Natural command variation, scheduled calibration intervals, or stored-data methods can sometimes provide excitation with less disruption.

Useful online metrics include

\[ \eta_D(t)=\left\| \hat{\mathbf{B} }(t)\mathbf{D}_{\rho}(t)-\mathbf{I} \right\|_F, \qquad \kappa_{\rho}(t)=\frac{\sigma_{\max}(\hat{\mathbf{B} })} {\max\{\sigma_{\min}(\hat{\mathbf{B} }),\sqrt{\rho}\} }. \]

A large decoupling residual or a rapidly increasing condition indicator is a reason to slow adaptation, increase regularization, freeze the inverse, or fall back to a conservative controller.

flowchart TD
  A["Read r, y, and actuator data"] --> B["Filter signals and form regression"]
  B --> C["Update coupling estimate"]
  C --> D["Project estimate into allowed set"]
  D --> E["Check singular values and condition indicator"]
  E -->|"acceptable"| F["Compute regularized decoupler"]
  E -->|"unsafe"| G["Freeze adaptation or use fallback map"]
  F --> H["Compute virtual command and physical input"]
  G --> H
  H --> I["Apply saturation and rate limits"]
  I --> J["Monitor tracking and interaction residuals"]
  J --> A
        

Important implementation rules are:

• Filter derivatives: do not differentiate noisy measurements directly.
• Use projection: preserve known signs, magnitude bounds, and nonsingularity margins.
• Regularize every inverse: an adaptive estimate can cross a singular set transiently.
• Coordinate with saturation: estimator updates based on commanded rather than applied input create bias.
• Limit adaptation rate: fast matrix updates amplify noise and can excite neglected dynamics.
• Provide fallback control: freeze or blend out the adaptive decoupler when monitors fail.
• Log matrix diagnostics: singular values, parameter rates, residual interaction, and actuator utilization are essential evidence.

The executable example uses

\[ \dot{\mathbf{x} }= -\begin{bmatrix}1&0\\0&1.3\end{bmatrix}\mathbf{x} +\begin{bmatrix}1&0.65\\0.45&1.10\end{bmatrix}\mathbf{u}, \qquad \mathbf{y}=\mathbf{x}. \]

Both cross gains are substantial. The estimate starts from a diagonal matrix and is updated from the normalized regression. The virtual command is

\[ \mathbf{v}=\mathbf{A}\mathbf{r} +2.5\mathbf{I}(\mathbf{r}-\mathbf{x})+\mathbf{p}(t), \]

where \(\mathbf{p}(t)\) is a small probing vector. The physical input is the regularized inverse command, followed by saturation. The programs compare interaction and tracking integrals and log the online matrix estimate.

With the supplied parameters, the tested implementations converge to an estimate near \(\begin{bmatrix}0.927&0.603\\0.407&1.072\end{bmatrix}\). In the first command interval, the interaction integral for the nominally uncommanded second output is approximately 0.107, compared with approximately 0.688 for the non-decoupled baseline. These values illustrate the example only; they are not general performance guarantees.

The Python implementation uses NumPy for matrix operations and Matplotlib for output and parameter trajectories. Install dependencies with python -m pip install numpy matplotlib.

Chapter23_Lesson4.py

"""Chapter23_Lesson4.py
Adaptive decoupling of a 2x2 first-order MIMO plant.

The estimator learns the unknown input matrix B from
    z = x_dot + A x = B u
and the controller uses a regularized right inverse of B_hat.
"""

from __future__ import annotations

import numpy as np
import matplotlib.pyplot as plt


def reference(t: float) -> np.ndarray:
    if t < 5.0:
        return np.array([1.0, 0.0])
    if t < 10.0:
        return np.array([1.0, 1.0])
    return np.array([0.0, 1.0])


def regularized_right_inverse(matrix: np.ndarray, rho: float) -> np.ndarray:
    rows = matrix.shape[0]
    return matrix.T @ np.linalg.inv(matrix @ matrix.T + rho * np.eye(rows))


def simulate(adaptive: bool = True) -> dict[str, np.ndarray | float]:
    dt = 0.002
    final_time = 15.0
    time = np.arange(0.0, final_time + dt, dt)

    # x_dot = -A x + B u, y = x
    a_matrix = np.diag([1.0, 1.3])
    b_true = np.array([[1.0, 0.65], [0.45, 1.10]])
    b_hat = np.array([[0.80, 0.0], [0.0, 0.80]], dtype=float)

    k_matrix = np.diag([2.5, 2.5])
    gamma = 8.0
    rho = 0.04
    epsilon = 0.02
    u_limit = 8.0

    x = np.zeros(2)
    y_log = np.zeros((time.size, 2))
    r_log = np.zeros((time.size, 2))
    u_log = np.zeros((time.size, 2))
    b_log = np.zeros((time.size, 4))

    interaction_iae = 0.0
    tracking_iae = 0.0

    for index, t in enumerate(time):
        r = reference(t)
        error = r - x

        # A small bounded probe improves identifiability. In practice it may be
        # replaced by naturally rich commands or a scheduled test signal.
        probe = (
            np.array([0.08 * np.sin(2.3 * t), 0.08 * np.cos(1.7 * t)])
            if adaptive
            else np.zeros(2)
        )
        virtual_input = a_matrix @ r + k_matrix @ error + probe

        if adaptive:
            decoupler = regularized_right_inverse(b_hat, rho)
            u = decoupler @ virtual_input
        else:
            # Baseline: two nominal diagonal loops with no interaction compensation.
            u = virtual_input

        u = np.clip(u, -u_limit, u_limit)
        x_dot = -a_matrix @ x + b_true @ u

        if adaptive:
            measured_regression_output = x_dot + a_matrix @ x
            prediction_error = measured_regression_output - b_hat @ u
            normalization = epsilon + float(u @ u)
            b_hat += dt * gamma * np.outer(prediction_error, u) / normalization
            b_hat = np.clip(b_hat, -2.0, 2.0)

        x += dt * x_dot

        y_log[index] = x
        r_log[index] = r
        u_log[index] = u
        b_log[index] = b_hat.reshape(-1)

        if t < 5.0:
            interaction_iae += abs(x[1]) * dt
        tracking_iae += np.sum(np.abs(error)) * dt

    return {
        "time": time,
        "y": y_log,
        "r": r_log,
        "u": u_log,
        "b_hat": b_log,
        "interaction_iae": interaction_iae,
        "tracking_iae": tracking_iae,
    }


def main() -> None:
    adaptive_result = simulate(adaptive=True)
    baseline_result = simulate(adaptive=False)

    print("Adaptive final B_hat:")
    print(adaptive_result["b_hat"][-1].reshape(2, 2))
    print(f"Adaptive interaction IAE: {adaptive_result['interaction_iae']:.6f}")
    print(f"Baseline interaction IAE: {baseline_result['interaction_iae']:.6f}")
    print(f"Adaptive total tracking IAE: {adaptive_result['tracking_iae']:.6f}")
    print(f"Baseline total tracking IAE: {baseline_result['tracking_iae']:.6f}")

    time = adaptive_result["time"]
    y = adaptive_result["y"]
    r = adaptive_result["r"]
    b_hat = adaptive_result["b_hat"]

    plt.figure()
    plt.plot(time, r[:, 0], "--", label="r1")
    plt.plot(time, y[:, 0], label="y1")
    plt.plot(time, r[:, 1], "--", label="r2")
    plt.plot(time, y[:, 1], label="y2")
    plt.xlabel("Time (s)")
    plt.ylabel("Output")
    plt.title("Adaptive decoupling: reference tracking")
    plt.grid(True)
    plt.legend()

    plt.figure()
    labels = ["Bhat11", "Bhat12", "Bhat21", "Bhat22"]
    for column, label in enumerate(labels):
        plt.plot(time, b_hat[:, column], label=label)
    plt.xlabel("Time (s)")
    plt.ylabel("Estimated coefficient")
    plt.title("Online estimate of the input-coupling matrix")
    plt.grid(True)
    plt.legend()
    plt.show()


if __name__ == "__main__":
    main()

The C++17 implementation performs all 2×2 operations explicitly and writes a CSV file for plotting in another tool. It requires only the standard library.

Chapter23_Lesson4.cpp

// Chapter23_Lesson4.cpp
// Adaptive decoupling of a 2x2 first-order MIMO plant.
#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>

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

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

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

Vec2 scale(double value, const Vec2& vector) {
    return {value * vector[0], value * vector[1]};
}

Vec2 multiply(const Mat2& matrix, const Vec2& vector) {
    return {
        matrix[0][0] * vector[0] + matrix[0][1] * vector[1],
        matrix[1][0] * vector[0] + matrix[1][1] * vector[1]
    };
}

Mat2 transpose(const Mat2& matrix) {
    return { { {matrix[0][0], matrix[1][0]}, {matrix[0][1], matrix[1][1]} } };
}

Mat2 multiply(const Mat2& left, const Mat2& right) {
    Mat2 result{};
    for (int i = 0; i < 2; ++i) {
        for (int j = 0; j < 2; ++j) {
            for (int k = 0; k < 2; ++k) {
                result[i][j] += left[i][k] * right[k][j];
            }
        }
    }
    return result;
}

Mat2 inverse(const Mat2& matrix) {
    const double determinant = matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
    if (std::abs(determinant) < 1.0e-12) {
        throw std::runtime_error("Singular 2x2 matrix");
    }
    return { { {matrix[1][1] / determinant, -matrix[0][1] / determinant},
             {-matrix[1][0] / determinant, matrix[0][0] / determinant} } };
}

Mat2 regularizedRightInverse(const Mat2& matrix, double rho) {
    const Mat2 matrixTranspose = transpose(matrix);
    Mat2 gram = multiply(matrix, matrixTranspose);
    gram[0][0] += rho;
    gram[1][1] += rho;
    return multiply(matrixTranspose, inverse(gram));
}

Vec2 reference(double time) {
    if (time < 5.0) return {1.0, 0.0};
    if (time < 10.0) return {1.0, 1.0};
    return {0.0, 1.0};
}

int main() {
    const double dt = 0.002;
    const double finalTime = 15.0;
    const Mat2 A{ { {1.0, 0.0}, {0.0, 1.3} } };
    const Mat2 B{ { {1.0, 0.65}, {0.45, 1.10} } };
    Mat2 Bhat{ { {0.80, 0.0}, {0.0, 0.80} } };
    const Mat2 K{ { {2.5, 0.0}, {0.0, 2.5} } };
    const double gamma = 8.0;
    const double rho = 0.04;
    const double epsilon = 0.02;
    const double inputLimit = 8.0;

    Vec2 x{0.0, 0.0};
    double interactionIAE = 0.0;
    double trackingIAE = 0.0;

    std::ofstream csv("Chapter23_Lesson4_cpp_results.csv");
    if (!csv) {
        std::cerr << "Could not open output CSV file.\n";
        return 1;
    }
    csv << "time,r1,r2,y1,y2,u1,u2,Bhat11,Bhat12,Bhat21,Bhat22\n";

    const int steps = static_cast<int>(finalTime / dt);
    for (int step = 0; step <= steps; ++step) {
        const double time = step * dt;
        const Vec2 r = reference(time);
        const Vec2 error = subtract(r, x);
        const Vec2 probe{0.08 * std::sin(2.3 * time), 0.08 * std::cos(1.7 * time)};
        const Vec2 virtualInput = add(add(multiply(A, r), multiply(K, error)), probe);

        const Mat2 decoupler = regularizedRightInverse(Bhat, rho);
        Vec2 u = multiply(decoupler, virtualInput);
        u[0] = std::clamp(u[0], -inputLimit, inputLimit);
        u[1] = std::clamp(u[1], -inputLimit, inputLimit);

        const Vec2 xDot = add(scale(-1.0, multiply(A, x)), multiply(B, u));
        const Vec2 regressionOutput = add(xDot, multiply(A, x));
        const Vec2 predictionError = subtract(regressionOutput, multiply(Bhat, u));
        const double normalization = epsilon + u[0] * u[0] + u[1] * u[1];

        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                Bhat[i][j] += dt * gamma * predictionError[i] * u[j] / normalization;
                Bhat[i][j] = std::clamp(Bhat[i][j], -2.0, 2.0);
            }
        }

        x = add(x, scale(dt, xDot));
        if (time < 5.0) interactionIAE += std::abs(x[1]) * dt;
        trackingIAE += (std::abs(error[0]) + std::abs(error[1])) * dt;

        if (step % 5 == 0) {
            csv << std::setprecision(10) << time << ',' << r[0] << ',' << r[1] << ','
                << x[0] << ',' << x[1] << ',' << u[0] << ',' << u[1] << ','
                << Bhat[0][0] << ',' << Bhat[0][1] << ','
                << Bhat[1][0] << ',' << Bhat[1][1] << '\n';
        }
    }

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Final Bhat = [[" << Bhat[0][0] << ", " << Bhat[0][1]
              << "], [" << Bhat[1][0] << ", " << Bhat[1][1] << "]]\n";
    std::cout << "Interaction IAE = " << interactionIAE << '\n';
    std::cout << "Tracking IAE = " << trackingIAE << '\n';
    std::cout << "Wrote Chapter23_Lesson4_cpp_results.csv\n";
    return 0;
}

The Java implementation uses primitive arrays and writes a CSV file. It can be compiled with a standard JDK and has no third-party dependency.

Chapter23_Lesson4.java

// Chapter23_Lesson4.java
// Adaptive decoupling of a 2x2 first-order MIMO plant.
import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.Locale;

public final class Chapter23_Lesson4 {
    private Chapter23_Lesson4() {}

    private static double[] reference(double time) {
        if (time < 5.0) return new double[] {1.0, 0.0};
        if (time < 10.0) return new double[] {1.0, 1.0};
        return new double[] {0.0, 1.0};
    }

    private static double[] multiply(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[][] multiply(double[][] left, double[][] right) {
        double[][] result = new double[2][2];
        for (int i = 0; i < 2; i++) {
            for (int j = 0; j < 2; j++) {
                for (int k = 0; k < 2; k++) {
                    result[i][j] += left[i][k] * right[k][j];
                }
            }
        }
        return result;
    }

    private static double[][] transpose(double[][] matrix) {
        return new double[][] {
            {matrix[0][0], matrix[1][0]},
            {matrix[0][1], matrix[1][1]}
        };
    }

    private static double[][] inverse(double[][] matrix) {
        double determinant = matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
        if (Math.abs(determinant) < 1.0e-12) {
            throw new IllegalArgumentException("Singular 2x2 matrix");
        }
        return new double[][] {
            { matrix[1][1] / determinant, -matrix[0][1] / determinant },
            { -matrix[1][0] / determinant, matrix[0][0] / determinant }
        };
    }

    private static double[][] regularizedRightInverse(double[][] matrix, double rho) {
        double[][] matrixTranspose = transpose(matrix);
        double[][] gram = multiply(matrix, matrixTranspose);
        gram[0][0] += rho;
        gram[1][1] += rho;
        return multiply(matrixTranspose, inverse(gram));
    }

    private static double clamp(double value, double lower, double upper) {
        return Math.max(lower, Math.min(upper, value));
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        double dt = 0.002;
        double finalTime = 15.0;
        double[][] aMatrix = { {1.0, 0.0}, {0.0, 1.3} };
        double[][] bTrue = { {1.0, 0.65}, {0.45, 1.10} };
        double[][] bHat = { {0.80, 0.0}, {0.0, 0.80} };
        double[][] kMatrix = { {2.5, 0.0}, {0.0, 2.5} };
        double gamma = 8.0;
        double rho = 0.04;
        double epsilon = 0.02;
        double inputLimit = 8.0;

        double[] x = {0.0, 0.0};
        double interactionIAE = 0.0;
        double trackingIAE = 0.0;

        try (PrintWriter csv = new PrintWriter(new BufferedWriter(
                new FileWriter("Chapter23_Lesson4_java_results.csv")))) {
            csv.println("time,r1,r2,y1,y2,u1,u2,Bhat11,Bhat12,Bhat21,Bhat22");
            int steps = (int) Math.round(finalTime / dt);
            for (int step = 0; step <= steps; step++) {
                double time = step * dt;
                double[] r = reference(time);
                double[] error = {r[0] - x[0], r[1] - x[1]};
                double[] ar = multiply(aMatrix, r);
                double[] ke = multiply(kMatrix, error);
                double[] probe = {0.08 * Math.sin(2.3 * time), 0.08 * Math.cos(1.7 * time)};
                double[] virtualInput = {
                    ar[0] + ke[0] + probe[0],
                    ar[1] + ke[1] + probe[1]
                };

                double[][] decoupler = regularizedRightInverse(bHat, rho);
                double[] u = multiply(decoupler, virtualInput);
                u[0] = clamp(u[0], -inputLimit, inputLimit);
                u[1] = clamp(u[1], -inputLimit, inputLimit);

                double[] ax = multiply(aMatrix, x);
                double[] bu = multiply(bTrue, u);
                double[] xDot = {-ax[0] + bu[0], -ax[1] + bu[1]};
                double[] regressionOutput = {xDot[0] + ax[0], xDot[1] + ax[1]};
                double[] predicted = multiply(bHat, u);
                double[] predictionError = {
                    regressionOutput[0] - predicted[0],
                    regressionOutput[1] - predicted[1]
                };
                double normalization = epsilon + u[0] * u[0] + u[1] * u[1];

                for (int i = 0; i < 2; i++) {
                    for (int j = 0; j < 2; j++) {
                        bHat[i][j] += dt * gamma * predictionError[i] * u[j] / normalization;
                        bHat[i][j] = clamp(bHat[i][j], -2.0, 2.0);
                    }
                }

                x[0] += dt * xDot[0];
                x[1] += dt * xDot[1];
                if (time < 5.0) interactionIAE += Math.abs(x[1]) * dt;
                trackingIAE += (Math.abs(error[0]) + Math.abs(error[1])) * dt;

                if (step % 5 == 0) {
                    csv.printf(Locale.US,
                        "%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f%n",
                        time, r[0], r[1], x[0], x[1], u[0], u[1],
                        bHat[0][0], bHat[0][1], bHat[1][0], bHat[1][1]);
                }
            }
        }

        System.out.printf(Locale.US, "Final Bhat = [[%.6f, %.6f], [%.6f, %.6f]]%n",
            bHat[0][0], bHat[0][1], bHat[1][0], bHat[1][1]);
        System.out.printf(Locale.US, "Interaction IAE = %.6f%n", interactionIAE);
        System.out.printf(Locale.US, "Tracking IAE = %.6f%n", trackingIAE);
        System.out.println("Wrote Chapter23_Lesson4_java_results.csv");
    }
}

The MATLAB script performs the simulation and includes a block-level Simulink mapping. In Simulink, the estimator and regularized inverse are conveniently implemented with MATLAB Function blocks, while State-Space, Saturation, Sum, Gain, and logging blocks implement the remaining loop.

Chapter23_Lesson4.m

%% Chapter23_Lesson4.m
% Adaptive decoupling of a 2x2 first-order MIMO plant.
clear; clc; close all;

dt = 0.002;
T = 15;
time = 0:dt:T;

% x_dot = -A*x + B*u, y = x
A = diag([1.0, 1.3]);
B = [1.0, 0.65; 0.45, 1.10];
Bhat = [0.80, 0.0; 0.0, 0.80];
K = diag([2.5, 2.5]);

gamma = 8.0;
rho = 0.04;
epsilon = 0.02;
uLimit = 8.0;

x = zeros(2,1);
yLog = zeros(2,numel(time));
rLog = zeros(2,numel(time));
uLog = zeros(2,numel(time));
BLog = zeros(4,numel(time));
interactionIAE = 0;
trackingIAE = 0;

for k = 1:numel(time)
    t = time(k);
    if t < 5
        r = [1; 0];
    elseif t < 10
        r = [1; 1];
    else
        r = [0; 1];
    end

    e = r - x;
    probe = [0.08*sin(2.3*t); 0.08*cos(1.7*t)];
    v = A*r + K*e + probe;

    % Tikhonov-regularized right inverse: D = Bhat'/(Bhat*Bhat' + rho*I)
    D = Bhat' / (Bhat*Bhat' + rho*eye(2));
    u = D*v;
    u = min(max(u, -uLimit), uLimit);

    xDot = -A*x + B*u;

    % Regression z = x_dot + A*x = B*u.
    z = xDot + A*x;
    predictionError = z - Bhat*u;
    Bhat = Bhat + dt*gamma*(predictionError*u')/(epsilon + u'*u);
    Bhat = min(max(Bhat, -2.0), 2.0);

    x = x + dt*xDot;

    yLog(:,k) = x;
    rLog(:,k) = r;
    uLog(:,k) = u;
    BLog(:,k) = Bhat(:);

    if t < 5
        interactionIAE = interactionIAE + abs(x(2))*dt;
    end
    trackingIAE = trackingIAE + sum(abs(e))*dt;
end

fprintf('Final Bhat =\n');
disp(Bhat);
fprintf('Interaction IAE = %.6f\n', interactionIAE);
fprintf('Tracking IAE = %.6f\n', trackingIAE);

figure;
plot(time, rLog(1,:), '--', time, yLog(1,:), ...
     time, rLog(2,:), '--', time, yLog(2,:));
grid on;
xlabel('Time (s)'); ylabel('Output');
title('Adaptive decoupling: reference tracking');
legend('r_1','y_1','r_2','y_2','Location','best');

figure;
plot(time, BLog(1,:), time, BLog(2,:), time, BLog(3,:), time, BLog(4,:));
grid on;
xlabel('Time (s)'); ylabel('Estimated coefficient');
title('Online estimate of the input-coupling matrix');
legend('Bhat_{11}','Bhat_{21}','Bhat_{12}','Bhat_{22}','Location','best');

%% Simulink implementation mapping
% 1. Use a State-Space block for x_dot = -A*x + B*u.
% 2. Build z = x_dot + A*x using Sum and Gain blocks.
% 3. Implement the normalized estimator in a MATLAB Function block.
% 4. Compute D with a MATLAB Function block using the regularized inverse.
% 5. Place Saturation blocks after u = D*v.
% 6. Log y, r, u, and Bhat with To Workspace blocks.

The notebook contains one executable input cell that runs the matrix estimator, adaptive inverse, trajectory logging, and plotting workflow.

Chapter23_Lesson4.nb

Notebook[{
  Cell["Chapter 23, Lesson 4: Adaptive Decoupling and Interaction Compensation", "Title"],
  Cell["Run the input cell to simulate online interaction estimation and regularized adaptive decoupling.", "Text"],
  Cell[BoxData["ClearAll[\"Global`*\"];
dt = 0.002; tFinal = 15.0;
a = DiagonalMatrix[{1.0, 1.3}]; bTrue = { {1.0, 0.65}, {0.45, 1.10} };
bHat = { {0.80, 0.0}, {0.0, 0.80} }; kGain = DiagonalMatrix[{2.5, 2.5}];
gamma = 8.0; rho = 0.04; epsilon = 0.02; inputLimit = 8.0;

reference[t_] := Piecewise[{
  { {1.0, 0.0}, t < 5.0},
  { {1.0, 1.0}, t < 10.0}
}, {0.0, 1.0}];

regularizedRightInverse[m_, regularization_] :=
  Transpose[m].Inverse[m.Transpose[m] + regularization IdentityMatrix[2]];

x = {0.0, 0.0}; interactionIAE = 0.0; trackingIAE = 0.0;
history = Reap[
  Do[
    r = reference[t]; error = r - x;
    probe = {0.08 Sin[2.3 t], 0.08 Cos[1.7 t]};
    virtualInput = a.r + kGain.error + probe;
    decoupler = regularizedRightInverse[bHat, rho];
    u = Clip[decoupler.virtualInput, {-inputLimit, inputLimit}];
    xDot = -a.x + bTrue.u;

    regressionOutput = xDot + a.x;
    predictionError = regressionOutput - bHat.u;
    normalization = epsilon + u.u;
    bHat = bHat + dt gamma Outer[Times, predictionError, u]/normalization;
    bHat = Map[Clip[#, {-2.0, 2.0}] &, bHat, {2}];

    x = x + dt xDot;
    If[t < 5.0, interactionIAE += Abs[x[[2]]] dt];
    trackingIAE += Total[Abs[error]] dt;
    Sow[Join[{t}, r, x, u, Flatten[bHat]]],
    {t, 0.0, tFinal, dt}
  ]
][[2, 1]];

Print[\"Final Bhat = \", MatrixForm[bHat]];
Print[\"Interaction IAE = \", interactionIAE];
Print[\"Tracking IAE = \", trackingIAE];

trackingPlot = ListLinePlot[
  {history[[All, {1, 2}]], history[[All, {1, 4}]],
   history[[All, {1, 3}]], history[[All, {1, 5}]]},
  PlotLegends -> {\"r1\", \"y1\", \"r2\", \"y2\"}, Frame -> True,
  FrameLabel -> {\"Time (s)\", \"Output\"},
  PlotLabel -> \"Adaptive decoupling: reference tracking\"];

parameterPlot = ListLinePlot[
  Table[history[[All, {1, column}]], {column, 8, 11}],
  PlotLegends -> {\"Bhat11\", \"Bhat12\", \"Bhat21\", \"Bhat22\"},
  Frame -> True, FrameLabel -> {\"Time (s)\", \"Estimated coefficient\"},
  PlotLabel -> \"Online estimate of the input-coupling matrix\"];

Column[{trackingPlot, parameterPlot}]"], "Input"]
}, WindowSize -> {1100, 800}, StyleDefinitions -> "Default.nb"]        

Problem 1 — Static interaction analysis: Let \(\mathbf{K}=\begin{bmatrix}1&0.6\\0.4&1.2\end{bmatrix}\). Compute its relative gain array and the static decoupler that produces unit steady-state diagonal gain.

Solution: The determinant is

\[ \det(\mathbf{K})=1.2-(0.6)(0.4)=0.96. \]

Therefore,

\[ \mathbf{K}^{-1}=\frac{1}{0.96} \begin{bmatrix}1.2&-0.6\\-0.4&1\end{bmatrix} =\begin{bmatrix}1.25&-0.625\\-0.4167&1.0417\end{bmatrix}. \]

The relative gain array is

\[ \boldsymbol{\Lambda} =\mathbf{K}\circ\mathbf{K}^{-T} =\begin{bmatrix}1.25&-0.25\\-0.25&1.25\end{bmatrix}. \]

The negative off-diagonal relative gains indicate unfavorable cross pairings. For unit diagonal steady-state gain, choose \(\mathbf{D}_0=\mathbf{K}^{-1}\), which gives \(\mathbf{K}\mathbf{D}_0=\mathbf{I}\).

Problem 2 — Bound on a regularized inverse: Prove that \(\|\mathbf{D}_{\rho}\|_2\leq 1/(2\sqrt{\rho})\) for \(\rho>0\).

Solution: Under singular value decomposition, the singular gains of the regularized inverse are

\[ f(\sigma)=\frac{\sigma}{\sigma^2+\rho}. \]

Differentiating gives

\[ f'(\sigma)=\frac{\rho-\sigma^2}{(\sigma^2+\rho)^2}. \]

The positive maximum occurs at \(\sigma=\sqrt{\rho}\). Substitution gives \(f_{\max}=1/(2\sqrt{\rho})\). The spectral norm equals the largest singular gain, proving the result.

Problem 3 — Nominal decoupled error dynamics: For \(\dot{\mathbf{x} }=-\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u}\), use \(\mathbf{u}=\mathbf{B}^{-1}[\mathbf{A}\mathbf{r}+ \mathbf{K}(\mathbf{r}-\mathbf{x})]\). Assume a constant reference. Derive the tracking-error dynamics and prove exponential stability when \(\mathbf{A}+\mathbf{K}>0\).

Solution: Substitution gives

\[ \dot{\mathbf{x} }=-\mathbf{A}\mathbf{x} +\mathbf{A}\mathbf{r}+\mathbf{K}(\mathbf{r}-\mathbf{x}). \]

For \(\mathbf{e}=\mathbf{x}-\mathbf{r}\) and constant \(\mathbf{r}\),

\[ \dot{\mathbf{e} }=-(\mathbf{A}+\mathbf{K})\mathbf{e}. \]

With \(V=\tfrac{1}{2}\mathbf{e}^T\mathbf{e}\),

\[ \dot V=-\mathbf{e}^T(\mathbf{A}+\mathbf{K})\mathbf{e} \leq -\lambda_{\min}(\mathbf{A}+\mathbf{K})\|\mathbf{e}\|^2. \]

Positive definiteness yields global exponential convergence. If the matrix is diagonal, each error channel evolves independently.

Problem 4 — Residual-interaction bound: Suppose \(\dot{\mathbf{e} }=-\mathbf{H}\mathbf{e}+ \boldsymbol{\Delta}\mathbf{v}+\mathbf{d}\), where \(\mathbf{H}=\mathbf{H}^T>0\), \(\|\mathbf{v}\|\leq\bar v\), and \(\|\mathbf{d}\|\leq\bar d\). Find an ultimate bound.

Solution: With \(V=\tfrac{1}{2}\|\mathbf{e}\|^2\),

\[ \dot V\leq -\lambda_{\min}(\mathbf{H})\|\mathbf{e}\|^2 +\|\mathbf{e}\|(\|\boldsymbol{\Delta}\|\bar v+\bar d). \]

Hence \(\dot V<0\) outside the ball

\[ \|\mathbf{e}\|\leq \frac{\|\boldsymbol{\Delta}\|\bar v+\bar d} {\lambda_{\min}(\mathbf{H})}. \]

The bound explicitly separates residual coupling, command magnitude, disturbance magnitude, and nominal feedback strength.

Problem 5 — Excitation and saturation: Explain why repeatedly commanding \(\mathbf{u}(t)=a(t)[1\;1]^T\) cannot identify all four entries of an unknown 2×2 matrix \(\mathbf{B}\). Then explain why estimator updates must use the applied saturated input rather than the unconstrained command.

Solution: The regression is

\[ \mathbf{z}=\mathbf{B}\mathbf{u} =a(t)\mathbf{B}\begin{bmatrix}1\\1\end{bmatrix}. \]

Measurements reveal only the sum of the two columns of \(\mathbf{B}\). The excitation Gramian has rank one, so infinitely many matrices with the same column sum produce identical data. Independent actuator directions are required for full matrix identification.

Under saturation, the plant receives \(\mathbf{u}_{\mathrm{applied} }\), not the larger requested vector. Using the unconstrained command in the regression attributes the saturation discrepancy to plant parameters, causing biased estimates and possible parameter drift. The actual applied input must therefore be fed back to the estimator.

Adaptive decoupling treats interaction as an uncertain multivariable map that must be estimated or compensated online. Static inverse decoupling is simple but only removes steady-state interaction and can become dangerous near singularity. A normalized matrix estimator combined with a regularized inverse provides an implementable conceptual architecture. Perfect inversion yields diagonal nominal error dynamics; estimation error, disturbances, and regularization appear as residual interaction terms and lead naturally to ultimate-bound analysis. Direct MIMO MRAC provides a second viewpoint in which off-diagonal controller gains adapt directly, but it requires matching and control-direction assumptions. Practical success depends on excitation, projection, filtered regressions, saturation-aware updates, conditioning monitors, and a safe fallback controller.