Consider the continuous-time LTI system \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \), \( \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u} \), where \( \mathbf{x}\in\mathbb{R}^n \), \( \mathbf{u}\in\mathbb{R}^m \), and \( \mathbf{y}\in\mathbb{R}^p \). Observability asks whether the initial state can be inferred from the measured output over a time interval, assuming the input and model are known.

Since the known input contribution can be subtracted from the output, the essence of observability is the zero-input response \( \mathbf{y}_0(t)=\mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0 \). Two initial states \( \mathbf{x}_1 \) and \( \mathbf{x}_2 \) are output-indistinguishable when their difference generates zero output for all relevant times:

\[ \mathbf{C}e^{\mathbf{A}t}(\mathbf{x}_1-\mathbf{x}_2)=\mathbf{0}, \qquad t\ge 0. \]

Therefore, the central object in this lesson is the set of initial-state directions that produce no zero-input output. This is the unobservable subspace. Its complement, more precisely the quotient of the state space by this hidden subspace, contains the state information that is theoretically reconstructible from the output.

flowchart TD
  X["Initial state x0"] --> E["Internal evolution exp(A t) x0"]
  E --> C["Sensor map C"]
  C --> Y["Measured output y0(t)"]
  X --> H["Hidden directions"]
  H --> Z["C exp(A t) x_hidden = 0 for all t"]
  Z --> U["Unobservable subspace"]
  Y --> O["Observable information: state modulo hidden directions"]
        

The unobservable subspace of the pair \( (\mathbf{A},\mathbf{C}) \) is

\[ \mathscr{N}_{\bar{o} } = \left\{\mathbf{x}\in\mathbb{R}^n: \mathbf{C}e^{\mathbf{A}t}\mathbf{x}=\mathbf{0}\;\text{for all}\;t\ge 0 \right\}. \]

Equivalently, a vector belongs to \( \mathscr{N}_{\bar{o} } \) if it cannot be detected by the output, even though it may generate nontrivial internal state motion. The finite-dimensional rank test follows from the observability matrix

\[ \mathcal{O}_n = \begin{bmatrix} \mathbf{C} \\ \mathbf{C}\mathbf{A} \\ \mathbf{C}\mathbf{A}^2 \\ \vdots \\ \mathbf{C}\mathbf{A}^{n-1} \end{bmatrix}\in\mathbb{R}^{pn\times n}. \]

The unobservable subspace is exactly the null space of this matrix:

\[ \boxed{\mathscr{N}_{\bar{o} }=\ker(\mathcal{O}_n)}. \]

The system pair \( (\mathbf{A},\mathbf{C}) \) is observable when \( \mathscr{N}_{\bar{o} }=\{\mathbf{0}\} \), or equivalently

\[ \operatorname{rank}(\mathcal{O}_n)=n. \]

If \( \operatorname{rank}(\mathcal{O}_n)=r<n \), then there are \( n-r \) linearly independent hidden initial-state directions. These directions are not numerical artifacts; they represent structural loss of information through the sensor matrix \( \mathbf{C} \) and its interaction with the dynamics \( \mathbf{A} \).

The definition of \( \mathscr{N}_{\bar{o} } \) involves all \( t\ge 0 \), but finite-dimensional linear algebra reduces the test to \( n \) block rows. Expanding the matrix exponential gives

\[ \mathbf{C}e^{\mathbf{A}t}\mathbf{x}=\sum_{k=0}^\infty \frac{t^k}{k!}\mathbf{C}\mathbf{A}^k\mathbf{x}. \]

If this analytic function is identically zero, then all its Taylor coefficients vanish:

\[ \mathbf{C}\mathbf{A}^k\mathbf{x}=\mathbf{0},\qquad k=0,1,2,\ldots. \]

Conversely, if the first \( n \) coefficients vanish, then all higher ones vanish by the Cayley-Hamilton theorem. If the characteristic polynomial of \( \mathbf{A} \) is

\[ q(s)=s^n+a_{n-1}s^{n-1}+\cdots+a_1s+a_0, \]

then \( q(\mathbf{A})=\mathbf{0} \), so

\[ \mathbf{A}^n=-a_{n-1}\mathbf{A}^{n-1}-\cdots-a_1\mathbf{A}-a_0\mathbf{I}. \]

Multiplying by \( \mathbf{C} \) and applying the same relation repeatedly shows that \( \mathbf{C}\mathbf{A}^k\mathbf{x}=\mathbf{0} \) for every \( k\ge n \). Hence \( \ker(\mathcal{O}_n) \) captures all hidden directions.

In Euclidean coordinates, it is common to call \( \mathscr{X}_o=\operatorname{range}(\mathcal{O}_n^T) \) the observable subspace. Since the fundamental theorem of linear algebra gives

\[ \operatorname{range}(\mathcal{O}_n^T)=\ker(\mathcal{O}_n)^\perp, \]

we obtain the orthogonal decomposition

\[ \mathbb{R}^n=\mathscr{X}_o\oplus\mathscr{N}_{\bar{o} }, \qquad \dim(\mathscr{X}_o)=\operatorname{rank}(\mathcal{O}_n). \]

The phrase “observable subspace” must be interpreted carefully. The reconstructible object is not always a coordinate subset of \( \mathbf{x} \); it is the equivalence class of \( \mathbf{x} \) modulo the unobservable subspace:

\[ \mathbf{x}\sim\mathbf{z} \quad \Longleftrightarrow \quad \mathbf{x}-\mathbf{z}\in\mathscr{N}_{\bar{o} }. \]

Thus the output determines the coset \( \mathbf{x}+\mathscr{N}_{\bar{o} } \), not the hidden component inside \( \mathscr{N}_{\bar{o} } \). This distinction becomes crucial when forming minimal realizations: a nonminimal realization can contain internal modes that never appear in the transfer function.

A central property is that the unobservable subspace is invariant under the autonomous dynamics:

\[ \boxed{\mathbf{A}\mathscr{N}_{\bar{o} }\subseteq \mathscr{N}_{\bar{o} } }. \]

Proof. Let \( \mathbf{v}\in\mathscr{N}_{\bar{o} }=\ker(\mathcal{O}_n) \). Then \( \mathbf{C}\mathbf{A}^k\mathbf{v}=\mathbf{0} \) for \( k=0,1,\ldots,n-1 \). For \( k=0,1,\ldots,n-2 \),

\[ \mathbf{C}\mathbf{A}^k(\mathbf{A}\mathbf{v})= \mathbf{C}\mathbf{A}^{k+1}\mathbf{v}=\mathbf{0}. \]

The remaining condition \( \mathbf{C}\mathbf{A}^{n-1}(\mathbf{A}\mathbf{v})=\mathbf{C}\mathbf{A}^n\mathbf{v} \) is also zero by Cayley-Hamilton because \( \mathbf{A}^n \) is a linear combination of \( \mathbf{I},\mathbf{A},\ldots,\mathbf{A}^{n-1} \). Therefore \( \mathbf{A}\mathbf{v}\in\ker(\mathcal{O}_n) \).

This invariance is the algebraic reason a coordinate transformation can isolate hidden dynamics in a lower block of the transformed state matrix.

Choose columns of \( \mathbf{T}_u \) as a basis for \( \mathscr{N}_{\bar{o} } \), and choose \( \mathbf{T}_o \) so that \( \mathbf{T}=[\mathbf{T}_o\;\mathbf{T}_u] \) is nonsingular. With \( \mathbf{x}=\mathbf{T}\mathbf{z} \) and \( \mathbf{z}=[\mathbf{z}_o^T\;\mathbf{z}_u^T]^T \), define

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

Since \( \mathbf{A}\mathscr{N}_{\bar{o} }\subseteq \mathscr{N}_{\bar{o} } \) and \( \mathbf{C}\mathscr{N}_{\bar{o} }=\{\mathbf{0}\} \), the pair has the block structure

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

The zero top-right block means hidden coordinates \( \mathbf{z}_u \) do not drive observable coordinates \( \mathbf{z}_o \). The zero right block of \( \bar{\mathbf{C} } \) means hidden coordinates do not appear directly in the output. The lower-left block \( \mathbf{A}_{uo} \) may be nonzero: observable dynamics can excite hidden coordinates, but those hidden coordinates remain invisible in the measured output.

The PBH test gives a modal interpretation. The pair \( (\mathbf{A},\mathbf{C}) \) is observable if and only if

\[ \operatorname{rank}\!\begin{bmatrix} \lambda\mathbf{I}-\mathbf{A} \\ \mathbf{C} \end{bmatrix}=n, \qquad \text{for every }\lambda\in\sigma(\mathbf{A}). \]

If there exists a nonzero eigenvector \( \mathbf{w} \) such that \( \mathbf{A}\mathbf{w}=\lambda\mathbf{w} \) and \( \mathbf{C}\mathbf{w}=\mathbf{0} \), then the mode \( e^{\lambda t}\mathbf{w} \) is hidden from the output. For repeated eigenvalues and Jordan chains, the generalized eigenvector conditions are captured by the same PBH rank test.

In numerical work, the observability matrix can become ill-conditioned. Rank decisions should therefore be made using singular values, not exact symbolic equality to zero unless the matrices are rational or symbolic.

flowchart TD
  A["Given matrices A and C"] --> O["Build O_n = stack(C A^k), k=0,...,n-1"]
  O --> SVD["Compute SVD or rank-revealing factorization"]
  SVD --> R["Rank r and nullity n-r"]
  R --> U["Null space gives hidden basis T_u"]
  R --> V["Row space gives observable basis T_o"]
  U --> T["Form T = [T_o T_u]"]
  V --> T
  T --> B["Transform A,B,C and inspect block zeros"]
        

Useful software tools include NumPy/SciPy and python-control in Python, Eigen or Armadillo in C++, EJML or Apache Commons Math in Java, obsv, null, ss, and Simulink State-Space blocks in MATLAB, and MatrixPower, NullSpace, and MatrixExp in Wolfram Mathematica. The downloadable files for this lesson include both library-based and from-scratch implementations where appropriate.

The following script computes \( \mathcal{O}_n \), estimates its numerical rank using SVD, extracts the unobservable subspace, checks invariance, and verifies that adding a hidden initial state does not change the zero-input output.

Chapter19_Lesson2.py

#!/usr/bin/env python3
"""
Chapter19_Lesson2.py
Observable/unobservable subspaces for a continuous-time LTI system.

Libraries covered:
  - NumPy: matrix powers, SVD rank, numerical null space.
  - SciPy / python-control can be used for larger control workflows, but this
    file intentionally implements the core observability operations from scratch.
"""
import numpy as np


def observability_matrix(A: np.ndarray, C: np.ndarray) -> np.ndarray:
    """Return O_n = [C; C A; ...; C A^(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 numerical_rank(M: np.ndarray, tol: float | None = None) -> int:
    """SVD numerical rank with MATLAB-like default tolerance."""
    singular_values = np.linalg.svd(M, compute_uv=False)
    if tol is None:
        tol = max(M.shape) * np.finfo(float).eps * (singular_values[0] if singular_values.size else 1.0)
    return int(np.sum(singular_values > tol))


def nullspace(M: np.ndarray, tol: float | None = None) -> np.ndarray:
    """Orthonormal basis for ker(M), returned as columns."""
    U, S, Vt = np.linalg.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)
    rank = int(np.sum(S > tol))
    return Vt[rank:, :].T.copy()


def row_space_basis(M: np.ndarray, tol: float | None = None) -> np.ndarray:
    """Orthonormal basis for range(M^T), returned as columns."""
    U, S, Vt = np.linalg.svd(M, full_matrices=False)
    if tol is None:
        tol = max(M.shape) * np.finfo(float).eps * (S[0] if S.size else 1.0)
    rank = int(np.sum(S > tol))
    return Vt[:rank, :].T.copy()


def a_invariance_residual(A: np.ndarray, N: np.ndarray) -> float:
    """
    For an orthonormal basis N of the unobservable subspace, compute
    ||(I - N N^T) A N||_F.  It should be near zero when ker(O_n) is A-invariant.
    """
    if N.size == 0:
        return 0.0
    n = A.shape[0]
    projector_perp = np.eye(n) - N @ N.T
    return float(np.linalg.norm(projector_perp @ A @ N, ord="fro"))


def observable_decomposition(A: np.ndarray, C: np.ndarray):
    """
    Build an orthonormal coordinate matrix Q = [Q_o Q_u], where Q_u spans the
    unobservable subspace and Q_o spans its orthogonal complement.  With x = Q z,
    the transformed pair has Cbar = [C_o 0] and Abar has a zero top-right block.
    """
    O = observability_matrix(A, C)
    Qu = nullspace(O)
    Qo = row_space_basis(O)

    # Complete the basis robustly if roundoff or rank decisions leave a gap.
    Q = np.hstack([Qo, Qu]) if Qu.size else Qo
    if Q.shape[1] < A.shape[0]:
        # QR completion from random directions projected against existing Q.
        rng = np.random.default_rng(0)
        R = rng.normal(size=(A.shape[0], A.shape[0] - Q.shape[1]))
        R = R - Q @ (Q.T @ R)
        extra, _ = np.linalg.qr(R)
        Q = np.hstack([Q, extra[:, : A.shape[0] - Q.shape[1]]])

    Abar = Q.T @ A @ Q
    Cbar = C @ Q
    return O, Qo, Qu, Q, Abar, Cbar


def main() -> None:
    # Two decoupled second-order modes.  The sensor measures only the first mode.
    A = np.array([
        [0.0, 1.0, 0.0, 0.0],
        [-2.0, -3.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 1.0],
        [0.0, 0.0, -5.0, -1.0],
    ])
    B = np.array([[0.0], [1.0], [0.0], [1.0]])
    C = np.array([[1.0, 0.0, 0.0, 0.0]])
    D = np.array([[0.0]])

    O, Qo, Qu, Q, Abar, Cbar = observable_decomposition(A, C)

    print("Observability matrix O_n:")
    print(O)
    print("rank(O_n) =", numerical_rank(O), "out of n =", A.shape[0])
    print("\nBasis for unobservable subspace ker(O_n), columns of Q_u:")
    print(Qu)
    print("A-invariance residual ||(I-Q_u Q_u^T) A Q_u||_F =", a_invariance_residual(A, Qu))
    print("\nTransformed Abar = Q^T A Q:")
    print(np.round(Abar, 10))
    print("\nTransformed Cbar = C Q:")
    print(np.round(Cbar, 10))

    # Demonstrate output indistinguishability: adding an unobservable initial state
    # produces no change in zero-input output y(t) = C exp(A t) x0.
    try:
        from scipy.linalg import expm
        x0_observable = np.array([1.0, 0.0, 0.0, 0.0])
        x0_hidden = np.array([0.0, 0.0, 1.0, -1.0])
        times = np.linspace(0.0, 5.0, 6)
        print("\nZero-input output comparison y(t):")
        for t in times:
            y1 = C @ expm(A * t) @ x0_observable
            y2 = C @ expm(A * t) @ (x0_observable + x0_hidden)
            print(f"t={t:4.1f}: y1={y1[0]: .8f}, y2={y2[0]: .8f}, difference={abs(y1[0]-y2[0]):.2e}")
    except ImportError:
        print("\nSciPy not installed; skipping expm-based output comparison.")


if __name__ == "__main__":
    main()

The C++ version uses Eigen's SVD routines. This is the preferred approach for numerical rank and null-space computations because exact Gaussian elimination is unreliable for nearly unobservable systems.

Chapter19_Lesson2.cpp

// Chapter19_Lesson2.cpp
// Observable/unobservable subspaces for an LTI pair (A,C).
// Dependency: Eigen 3 (https://eigen.tuxfamily.org). Compile example:
//   g++ -std=c++17 Chapter19_Lesson2.cpp -I /path/to/eigen -O2 -o Chapter19_Lesson2

#include <Eigen/Dense>
#include <iostream>
#include <iomanip>

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

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

int numericalRank(const MatrixXd& M, double tol = -1.0) {
    Eigen::JacobiSVD<MatrixXd> svd(M, Eigen::ComputeThinU | Eigen::ComputeThinV);
    VectorXd s = svd.singularValues();
    if (tol < 0.0) {
        tol = std::max(M.rows(), M.cols()) * std::numeric_limits<double>::epsilon() * (s.size() ? s(0) : 1.0);
    }
    int r = 0;
    for (int i = 0; i < s.size(); ++i) {
        if (s(i) > tol) { ++r; }
    }
    return r;
}

MatrixXd nullspace(const MatrixXd& M, double tol = -1.0) {
    Eigen::JacobiSVD<MatrixXd> svd(M, Eigen::ComputeFullU | Eigen::ComputeFullV);
    VectorXd s = svd.singularValues();
    if (tol < 0.0) {
        tol = std::max(M.rows(), M.cols()) * std::numeric_limits<double>::epsilon() * (s.size() ? s(0) : 1.0);
    }
    int r = 0;
    for (int i = 0; i < s.size(); ++i) {
        if (s(i) > tol) { ++r; }
    }
    return svd.matrixV().rightCols(M.cols() - r);
}

MatrixXd rowSpaceBasis(const MatrixXd& M, double tol = -1.0) {
    Eigen::JacobiSVD<MatrixXd> svd(M, Eigen::ComputeFullU | Eigen::ComputeFullV);
    VectorXd s = svd.singularValues();
    if (tol < 0.0) {
        tol = std::max(M.rows(), M.cols()) * std::numeric_limits<double>::epsilon() * (s.size() ? s(0) : 1.0);
    }
    int r = 0;
    for (int i = 0; i < s.size(); ++i) {
        if (s(i) > tol) { ++r; }
    }
    return svd.matrixV().leftCols(r);
}

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, 0.0, 1.0,
         0.0, 0.0, -5.0, -1.0;

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

    MatrixXd O = observabilityMatrix(A, C);
    MatrixXd Qu = nullspace(O);
    MatrixXd Qo = rowSpaceBasis(O);

    MatrixXd Q(A.rows(), A.cols());
    Q << Qo, Qu;  // For this example, dimensions add to n.

    MatrixXd Abar = Q.transpose() * A * Q;
    MatrixXd Cbar = C * Q;

    MatrixXd I = MatrixXd::Identity(A.rows(), A.cols());
    double invarianceResidual = ((I - Qu * Qu.transpose()) * A * Qu).norm();

    std::cout << std::setprecision(6) << std::fixed;
    std::cout << "Observability matrix O_n:\n" << O << "\n\n";
    std::cout << "rank(O_n) = " << numericalRank(O) << " out of n = " << A.rows() << "\n\n";
    std::cout << "Basis for unobservable subspace ker(O_n):\n" << Qu << "\n\n";
    std::cout << "A-invariance residual = " << invarianceResidual << "\n\n";
    std::cout << "Abar = Q^T A Q:\n" << Abar << "\n\n";
    std::cout << "Cbar = C Q:\n" << Cbar << "\n";

    return 0;
}

For teaching purposes, the Java version implements reduced row-echelon form directly. In production numerical work, replace this with EJML or Apache Commons Math SVD/QR routines.

Chapter19_Lesson2.java

// Chapter19_Lesson2.java
// Observable/unobservable subspaces using from-scratch RREF operations.
// Compile/run:
//   javac Chapter19_Lesson2.java
//   java Chapter19_Lesson2

public class Chapter19_Lesson2 {
    static final double EPS = 1e-10;

    static double[][] multiply(double[][] A, double[][] B) {
        int m = A.length, n = B[0].length, inner = B.length;
        double[][] C = new double[m][n];
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                double sum = 0.0;
                for (int k = 0; k < inner; k++) sum += A[i][k] * B[k][j];
                C[i][j] = sum;
            }
        }
        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[][] observabilityMatrix(double[][] A, double[][] C) {
        int n = A.length, p = C.length;
        double[][] O = new double[p * n][n];
        double[][] Ak = eye(n);
        for (int block = 0; block < n; block++) {
            double[][] CAk = multiply(C, Ak);
            for (int i = 0; i < p; i++) {
                System.arraycopy(CAk[i], 0, O[block * p + i], 0, n);
            }
            Ak = multiply(Ak, A);
        }
        return O;
    }

    static double[][] rref(double[][] M) {
        int rows = M.length, cols = M[0].length;
        double[][] R = new double[rows][cols];
        for (int i = 0; i < rows; i++) System.arraycopy(M[i], 0, R[i], 0, cols);

        int lead = 0;
        for (int r = 0; r < rows && lead < cols; r++) {
            int i = r;
            while (i < rows && Math.abs(R[i][lead]) < EPS) i++;
            if (i == rows) { lead++; r--; continue; }
            double[] temp = R[r]; R[r] = R[i]; R[i] = temp;
            double pivot = R[r][lead];
            for (int j = 0; j < cols; j++) R[r][j] /= pivot;
            for (int rr = 0; rr < rows; rr++) {
                if (rr == r) continue;
                double factor = R[rr][lead];
                for (int j = 0; j < cols; j++) R[rr][j] -= factor * R[r][j];
            }
            lead++;
        }
        return R;
    }

    static int rank(double[][] M) {
        double[][] R = rref(M);
        int rank = 0;
        for (double[] row : R) {
            boolean nonzero = false;
            for (double v : row) if (Math.abs(v) > EPS) { nonzero = true; break; }
            if (nonzero) rank++;
        }
        return rank;
    }

    static double[][] nullspace(double[][] M) {
        double[][] R = rref(M);
        int rows = R.length, cols = R[0].length;
        boolean[] pivot = new boolean[cols];
        int[] pivotColForRow = new int[rows];
        for (int i = 0; i < rows; i++) pivotColForRow[i] = -1;

        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                if (Math.abs(R[i][j] - 1.0) < EPS) {
                    boolean isPivot = true;
                    for (int k = 0; k < rows; k++) {
                        if (k != i && Math.abs(R[k][j]) > EPS) { isPivot = false; break; }
                    }
                    if (isPivot) { pivot[j] = true; pivotColForRow[i] = j; break; }
                }
            }
        }

        int freeCount = 0;
        for (boolean b : pivot) if (!b) freeCount++;
        double[][] N = new double[cols][freeCount];
        int colIndex = 0;
        for (int free = 0; free < cols; free++) {
            if (pivot[free]) continue;
            N[free][colIndex] = 1.0;
            for (int i = 0; i < rows; i++) {
                int pc = pivotColForRow[i];
                if (pc >= 0) N[pc][colIndex] = -R[i][free];
            }
            colIndex++;
        }
        return N;
    }

    static double frobeniusNorm(double[][] M) {
        double s = 0.0;
        for (double[] row : M) for (double v : row) s += v * v;
        return Math.sqrt(s);
    }

    static void printMatrix(String name, double[][] M) {
        System.out.println(name + ":");
        for (double[] row : M) {
            for (double v : row) System.out.printf("%10.5f ", v);
            System.out.println();
        }
        System.out.println();
    }

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

        double[][] O = observabilityMatrix(A, C);
        double[][] N = nullspace(O);
        double[][] OAN = multiply(O, multiply(A, N));

        printMatrix("Observability matrix O_n", O);
        System.out.println("rank(O_n) = " + rank(O) + " out of n = " + A.length + "\n");
        printMatrix("Basis for unobservable subspace ker(O_n)", N);
        System.out.printf("Invariance check ||O_n A N||_F = %.5e\n", frobeniusNorm(OAN));
    }
}

MATLAB provides both high-level control functions and base linear algebra tools. The script below constructs the observability matrix from scratch, compares with obsv when available, and optionally creates a Simulink model containing a continuous State-Space block.

Chapter19_Lesson2.m

% Chapter19_Lesson2.m
% Observable/unobservable subspaces for a continuous-time LTI system.
% Toolboxes/libraries covered:
%   - Control System Toolbox: obsv, ss, initial.
%   - Base MATLAB: matrix powers, rank, null, orth.
%   - Simulink: optional State-Space block model construction.

clear; clc;

A = [ 0  1  0  0;
     -2 -3  0  0;
      0  0  0  1;
      0  0 -5 -1];
B = [0; 1; 0; 1];
C = [1 0 0 0];
D = 0;

n = size(A,1);

% From-scratch observability matrix
O = [];
Ak = eye(n);
for k = 0:n-1
    O = [O; C*Ak]; %#ok<AGROW>
    Ak = Ak*A;
end

fprintf('Observability matrix O_n =\n'); disp(O);
fprintf('rank(O_n) = %d out of n = %d\n\n', rank(O), n);

% Control System Toolbox equivalent, when available:
if exist('obsv','file') == 2
    fprintf('Control System Toolbox obsv(A,C) matches from-scratch O: %d\n\n', norm(obsv(A,C)-O,'fro') < 1e-12);
end

Qu = null(O, 'r');        % rational-looking basis for ker(O)
Qo = orth(O');            % basis for row space of O
Q  = [Qo Qu];             % x = Q z, for this example Q is nonsingular/orthonormal enough

fprintf('Basis Qu for unobservable subspace ker(O_n):\n'); disp(Qu);
fprintf('A-invariance residual ||(I-Qu*pinv(Qu))*A*Qu||_F = %.3e\n\n', norm((eye(n)-Qu*pinv(Qu))*A*Qu, 'fro'));

Abar = Q\A*Q;             % equivalent to inv(Q)*A*Q, better numerically
Bbar = Q\B;
Cbar = C*Q;

fprintf('Abar = inv(Q)*A*Q =\n'); disp(Abar);
fprintf('Bbar = inv(Q)*B =\n'); disp(Bbar);
fprintf('Cbar = C*Q =\n'); disp(Cbar);

% Demonstrate output indistinguishability using initial response.
if exist('ss','file') == 2
    sys = ss(A,B,C,D);
    t = linspace(0,5,200);
    x_observable = [1;0;0;0];
    x_hidden = [0;0;1;-1];
    y1 = initial(sys, x_observable, t);
    y2 = initial(sys, x_observable + x_hidden, t);
    fprintf('max |y1(t)-y2(t)| over t-grid = %.3e\n\n', max(abs(y1-y2)));
end

% Optional Simulink model construction.  Requires Simulink license.
if exist('simulink','file') == 2
    model = 'Chapter19_Lesson2_Simulink';
    if bdIsLoaded(model), close_system(model, 0); end
    new_system(model);
    add_block('simulink/Continuous/State-Space', [model '/State-Space Pair']);
    set_param([model '/State-Space Pair'], 'A', 'A', 'B', 'B', 'C', 'C', 'D', 'D');
    add_block('simulink/Sources/Step', [model '/Input Step']);
    add_block('simulink/Sinks/Scope', [model '/Output Scope']);
    add_line(model, 'Input Step/1', 'State-Space Pair/1');
    add_line(model, 'State-Space Pair/1', 'Output Scope/1');
    save_system(model);
    fprintf('Created optional Simulink model: %s.slx\n', model);
end

Mathematica is particularly useful when the matrices have exact rational entries and symbolic verification is desired. The notebook-style file computes the null space and simplifies the output difference exactly.

Chapter19_Lesson2.nb

Notebook[{
Cell["Chapter19_Lesson2.nb", "Title"],
Cell["Observable/unobservable subspaces for an LTI system.", "Text"],
Cell[BoxData['(* Chapter19_Lesson2.nb *)\n(* Observable/unobservable subspaces for a continuous-time LTI system. *)\n\nClearAll[ObservabilityMatrix, NumericalRank, UnobservableBasis];\n\nObservabilityMatrix[A_, C_] := Module[{n = Length[A]},\n  Join @@ Table[C . MatrixPower[A, k], {k, 0, n - 1}]\n];\n\nNumericalRank[M_, tol_: 10^-10] := MatrixRank[N[M], Tolerance -> tol];\n\nUnobservableBasis[A_, C_] := NullSpace[ObservabilityMatrix[A, C]];\n\nA = { {0, 1, 0, 0}, {-2, -3, 0, 0}, {0, 0, 0, 1}, {0, 0, -5, -1} };\nB = { {0}, {1}, {0}, {1} };\nCmat = { {1, 0, 0, 0} };\nDmat = { {0} };\n\nOmat = ObservabilityMatrix[A, Cmat];\nrank = NumericalRank[Omat];\nQuRows = UnobservableBasis[A, Cmat];\nQu = Transpose[QuRows];\nQo = Orthogonalize[Transpose[Omat]];\nQ = Transpose[Join[Qo, QuRows]];\n\nAbar = Inverse[Q].A.Q;\nBbar = Inverse[Q].B;\nCbar = Cmat.Q;\n\nPrint["Observability matrix O_n:"];\nPrint[MatrixForm[Omat]];\nPrint["rank(O_n) = ", rank, " out of n = ", Length[A]];\nPrint["Basis for unobservable subspace ker(O_n):"];\nPrint[MatrixForm[Qu]];\nPrint["Abar = inv(Q) A Q:"];\nPrint[MatrixForm[Simplify[Abar]]];\nPrint["Cbar = C Q:"];\nPrint[MatrixForm[Simplify[Cbar]]];\n\n(* Output indistinguishability check: C exp(A t) x0 is unchanged by adding an\n   initial condition from the unobservable subspace. *)\nxObservable = { {1}, {0}, {0}, {0} };\nxHidden = { {0}, {0}, {1}, {-1} };\ny1[t_] := Cmat.MatrixExp[A t].xObservable;\ny2[t_] := Cmat.MatrixExp[A t].(xObservable + xHidden);\nPrint["Simplified output difference:"];\nPrint[Simplify[y1[t] - y2[t]]];\n'], "Input"]
}]

Problem 1 (Computing the Unobservable Subspace): For

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

compute \( \mathcal{O}_4 \), its rank, and \( \mathscr{N}_{\bar{o} } \).

Solution: The observability matrix is

\[ \mathcal{O}_4=\begin{bmatrix} 1&0&0&0\\ 0&1&0&0\\ -2&-3&0&0\\ 6&7&0&0 \end{bmatrix}. \]

Its first two columns are independent and the last two columns are zero, so \( \operatorname{rank}(\mathcal{O}_4)=2 \). Hence

\[ \mathscr{N}_{\bar{o} }=\ker(\mathcal{O}_4)= \operatorname{span}\left\{ \begin{bmatrix}0\\0\\1\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\1\end{bmatrix} \right\}. \]

Problem 2 (A-Invariance): Prove that the subspace found in Problem 1 is invariant under \( \mathbf{A} \).

Solution: Let

\[ \mathbf{v}=\begin{bmatrix}0\\0\\\alpha\\\beta\end{bmatrix}. \]

Then

\[ \mathbf{A}\mathbf{v}=\begin{bmatrix}0\\0\\\beta\\-5\alpha-\beta\end{bmatrix}, \]

which is again in the span of the third and fourth coordinate vectors. Therefore \( \mathbf{A}\mathscr{N}_{\bar{o} }\subseteq \mathscr{N}_{\bar{o} } \).

Problem 3 (Coordinate Block Structure): Use the standard basis ordering \( \mathbf{T}=[\mathbf{e}_1\;\mathbf{e}_2\;\mathbf{e}_3\;\mathbf{e}_4] \) for the same system. Identify the observable and unobservable blocks.

Solution: In this example the original coordinates already separate the measured and hidden modes:

\[ \mathbf{A}=\begin{bmatrix} \mathbf{A}_{oo}&\mathbf{0}\\ \mathbf{0}&\mathbf{A}_{uu} \end{bmatrix},\qquad \mathbf{A}_{oo}=\begin{bmatrix}0&1\\-2&-3\end{bmatrix},\qquad \mathbf{A}_{uu}=\begin{bmatrix}0&1\\-5&-1\end{bmatrix}, \]

\[ \mathbf{C}=\begin{bmatrix}\mathbf{C}_o&\mathbf{0}\end{bmatrix}, \qquad \mathbf{C}_o=\begin{bmatrix}1&0\end{bmatrix}. \]

The unobservable block \( \mathbf{A}_{uu} \) has internal dynamics but no output signature.

Problem 4 (PBH Hidden Mode): Suppose \( \mathbf{A}\mathbf{w}=\lambda\mathbf{w} \) and \( \mathbf{C}\mathbf{w}=\mathbf{0} \). Show that the initial condition \( \mathbf{x}_0=\mathbf{w} \) is unobservable.

Solution: The zero-input response is

\[ \mathbf{y}_0(t)=\mathbf{C}e^{\mathbf{A}t}\mathbf{w} =\mathbf{C}e^{\lambda t}\mathbf{w} =e^{\lambda t}\mathbf{C}\mathbf{w}=\mathbf{0}. \]

Thus the eigenmode is hidden from the output. This proves the intuitive direction of the PBH observability test.

Problem 5 (Detectability Preview): A pair has one unobservable eigenvalue \( \lambda=-4 \) and all observable modes are arbitrary. Is the pair detectable? What if the hidden eigenvalue is \( \lambda=2 \)?

Solution: A pair is detectable when every unobservable mode is asymptotically stable. Since \( \operatorname{Re}(-4)<0 \), the first pair is detectable. Since \( \operatorname{Re}(2)>0 \), the second pair is not detectable. This preview connects the present lesson to the later observer-design requirement that unstable hidden modes cannot be tolerated.

The unobservable subspace is the set of initial states that generate zero output under zero input. Algebraically, it is \( \ker(\mathcal{O}_n) \). Its orthogonal complement is \( \operatorname{range}(\mathcal{O}_n^T) \), and the truly reconstructible object is the state modulo the hidden subspace. Because the hidden subspace is \( \mathbf{A} \)-invariant, a similarity transformation can separate the observable and unobservable dynamics into block form. This prepares the exact block structure of the Kalman decomposition in Lesson 3.