In exact algebra, a finite-dimensional continuous-time LTI realization \( (\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}) \) is minimal exactly when it is both controllable and observable. The transfer matrix is

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

In finite precision arithmetic, however, a mode may be formally controllable and observable while being so weakly excited or sensed that its contribution is below the numerical noise floor. This creates a practical distinction between exact minimality and numerical minimality.

\[ \operatorname{rank}(\mathbf{R}_c)=n,\quad \operatorname{rank}(\mathbf{R}_o)=n \quad \text{does not imply} \quad \kappa(\mathbf{R}_c),\kappa(\mathbf{R}_o) \text{ are small}. \]

The controllability and observability matrices are

\[ \mathbf{R}_c= \begin{bmatrix} \mathbf{B} & \mathbf{A}\mathbf{B} & \cdots & \mathbf{A}^{n-1}\mathbf{B} \end{bmatrix}, \quad \mathbf{R}_o= \begin{bmatrix} \mathbf{C} \\ \mathbf{C}\mathbf{A} \\ \vdots \\ \mathbf{C}\mathbf{A}^{n-1} \end{bmatrix}. \]

A numerical algorithm should therefore report singular values, condition numbers, and chosen tolerances—not only a binary rank statement.

flowchart TD
  A["State-space model"] --> B["Scale states, inputs, outputs"]
  B --> C["Compute rank with QR or SVD"]
  C --> D["Inspect singular values and tolerances"]
  D --> E["Separate strong and weak directions"]
  E --> F["Use balanced or orthogonal coordinates"]
  F --> G["Choose reduced order from HSV decay"]
  G --> H["Validate frequency and time responses"]
        

Let \( \mathbf{M}\in\mathbb{R}^{p\times q} \) be a matrix whose rank is needed in realization reduction. Its singular value decomposition is

\[ \mathbf{M}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^T, \quad \boldsymbol{\Sigma}= \operatorname{diag}(\sigma_1,\sigma_2,\ldots,\sigma_m), \quad \sigma_1\ge \sigma_2\ge \cdots \ge \sigma_m\ge 0. \]

The exact rank counts nonzero singular values. The numerical rank counts singular values larger than a tolerance:

\[ \operatorname{rank}_\varepsilon(\mathbf{M}) = \#\{i:\sigma_i > \varepsilon\}, \quad \varepsilon = \gamma\,\max(p,q)\,\epsilon_{mach}\,\sigma_1. \]

Here \( \epsilon_{mach} \) is machine precision and \( \gamma \) is a user-selected safety factor. When there is a clear gap \( \sigma_r \gg \sigma_{r+1} \), the numerical rank is reliable. When the singular values decay gradually, the reduced order is a modeling decision rather than a purely algebraic fact.

Proposition 1 (rank instability under perturbation). If \( \sigma_r(\mathbf{M}) > \delta \) and \( \sigma_{r+1}(\mathbf{M}) < \delta \), then any perturbation \( \Delta\mathbf{M} \) satisfying \( \|\Delta\mathbf{M}\|_2 < \frac{1}{2}(\sigma_r-\sigma_{r+1}) \) cannot move both singular values across the same threshold \( \delta \).

Proof. Weyl's singular value perturbation inequality gives

\[ |\sigma_i(\mathbf{M}+\Delta\mathbf{M})-\sigma_i(\mathbf{M})| \le \|\Delta\mathbf{M}\|_2. \]

Therefore \( \sigma_r \) remains above the midpoint between \( \sigma_r \) and \( \sigma_{r+1} \), while \( \sigma_{r+1} \) remains below it. A large singular-value gap produces a stable rank decision. A small gap makes rank tolerance-dependent.

Two internally equivalent realizations are related by a nonsingular similarity transformation:

\[ \bar{\mathbf{A} }=\mathbf{T}^{-1}\mathbf{A}\mathbf{T}, \quad \bar{\mathbf{B} }=\mathbf{T}^{-1}\mathbf{B}, \quad \bar{\mathbf{C} }=\mathbf{C}\mathbf{T}, \quad \bar{\mathbf{D} }=\mathbf{D}. \]

Algebraically, the transfer function is unchanged. Numerically, a poorly conditioned \( \mathbf{T} \) magnifies roundoff:

\[ \frac{\|\Delta(\mathbf{T}^{-1}\mathbf{A}\mathbf{T})\|_2} {\|\mathbf{T}^{-1}\mathbf{A}\mathbf{T}\|_2} \lesssim \kappa_2(\mathbf{T}) \frac{\|\Delta\mathbf{A}\|_2}{\|\mathbf{A}\|_2} +O(\epsilon_{mach}). \]

This is why canonical forms constructed from monomial bases, companion matrices, or nearly dependent Krylov vectors may be poor numerical choices even though they are excellent conceptual tools.

Scaling rule. If state components have very different magnitudes, choose a diagonal scaling \( \mathbf{S}=\operatorname{diag}(s_1,\ldots,s_n) \) and replace \( \mathbf{x}=\mathbf{S}\mathbf{z} \). Then

\[ \dot{\mathbf{z} }=\mathbf{S}^{-1}\mathbf{A}\mathbf{S}\mathbf{z} +\mathbf{S}^{-1}\mathbf{B}\mathbf{u}, \quad \mathbf{y}=\mathbf{C}\mathbf{S}\mathbf{z}+\mathbf{D}\mathbf{u}. \]

Good scaling makes the rows and columns of important matrices have comparable norms. It does not change the system, but it can greatly improve numerical reliability.

A nonminimal realization may contain an exact uncontrollable or unobservable mode that cancels from the transfer function. In computation, cancellations are often near cancellations. Consider the scalar transfer function

\[ G_\eta(s)=\frac{s+a+\eta}{(s+a)(s+b)}, \quad 0<|\eta|\ll 1. \]

If \( \eta=0 \), the pole at \( s=-a \) cancels exactly. If \( \eta\ne 0 \), the system is mathematically second order, but for small \( |\eta| \) the residue of the nearly cancelled pole is small. Numerically, the decision “remove or retain” depends on the tolerance and on the intended frequency band.

\[ G_\eta(s)= \frac{\alpha(\eta)}{s+a}+ \frac{\beta(\eta)}{s+b}, \quad |\alpha(\eta)|=O(|\eta|). \]

Thus a small pole residue indicates that the mode has weak external effect, even if it is exactly present in the internal state description.

For a stable continuous-time realization, the controllability and observability Gramians solve

\[ \mathbf{A}\mathbf{W}_c+\mathbf{W}_c\mathbf{A}^T+ \mathbf{B}\mathbf{B}^T=\mathbf{0}, \quad \mathbf{A}^T\mathbf{W}_o+\mathbf{W}_o\mathbf{A}+ \mathbf{C}^T\mathbf{C}=\mathbf{0}. \]

A state direction is important for input-output behavior only if it is both easily reachable and easily observable. This joint importance is measured by the Hankel singular values:

\[ \sigma_i^H = \sqrt{\lambda_i(\mathbf{W}_c\mathbf{W}_o)}, \quad \sigma_1^H\ge \sigma_2^H\ge \cdots \ge \sigma_n^H\ge 0. \]

A balanced realization satisfies

\[ \bar{\mathbf{W} }_c=\bar{\mathbf{W} }_o= \boldsymbol{\Sigma}= \operatorname{diag}(\sigma_1^H,\ldots,\sigma_n^H). \]

In balanced coordinates, small diagonal entries identify states that are simultaneously hard to reach and hard to observe. These are natural candidates for truncation.

flowchart TD
  A["Stable realization"] --> B["Solve controllability Gramian"]
  A --> C["Solve observability Gramian"]
  B --> D["Square-root factors"]
  C --> D
  D --> E["SVD of product of factors"]
  E --> F["Balanced coordinates"]
  F --> G["Partition strong and weak states"]
  G --> H["Truncate weak states and validate"]
        

Directly forming \( \mathbf{W}_c\mathbf{W}_o \) can square condition numbers and destroy accuracy. The square-root method works with Cholesky-like factors:

\[ \mathbf{W}_c=\mathbf{S}\mathbf{S}^T, \quad \mathbf{W}_o=\mathbf{R}\mathbf{R}^T. \]

Compute the SVD

\[ \mathbf{R}^T\mathbf{S} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^T. \]

Then one balancing transformation is

\[ \mathbf{T}=\mathbf{S}\mathbf{V}\boldsymbol{\Sigma}^{-1/2}, \quad \mathbf{T}^{-1}= \boldsymbol{\Sigma}^{-1/2}\mathbf{U}^T\mathbf{R}^T. \]

Transform the realization and partition:

\[ \bar{\mathbf{A} }= \begin{bmatrix} \mathbf{A}_{11} & \mathbf{A}_{12} \\ \mathbf{A}_{21} & \mathbf{A}_{22} \end{bmatrix}, \quad \bar{\mathbf{B} }= \begin{bmatrix} \mathbf{B}_1 \\ \mathbf{B}_2 \end{bmatrix}, \quad \bar{\mathbf{C} }= \begin{bmatrix} \mathbf{C}_1 & \mathbf{C}_2 \end{bmatrix}. \]

The order-\( r \) reduced model is

\[ \mathbf{A}_r=\mathbf{A}_{11}, \quad \mathbf{B}_r=\mathbf{B}_1, \quad \mathbf{C}_r=\mathbf{C}_1, \quad \mathbf{D}_r=\mathbf{D}. \]

For stable balanced truncation, a fundamental error estimate is

\[ \|\mathbf{G}-\mathbf{G}_r\|_\infty \le 2\sum_{i=r+1}^n \sigma_i^H. \]

The bound is conservative but extremely useful: discarded Hankel singular values quantify the worst-case input-output error.

Proposition 2 (similarity invariance of Hankel singular values). Under a nonsingular similarity transformation \( \bar{\mathbf{A} }=\mathbf{T}^{-1}\mathbf{A}\mathbf{T} \), \( \bar{\mathbf{B} }=\mathbf{T}^{-1}\mathbf{B} \), and \( \bar{\mathbf{C} }=\mathbf{C}\mathbf{T} \), the Hankel singular values are unchanged.

Proof. The transformed Gramians are

\[ \bar{\mathbf{W} }_c=\mathbf{T}^{-1}\mathbf{W}_c\mathbf{T}^{-T}, \quad \bar{\mathbf{W} }_o=\mathbf{T}^T\mathbf{W}_o\mathbf{T}. \]

Hence

\[ \bar{\mathbf{W} }_c\bar{\mathbf{W} }_o = \mathbf{T}^{-1}\mathbf{W}_c\mathbf{W}_o\mathbf{T}, \]

which is similar to \( \mathbf{W}_c\mathbf{W}_o \). Similar matrices have the same eigenvalues, so the square roots of these eigenvalues, the Hankel singular values, are invariant.

Proposition 3 (balanced truncation as an energy argument). If a balanced system has \( \sigma_i^H \) small for \( i>r \), then the truncated states have small joint reachability-observability energy.

Proof sketch. In balanced coordinates, both the minimum input energy interpretation from the controllability Gramian and the output energy interpretation from the observability Gramian assign the same diagonal weight \( \sigma_i^H \) to state coordinate \( i \). A small diagonal entry means that the coordinate is neither strongly produced by inputs nor strongly visible at outputs. Truncation therefore removes directions that are weak in the combined input-output map rather than merely small in a coordinate-dependent state norm.

The following table summarizes frequent numerical issues in minimal realization and model reduction.

A reliable workflow is: scale first, use orthogonal transformations, compute SVD diagnostics, select model order from gaps or error goals, and validate against original input-output behavior.

The Python script uses NumPy and SciPy to compute numerical rank, Gramians, Hankel singular values, square-root balancing, and reduced frequency response error.

# Chapter20_Lesson5.py
"""
Numerical issues in realization and model reduction.

This script demonstrates:
1. numerical rank tests for controllability/observability,
2. balancing transformation from Lyapunov Gramians,
3. balanced truncation and a frequency-domain error check.

Dependencies:
    numpy scipy
Optional:
    matplotlib (only if you want plots)
"""

import numpy as np
from numpy.linalg import svd, cond, norm
from scipy.linalg import solve_continuous_lyapunov, cholesky


def ctrb(A, B):
    """Kalman controllability matrix [B, AB, ..., A^(n-1)B]."""
    n = A.shape[0]
    blocks = [B]
    Ap = np.eye(n)
    for _ in range(1, n):
        Ap = Ap @ A
        blocks.append(Ap @ B)
    return np.hstack(blocks)


def obsv(A, C):
    """Kalman observability matrix [C; CA; ...; CA^(n-1)]."""
    n = A.shape[0]
    blocks = [C]
    Ap = np.eye(n)
    for _ in range(1, n):
        Ap = Ap @ A
        blocks.append(C @ Ap)
    return np.vstack(blocks)


def numerical_rank(M, rtol=None):
    """SVD-based numerical rank."""
    s = svd(M, compute_uv=False)
    if rtol is None:
        rtol = max(M.shape) * np.finfo(float).eps
    threshold = rtol * s[0] if s.size else 0.0
    return int(np.sum(s > threshold)), s, threshold


def continuous_gramians(A, B, C):
    """
    Solve continuous-time Lyapunov equations:
        A Wc + Wc A.T + B B.T = 0
        A.T Wo + Wo A + C.T C = 0
    Requires A Hurwitz.
    """
    Wc = solve_continuous_lyapunov(A, -(B @ B.T))
    Wo = solve_continuous_lyapunov(A.T, -(C.T @ C))
    Wc = 0.5 * (Wc + Wc.T)
    Wo = 0.5 * (Wo + Wo.T)
    return Wc, Wo


def low_rank_cholesky_psd(W, jitter=1e-12):
    """
    Cholesky factor for a positive semidefinite Gramian.
    Adds small diagonal jitter if roundoff produces tiny negative eigenvalues.
    """
    n = W.shape[0]
    try:
        return cholesky(W, lower=True)
    except Exception:
        return cholesky(W + jitter * np.eye(n), lower=True)


def balanced_realization(A, B, C, D):
    """
    Square-root balanced realization.
    If Wc = S S.T and Wo = R R.T, compute SVD:
        R.T S = U Sigma V.T
    Then:
        T    = S V Sigma^(-1/2)
        Tinv = Sigma^(-1/2) U.T R.T
    """
    Wc, Wo = continuous_gramians(A, B, C)
    S = low_rank_cholesky_psd(Wc)
    R = low_rank_cholesky_psd(Wo)

    U, hsv, Vt = svd(R.T @ S)
    inv_sqrt = np.diag(1.0 / np.sqrt(hsv))

    T = S @ Vt.T @ inv_sqrt
    Tinv = inv_sqrt @ U.T @ R.T

    Ab = Tinv @ A @ T
    Bb = Tinv @ B
    Cb = C @ T
    Db = D.copy()

    return Ab, Bb, Cb, Db, hsv, Wc, Wo, T, Tinv


def balanced_truncate(A, B, C, D, r):
    """Keep the first r balanced states."""
    Ab, Bb, Cb, Db, hsv, Wc, Wo, T, Tinv = balanced_realization(A, B, C, D)
    Ar = Ab[:r, :r]
    Br = Bb[:r, :]
    Cr = Cb[:, :r]
    Dr = Db
    return Ar, Br, Cr, Dr, hsv


def freq_response(A, B, C, D, w):
    """Evaluate G(jw)=C(jwI-A)^(-1)B + D for SISO or MIMO systems."""
    n = A.shape[0]
    I = np.eye(n)
    values = []
    for wi in w:
        G = C @ np.linalg.solve(1j * wi * I - A, B) + D
        values.append(G)
    return np.array(values)


def main():
    # A stable fourth-order example with weakly coupled fast states.
    A = np.array([
        [-0.20,  0.05,  0.00,  0.00],
        [ 0.00, -1.00,  0.10,  0.00],
        [ 0.00,  0.00, -8.00,  0.20],
        [ 0.00,  0.00,  0.00, -20.0],
    ], dtype=float)
    B = np.array([[1.0], [0.4], [0.05], [0.01]])
    C = np.array([[1.0, 0.3, 0.02, 0.005]])
    D = np.array([[0.0]])

    print("Condition number of A:", cond(A))

    Rc = ctrb(A, B)
    Ro = obsv(A, C)
    rank_c, sc, thc = numerical_rank(Rc)
    rank_o, so, tho = numerical_rank(Ro)

    print("\nControllability singular values:", sc)
    print("Numerical controllability rank:", rank_c, "threshold:", thc)
    print("\nObservability singular values:", so)
    print("Numerical observability rank:", rank_o, "threshold:", tho)

    Wc, Wo = continuous_gramians(A, B, C)
    print("\ncond(Wc):", cond(Wc))
    print("cond(Wo):", cond(Wo))

    Ar, Br, Cr, Dr, hsv = balanced_truncate(A, B, C, D, r=2)
    print("\nHankel singular values:", hsv)
    print("Reduced A:\n", Ar)
    print("Reduced B:\n", Br)
    print("Reduced C:\n", Cr)

    # Balanced truncation H-infinity error upper bound:
    # ||G-Gr||_inf <= 2 * sum(discarded Hankel singular values)
    bound = 2.0 * np.sum(hsv[2:])
    print("\nBalanced-truncation error bound:", bound)

    # Coarse numerical frequency-response check.
    w = np.logspace(-3, 3, 300)
    G = freq_response(A, B, C, D, w)
    Gr = freq_response(Ar, Br, Cr, Dr, w)
    err = np.max(np.abs(G[:, 0, 0] - Gr[:, 0, 0]))
    print("Coarse sampled max frequency error:", err)


if __name__ == "__main__":
    main()
      

The C++ program uses the Eigen library. It focuses on robust SVD rank diagnostics, which are the first step before performing any realization reduction.

// Chapter20_Lesson5.cpp
/*
Numerical rank diagnostics for realization reduction.

Dependencies:
    Eigen 3
Compile example:
    g++ -std=c++17 Chapter20_Lesson5.cpp -I /path/to/eigen -O2 -o Chapter20_Lesson5
*/

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

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

MatrixXd matrixPower(const MatrixXd& A, int k) {
    MatrixXd P = MatrixXd::Identity(A.rows(), A.cols());
    for (int i = 0; i < k; ++i) {
        P = P * A;
    }
    return P;
}

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

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

int numericalRank(const MatrixXd& M, double rtol) {
    Eigen::JacobiSVD<MatrixXd> svd(M);
    VectorXd s = svd.singularValues();
    if (s.size() == 0) return 0;
    double threshold = rtol * s(0);
    int r = 0;
    for (int i = 0; i < s.size(); ++i) {
        if (s(i) > threshold) ++r;
    }
    return r;
}

void printSingularValues(const std::string& name, const MatrixXd& M) {
    Eigen::JacobiSVD<MatrixXd> svd(M);
    std::cout << name << " singular values:\n"
              << svd.singularValues().transpose() << "\n\n";
}

int main() {
    MatrixXd A(4, 4);
    A << -0.20,  0.05,  0.00,  0.00,
          0.00, -1.00,  0.10,  0.00,
          0.00,  0.00, -8.00,  0.20,
          0.00,  0.00,  0.00, -20.0;

    MatrixXd B(4, 1);
    B << 1.0, 0.4, 0.05, 0.01;

    MatrixXd C(1, 4);
    C << 1.0, 0.3, 0.02, 0.005;

    MatrixXd Rc = controllabilityMatrix(A, B);
    MatrixXd Ro = observabilityMatrix(A, C);

    double eps = std::numeric_limits<double>::epsilon();
    double rtolC = std::max(Rc.rows(), Rc.cols()) * eps;
    double rtolO = std::max(Ro.rows(), Ro.cols()) * eps;

    printSingularValues("Controllability matrix", Rc);
    printSingularValues("Observability matrix", Ro);

    std::cout << "Numerical controllability rank = "
              << numericalRank(Rc, rtolC) << "\n";
    std::cout << "Numerical observability rank = "
              << numericalRank(Ro, rtolO) << "\n";

    Eigen::JacobiSVD<MatrixXd> svdC(Rc);
    double conditionEstimate = svdC.singularValues()(0) /
                               svdC.singularValues()(svdC.singularValues().size() - 1);

    std::cout << "Estimated condition number of controllability matrix = "
              << conditionEstimate << "\n";

    std::cout << "\nInterpretation:\n";
    std::cout << "A formal full rank result can be misleading when the smallest\n";
    std::cout << "singular values are near the tolerance. In that case, reduction\n";
    std::cout << "or controller synthesis should be based on SVD/QR diagnostics\n";
    std::cout << "rather than exact symbolic rank.\n";

    return 0;
}
      

The Java example uses Apache Commons Math for singular value decomposition. It computes controllability and observability matrices and reports numerical ranks.

// Chapter20_Lesson5.java
/*
Numerical realization diagnostics using Apache Commons Math.

Dependency:
    commons-math3-3.6.1.jar

Compile:
    javac -cp commons-math3-3.6.1.jar Chapter20_Lesson5.java

Run:
    java -cp .:commons-math3-3.6.1.jar Chapter20_Lesson5
On Windows use:
    java -cp .;commons-math3-3.6.1.jar Chapter20_Lesson5
*/

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

public class Chapter20_Lesson5 {
    static RealMatrix multiplyPower(RealMatrix A, int k) {
        RealMatrix P = MatrixUtils.createRealIdentityMatrix(A.getRowDimension());
        for (int i = 0; i < k; i++) {
            P = P.multiply(A);
        }
        return P;
    }

    static RealMatrix horizontalStack(RealMatrix[] blocks) {
        int rows = blocks[0].getRowDimension();
        int cols = 0;
        for (RealMatrix b : blocks) {
            cols += b.getColumnDimension();
        }
        RealMatrix out = MatrixUtils.createRealMatrix(rows, cols);
        int c0 = 0;
        for (RealMatrix b : blocks) {
            out.setSubMatrix(b.getData(), 0, c0);
            c0 += b.getColumnDimension();
        }
        return out;
    }

    static RealMatrix verticalStack(RealMatrix[] blocks) {
        int cols = blocks[0].getColumnDimension();
        int rows = 0;
        for (RealMatrix b : blocks) {
            rows += b.getRowDimension();
        }
        RealMatrix out = MatrixUtils.createRealMatrix(rows, cols);
        int r0 = 0;
        for (RealMatrix b : blocks) {
            out.setSubMatrix(b.getData(), r0, 0);
            r0 += b.getRowDimension();
        }
        return out;
    }

    static RealMatrix ctrb(RealMatrix A, RealMatrix B) {
        int n = A.getRowDimension();
        RealMatrix[] blocks = new RealMatrix[n];
        RealMatrix Ap = MatrixUtils.createRealIdentityMatrix(n);
        for (int k = 0; k < n; k++) {
            blocks[k] = Ap.multiply(B);
            Ap = Ap.multiply(A);
        }
        return horizontalStack(blocks);
    }

    static RealMatrix obsv(RealMatrix A, RealMatrix C) {
        int n = A.getRowDimension();
        RealMatrix[] blocks = new RealMatrix[n];
        RealMatrix Ap = MatrixUtils.createRealIdentityMatrix(n);
        for (int k = 0; k < n; k++) {
            blocks[k] = C.multiply(Ap);
            Ap = Ap.multiply(A);
        }
        return verticalStack(blocks);
    }

    static int numericalRank(RealMatrix M) {
        SingularValueDecomposition svd = new SingularValueDecomposition(M);
        double[] s = svd.getSingularValues();
        if (s.length == 0) return 0;
        double eps = Math.ulp(1.0);
        double tol = Math.max(M.getRowDimension(), M.getColumnDimension()) * eps * s[0];
        int rank = 0;
        for (double value : s) {
            if (value > tol) rank++;
        }
        return rank;
    }

    static void printSingularValues(String name, RealMatrix M) {
        SingularValueDecomposition svd = new SingularValueDecomposition(M);
        System.out.println(name + " singular values:");
        for (double s : svd.getSingularValues()) {
            System.out.printf("% .6e ", s);
        }
        System.out.println("\n");
    }

    public static void main(String[] args) {
        double[][] Adata = {
            {-0.20,  0.05,  0.00,  0.00},
            { 0.00, -1.00,  0.10,  0.00},
            { 0.00,  0.00, -8.00,  0.20},
            { 0.00,  0.00,  0.00, -20.0}
        };
        double[][] Bdata = { {1.0}, {0.4}, {0.05}, {0.01} };
        double[][] Cdata = { {1.0, 0.3, 0.02, 0.005} };

        RealMatrix A = MatrixUtils.createRealMatrix(Adata);
        RealMatrix B = MatrixUtils.createRealMatrix(Bdata);
        RealMatrix C = MatrixUtils.createRealMatrix(Cdata);

        RealMatrix Rc = ctrb(A, B);
        RealMatrix Ro = obsv(A, C);

        printSingularValues("Controllability matrix", Rc);
        printSingularValues("Observability matrix", Ro);

        System.out.println("Numerical controllability rank = " + numericalRank(Rc));
        System.out.println("Numerical observability rank = " + numericalRank(Ro));

        SingularValueDecomposition svdC = new SingularValueDecomposition(Rc);
        double[] s = svdC.getSingularValues();
        double conditionEstimate = s[0] / s[s.length - 1];

        System.out.printf("Estimated condition number of controllability matrix = %.6e%n",
                          conditionEstimate);

        System.out.println("\nLesson point:");
        System.out.println("Near-zero singular values indicate weakly reachable or weakly");
        System.out.println("observable directions. Such directions are not safely detected");
        System.out.println("by exact rank logic in finite precision arithmetic.");
    }
}
      

MATLAB directly exposes realization reduction tools through the Control System Toolbox. The script also optionally creates a simple Simulink model with original and reduced state-space blocks.

% Chapter20_Lesson5.m
% Numerical issues in realization and model reduction.
% Requires Control System Toolbox for ss, gram, balreal, modred, minreal.
% Simulink section requires Simulink.

clear; clc; close all;

A = [-0.20  0.05  0.00  0.00;
      0.00 -1.00  0.10  0.00;
      0.00  0.00 -8.00  0.20;
      0.00  0.00  0.00 -20.0];

B = [1.0; 0.4; 0.05; 0.01];
C = [1.0 0.3 0.02 0.005];
D = 0;

sys = ss(A,B,C,D);

fprintf('Condition number of controllability matrix: %.6e\n', cond(ctrb(A,B)));
fprintf('Condition number of observability matrix: %.6e\n', cond(obsv(A,C)));

svC = svd(ctrb(A,B));
svO = svd(obsv(A,C));

disp('Controllability singular values:');
disp(svC.');

disp('Observability singular values:');
disp(svO.');

% MATLAB minreal uses tolerance logic to cancel weak pole-zero pairs.
sysMinDefault = minreal(sys);
sysMinTight   = minreal(sys, 1e-10);

disp('Default minreal order:');
disp(order(sysMinDefault));

disp('Tight tolerance minreal order:');
disp(order(sysMinTight));

% Balanced realization and model reduction.
[sysBal, hsv] = balreal(sys);
disp('Hankel singular values:');
disp(hsv.');

r = 2;
elim = (r+1):order(sysBal);
sysRed = modred(sysBal, elim, 'truncate');

fprintf('Reduced model order: %d\n', order(sysRed));
fprintf('Balanced truncation error upper bound: %.6e\n', 2*sum(hsv(r+1:end)));

figure;
bodemag(sys, sysRed);
grid on;
legend('Original','Balanced truncation r=2');
title('Chapter 20 Lesson 5: Numerical model reduction comparison');

% Optional Simulink block creation for side-by-side State-Space blocks.
% This programmatically creates a simple model with original and reduced systems.
modelName = 'Chapter20_Lesson5_Simulink';
if exist('simulink', 'file') == 4
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    new_system(modelName);
    open_system(modelName);

    add_block('simulink/Sources/Step', [modelName '/Step'], ...
        'Position', [40 100 70 130]);

    add_block('simulink/Continuous/State-Space', [modelName '/Original State-Space'], ...
        'A', 'A', 'B', 'B', 'C', 'C', 'D', 'D', ...
        'Position', [140 60 280 120]);

    Ar = sysRed.A; Br = sysRed.B; Cr = sysRed.C; Dr = sysRed.D; %#ok<NASGU>
    assignin('base','Ar',Ar);
    assignin('base','Br',Br);
    assignin('base','Cr',Cr);
    assignin('base','Dr',Dr);

    add_block('simulink/Continuous/State-Space', [modelName '/Reduced State-Space'], ...
        'A', 'Ar', 'B', 'Br', 'C', 'Cr', 'D', 'Dr', ...
        'Position', [140 160 280 220]);

    add_block('simulink/Sinks/Scope', [modelName '/Scope'], ...
        'Position', [370 90 400 180]);

    add_line(modelName, 'Step/1', 'Original State-Space/1');
    add_line(modelName, 'Step/1', 'Reduced State-Space/1');
    add_line(modelName, 'Original State-Space/1', 'Scope/1');
    add_line(modelName, 'Reduced State-Space/1', 'Scope/2');

    save_system(modelName);
    fprintf('Created Simulink model: %s.slx\n', modelName);
end
      

Mathematica is useful for mixing exact matrix construction with numerical singular-value diagnostics. The downloadable file is named Chapter20_Lesson5.nb.

(* Chapter20_Lesson5.nb *)
(* Wolfram Mathematica / Wolfram Language code for numerical realization diagnostics *)

ClearAll["Global`*"];

A = { {-0.20, 0.05, 0.00, 0.00},
     { 0.00,-1.00, 0.10, 0.00},
     { 0.00, 0.00,-8.00, 0.20},
     { 0.00, 0.00, 0.00,-20.0} };

B = { {1.0}, {0.4}, {0.05}, {0.01} };
Cmat = { {1.0, 0.3, 0.02, 0.005} };
Dmat = { {0.0} };

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

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

numericalRank[M_] := Module[{s, tol},
  s = SingularValueList[N[M]];
  tol = Max[Dimensions[M]]*$MachineEpsilon*First[s];
  Count[s, x_ /; x > tol]
];

Rc = ctrb[A, B];
Ro = obsv[A, Cmat];

Print["Controllability singular values = ", SingularValueList[N[Rc]]];
Print["Observability singular values = ", SingularValueList[N[Ro]]];

Print["Numerical controllability rank = ", numericalRank[Rc]];
Print["Numerical observability rank = ", numericalRank[Ro]];

Print["Condition number of controllability matrix = ", ConditionNumber[N[Rc]]];
Print["Condition number of observability matrix = ", ConditionNumber[N[Ro]]];

(* Continuous-time Gramians using matrix equations:
   A.Wc + Wc.Transpose[A] + B.Transpose[B] == 0
   Transpose[A].Wo + Wo.A + Transpose[C].C == 0
*)
varsWc = Array[wc, {4, 4}];
varsWo = Array[wo, {4, 4}];

eqWc = Thread[Flatten[A.varsWc + varsWc.Transpose[A] + B.Transpose[B]] == 0];
eqWo = Thread[Flatten[Transpose[A].varsWo + varsWo.A + Transpose[Cmat].Cmat] == 0];

solWc = Solve[eqWc, Flatten[varsWc]][[1]];
solWo = Solve[eqWo, Flatten[varsWo]][[1]];

Wc = N[varsWc /. solWc];
Wo = N[varsWo /. solWo];

Print["Eigenvalues of Wc = ", Eigenvalues[Wc]];
Print["Eigenvalues of Wo = ", Eigenvalues[Wo]];

hsv = Sqrt[Eigenvalues[Wc.Wo]];
Print["Hankel singular values (unsorted) = ", hsv];

sortedHsv = Reverse[Sort[Abs[hsv]]];
Print["Hankel singular values (sorted) = ", sortedHsv];

r = 2;
Print["Balanced truncation theoretical error bound for r=2 = ",
      2 Total[sortedHsv[[r + 1 ;;]]]];
      

Problem 1 (Numerical Rank Decision): Suppose the singular values of a controllability matrix are \( 10^2, 10^0, 10^{-6}, 10^{-13} \). In double precision, estimate the numerical rank using \( \varepsilon=\max(p,q)\epsilon_{mach}\sigma_1 \) with \( \max(p,q)=8 \).

Solution: Since \( \epsilon_{mach}\approx 2.22\times 10^{-16} \),

\[ \varepsilon = 8(2.22\times 10^{-16})(10^2) \approx 1.78\times 10^{-13}. \]

The first three singular values are above this threshold. The fourth singular value \( 10^{-13} \) is below it. Thus the numerical rank is \( 3 \). If the model has four states, the fourth state is numerically weakly reachable.

Problem 2 (Effect of Bad Scaling): Let \( \mathbf{x}=\mathbf{S}\mathbf{z} \) with \( \mathbf{S}=\operatorname{diag}(10^6,1) \). Explain why exact transfer equivalence does not guarantee numerical equivalence.

Solution: The transformed matrices are \( \bar{\mathbf{A} }=\mathbf{S}^{-1}\mathbf{A}\mathbf{S} \), \( \bar{\mathbf{B} }=\mathbf{S}^{-1}\mathbf{B} \), and \( \bar{\mathbf{C} }=\mathbf{C}\mathbf{S} \). Since \( \kappa_2(\mathbf{S})=10^6 \), roundoff errors may be amplified by roughly six orders of magnitude in computations involving the transformation. The transfer function is algebraically unchanged, but rank tests, matrix inversions, and Gramian computations may become much less accurate.

Problem 3 (Near Cancellation): For \( G_\eta(s)=\frac{s+2+\eta}{(s+2)(s+5)} \), compute the residue of the pole at \( s=-2 \).

Solution: Write

\[ G_\eta(s)=\frac{\alpha}{s+2}+\frac{\beta}{s+5}. \]

Multiplying by \( (s+2)(s+5) \) gives \( s+2+\eta=\alpha(s+5)+\beta(s+2) \). Setting \( s=-2 \) yields

\[ \eta=3\alpha \quad \Rightarrow \quad \alpha=\frac{\eta}{3}. \]

Thus the pole at \( -2 \) has residue proportional to \( \eta \). When \( |\eta| \) is very small, this mode has a weak input-output contribution.

Problem 4 (Balanced-Truncation Error Bound): A stable balanced realization has Hankel singular values \( 4,1,0.05,0.01 \). If an order-2 model is retained, bound the approximation error in the \( \mathcal{H}_\infty \) norm.

Solution: The discarded Hankel singular values are \( 0.05 \) and \( 0.01 \). Hence

\[ \|\mathbf{G}-\mathbf{G}_2\|_\infty \le 2(0.05+0.01)=0.12. \]

This means the worst-case gain error over all frequencies is bounded by \( 0.12 \).

Problem 5 (Why Square-Root Balancing?): Explain why computing the eigenvalues of \( \mathbf{W}_c\mathbf{W}_o \) directly may be less reliable than computing the SVD of \( \mathbf{R}^T\mathbf{S} \), where \( \mathbf{W}_c=\mathbf{S}\mathbf{S}^T \) and \( \mathbf{W}_o=\mathbf{R}\mathbf{R}^T \).

Solution: Forming \( \mathbf{W}_c\mathbf{W}_o \) can multiply two ill-conditioned matrices and can also produce a nonsymmetric matrix whose eigenvalues are sensitive to roundoff. The square-root method uses Gramian factors and an SVD, which is backward stable and avoids directly squaring the conditioning. Therefore it is preferred in numerical model reduction.

Minimal realization is an exact concept, but practical realization reduction is a numerical decision. Reliable computation requires scaling, SVD or QR rank tests, condition-number inspection, tolerance selection, stable Gramian computations, and validation of the reduced input-output behavior. Balanced truncation provides a principled way to reduce stable systems because it ranks states by joint reachability and observability through Hankel singular values.