Consider a continuous-time finite-dimensional LTI realization \( \Sigma=(A,B,C,D) \):

\[ \dot{x}(t)=Ax(t)+Bu(t),\qquad y(t)=Cx(t)+Du(t), \]

where \( x(t)\in\mathbb{R}^n \), \( u(t)\in\mathbb{R}^m \), and \( y(t)\in\mathbb{R}^p \). The associated transfer matrix is

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

A realization is minimal when no lower-dimensional realization has the same transfer matrix. In this lesson, minimality is identified through the decomposition tools introduced in Lessons 1--3: reachable states, unobservable states, and the Kalman block decomposition.

The main theorem is: a finite-dimensional realization is minimal if and only if it is controllable/reachable and observable. Therefore, the practical task is to isolate the controllable-observable part of the model and discard the rest without changing \( G(s) \).

flowchart TD
  A["Given realization (A,B,C,D)"] --> B["Compute reachable subspace"]
  B --> C["Remove unreachable coordinates"]
  C --> D["Compute unobservable subspace of reachable part"]
  D --> E["Remove unobservable reachable coordinates"]
  E --> F["Minimal realization (Amin,Bmin,Cmin,D)"]
  F --> G["Same transfer matrix G(s), smaller state dimension"]
        

The reachable subspace is generated by the columns of the controllability matrix

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

Its dimension is \( r_c=\operatorname{rank}\mathcal{C}(A,B) \). If \( r_c=n \), every state direction is reachable. Otherwise, the model contains input-inaccessible internal modes.

The unobservable subspace is the null space of the observability matrix

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

Its codimension is \( r_o=\operatorname{rank}\mathcal{O}(A,C) \). If \( r_o=n \), every state direction affects the measured output over some time interval. If \( r_o<n \), there are internal directions that are invisible at the output.

Both subspaces are invariant in the correct sense:

\[ A\mathcal{R}\subseteq\mathcal{R},\qquad A\mathcal{N}\subseteq\mathcal{N}. \]

The first inclusion follows because multiplication by \( A \) maps \( A^kB \) into \( A^{k+1}B \), and higher powers reduce by the Cayley-Hamilton theorem. The second follows because if \( x\in\ker\mathcal{O} \), then \( CA^k x=0 \) for all \( k=0,\dots,n-1 \), hence \( CA^k(Ax)=CA^{k+1}x=0 \) for the required powers, again using Cayley-Hamilton for closure.

The full Kalman decomposition separates the state space into four conceptual components:

\[ \mathbb{R}^n=\mathcal{X}_{co}\oplus\mathcal{X}_{c\bar o} \oplus\mathcal{X}_{\bar c o}\oplus\mathcal{X}_{\bar c\bar o}, \]

where the subscript \( c \) means controllable, \( o \) means observable, and a bar denotes the absence of that property. Only \( \mathcal{X}_{co} \) contributes to the zero-state transfer matrix.

flowchart LR
  U["Input u"] --> CO["controllable + observable"]
  U --> CU["controllable + unobservable"]
  U -. no access .-> UO["uncontrollable + observable"]
  U -. no access .-> UU["uncontrollable + unobservable"]
  CO --> Y["Output y"]
  CU -. invisible .-> Y
  UO --> Y
  UU -. invisible .-> Y
  CO --> MIN["minimal realization"]
        

In a basis compatible with the decomposition, the matrices acquire a block form. The exact ordering of blocks depends on the chosen basis, but the transfer matrix depends only on the controllable-observable block:

\[ G(s)=C_{co}(sI-A_{co})^{-1}B_{co}+D. \]

The other blocks may influence nonzero-initial-condition responses, but they do not define a smaller or larger zero-state input-output transfer matrix. This is why a minimal realization is an input-output model, not necessarily a complete physical model of every internal mode.

A numerically convenient way to identify a minimal realization is to use two reductions. First remove unreachable states; then remove unobservable states of the reachable subsystem.

Step 1: Reachable reduction. Let \( T_c=\begin{bmatrix}T_r&T_u\end{bmatrix} \), where the columns of \( T_r \) form a basis of \( \mathcal{R} \). In coordinates \( x=T_c z \),

\[ \bar{A}=T_c^{-1}AT_c=\begin{bmatrix}A_r&A_{ru}\\0&A_u\end{bmatrix},\qquad \bar{B}=T_c^{-1}B=\begin{bmatrix}B_r\\0\end{bmatrix},\qquad \bar{C}=CT_c=\begin{bmatrix}C_r&C_u\end{bmatrix}. \]

Since the lower block receives no input from the zero initial state, the transfer matrix is unchanged by retaining only \( (A_r,B_r,C_r,D) \).

Step 2: Observable reduction of the reachable part. For \( (A_r,B_r,C_r,D) \), compute \( \mathcal{N}_r=\ker\mathcal{O}(A_r,C_r) \). Let \( T_o=\begin{bmatrix}T_n&T_q\end{bmatrix} \), where \( T_n \) spans the unobservable subspace. In coordinates \( z_r=T_o \xi \),

\[ \tilde{A}=T_o^{-1}A_rT_o=\begin{bmatrix}A_n&A_{nq}\\0&A_m\end{bmatrix},\qquad \tilde{B}=T_o^{-1}B_r=\begin{bmatrix}B_n\\B_m\end{bmatrix},\qquad \tilde{C}=C_rT_o=\begin{bmatrix}0&C_m\end{bmatrix}. \]

The subsystem \( (A_n,B_n,0,D) \) is invisible at the output and cannot influence the observable quotient dynamics because of the lower-left zero block. Hence the minimal realization is

\[ \Sigma_{min}=(A_m,B_m,C_m,D). \]

The proof relies on block inversion and the zero-state assumption. For the reachable reduction, the transformed matrices have the form

\[ \bar{A}=\begin{bmatrix}A_r&A_{ru}\\0&A_u\end{bmatrix},\quad \bar{B}=\begin{bmatrix}B_r\\0\end{bmatrix},\quad \bar{C}=\begin{bmatrix}C_r&C_u\end{bmatrix}. \]

Since

\[ sI-\bar{A}=\begin{bmatrix}sI-A_r&-A_{ru}\\0&sI-A_u\end{bmatrix}, \]

solving \( (sI-\bar{A})v=\bar{B} \) gives \( v_2=0 \) and \( v_1=(sI-A_r)^{-1}B_r \). Thus

\[ \bar{C}(sI-\bar{A})^{-1}\bar{B}+D =C_r(sI-A_r)^{-1}B_r+D. \]

For the observable reduction, use

\[ \tilde{A}=\begin{bmatrix}A_n&A_{nq}\\0&A_m\end{bmatrix},\quad \tilde{B}=\begin{bmatrix}B_n\\B_m\end{bmatrix},\quad \tilde{C}=\begin{bmatrix}0&C_m\end{bmatrix}. \]

Solving \( (sI-\tilde{A})w=\tilde{B} \) gives \( w_2=(sI-A_m)^{-1}B_m \). The output ignores \( w_1 \), so

\[ \tilde{C}(sI-\tilde{A})^{-1}\tilde{B}+D =C_m(sI-A_m)^{-1}B_m+D. \]

Therefore the two-step decomposition produces a realization of the same transfer matrix. Since the remaining realization is both reachable and observable, it is minimal.

Consider the four-state realization

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

The first two states are both controllable and observable. The third state is controllable but unobservable. The fourth state is observable from a nonzero initial condition but uncontrollable from the input.

The controllability matrix is

\[ \mathcal{C}=\begin{bmatrix} 0&1&-3&7\\1&-3&7&-15\\1&-4&16&-64\\0&0&0&0 \end{bmatrix},\qquad \operatorname{rank}\mathcal{C}=3. \]

The transfer function is

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

Thus a minimal realization is

\[ A_m=\begin{bmatrix}0&1\\-2&-3\end{bmatrix},\qquad B_m=\begin{bmatrix}0\\1\end{bmatrix},\qquad C_m=\begin{bmatrix}1&0\end{bmatrix},\qquad D_m=0. \]

The removed modes \( -4 \) and \( -5 \) are internal nonminimal modes. They may appear in the state trajectory under certain initial conditions or internal coordinates, but not in the zero-state input-output transfer function.

In exact mathematics, ranks are unambiguous. In numerical computation, the singular values of \( \mathcal{C} \) and \( \mathcal{O} \) determine whether a direction is kept or discarded. If

\[ \sigma_1\geq\sigma_2\geq\cdots\geq\sigma_q\geq0, \]

a practical rank rule is

\[ r=\#\left\{i:\sigma_i>\epsilon\max(n,m)\sigma_1\right\}. \]

Small nonzero singular values may represent weakly controllable or weakly observable states rather than exactly removable states. In such cases, the exact Kalman decomposition should be treated as a structural diagnostic, while balanced truncation or model reduction methods should be used for approximate reduction. Those approximate methods are not the topic of this lesson, but this distinction is crucial in engineering computation.

This Python implementation uses NumPy and SciPy. The decomposition is performed from scratch using controllability/observability matrices and singular-value decompositions.

# Chapter19_Lesson4.py
# Identification of a minimal realization via sequential Kalman decomposition.
# Requires: numpy, scipy

import numpy as np
from scipy.linalg import block_diag


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


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


def svd_rank(M, tol=1e-10):
    s = np.linalg.svd(M, compute_uv=False)
    if s.size == 0:
        return 0
    return int(np.sum(s > tol * max(M.shape) * max(s[0], 1.0)))


def kalman_minimal_realization(A, B, C, D, tol=1e-10):
    """
    Reduce a continuous-time LTI realization to a minimal one by:
      1. retaining only the reachable subspace,
      2. quotienting/removing the unobservable subspace of that reachable part.

    The construction uses orthonormal bases from SVD, so transformations are
    similarity transformations with T^{-1} = T.T whenever T is square.
    """
    n = A.shape[0]

    # Step 1: reachable reduction.
    Wc = controllability_matrix(A, B)
    Uc, sc, _ = np.linalg.svd(Wc, full_matrices=True)
    rc = svd_rank(Wc, tol)
    Tc = Uc  # first rc columns span reachable subspace
    Ac = Tc.T @ A @ Tc
    Bc = Tc.T @ B
    Cc = C @ Tc
    Ar = Ac[:rc, :rc]
    Br = Bc[:rc, :]
    Cr = Cc[:, :rc]

    # Step 2: observable reduction of the reachable subsystem.
    Wo = observability_matrix(Ar, Cr)
    Uo, so, Vho = np.linalg.svd(Wo, full_matrices=True)
    ro = svd_rank(Wo, tol)
    V = Vho.T

    # V[:, ro:] spans the unobservable subspace; V[:, :ro] is an orthonormal
    # complement. Put unobservable coordinates first, observable coordinates last.
    To = np.hstack([V[:, ro:], V[:, :ro]]) if ro < rc else V[:, :ro]
    Ao = To.T @ Ar @ To
    Bo = To.T @ Br
    Co = Cr @ To

    n_unobs = rc - ro
    Amin = Ao[n_unobs:, n_unobs:]
    Bmin = Bo[n_unobs:, :]
    Cmin = Co[:, n_unobs:]

    return {
        "Amin": Amin,
        "Bmin": Bmin,
        "Cmin": Cmin,
        "Dmin": D.copy(),
        "rank_controllability": rc,
        "rank_observability_after_reachable": ro,
        "Tc": Tc,
        "To": To,
        "reachable_realization": (Ar, Br, Cr, D.copy()),
    }


def transfer_value(A, B, C, D, s):
    n = A.shape[0]
    return C @ np.linalg.solve(s * np.eye(n) - A, B) + D


def demo():
    np.set_printoptions(precision=6, suppress=True)

    # Transparent nonminimal realization: first two states are controllable and observable.
    Aco = np.array([[0.0, 1.0], [-2.0, -3.0]])
    Bco = np.array([[0.0], [1.0]])
    Cco = np.array([[1.0, 0.0]])

    # Hidden modes: controllable-unobservable and uncontrollable-observable.
    A0 = block_diag(Aco, np.array([[-4.0]]), np.array([[-5.0]]))
    B0 = np.array([[0.0], [1.0], [1.0], [0.0]])
    C0 = np.array([[1.0, 0.0, 0.0, 2.0]])
    D0 = np.array([[0.0]])

    # Hide the structure with an orthogonal similarity transformation.
    rng = np.random.default_rng(19)
    Q, _ = np.linalg.qr(rng.normal(size=(4, 4)))
    A = Q @ A0 @ Q.T
    B = Q @ B0
    C = C0 @ Q.T
    D = D0

    result = kalman_minimal_realization(A, B, C, D)
    Amin, Bmin, Cmin, Dmin = result["Amin"], result["Bmin"], result["Cmin"], result["Dmin"]

    print("rank controllability =", result["rank_controllability"])
    print("rank observability after reachable reduction =", result["rank_observability_after_reachable"])
    print("Amin =\n", Amin)
    print("Bmin =\n", Bmin)
    print("Cmin =\n", Cmin)

    for s in [0.0, 1.0, 2.0, 1.0 + 2.0j]:
        G_full = transfer_value(A, B, C, D, s)
        G_min = transfer_value(Amin, Bmin, Cmin, Dmin, s)
        print(f"s={s:>8}: G_full={G_full[0,0]: .8f}, G_min={G_min[0,0]: .8f}, error={abs(G_full[0,0]-G_min[0,0]):.2e}")


if __name__ == "__main__":
    demo()
      

This C++ implementation uses Eigen for matrix algebra and SVD. It is suitable for students who want to see the matrix construction explicitly rather than calling a high-level control toolbox.

// Chapter19_Lesson4.cpp
// Identification of a minimal realization via sequential Kalman decomposition.
// Requires Eigen 3: g++ -std=c++17 Chapter19_Lesson4.cpp -I /path/to/eigen -O2

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

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

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

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

int svdRank(const MatrixXd& M, double tol = 1e-10) {
    JacobiSVD<MatrixXd> svd(M);
    const VectorXd s = svd.singularValues();
    if (s.size() == 0) return 0;
    const double threshold = tol * std::max(M.rows(), M.cols()) * std::max(s(0), 1.0);
    int r = 0;
    for (int i = 0; i < s.size(); ++i) {
        if (s(i) > threshold) ++r;
    }
    return r;
}

struct MinimalResult {
    MatrixXd A, B, C, D;
    int reachableRank;
    int observableRankAfterReachable;
};

MinimalResult kalmanMinimalRealization(const MatrixXd& A, const MatrixXd& B,
                                       const MatrixXd& C, const MatrixXd& D) {
    const int n = A.rows();

    MatrixXd Wc = controllabilityMatrix(A, B);
    JacobiSVD<MatrixXd> svdc(Wc, Eigen::ComputeFullU | Eigen::ComputeFullV);
    MatrixXd Tc = svdc.matrixU();
    int rc = svdRank(Wc);

    MatrixXd Ac = Tc.transpose() * A * Tc;
    MatrixXd Bc = Tc.transpose() * B;
    MatrixXd Cc = C * Tc;

    MatrixXd Ar = Ac.block(0, 0, rc, rc);
    MatrixXd Br = Bc.block(0, 0, rc, B.cols());
    MatrixXd Cr = Cc.block(0, 0, C.rows(), rc);

    MatrixXd Wo = observabilityMatrix(Ar, Cr);
    JacobiSVD<MatrixXd> svdo(Wo, Eigen::ComputeFullU | Eigen::ComputeFullV);
    MatrixXd V = svdo.matrixV();
    int ro = svdRank(Wo);

    MatrixXd To(rc, rc);
    if (ro < rc) {
        To << V.rightCols(rc - ro), V.leftCols(ro);
    } else {
        To = V.leftCols(ro);
    }

    MatrixXd Ao = To.transpose() * Ar * To;
    MatrixXd Bo = To.transpose() * Br;
    MatrixXd Co = Cr * To;

    int nUnobs = rc - ro;
    MinimalResult result;
    result.A = Ao.block(nUnobs, nUnobs, ro, ro);
    result.B = Bo.block(nUnobs, 0, ro, B.cols());
    result.C = Co.block(0, nUnobs, C.rows(), ro);
    result.D = D;
    result.reachableRank = rc;
    result.observableRankAfterReachable = ro;
    return result;
}

std::complex<double> transferValue(const MatrixXd& A, const MatrixXd& B,
                                   const MatrixXd& C, const MatrixXd& D,
                                   std::complex<double> s) {
    using MatrixXcd = Eigen::MatrixXcd;
    const int n = A.rows();
    MatrixXcd Ac = A.cast<std::complex<double>>();
    MatrixXcd Bc = B.cast<std::complex<double>>();
    MatrixXcd Cc = C.cast<std::complex<double>>();
    MatrixXcd Dc = D.cast<std::complex<double>>();
    MatrixXcd SIminusA = s * MatrixXcd::Identity(n, n) - Ac;
    MatrixXcd G = Cc * SIminusA.fullPivLu().solve(Bc) + Dc;
    return G(0, 0);
}

int main() {
    MatrixXd A(4, 4);
    A << 0.0, 1.0, 0.0, 0.0,
        -2.0, -3.0, 0.0, 0.0,
         0.0, 0.0, -4.0, 0.0,
         0.0, 0.0, 0.0, -5.0;

    MatrixXd B(4, 1);
    B << 0.0, 1.0, 1.0, 0.0;

    MatrixXd C(1, 4);
    C << 1.0, 0.0, 0.0, 2.0;

    MatrixXd D(1, 1);
    D << 0.0;

    MinimalResult minsys = kalmanMinimalRealization(A, B, C, D);

    std::cout << "rank controllability = " << minsys.reachableRank << "\n";
    std::cout << "rank observability after reachable reduction = "
              << minsys.observableRankAfterReachable << "\n";
    std::cout << "Amin:\n" << minsys.A << "\n";
    std::cout << "Bmin:\n" << minsys.B << "\n";
    std::cout << "Cmin:\n" << minsys.C << "\n";

    std::vector<std::complex<double>> testPoints = { {0.0, 0.0}, {1.0, 0.0}, {2.0, 0.0}, {1.0, 2.0} };
    for (const auto& s : testPoints) {
        auto gf = transferValue(A, B, C, D, s);
        auto gm = transferValue(minsys.A, minsys.B, minsys.C, minsys.D, s);
        std::cout << "s = " << s << ", Gfull = " << gf
                  << ", Gmin = " << gm
                  << ", error = " << std::abs(gf - gm) << "\n";
    }

    return 0;
}
      

Java has no standard built-in numerical linear algebra library, so this example uses Apache Commons Math. EJML or ojAlgo can be used similarly.

// Chapter19_Lesson4.java
// Identification of a minimal realization via sequential Kalman decomposition.
// Requires Apache Commons Math 3.6.1 on the classpath.
// Compile: javac -cp commons-math3-3.6.1.jar Chapter19_Lesson4.java
// Run:     java  -cp .:commons-math3-3.6.1.jar Chapter19_Lesson4

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.SingularValueDecomposition;

public class Chapter19_Lesson4 {
    static RealMatrix controllabilityMatrix(RealMatrix A, RealMatrix B) {
        int n = A.getRowDimension();
        int m = B.getColumnDimension();
        RealMatrix W = new Array2DRowRealMatrix(n, n * m);
        RealMatrix Ak = MatrixUtils.createRealIdentityMatrix(n);
        for (int k = 0; k < n; k++) {
            W.setSubMatrix(Ak.multiply(B).getData(), 0, k * m);
            Ak = Ak.multiply(A);
        }
        return W;
    }

    static RealMatrix observabilityMatrix(RealMatrix A, RealMatrix C) {
        int n = A.getRowDimension();
        int p = C.getRowDimension();
        RealMatrix W = new Array2DRowRealMatrix(n * p, n);
        RealMatrix Ak = MatrixUtils.createRealIdentityMatrix(n);
        for (int k = 0; k < n; k++) {
            W.setSubMatrix(C.multiply(Ak).getData(), k * p, 0);
            Ak = Ak.multiply(A);
        }
        return W;
    }

    static RealMatrix columns(RealMatrix M, int firstInclusive, int lastExclusive) {
        int rows = M.getRowDimension();
        int cols = lastExclusive - firstInclusive;
        RealMatrix out = new Array2DRowRealMatrix(rows, cols);
        for (int j = 0; j < cols; j++) {
            out.setColumnVector(j, M.getColumnVector(firstInclusive + j));
        }
        return out;
    }

    static RealMatrix hstack(RealMatrix A, RealMatrix B) {
        RealMatrix out = new Array2DRowRealMatrix(A.getRowDimension(), A.getColumnDimension() + B.getColumnDimension());
        out.setSubMatrix(A.getData(), 0, 0);
        out.setSubMatrix(B.getData(), 0, A.getColumnDimension());
        return out;
    }

    static RealMatrix sub(RealMatrix M, int r0, int r1, int c0, int c1) {
        return M.getSubMatrix(r0, r1 - 1, c0, c1 - 1);
    }

    static class MinimalResult {
        RealMatrix A, B, C, D;
        int reachableRank;
        int observableRankAfterReachable;
    }

    static MinimalResult kalmanMinimalRealization(RealMatrix A, RealMatrix B, RealMatrix C, RealMatrix D) {
        RealMatrix Wc = controllabilityMatrix(A, B);
        SingularValueDecomposition svdc = new SingularValueDecomposition(Wc);
        int rc = svdc.getRank();
        RealMatrix Tc = svdc.getU();

        RealMatrix Ac = Tc.transpose().multiply(A).multiply(Tc);
        RealMatrix Bc = Tc.transpose().multiply(B);
        RealMatrix Cc = C.multiply(Tc);

        RealMatrix Ar = sub(Ac, 0, rc, 0, rc);
        RealMatrix Br = sub(Bc, 0, rc, 0, B.getColumnDimension());
        RealMatrix Cr = sub(Cc, 0, C.getRowDimension(), 0, rc);

        RealMatrix Wo = observabilityMatrix(Ar, Cr);
        SingularValueDecomposition svdo = new SingularValueDecomposition(Wo);
        int ro = svdo.getRank();
        RealMatrix V = svdo.getV();

        RealMatrix To;
        if (ro < rc) {
            RealMatrix unobservable = columns(V, ro, rc);
            RealMatrix observableComplement = columns(V, 0, ro);
            To = hstack(unobservable, observableComplement);
        } else {
            To = columns(V, 0, ro);
        }

        RealMatrix Ao = To.transpose().multiply(Ar).multiply(To);
        RealMatrix Bo = To.transpose().multiply(Br);
        RealMatrix Co = Cr.multiply(To);

        int nUnobs = rc - ro;
        MinimalResult result = new MinimalResult();
        result.A = sub(Ao, nUnobs, nUnobs + ro, nUnobs, nUnobs + ro);
        result.B = sub(Bo, nUnobs, nUnobs + ro, 0, B.getColumnDimension());
        result.C = sub(Co, 0, C.getRowDimension(), nUnobs, nUnobs + ro);
        result.D = D.copy();
        result.reachableRank = rc;
        result.observableRankAfterReachable = ro;
        return result;
    }

    static double transferAtRealS(RealMatrix A, RealMatrix B, RealMatrix C, RealMatrix D, double s) {
        int n = A.getRowDimension();
        RealMatrix sIminusA = MatrixUtils.createRealIdentityMatrix(n).scalarMultiply(s).subtract(A);
        RealMatrix X = new LUDecomposition(sIminusA).getSolver().solve(B);
        return C.multiply(X).add(D).getEntry(0, 0);
    }

    static void printMatrix(String name, RealMatrix M) {
        System.out.println(name + ":");
        for (int i = 0; i < M.getRowDimension(); i++) {
            for (int j = 0; j < M.getColumnDimension(); j++) {
                System.out.printf("% .6f ", M.getEntry(i, j));
            }
            System.out.println();
        }
    }

    public static void main(String[] args) {
        RealMatrix A = MatrixUtils.createRealMatrix(new double[][]{
            {0.0, 1.0, 0.0, 0.0},
            {-2.0, -3.0, 0.0, 0.0},
            {0.0, 0.0, -4.0, 0.0},
            {0.0, 0.0, 0.0, -5.0}
        });
        RealMatrix B = MatrixUtils.createColumnRealMatrix(new double[]{0.0, 1.0, 1.0, 0.0});
        RealMatrix C = MatrixUtils.createRowRealMatrix(new double[]{1.0, 0.0, 0.0, 2.0});
        RealMatrix D = MatrixUtils.createRealMatrix(new double[][]{ {0.0} });

        MinimalResult minsys = kalmanMinimalRealization(A, B, C, D);
        System.out.println("rank controllability = " + minsys.reachableRank);
        System.out.println("rank observability after reachable reduction = " + minsys.observableRankAfterReachable);
        printMatrix("Amin", minsys.A);
        printMatrix("Bmin", minsys.B);
        printMatrix("Cmin", minsys.C);

        double[] sValues = {0.0, 1.0, 2.0};
        for (double s : sValues) {
            double gf = transferAtRealS(A, B, C, D, s);
            double gm = transferAtRealS(minsys.A, minsys.B, minsys.C, minsys.D, s);
            System.out.printf("s = %.1f, Gfull = %.8f, Gmin = %.8f, error = %.2e%n", s, gf, gm, Math.abs(gf - gm));
        }
    }
}
      

MATLAB provides ctrb, obsv, ss, tf, and minreal. The script below still shows the decomposition steps explicitly. In Simulink, the full and minimal systems can be compared by placing their matrices in two State-Space blocks driven by the same input.

% Chapter19_Lesson4.m
% Identification of a minimal realization via sequential Kalman decomposition.
% Requires Control System Toolbox for ctrb, obsv, ss, tf.

clear; clc;

% Transparent nonminimal realization.
Aco = [0 1; -2 -3];
Bco = [0; 1];
Cco = [1 0];

A = blkdiag(Aco, -4, -5);
B = [0; 1; 1; 0];
C = [1 0 0 2];
D = 0;

[Amin, Bmin, Cmin, Dmin, info] = kalman_minimal_decomposition(A, B, C, D);

disp('rank controllability ='); disp(info.reachable_rank);
disp('rank observability after reachable reduction ='); disp(info.observable_rank_after_reachable);
disp('Amin ='); disp(Amin);
disp('Bmin ='); disp(Bmin);
disp('Cmin ='); disp(Cmin);

sys_full = ss(A, B, C, D);
sys_min  = ss(Amin, Bmin, Cmin, Dmin);

disp('Full transfer function:');
tf(sys_full)
disp('Minimal transfer function:');
tf(sys_min)

% Optional comparison with MATLAB minreal.
disp('MATLAB minreal result:');
minreal(sys_full)

% Simulink note:
% The same matrices can be placed in a State-Space block. The reduced matrices
% Amin, Bmin, Cmin, Dmin define an equivalent State-Space block with fewer states.

function [Amin, Bmin, Cmin, Dmin, info] = kalman_minimal_decomposition(A, B, C, D)
    n = size(A, 1);

    % Step 1: reachable reduction.
    Wc = ctrb(A, B);
    R = orth(Wc);
    rc = size(R, 2);
    N = null(R');
    Tc = [R N];

    Ac = Tc' * A * Tc;
    Bc = Tc' * B;
    Cc = C * Tc;

    Ar = Ac(1:rc, 1:rc);
    Br = Bc(1:rc, :);
    Cr = Cc(:, 1:rc);

    % Step 2: observable reduction of the reachable subsystem.
    Wo = obsv(Ar, Cr);
    ro = rank(Wo);
    Nobs = null(Wo);
    Ocomp = orth(Wo');
    To = [Nobs Ocomp];

    Ao = To' * Ar * To;
    Bo = To' * Br;
    Co = Cr * To;

    n_unobs = rc - ro;
    Amin = Ao(n_unobs+1:end, n_unobs+1:end);
    Bmin = Bo(n_unobs+1:end, :);
    Cmin = Co(:, n_unobs+1:end);
    Dmin = D;

    info.reachable_rank = rc;
    info.observable_rank_after_reachable = ro;
    info.T_reachable = Tc;
    info.T_observable = To;
    info.A_reachable = Ar;
    info.B_reachable = Br;
    info.C_reachable = Cr;
end
      

Mathematica is especially useful for symbolic verification of transfer equivalence. The notebook code below computes ranks and proves that the full and reduced transfer functions are identical for the worked example.

Notebook[{
Cell["Chapter 19 - Lesson 4: Minimal Realization via Decomposition", "Title"],
Cell["This notebook computes controllability and observability matrices, compares full and minimal transfer functions, and verifies transfer equivalence for the worked example.", "Text"],
Cell["ClearAll[controllabilityMatrix, observabilityMatrix, transferFunction];

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

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

transferFunction[A_, B_, C_, D_, s_] :=
  Simplify[C.Inverse[s IdentityMatrix[Length[A]] - A].B + D];

Aco = { {0, 1}, {-2, -3} };
Bco = { {0}, {1} };
Cco = { {1, 0} };

A = ArrayFlatten[{ {Aco, ConstantArray[0, {2, 1}], ConstantArray[0, {2, 1}]},
                  {ConstantArray[0, {1, 2}], { {-4} }, { {0} } },
                  {ConstantArray[0, {1, 2}], { {0} }, { {-5} } } }];
B = { {0}, {1}, {1}, {0} };
C = { {1, 0, 0, 2} };
D = { {0} };

Wc = controllabilityMatrix[A, B];
Wo = observabilityMatrix[A, C];

{MatrixRank[Wc], MatrixRank[Wo]}

Gfull = transferFunction[A, B, C, D, s]

Amin = A[[1 ;; 2, 1 ;; 2]];
Bmin = B[[1 ;; 2, All]];
Cmin = C[[All, 1 ;; 2]];
Dmin = D;

Gmin = transferFunction[Amin, Bmin, Cmin, Dmin, s]

FullSimplify[Gfull == Gmin]
", "Input"]
},
WindowSize -> {900, 700},
StyleDefinitions -> "Default.nb"]
      

Problem 1 (Reachable Reduction): Let \( \mathcal{R}=\operatorname{im}\mathcal{C}(A,B) \) have dimension \( r_c<n \). Show that in a basis whose first \( r_c \) vectors span \( \mathcal{R} \), the transformed input matrix has the form \( \bar{B}=\begin{bmatrix}B_r^T&0^T\end{bmatrix}^T \) and the lower-left block of \( \bar{A} \) is zero.

Solution: Since every column of \( B \) belongs to \( \mathcal{R} \), it has zero coordinate component in any complement of \( \mathcal{R} \). Therefore the lower block of \( \bar{B} \) is zero. Also, \( A\mathcal{R}\subseteq\mathcal{R} \), so if the state starts in \( \mathcal{R} \), multiplication by \( A \) cannot generate a component outside \( \mathcal{R} \). Hence the block mapping reachable coordinates into unreachable coordinates is zero.

Problem 2 (Observable Quotient): Suppose \( \mathcal{N}=\ker\mathcal{O}(A,C) \). Prove that in a basis whose first vectors span \( \mathcal{N} \), the output matrix has the form \( \bar{C}=\begin{bmatrix}0&C_q\end{bmatrix} \).

Solution: If \( x\in\mathcal{N} \), then the first block row of the observability equations gives \( Cx=0 \). Thus every vector in the unobservable basis is annihilated by \( C \). Therefore the columns of \( \bar{C} \) associated with the unobservable coordinates are zero.

Problem 3 (Minimal Order of a Given System): For

\[ A=\begin{bmatrix}0&1&0\\-6&-5&0\\0&0&-8\end{bmatrix},\quad B=\begin{bmatrix}0\\1\\1\end{bmatrix},\quad C=\begin{bmatrix}1&0&0\end{bmatrix},\quad D=0, \]

find a minimal realization.

Solution: The third state is controllable because the input enters it, but it is unobservable because it never appears in \( Cx \) and it is dynamically decoupled from the first two states. The first two-state subsystem is controllable and observable. Thus a minimal realization is

\[ A_m=\begin{bmatrix}0&1\\-6&-5\end{bmatrix},\quad B_m=\begin{bmatrix}0\\1\end{bmatrix},\quad C_m=\begin{bmatrix}1&0\end{bmatrix},\quad D_m=0. \]

Problem 4 (Transfer Equivalence): For the realization in Problem 3, verify by direct calculation that the full and minimal transfer functions are equal.

Solution: Since the full matrix is block diagonal,

\[ (sI-A)^{-1}B=\begin{bmatrix} (sI-A_m)^{-1}B_m\\\frac{1}{s+8} \end{bmatrix}. \]

Multiplication by \( C=\begin{bmatrix}C_m&0\end{bmatrix} \) eliminates the third component. Hence

\[ C(sI-A)^{-1}B=C_m(sI-A_m)^{-1}B_m=\frac{1}{s^2+5s+6}. \]

Problem 5 (Diagnosing Nearly Nonminimal Models): A numerical model gives singular values \( 10^1,10^0,10^{-2},10^{-9} \) for its controllability matrix. Should the fourth direction always be deleted?

Solution: Not automatically. If the model is exact and the fourth singular value is zero within numerical roundoff, deletion is justified. If the value reflects a weak but real actuator path, deletion changes the transfer matrix. In engineering computation, rank decisions must be tied to scale, modeling uncertainty, and the intended frequency range. Exact minimal realization removes only exactly uncontrollable and exactly unobservable directions; approximate reduction is a separate model-reduction problem.

Minimal realization via decomposition is based on a simple structural principle: keep precisely the states that are both reachable from the input and observable at the output. The reachable reduction removes input-inaccessible states, and the observable reduction removes states that cannot affect measured outputs. The resulting realization preserves \( G(s) \), has smaller dimension, and is minimal when it is both controllable and observable.