Let a continuous-time LTI realization be \( \Sigma=(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}) \), where \( \mathbf{A}\in\mathbb{R}^{n\times n} \), \( \mathbf{B}\in\mathbb{R}^{n\times m} \), \( \mathbf{C}\in\mathbb{R}^{p\times n} \), and \( \mathbf{D}\in\mathbb{R}^{p\times m} \):

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

The transfer matrix associated with this realization is

\[ \mathbf{G}(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+ \mathbf{D}. \]

A realization reduction algorithm constructs another realization \( \Sigma_r=(\mathbf{A}_r,\mathbf{B}_r,\mathbf{C}_r,\mathbf{D}) \) of order \( r\le n \) such that the transfer matrix is preserved exactly:

\[ \mathbf{C}_r(s\mathbf{I}_r-\mathbf{A}_r)^{-1}\mathbf{B}_r+ \mathbf{D} = \mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+ \mathbf{D}. \]

In this chapter, reduction means eliminating dynamically redundant states, not approximating the input-output behavior. Thus the reduced order is the McMillan degree of \( \mathbf{G}(s) \) when the realization is reduced to minimal order.

Reduction algorithms use the reachable subspace and the unobservable subspace. The finite-dimensional reachability matrix is

\[ \mathcal{R}=\begin{bmatrix}\mathbf{B}&\mathbf{AB}&\mathbf{A}^2\mathbf{B}&\cdots&\mathbf{A}^{n-1}\mathbf{B}\end{bmatrix},\qquad \mathscr{R}=\operatorname{range}(\mathcal{R}). \]

The finite-dimensional observability matrix is

\[ \mathcal{O}=\begin{bmatrix}\mathbf{C}\\ \mathbf{CA}\\ \mathbf{CA}^2\\ \vdots\\ \mathbf{CA}^{n-1} \end{bmatrix},\qquad \mathscr{N}=\ker(\mathcal{O}). \]

The subspace \( \mathscr{R} \) contains all states that can be generated from zero initial condition by some input. The subspace \( \mathscr{N} \) contains all initial states that produce identically zero free output. Therefore, unreachable states cannot be excited by the input, while unobservable states cannot be detected at the output.

\[ \operatorname{rank}(\mathcal{R})=n \quad\Longleftrightarrow\quad \text{reachable},\qquad \operatorname{rank}(\mathcal{O})=n \quad\Longleftrightarrow\quad \text{observable}. \]

A finite-dimensional LTI realization is minimal exactly when both rank conditions hold. Hence every exact realization reduction algorithm is a constructive method for discarding directions outside these two properties.

The first reduction removes states that cannot be reached from the input. Choose a full-column-rank matrix \( \mathbf{T}_r\in\mathbb{R}^{n\times r} \) whose columns form a basis for \( \mathscr{R} \):

\[ \operatorname{range}(\mathbf{T}_r)=\mathscr{R},\qquad r=\operatorname{rank}(\mathcal{R}). \]

Because \( \mathscr{R} \) is invariant under \( \mathbf{A} \) and contains \( \operatorname{range}(\mathbf{B}) \), there exist matrices \( \mathbf{A}_r \) and \( \mathbf{B}_r \) satisfying

\[ \mathbf{A}\mathbf{T}_r=\mathbf{T}_r\mathbf{A}_r, \qquad \mathbf{B}=\mathbf{T}_r\mathbf{B}_r. \]

With any left inverse \( \mathbf{T}_r^\dagger \) satisfying \( \mathbf{T}_r^\dagger\mathbf{T}_r=\mathbf{I}_r \), the reachable restricted realization is

\[ \mathbf{A}_r=\mathbf{T}_r^\dagger\mathbf{A}\mathbf{T}_r, \qquad \mathbf{B}_r=\mathbf{T}_r^\dagger\mathbf{B}, \qquad \mathbf{C}_r=\mathbf{C}\mathbf{T}_r, \qquad \mathbf{D}_r=\mathbf{D}. \]

If \( \mathbf{T}_r \) is chosen orthonormally, then \( \mathbf{T}_r^\dagger=\mathbf{T}_r^T \). In exact arithmetic, reachability trimming changes internal coordinates but does not alter the transfer matrix because all terms \( \mathbf{A}^k\mathbf{B} \) already lie in \( \mathscr{R} \).

After restricting to reachable states, the remaining realization may still contain unobservable internal directions. Let \( \mathscr{N}=\ker(\mathcal{O}) \) for the current reachable realization. Choose a coordinate transformation \( \mathbf{T}=[\mathbf{T}_u\;\mathbf{T}_o] \) such that \( \operatorname{range}(\mathbf{T}_u)=\mathscr{N} \). In transformed coordinates \( \mathbf{x}=\mathbf{T}\begin{bmatrix}\mathbf{x}_u^T&\mathbf{x}_o^T\end{bmatrix}^T \), the realization has the block form

\[ \bar{\mathbf{A} }=\mathbf{T}^{-1}\mathbf{A}\mathbf{T} =\begin{bmatrix} \mathbf{A}_{uu} & \mathbf{A}_{uo}\\ \mathbf{0} & \mathbf{A}_{oo} \end{bmatrix},\qquad \bar{\mathbf{B} }=\mathbf{T}^{-1}\mathbf{B} =\begin{bmatrix}\mathbf{B}_u\\\mathbf{B}_o\end{bmatrix},\qquad \bar{\mathbf{C} }=\mathbf{C}\mathbf{T} =\begin{bmatrix}\mathbf{0}&\mathbf{C}_o\end{bmatrix}. \]

The unobservable coordinate \( \mathbf{x}_u \) has no direct output contribution. Since \( \mathscr{N} \) is \( \mathbf{A} \)-invariant, there is no coupling from \( \mathbf{x}_u \) into the observable quotient dynamics. The reduced observable realization is

\[ \mathbf{A}_{min}=\mathbf{A}_{oo},\qquad \mathbf{B}_{min}=\mathbf{B}_o, \qquad \mathbf{C}_{min}=\mathbf{C}_o, \qquad \mathbf{D}_{min}=\mathbf{D}. \]

The order after both trimming steps is the dimension of the reachable part modulo the unobservable directions that remain inside it.

\[ n_{min}=\dim(\mathscr{R})-\dim(\mathscr{R}\cap\mathscr{N}). \]

The complete conceptual procedure is a two-pass trimming algorithm. In exact symbolic work the ranks are exact. In floating-point work the same steps are used, but singular-value thresholds become decisive; that numerical issue is the subject of the next lesson.

flowchart TD
  A["Start with realization A, B, C, D"] --> B["Build reachability matrix R"]
  B --> C["Find basis for range(R)"]
  C --> D["Project onto reachable coordinates"]
  D --> E["Build observability matrix of reachable subsystem"]
  E --> F["Find basis for null(O)"]
  F --> G["Transform so unobservable states come first"]
  G --> H["Discard unobservable block"]
  H --> I["Return minimal realization with same transfer matrix"]
        

A second basic route starts not from an existing state-space realization but from input-output data such as Markov parameters. For a strictly proper discrete-time system, the Markov parameters are

\[ \mathbf{H}_0=\mathbf{D},\qquad \mathbf{H}_k=\mathbf{C}\mathbf{A}^{k-1}\mathbf{B},\quad k=1,2,\dots. \]

A block Hankel matrix built from these parameters factors as an observability factor times a reachability factor:

\[ \mathcal{H}_{q,r}=\begin{bmatrix} \mathbf{H}_1 & \mathbf{H}_2 & \cdots & \mathbf{H}_r\\ \mathbf{H}_2 & \mathbf{H}_3 & \cdots & \mathbf{H}_{r+1}\\ \vdots & \vdots & \ddots & \vdots\\ \mathbf{H}_q & \mathbf{H}_{q+1} & \cdots & \mathbf{H}_{q+r-1} \end{bmatrix} = \begin{bmatrix}\mathbf{C}\\\mathbf{CA}\\\vdots\\\mathbf{CA}^{q-1}\end{bmatrix} \begin{bmatrix}\mathbf{B}&\mathbf{AB}&\cdots&\mathbf{A}^{r-1}\mathbf{B}\end{bmatrix}. \]

Therefore the rank of a sufficiently large Hankel matrix gives the minimal state dimension:

\[ n_{min}=\operatorname{rank}(\mathcal{H}_{q,r}) \quad \text{for sufficiently large } q,r. \]

Conceptually, the Ho-Kalman construction factorizes \( \mathcal{H}_{q,r} \) into an observability-like part and a reachability-like part, then extracts \( \mathbf{A} \) from the one-step shifted Hankel matrix. This is realization construction and reduction at the same time: only the rank-revealed dynamic directions are retained.

flowchart TD
  M["Markov parameters H0, H1, H2, ..."] --> H["Build block Hankel matrix"]
  H --> R["Rank reveals minimal order"]
  R --> F["Factor Hankel into observability factor and reachability factor"]
  F --> S["Use shifted Hankel matrix to recover A"]
  S --> G["Recover B, C, and D"]
        

The exactness of the trimming algorithm can be proved by comparing Markov parameters. For the original realization, the strictly proper Markov sequence is \( \mathbf{M}_k=\mathbf{C}\mathbf{A}^k\mathbf{B} \) for \( k=0,1,2,\dots \). Since \( \mathbf{B}=\mathbf{T}_r\mathbf{B}_r \) and \( \mathbf{A}\mathbf{T}_r=\mathbf{T}_r\mathbf{A}_r \), induction gives

\[ \mathbf{A}^k\mathbf{B}=\mathbf{T}_r\mathbf{A}_r^k\mathbf{B}_r, \qquad k=0,1,2,\dots. \]

Thus reachability trimming preserves all Markov parameters:

\[ \mathbf{C}\mathbf{A}^k\mathbf{B} =\mathbf{C}\mathbf{T}_r\mathbf{A}_r^k\mathbf{B}_r =\mathbf{C}_r\mathbf{A}_r^k\mathbf{B}_r. \]

For observability trimming, use the transformed block realization from Section 4. Because \( \bar{\mathbf{C} }=[\mathbf{0}\;\mathbf{C}_o] \) and the lower-left block of \( \bar{\mathbf{A} } \) is zero, the observable quotient gives

\[ \bar{\mathbf{C} }\bar{\mathbf{A} }^k\bar{\mathbf{B} } =\mathbf{C}_o\mathbf{A}_{oo}^k\mathbf{B}_o, \qquad k=0,1,2,\dots. \]

Hence all Markov parameters, and therefore the transfer matrix, remain unchanged. The result is minimal because the final realization is both reachable and observable; by the minimal-realization theorem from Lesson 1, no lower-order realization can realize the same transfer matrix.

Consider the nonminimal realization

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

The reachability matrix and observability matrix are

\[ \mathcal{R}=\begin{bmatrix}1&-1&1\\0&0&0\\1&-3&9\end{bmatrix}, \qquad \mathcal{O}=\begin{bmatrix}1&0&0\\-1&0&0\\1&0&0\end{bmatrix}. \]

Hence \( \operatorname{rank}(\mathcal{R})=2 \) and \( \operatorname{rank}(\mathcal{O})=1 \). The state associated with eigenvalue \( -2 \) is not reachable, and the state associated with eigenvalue \( -3 \) is reachable but unobservable. The transfer function is

\[ G(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B} =\begin{bmatrix}1&0&0\end{bmatrix} \begin{bmatrix}\frac{1}{s+1}&0&0\\0&\frac{1}{s+2}&0\\0&0&\frac{1}{s+3}\end{bmatrix} \begin{bmatrix}1\\0\\1\end{bmatrix} =\frac{1}{s+1}. \]

A minimal realization is therefore

\[ A_{min}=-1, \qquad B_{min}=1, \qquad C_{min}=1, \qquad D_{min}=0. \]

The Python version uses NumPy SVD routines to compute bases for the reachable column space and the unobservable null space. Related control libraries include python-control for state-space models and SciPy for numerical linear algebra.

Chapter20_Lesson4.py

# Chapter20_Lesson4.py
"""
Basic algorithms for exact realization reduction.

This educational script demonstrates the two-stage reduction:
1. restrict the realization to the reachable subspace;
2. quotient out the unobservable subspace of the reachable realization.

Required packages:
    pip install numpy scipy
For production control work, compare the result with python-control:
    pip install control
"""

from __future__ import annotations
import numpy as np
from numpy.linalg import svd


def controllability_matrix(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """Return C_R = [B, AB, ..., A^(n-1)B]."""
    n = A.shape[0]
    blocks = []
    Apow = np.eye(n)
    for _ in range(n):
        blocks.append(Apow @ B)
        Apow = Apow @ A
    return np.hstack(blocks)


def observability_matrix(A: np.ndarray, C: np.ndarray) -> np.ndarray:
    """Return O = [C; CA; ...; CA^(n-1)]."""
    n = A.shape[0]
    blocks = []
    Apow = np.eye(n)
    for _ in range(n):
        blocks.append(C @ Apow)
        Apow = Apow @ A
    return np.vstack(blocks)


def numerical_rank(M: np.ndarray, tol: float | None = 1e-10) -> int:
    s = svd(M, compute_uv=False)
    if tol is None:
        tol = max(M.shape) * np.finfo(float).eps * (s[0] if s.size else 1.0)
    return int(np.sum(s > tol))


def column_space_basis(M: np.ndarray, tol: float | None = 1e-10) -> np.ndarray:
    """Orthonormal basis for range(M), returned as columns."""
    U, s, _ = svd(M, full_matrices=True)
    if tol is None:
        tol = max(M.shape) * np.finfo(float).eps * (s[0] if s.size else 1.0)
    r = int(np.sum(s > tol))
    return U[:, :r]


def null_space_basis(M: np.ndarray, tol: float | None = 1e-10) -> np.ndarray:
    """Orthonormal basis for null(M), returned as columns."""
    _, s, Vt = svd(M, full_matrices=True)
    if tol is None:
        tol = max(M.shape) * np.finfo(float).eps * (s[0] if s.size else 1.0)
    r = int(np.sum(s > tol))
    return Vt.T[:, r:]


def exact_minimal_reduction(A: np.ndarray, B: np.ndarray, C: np.ndarray, D: np.ndarray):
    """
    Return a minimal realization externally equivalent to (A,B,C,D),
    assuming exact rank decisions are reliable for the given data.
    """
    n = A.shape[0]

    # Stage 1: keep only the reachable subspace.
    R = controllability_matrix(A, B)
    Qr = column_space_basis(R)
    A1 = Qr.T @ A @ Qr
    B1 = Qr.T @ B
    C1 = C @ Qr

    # Stage 2: remove the unobservable part of the reachable realization.
    O1 = observability_matrix(A1, C1)
    Q_unobs = null_space_basis(O1)
    q = Q_unobs.shape[1]
    if q == 0:
        return A1, B1, C1, D, Qr, np.eye(A1.shape[0])

    # Orthogonal complement completes the coordinate basis.
    Q_obs = null_space_basis(Q_unobs.T)
    T = np.hstack([Q_unobs, Q_obs])

    Ahat = T.T @ A1 @ T
    Bhat = T.T @ B1
    Chat = C1 @ T

    Amin = Ahat[q:, q:]
    Bmin = Bhat[q:, :]
    Cmin = Chat[:, q:]
    return Amin, Bmin, Cmin, D, Qr, T


def transfer_value(A: np.ndarray, B: np.ndarray, C: np.ndarray, D: np.ndarray, s: complex) -> np.ndarray:
    """Evaluate G(s)=C(sI-A)^(-1)B+D."""
    n = A.shape[0]
    return C @ np.linalg.solve(s * np.eye(n) - A, B) + D


def main() -> None:
    # Nonminimal realization: states 2 and 3 do not affect the input-output map.
    A = np.diag([-1.0, -2.0, -3.0])
    B = np.array([[1.0], [0.0], [1.0]])
    C = np.array([[1.0, 0.0, 0.0]])
    D = np.array([[0.0]])

    R = controllability_matrix(A, B)
    O = observability_matrix(A, C)
    print("rank controllability matrix:", numerical_rank(R), "of", A.shape[0])
    print("rank observability matrix:", numerical_rank(O), "of", A.shape[0])

    Amin, Bmin, Cmin, Dmin, _, _ = exact_minimal_reduction(A, B, C, D)
    print("Amin =\n", Amin)
    print("Bmin =\n", Bmin)
    print("Cmin =\n", Cmin)
    print("Dmin =\n", Dmin)

    for s in [0.2, 1.0, 2.5]:
        g_old = transfer_value(A, B, C, D, s)
        g_new = transfer_value(Amin, Bmin, Cmin, Dmin, s)
        print(f"s={s:3.1f}: original {g_old.ravel()[0]: .8f}, reduced {g_new.ravel()[0]: .8f}")


if __name__ == "__main__":
    main()

The C++ version uses Eigen for matrix and SVD operations. In industrial C++ control software, Eigen is commonly used because it provides efficient dense linear algebra without requiring a large runtime dependency.

Chapter20_Lesson4.cpp

// Chapter20_Lesson4.cpp
// Basic rank-based realization-reduction diagnostics using Eigen.
// Compile example:
//   g++ -std=c++17 Chapter20_Lesson4.cpp -I /path/to/eigen -O2 -o Chapter20_Lesson4

#include <Eigen/Dense>
#include <iostream>
#include <vector>
#include <cmath>

using Eigen::MatrixXd;
using Eigen::VectorXd;

MatrixXd controllabilityMatrix(const MatrixXd& A, const MatrixXd& B) {
    const int n = A.rows();
    MatrixXd R(n, n * B.cols());
    MatrixXd Apow = MatrixXd::Identity(n, n);
    for (int k = 0; k < n; ++k) {
        R.block(0, k * B.cols(), n, B.cols()) = Apow * B;
        Apow = Apow * A;
    }
    return R;
}

MatrixXd observabilityMatrix(const MatrixXd& A, const MatrixXd& C) {
    const int n = A.rows();
    MatrixXd O(n * C.rows(), n);
    MatrixXd Apow = MatrixXd::Identity(n, n);
    for (int k = 0; k < n; ++k) {
        O.block(k * C.rows(), 0, C.rows(), n) = C * Apow;
        Apow = Apow * A;
    }
    return O;
}

int svdRank(const MatrixXd& M, double tol = 1e-10) {
    Eigen::JacobiSVD<MatrixXd> svd(M, Eigen::ComputeFullU | Eigen::ComputeFullV);
    VectorXd s = svd.singularValues();
    if (s.size() == 0) return 0;
    int r = 0;
    for (int i = 0; i < s.size(); ++i) if (s(i) > tol) ++r;
    return r;
}

MatrixXd columnSpaceBasis(const MatrixXd& M) {
    Eigen::JacobiSVD<MatrixXd> svd(M, Eigen::ComputeFullU | Eigen::ComputeFullV);
    VectorXd s = svd.singularValues();
    double tol = 1e-10;
    int r = 0;
    for (int i = 0; i < s.size(); ++i) if (s(i) > tol) ++r;
    return svd.matrixU().leftCols(r);
}

MatrixXd nullSpaceBasis(const MatrixXd& M) {
    Eigen::JacobiSVD<MatrixXd> svd(M, Eigen::ComputeFullU | Eigen::ComputeFullV);
    VectorXd s = svd.singularValues();
    double tol = 1e-10;
    int r = 0;
    for (int i = 0; i < s.size(); ++i) if (s(i) > tol) ++r;
    return svd.matrixV().rightCols(M.cols() - r);
}

int main() {
    MatrixXd A(3, 3);
    A << -1.0, 0.0, 0.0,
          0.0,-2.0, 0.0,
          0.0, 0.0,-3.0;
    MatrixXd B(3, 1);
    B << 1.0, 0.0, 1.0;
    MatrixXd C(1, 3);
    C << 1.0, 0.0, 0.0;
    MatrixXd D(1, 1);
    D << 0.0;

    MatrixXd R = controllabilityMatrix(A, B);
    MatrixXd O = observabilityMatrix(A, C);

    std::cout << "rank(R) = " << svdRank(R) << " of " << A.rows() << "\n";
    std::cout << "rank(O) = " << svdRank(O) << " of " << A.rows() << "\n\n";

    // Stage 1: reachable projection.
    MatrixXd Qr = columnSpaceBasis(R);
    MatrixXd A1 = Qr.transpose() * A * Qr;
    MatrixXd B1 = Qr.transpose() * B;
    MatrixXd C1 = C * Qr;

    // Stage 2: unobservable quotient within the reachable subsystem.
    MatrixXd O1 = observabilityMatrix(A1, C1);
    MatrixXd Qun = nullSpaceBasis(O1);
    int q = Qun.cols();

    MatrixXd Amin, Bmin, Cmin;
    if (q == 0) {
        Amin = A1; Bmin = B1; Cmin = C1;
    } else {
        MatrixXd Qobs = nullSpaceBasis(Qun.transpose());
        MatrixXd T(A1.rows(), A1.rows());
        T << Qun, Qobs;
        MatrixXd Ahat = T.transpose() * A1 * T;
        MatrixXd Bhat = T.transpose() * B1;
        MatrixXd Chat = C1 * T;
        int m = A1.rows() - q;
        Amin = Ahat.block(q, q, m, m);
        Bmin = Bhat.block(q, 0, m, Bhat.cols());
        Cmin = Chat.block(0, q, Chat.rows(), m);
    }

    std::cout << "Amin =\n" << Amin << "\n\n";
    std::cout << "Bmin =\n" << Bmin << "\n\n";
    std::cout << "Cmin =\n" << Cmin << "\n";
    return 0;
}

The Java version is deliberately library-free and demonstrates the rank tests using reduced row-echelon elimination. For production Java numerical work, libraries such as EJML, ojAlgo, or Apache Commons Math should be preferred.

Chapter20_Lesson4.java

// Chapter20_Lesson4.java
// Library-free educational implementation of controllability and observability ranks.
// For production Java control computation, use EJML, ojAlgo, or Apache Commons Math.

import java.util.Arrays;

public class Chapter20_Lesson4 {
    static double[][] multiply(double[][] A, double[][] B) {
        int m = A.length, n = B[0].length, p = B.length;
        double[][] C = new double[m][n];
        for (int i = 0; i < m; i++) {
            for (int k = 0; k < p; k++) {
                for (int j = 0; j < n; j++) C[i][j] += A[i][k] * B[k][j];
            }
        }
        return C;
    }

    static double[][] eye(int n) {
        double[][] I = new double[n][n];
        for (int i = 0; i < n; i++) I[i][i] = 1.0;
        return I;
    }

    static double[][] controllabilityMatrix(double[][] A, double[][] B) {
        int n = A.length, m = B[0].length;
        double[][] R = new double[n][n * m];
        double[][] Apow = eye(n);
        for (int k = 0; k < n; k++) {
            double[][] block = multiply(Apow, B);
            for (int i = 0; i < n; i++) {
                for (int j = 0; j < m; j++) R[i][k * m + j] = block[i][j];
            }
            Apow = multiply(Apow, A);
        }
        return R;
    }

    static double[][] observabilityMatrix(double[][] A, double[][] C) {
        int n = A.length, p = C.length;
        double[][] O = new double[n * p][n];
        double[][] Apow = eye(n);
        for (int k = 0; k < n; k++) {
            double[][] block = multiply(C, Apow);
            for (int i = 0; i < p; i++) {
                for (int j = 0; j < n; j++) O[k * p + i][j] = block[i][j];
            }
            Apow = multiply(Apow, A);
        }
        return O;
    }

    static int rankRref(double[][] M, double tol) {
        int rows = M.length, cols = M[0].length;
        double[][] A = new double[rows][cols];
        for (int i = 0; i < rows; i++) A[i] = Arrays.copyOf(M[i], cols);
        int r = 0;
        for (int c = 0; c < cols && r < rows; c++) {
            int pivot = r;
            for (int i = r + 1; i < rows; i++) {
                if (Math.abs(A[i][c]) > Math.abs(A[pivot][c])) pivot = i;
            }
            if (Math.abs(A[pivot][c]) <= tol) continue;
            double[] tmp = A[r]; A[r] = A[pivot]; A[pivot] = tmp;
            double piv = A[r][c];
            for (int j = c; j < cols; j++) A[r][j] /= piv;
            for (int i = 0; i < rows; i++) {
                if (i == r) continue;
                double factor = A[i][c];
                for (int j = c; j < cols; j++) A[i][j] -= factor * A[r][j];
            }
            r++;
        }
        return r;
    }

    static void printMatrix(String name, double[][] M) {
        System.out.println(name + " =");
        for (double[] row : M) System.out.println(Arrays.toString(row));
        System.out.println();
    }

    public static void main(String[] args) {
        double[][] A = {
            {-1.0, 0.0, 0.0},
            { 0.0,-2.0, 0.0},
            { 0.0, 0.0,-3.0}
        };
        double[][] B = { {1.0}, {0.0}, {1.0} };
        double[][] C = { {1.0, 0.0, 0.0} };

        double[][] R = controllabilityMatrix(A, B);
        double[][] O = observabilityMatrix(A, C);
        printMatrix("Controllability matrix", R);
        printMatrix("Observability matrix", O);

        int n = A.length;
        int rankR = rankRref(R, 1e-10);
        int rankO = rankRref(O, 1e-10);
        System.out.println("rank(R) = " + rankR + " of " + n);
        System.out.println("rank(O) = " + rankO + " of " + n);
        System.out.println("Conclusion: this realization is not minimal.");
        System.out.println("For this diagonal example, the input-output map is G(s)=1/(s+1), so a minimal realization is A=[-1], B=[1], C=[1], D=[0].");
    }
}

MATLAB directly supports state-space objects through ss, reachability and observability matrices through ctrb and obsv, and exact/minimal simplification through minreal. The following script also includes local implementations of the matrices so that the algorithmic idea is visible. In Simulink, the reduced state-space matrices can be placed directly in a State-Space block after verification of the transfer function.

Chapter20_Lesson4.m

% Chapter20_Lesson4.m
% Basic algorithms for exact realization reduction.
% Uses MATLAB functions orth, null, rank. Control System Toolbox functions
% ctrb, obsv, ss, and minreal are also shown for comparison.

clear; clc;

A = diag([-1 -2 -3]);
B = [1; 0; 1];
C = [1 0 0];
D = 0;

R = local_ctrb(A,B);
O = local_obsv(A,C);
fprintf('rank(R) = %d of %d\n', rank(R), size(A,1));
fprintf('rank(O) = %d of %d\n', rank(O), size(A,1));

[Amin,Bmin,Cmin,Dmin] = exact_minimal_reduction(A,B,C,D);
disp('Amin ='); disp(Amin);
disp('Bmin ='); disp(Bmin);
disp('Cmin ='); disp(Cmin);
disp('Dmin ='); disp(Dmin);

% Optional Control System Toolbox comparison:
% sys = ss(A,B,C,D);
% sys_min = minreal(sys);
% disp(sys_min);

function R = local_ctrb(A,B)
    n = size(A,1);
    R = [];
    Apow = eye(n);
    for k = 1:n
        R = [R, Apow*B]; %#ok<AGROW>
        Apow = Apow*A;
    end
end

function O = local_obsv(A,C)
    n = size(A,1);
    O = [];
    Apow = eye(n);
    for k = 1:n
        O = [O; C*Apow]; %#ok<AGROW>
        Apow = Apow*A;
    end
end

function [Amin,Bmin,Cmin,Dmin] = exact_minimal_reduction(A,B,C,D)
    % Stage 1: reachable restriction.
    R = local_ctrb(A,B);
    Qr = orth(R);
    A1 = Qr'*A*Qr;
    B1 = Qr'*B;
    C1 = C*Qr;

    % Stage 2: observable quotient inside reachable part.
    O1 = local_obsv(A1,C1);
    Qun = null(O1);
    q = size(Qun,2);
    if q == 0
        Amin = A1; Bmin = B1; Cmin = C1; Dmin = D;
        return;
    end

    Qobs = null(Qun');
    T = [Qun, Qobs];
    Ahat = T'*A1*T;
    Bhat = T'*B1;
    Chat = C1*T;

    Amin = Ahat(q+1:end, q+1:end);
    Bmin = Bhat(q+1:end, :);
    Cmin = Chat(:, q+1:end);
    Dmin = D;
end

Mathematica is useful for exact symbolic ranks and for checking whether a cancellation is structural rather than merely numerical. The notebook source below constructs reachability and observability matrices, performs the two trimming stages, and compares transfer-function models.

Chapter20_Lesson4.nb

(* Chapter20_Lesson4.nb *)
(* Wolfram Mathematica educational notebook source. Open as a notebook or paste into Mathematica. *)

ClearAll[controllabilityMatrix, observabilityMatrix, exactMinimalReduction];

controllabilityMatrix[A_, B_] := Module[{n = Length[A]},
  ArrayFlatten[{Table[MatrixPower[A, k].B, {k, 0, n - 1}]}]
];

observabilityMatrix[A_, C_] := Module[{n = Length[A]},
  Join @@ Table[C.MatrixPower[A, k], {k, 0, n - 1}]
];

exactMinimalReduction[A_, B_, C_, D_] := Module[
  {R, Qr, A1, B1, C1, O1, Qun, q, Qobs, T, Ahat, Bhat, Chat},
  R = controllabilityMatrix[A, B];
  Qr = Orthogonalize[Transpose[RowReduce[Transpose[R]]]];
  Qr = Select[Qr, Norm[#] > 10^-10 &] // Transpose;
  A1 = Transpose[Qr].A.Qr;
  B1 = Transpose[Qr].B;
  C1 = C.Qr;
  O1 = observabilityMatrix[A1, C1];
  Qun = Transpose[NullSpace[O1]];
  q = If[MatrixQ[Qun], Dimensions[Qun][[2]], 0];
  If[q == 0, Return[{A1, B1, C1, D}]];
  Qobs = Transpose[NullSpace[Transpose[Qun]]];
  T = ArrayFlatten[{ {Qun, Qobs} }];
  Ahat = Transpose[T].A1.T;
  Bhat = Transpose[T].B1;
  Chat = C1.T;
  {Ahat[[q + 1 ;;, q + 1 ;;]], Bhat[[q + 1 ;;, All]], Chat[[All, q + 1 ;;]], D}
];

A = DiagonalMatrix[{-1, -2, -3}];
B = { {1}, {0}, {1} };
Cmat = { {1, 0, 0} };
Dmat = { {0} };

R = controllabilityMatrix[A, B];
O = observabilityMatrix[A, Cmat];
Print["rank(R) = ", MatrixRank[R], " of ", Length[A]];
Print["rank(O) = ", MatrixRank[O], " of ", Length[A]];

{Amin, Bmin, Cmin, Dmin} = exactMinimalReduction[A, B, Cmat, Dmat];
Print["Amin = ", MatrixForm[Amin]];
Print["Bmin = ", MatrixForm[Bmin]];
Print["Cmin = ", MatrixForm[Cmin]];
Print["Dmin = ", MatrixForm[Dmin]];

TransferFunctionModel[StateSpaceModel[{A, B, Cmat, Dmat}], s]
TransferFunctionModel[StateSpaceModel[{Amin, Bmin, Cmin, Dmin}], s]

Problem 1 (Reachability Trimming): Let \( \mathbf{A}=\operatorname{diag}(-1,-2) \), \( \mathbf{B}=\begin{bmatrix}1&0\end{bmatrix}^T \), \( \mathbf{C}=\begin{bmatrix}1&1\end{bmatrix} \), and \( \mathbf{D}=0 \). Determine whether the realization is reachable and construct a reduced realization.

Solution: The reachability matrix is

\[ \mathcal{R}=\begin{bmatrix}\mathbf{B}&\mathbf{AB}\end{bmatrix} =\begin{bmatrix}1&-1\\0&0\end{bmatrix},\qquad \operatorname{rank}(\mathcal{R})=1. \]

Only the first state is reachable. Taking \( \mathbf{T}_r=\begin{bmatrix}1&0\end{bmatrix}^T \) gives

\[ A_r=-1, \qquad B_r=1, \qquad C_r=1, \qquad D_r=0. \]

The second state affects neither the response from zero initial condition nor the transfer function.

Problem 2 (Observability Trimming): Let \( \mathbf{A}=\operatorname{diag}(-1,-3) \), \( \mathbf{B}=\begin{bmatrix}1&1\end{bmatrix}^T \), \( \mathbf{C}=\begin{bmatrix}1&0\end{bmatrix} \), and \( \mathbf{D}=0 \). Is the realization minimal?

Solution: The reachability matrix is

\[ \mathcal{R}=\begin{bmatrix}1&-1\\1&-3\end{bmatrix}, \qquad \det(\mathcal{R})=-2, \]

so the realization is reachable. The observability matrix is

\[ \mathcal{O}=\begin{bmatrix}1&0\\-1&0\end{bmatrix}, \qquad \operatorname{rank}(\mathcal{O})=1. \]

Hence the second state is unobservable and the realization is not minimal. The transfer function is \( G(s)=1/(s+1) \), so a minimal realization is \( A_{min}=-1, B_{min}=1, C_{min}=1, D_{min}=0 \).

Problem 3 (Proof by Markov Parameters): Suppose \( \operatorname{range}(\mathbf{T}_r)=\mathscr{R} \) and \( \mathbf{B}=\mathbf{T}_r\mathbf{B}_r \), \( \mathbf{A}\mathbf{T}_r=\mathbf{T}_r\mathbf{A}_r \). Prove that reachability trimming preserves the transfer matrix.

Solution: For \( k=0 \), \( \mathbf{A}^0\mathbf{B}=\mathbf{B}=\mathbf{T}_r\mathbf{B}_r \). Assume \( \mathbf{A}^k\mathbf{B}=\mathbf{T}_r\mathbf{A}_r^k\mathbf{B}_r \). Then

\[ \mathbf{A}^{k+1}\mathbf{B} =\mathbf{A}\mathbf{T}_r\mathbf{A}_r^k\mathbf{B}_r =\mathbf{T}_r\mathbf{A}_r^{k+1}\mathbf{B}_r. \]

Therefore \( \mathbf{C}\mathbf{A}^k\mathbf{B}=\mathbf{C}_r\mathbf{A}_r^k\mathbf{B}_r \) for all \( k \). Since the direct term \( \mathbf{D} \) is unchanged, the transfer matrices are equal.

Problem 4 (Hankel Rank): A SISO discrete-time strictly proper system has Markov parameters \( H_k=2^{k-1}+3^{k-1} \) for \( k=1,2,\dots \). Use the Hankel idea to determine the minimal order.

Solution: Build a two-row Hankel matrix:

\[ \mathcal{H}_{2,2}=\begin{bmatrix} H_1&H_2\\H_2&H_3 \end{bmatrix} =\begin{bmatrix}2&5\\5&13\end{bmatrix}, \qquad \det(\mathcal{H}_{2,2})=26-25=1. \]

Thus \( \operatorname{rank}(\mathcal{H}_{2,2})=2 \), so the minimal order is at least 2. Since the sequence is the sum of two independent exponential modes, the minimal order is exactly 2.

Problem 5 (Dimension Formula): A realization has \( \dim(\mathscr{R})=4 \) and \( \dim(\mathscr{R}\cap\mathscr{N})=1 \). What is the order of the minimal realization obtained by exact trimming?

Solution: Exact trimming first restricts the system to the 4-dimensional reachable subspace. It then removes the 1-dimensional component that is reachable but unobservable. Hence

\[ n_{min}=\dim(\mathscr{R})-\dim(\mathscr{R}\cap\mathscr{N})=4-1=3. \]

Exact realization reduction is based on two algebraic eliminations: remove states outside the reachable subspace and remove states inside the unobservable subspace of the reachable part. Rank-revealing bases obtained from reachability and observability matrices produce the reduced coordinates. The Ho-Kalman/Hankel viewpoint reaches the same conclusion from input-output Markov parameters: the rank of the Hankel matrix reveals the minimal dynamic order. In exact arithmetic, these algorithms preserve the transfer matrix exactly; in floating-point arithmetic, rank tolerance and conditioning become central.