Consider a continuous-time linear time-invariant system

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

Here \( \mathbf{x}(t)\in\mathbb{R}^n \), \( \mathbf{u}(t)\in\mathbb{R}^m \), and \( \mathbf{y}(t)\in\mathbb{R}^p \). The observability question asks whether the unknown initial condition \( \mathbf{x}(0) \) can be uniquely determined from the measured output over a finite time interval, assuming the input \( \mathbf{u}(t) \) is known.

Since the known input contribution can be subtracted from the output, the core algebraic test depends only on the pair \( (\mathbf{A},\mathbf{C}) \). Therefore, for the observability test we examine the zero-input response

\[ \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}(0),\qquad \mathbf{y}(t)=\mathbf{C}e^{\mathbf{A}t}\mathbf{x}(0). \]

Observability means that no two distinct initial states can produce the same output signal. Equivalently, the only initial state that produces identically zero output is the zero state:

\[ \mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0=\mathbf{0} \text{ for all } t\ge 0 \quad \Longrightarrow \quad \mathbf{x}_0=\mathbf{0}. \]

flowchart TD
  A["State model: xdot = A x + B u"] --> B["Known input contribution can be removed"]
  B --> C["Study zero-input output: y(t) = C exp(A t) x0"]
  C --> D["Different x0 values produce output signals"]
  D --> E["Can output uniquely identify x0?"]
  E -->|yes| F["Observable pair (A, C)"]
  E -->|no| G["Nonzero hidden initial-state direction exists"]
  G --> H["Hidden direction lies in null space of O_n"]
        

The zero-input output is analytic because it is generated by the matrix exponential. Differentiating the output at \( t=0 \) gives

\[ \begin{aligned} \mathbf{y}(0) &= \mathbf{C}\mathbf{x}_0,\\ \dot{\mathbf{y} }(0) &= \mathbf{C}\mathbf{A}\mathbf{x}_0,\\ \mathbf{y}^{(2)}(0) &= \mathbf{C}\mathbf{A}^2\mathbf{x}_0,\\ &\,\vdots\\ \mathbf{y}^{(k)}(0) &= \mathbf{C}\mathbf{A}^k\mathbf{x}_0. \end{aligned} \]

Stacking the first \( n \) derivative equations produces the Kalman 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}. \]

If the stacked derivative vector is denoted by \( \boldsymbol{\eta}(0) \), then

\[ \boldsymbol{\eta}(0)= \begin{bmatrix} \mathbf{y}(0)\\ \dot{\mathbf{y} }(0)\\ \vdots\\ \mathbf{y}^{(n-1)}(0) \end{bmatrix} = \mathcal{O}_n\mathbf{x}_0. \]

Hence, recovering the initial state from ideal output derivatives is a linear algebra problem: solve \( \mathcal{O}_n\mathbf{x}_0=\boldsymbol{\eta}(0) \). Unique recovery is possible exactly when \( \mathcal{O}_n \) has full column rank.

The central finite-dimensional test is:

\[ \boxed{(\mathbf{A},\mathbf{C})\text{ is observable } \iff \operatorname{rank}(\mathcal{O}_n)=n.} \]

Since \( \mathcal{O}_n \) has exactly \( n \) columns, the condition \( \operatorname{rank}(\mathcal{O}_n)=n \) means that the columns of \( \mathcal{O}_n \) are linearly independent.

Proof: Full Rank Implies Observability

Suppose \( \operatorname{rank}(\mathcal{O}_n)=n \). If an initial state \( \mathbf{x}_0 \) produces zero output for all time, then all derivatives at \( t=0 \) are zero:

\[ \mathbf{C}\mathbf{A}^k\mathbf{x}_0=\mathbf{0}, \qquad k=0,1,\dots,n-1. \]

Stacking these equations gives \( \mathcal{O}_n\mathbf{x}_0=\mathbf{0} \). Since \( \mathcal{O}_n \) has full column rank, its null space is trivial:

\[ \mathcal{N}(\mathcal{O}_n)=\{\mathbf{0}\}. \]

Therefore \( \mathbf{x}_0=\mathbf{0} \). The only state producing zero output is the zero state, so the pair is observable.

Proof: Lack of Full Rank Implies Nonobservability

Suppose \( \operatorname{rank}(\mathcal{O}_n)<n \). Then there exists a nonzero vector \( \mathbf{v}\ne\mathbf{0} \) such that

\[ \mathcal{O}_n\mathbf{v}=\mathbf{0}, \qquad \mathbf{C}\mathbf{A}^k\mathbf{v}=\mathbf{0}, \quad k=0,1,\dots,n-1. \]

By the Cayley-Hamilton theorem, every higher power \( \mathbf{A}^k \) for \( k\ge n \) is a linear combination of \( \mathbf{I},\mathbf{A},\dots,\mathbf{A}^{n-1} \). Therefore

\[ \mathbf{C}\mathbf{A}^k\mathbf{v}=\mathbf{0}, \qquad \text{for all } k\ge 0. \]

Using the matrix exponential series,

\[ \mathbf{C}e^{\mathbf{A}t}\mathbf{v} = \mathbf{C}\left(\sum_{k=0}^{\infty} \frac{\mathbf{A}^k t^k}{k!}\right)\mathbf{v} = \sum_{k=0}^{\infty} \frac{t^k}{k!}\mathbf{C}\mathbf{A}^k\mathbf{v} = \mathbf{0}. \]

Thus the nonzero initial state \( \mathbf{x}_0=\mathbf{v} \) produces zero output for all time, so the system is not observable.

The null space of the observability matrix is the set of all initial state components invisible from the measured output:

\[ \mathcal{N}(\mathcal{O}_n)= \left\{\mathbf{x}_0\in\mathbb{R}^n: \mathbf{C}\mathbf{A}^k\mathbf{x}_0=\mathbf{0}, \; k=0,1,\dots,n-1\right\}. \]

The dimension of this unobservable subspace is obtained from the rank-nullity theorem:

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

Therefore, if \( \operatorname{rank}(\mathcal{O}_n)=r \), then \( r \) independent state directions are visible through the output derivative sequence, while \( n-r \) independent directions are hidden.

The unobservable subspace is invariant under \( \mathbf{A} \). If \( \mathbf{v}\in\mathcal{N}(\mathcal{O}_n) \), then

\[ \mathbf{C}\mathbf{A}^k(\mathbf{A}\mathbf{v}) =\mathbf{C}\mathbf{A}^{k+1}\mathbf{v}=\mathbf{0}, \qquad k=0,1,\dots,n-2. \]

The missing condition at power \( n \) follows again from Cayley-Hamilton. Thus \( \mathbf{A}\mathcal{N}(\mathcal{O}_n) \subseteq \mathcal{N}(\mathcal{O}_n) \). This is why unobservable modes can be treated as an internally evolving hidden subsystem in later decomposition results.

Observability is a property of the input-output/state description, not of the particular coordinate basis. Let \( \mathbf{x}=\mathbf{T}\mathbf{z} \) with nonsingular \( \mathbf{T} \). Then

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

The transformed pair is \( (\tilde{\mathbf{A} },\tilde{\mathbf{C} }) =(\mathbf{T}^{-1}\mathbf{A}\mathbf{T},\mathbf{C}\mathbf{T}) \). Its observability matrix is

\[ \tilde{\mathcal{O} }_n = \begin{bmatrix} \mathbf{C}\mathbf{T}\\ \mathbf{C}\mathbf{T}(\mathbf{T}^{-1}\mathbf{A}\mathbf{T})\\ \vdots\\ \mathbf{C}\mathbf{T}(\mathbf{T}^{-1}\mathbf{A}\mathbf{T})^{n-1} \end{bmatrix} = \begin{bmatrix} \mathbf{C}\\ \mathbf{C}\mathbf{A}\\ \vdots\\ \mathbf{C}\mathbf{A}^{n-1} \end{bmatrix}\mathbf{T} = \mathcal{O}_n\mathbf{T}. \]

Since multiplication by a nonsingular matrix does not change column rank,

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

In numerical work, exact rank is replaced by numerical rank. A standard singular-value criterion is

\[ \operatorname{rank}_{\epsilon}(\mathcal{O}_n) = \#\left\{\sigma_i(\mathcal{O}_n): \sigma_i(\mathcal{O}_n)>\epsilon\right\}, \qquad \epsilon = \max(pn,n)\,\varepsilon_{\text{mach} }\,\sigma_{\max}. \]

Poor scaling can make a mathematically observable system appear nearly unobservable. For this reason, rank should often be interpreted together with singular values, condition numbers, and physically meaningful state scaling.

flowchart TD
  A["Given A and C"] --> B["Build O_n = stack(C, C A, ..., C A^(n-1))"]
  B --> C["Compute singular values of O_n"]
  C --> D["Choose tolerance based on scale"]
  D --> E["Numerical rank r"]
  E --> F["r = n?"]
  F -->|yes| G["Observable in numerical test"]
  F -->|no| H["Hidden directions: null space of O_n"]
  H --> I["Inspect singular vectors and sensor placement"]
        

Example 1: Position Measurement Reveals Velocity Through Dynamics

Consider \( \mathbf{A}=\begin{bmatrix}0&1\\0&0\end{bmatrix} \) and \( \mathbf{C}=\begin{bmatrix}1&0\end{bmatrix} \). Then

\[ \mathcal{O}_2 = \begin{bmatrix} \mathbf{C}\\ \mathbf{C}\mathbf{A} \end{bmatrix} = \begin{bmatrix} 1&0\\ 0&1 \end{bmatrix}, \qquad \operatorname{rank}(\mathcal{O}_2)=2. \]

The system is observable. Physically, measuring position \( x_1(t) \) also reveals velocity \( x_2(t) \) through the derivative \( \dot{y}(0)=x_2(0) \).

Example 2: A Hidden Decoupled Mode

Now consider \( \mathbf{A}=\begin{bmatrix}0&0\\0&-2\end{bmatrix} \) and \( \mathbf{C}=\begin{bmatrix}1&0\end{bmatrix} \). Then

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

The second state is unobservable because neither \( y(t) \) nor its derivatives contain \( x_2(0) \). The unobservable subspace is

\[ \mathcal{N}(\mathcal{O}_2)= \operatorname{span}\left\{\begin{bmatrix}0\\1\end{bmatrix}\right\}. \]

For a finite-dimensional continuous-time LTI model, the practical algorithm is:

1. Verify that \( \mathbf{A}\in\mathbb{R}^{n\times n} \) and \( \mathbf{C}\in\mathbb{R}^{p\times n} \).
2. Construct \( \mathcal{O}_n=[\mathbf{C}^T,(\mathbf{C}\mathbf{A})^T,\dots,(\mathbf{C}\mathbf{A}^{n-1})^T]^T \).
3. Compute rank using exact symbolic algebra for small exact matrices, or SVD/QR for numerical matrices.
4. Declare the pair observable if \( \operatorname{rank}(\mathcal{O}_n)=n \).
5. If not observable, compute \( \mathcal{N}(\mathcal{O}_n) \) to identify hidden initial-state directions.

Useful libraries include NumPy/SciPy and python-control in Python, Eigen in C++, EJML or Apache Commons Math in Java, Control System Toolbox in MATLAB, Simulink linearization tools for models built as block diagrams, and Wolfram Mathematica for exact or symbolic rank calculations.

Chapter14_Lesson1.py

"""
Chapter14_Lesson1.py

Kalman observability matrix and rank condition for continuous-time LTI systems.

Model:
    x_dot = A x + B u
    y     = C x + D u

For observability of the initial state, only (A, C) matters.
"""

import numpy as np


def observability_matrix(A: np.ndarray, C: np.ndarray, n: int | None = None) -> np.ndarray:
    """Return O_n = [C; C A; ...; C A^(n-1)]."""
    A = np.asarray(A, dtype=float)
    C = np.asarray(C, dtype=float)

    if A.ndim != 2 or A.shape[0] != A.shape[1]:
        raise ValueError("A must be square.")
    if C.ndim != 2 or C.shape[1] != A.shape[0]:
        raise ValueError("C must have the same number of columns as A has states.")

    n_states = A.shape[0]
    n = n_states if n is None else int(n)
    blocks = []
    Ak = np.eye(n_states)
    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) -> tuple[int, np.ndarray, float]:
    """SVD-based numerical rank."""
    M = np.asarray(M, dtype=float)
    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 0.0)
    rank = int(np.sum(singular_values > tol))
    return rank, singular_values, float(tol)


def null_space_svd(M: np.ndarray, tol: float | None = None) -> np.ndarray:
    """Approximate null space basis using SVD."""
    M = np.asarray(M, dtype=float)
    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 0.0)
    rank = int(np.sum(S > tol))
    return Vt[rank:, :].T


def analyze_observability(A: np.ndarray, C: np.ndarray, name: str = "system") -> None:
    """Print observability matrix, rank, singular values, and unobservable directions."""
    O = observability_matrix(A, C)
    rank, singular_values, tol = numerical_rank(O)
    n = np.asarray(A).shape[0]
    N = null_space_svd(O, tol)

    print(f"\n{name}")
    print("-" * len(name))
    print("A =\n", np.asarray(A, dtype=float))
    print("C =\n", np.asarray(C, dtype=float))
    print("Observability matrix O_n =\n", O)
    print("Singular values =", singular_values)
    print("Tolerance =", tol)
    print(f"rank(O_n) = {rank} of n = {n}")
    print("Observable?", rank == n)
    if rank < n:
        print("Basis for unobservable initial-state directions:")
        print(N)


if __name__ == "__main__":
    # Example 1: position is measured in a double-integrator-like system.
    A1 = np.array([[0.0, 1.0],
                   [0.0, 0.0]])
    C1 = np.array([[1.0, 0.0]])
    analyze_observability(A1, C1, "Example 1: observable double integrator")

    # Example 2: second state does not affect the output or its derivatives.
    A2 = np.array([[0.0, 0.0],
                   [0.0, -2.0]])
    C2 = np.array([[1.0, 0.0]])
    analyze_observability(A2, C2, "Example 2: unobservable second mode")

    # Example 3: same input matrix B would not change observability.
    B = np.array([[0.0],
                  [1.0]])
    print("\nB does not enter the Kalman observability rank test for (A, C).")
    print("B =\n", B)

Chapter14_Lesson1.cpp

/*
Chapter14_Lesson1.cpp

Kalman observability matrix and rank condition using Eigen.

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

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

Eigen::MatrixXd observabilityMatrix(const Eigen::MatrixXd& A, const Eigen::MatrixXd& C) {
    const int n = static_cast<int>(A.rows());
    const int p = static_cast<int>(C.rows());

    Eigen::MatrixXd O(p * n, n);
    Eigen::MatrixXd Ak = Eigen::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 Eigen::MatrixXd& M, double* toleranceOut = nullptr) {
    Eigen::JacobiSVD<Eigen::MatrixXd> svd(M, Eigen::ComputeThinU | Eigen::ComputeThinV);
    Eigen::VectorXd s = svd.singularValues();

    double maxSingular = (s.size() > 0) ? s(0) : 0.0;
    double tol = std::max(M.rows(), M.cols()) * std::numeric_limits<double>::epsilon() * maxSingular;

    int rank = 0;
    for (int i = 0; i < s.size(); ++i) {
        if (s(i) > tol) {
            ++rank;
        }
    }

    if (toleranceOut != nullptr) {
        *toleranceOut = tol;
    }
    return rank;
}

Eigen::MatrixXd nullSpaceSVD(const Eigen::MatrixXd& M) {
    Eigen::JacobiSVD<Eigen::MatrixXd> svd(M, Eigen::ComputeFullV);
    Eigen::VectorXd s = svd.singularValues();
    double tol = std::max(M.rows(), M.cols()) * std::numeric_limits<double>::epsilon() * (s.size() > 0 ? s(0) : 0.0);

    int rank = 0;
    for (int i = 0; i < s.size(); ++i) {
        if (s(i) > tol) {
            ++rank;
        }
    }

    return svd.matrixV().rightCols(M.cols() - rank);
}

void analyzeObservability(const Eigen::MatrixXd& A, const Eigen::MatrixXd& C, const std::string& name) {
    Eigen::MatrixXd O = observabilityMatrix(A, C);
    double tol = 0.0;
    int rank = numericalRank(O, &tol);

    Eigen::JacobiSVD<Eigen::MatrixXd> svd(O);
    std::cout << "\n" << name << "\n";
    std::cout << std::string(name.size(), '-') << "\n";
    std::cout << "A =\n" << A << "\n";
    std::cout << "C =\n" << C << "\n";
    std::cout << "Observability matrix O_n =\n" << O << "\n";
    std::cout << "Singular values = " << svd.singularValues().transpose() << "\n";
    std::cout << "Tolerance = " << tol << "\n";
    std::cout << "rank(O_n) = " << rank << " of n = " << A.rows() << "\n";
    std::cout << "Observable? " << (rank == A.rows() ? "true" : "false") << "\n";

    if (rank < A.rows()) {
        std::cout << "Basis for unobservable initial-state directions:\n";
        std::cout << nullSpaceSVD(O) << "\n";
    }
}

int main() {
    Eigen::MatrixXd A1(2, 2);
    A1 << 0.0, 1.0,
          0.0, 0.0;
    Eigen::MatrixXd C1(1, 2);
    C1 << 1.0, 0.0;
    analyzeObservability(A1, C1, "Example 1: observable double integrator");

    Eigen::MatrixXd A2(2, 2);
    A2 << 0.0, 0.0,
          0.0, -2.0;
    Eigen::MatrixXd C2(1, 2);
    C2 << 1.0, 0.0;
    analyzeObservability(A2, C2, "Example 2: unobservable second mode");

    return 0;
}

Chapter14_Lesson1.java

/*
Chapter14_Lesson1.java

Kalman observability matrix and rank condition using EJML.

Maven dependency example:
  <dependency>
    <groupId>org.ejml</groupId>
    <artifactId>ejml-simple</artifactId>
    <version>0.43.1</version>
  </dependency>
*/

import org.ejml.simple.SimpleMatrix;
import org.ejml.simple.SimpleSVD;

public class Chapter14_Lesson1 {

    public static SimpleMatrix observabilityMatrix(SimpleMatrix A, SimpleMatrix C) {
        int n = A.numRows();
        int p = C.numRows();

        SimpleMatrix O = new SimpleMatrix(p * n, n);
        SimpleMatrix Ak = SimpleMatrix.identity(n);

        for (int k = 0; k < n; k++) {
            SimpleMatrix block = C.mult(Ak);
            O.insertIntoThis(k * p, 0, block);
            Ak = Ak.mult(A);
        }
        return O;
    }

    public static int numericalRank(SimpleMatrix M) {
        SimpleSVD<SimpleMatrix> svd = M.svd();
        double[] s = svd.getSingularValues();

        double maxSingular = s.length > 0 ? s[0] : 0.0;
        double tol = Math.max(M.numRows(), M.numCols()) * Math.ulp(1.0) * maxSingular;

        int rank = 0;
        for (double value : s) {
            if (value > tol) {
                rank++;
            }
        }
        return rank;
    }

    public static void analyzeObservability(SimpleMatrix A, SimpleMatrix C, String name) {
        SimpleMatrix O = observabilityMatrix(A, C);
        SimpleSVD<SimpleMatrix> svd = O.svd();
        int rank = numericalRank(O);

        System.out.println("\n" + name);
        System.out.println("-".repeat(name.length()));
        System.out.println("A =");
        A.print();
        System.out.println("C =");
        C.print();
        System.out.println("Observability matrix O_n =");
        O.print();
        System.out.print("Singular values = ");
        for (double value : svd.getSingularValues()) {
            System.out.printf("%.8f ", value);
        }
        System.out.println();
        System.out.println("rank(O_n) = " + rank + " of n = " + A.numRows());
        System.out.println("Observable? " + (rank == A.numRows()));
    }

    public static void main(String[] args) {
        SimpleMatrix A1 = new SimpleMatrix(new double[][] {
                {0.0, 1.0},
                {0.0, 0.0}
        });
        SimpleMatrix C1 = new SimpleMatrix(new double[][] {
                {1.0, 0.0}
        });
        analyzeObservability(A1, C1, "Example 1: observable double integrator");

        SimpleMatrix A2 = new SimpleMatrix(new double[][] {
                {0.0, 0.0},
                {0.0, -2.0}
        });
        SimpleMatrix C2 = new SimpleMatrix(new double[][] {
                {1.0, 0.0}
        });
        analyzeObservability(A2, C2, "Example 2: unobservable second mode");
    }
}

Chapter14_Lesson1.m

% Chapter14_Lesson1.m
%
% Kalman observability matrix and rank condition.
% This script uses only basic MATLAB. If Control System Toolbox is available,
% it also compares the result with obsv(A,C).

clear; clc;

A1 = [0 1; 0 0];
C1 = [1 0];

A2 = [0 0; 0 -2];
C2 = [1 0];

analyzeObservability(A1, C1, "Example 1: observable double integrator");
analyzeObservability(A2, C2, "Example 2: unobservable second mode");

% Control System Toolbox comparison:
if exist("obsv", "file") == 2
    fprintf("\nControl System Toolbox check for Example 1:\n");
    disp(obsv(A1, C1));
end

% Simulink workflow note:
% 1. Put A, B, C, D into a State-Space block.
% 2. Linearize the model or use the same matrices directly.
% 3. Run O = obsv(A,C); rank(O).
% Observability of the initial state depends on A and C, not on B or D.

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

function analyzeObservability(A, C, name)
    O = observabilityMatrix(A, C);
    s = svd(O);
    tol = max(size(O)) * eps(max(s));
    r = sum(s > tol);

    fprintf("\n%s\n", name);
    fprintf("%s\n", repmat("-", 1, strlength(name)));
    fprintf("A =\n"); disp(A);
    fprintf("C =\n"); disp(C);
    fprintf("Observability matrix O_n =\n"); disp(O);
    fprintf("Singular values =\n"); disp(s.');
    fprintf("Tolerance = %.4e\n", tol);
    fprintf("rank(O_n) = %d of n = %d\n", r, size(A,1));
    fprintf("Observable? %d\n", r == size(A,1));

    if r < size(A,1)
        fprintf("Basis for unobservable initial-state directions:\n");
        disp(null(O));
    end
end

Chapter14_Lesson1.nb

(* Chapter14_Lesson1.nb *)

ClearAll[ObservabilityMatrix, NumericalRank, AnalyzeObservability];

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

NumericalRank[M_, tol_: Automatic] := Module[
  {s = SingularValueList[N[M]], threshold},
  threshold = If[
    tol === Automatic,
    Max[Dimensions[M]] * $MachineEpsilon * Max[SingularValueList[N[M]]],
    tol
  ];
  Count[s, _?(# > threshold &)]
];

AnalyzeObservability[A_, C_, name_] := Module[
  {O, r, n, ns},
  O = ObservabilityMatrix[A, C];
  r = NumericalRank[O];
  n = Length[A];
  ns = NullSpace[O];

  Print["\n", name];
  Print["A = ", MatrixForm[A]];
  Print["C = ", MatrixForm[C]];
  Print["Observability matrix O_n = ", MatrixForm[O]];
  Print["Singular values = ", SingularValueList[N[O]]];
  Print["rank(O_n) = ", r, " of n = ", n];
  Print["Observable? ", r == n];

  If[r < n,
    Print["Basis for unobservable initial-state directions = ", ns]
  ];
];

A1 = { {0, 1}, {0, 0} };
C1 = { {1, 0} };
AnalyzeObservability[A1, C1, "Example 1: observable double integrator"];

A2 = { {0, 0}, {0, -2} };
C2 = { {1, 0} };
AnalyzeObservability[A2, C2, "Example 2: unobservable second mode"];

Problem 1 (Direct Rank Test): For \( \mathbf{A}=\begin{bmatrix}0&1\\-2&-3\end{bmatrix} \) and \( \mathbf{C}=\begin{bmatrix}1&0\end{bmatrix} \), determine whether the pair is observable.

Solution: Since \( n=2 \),

\[ \mathcal{O}_2= \begin{bmatrix} \mathbf{C}\\ \mathbf{C}\mathbf{A} \end{bmatrix} = \begin{bmatrix} 1&0\\ 0&1 \end{bmatrix}. \]

Hence \( \operatorname{rank}(\mathcal{O}_2)=2=n \), so the pair is observable.

Problem 2 (Output Choice Matters): For the same \( \mathbf{A}=\begin{bmatrix}0&1\\-2&-3\end{bmatrix} \), use \( \mathbf{C}=\begin{bmatrix}0&1\end{bmatrix} \). Is the pair observable?

Solution:

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

Since the determinant is nonzero, \( \operatorname{rank}(\mathcal{O}_2)=2 \). The pair is observable. In this example, measuring either state is enough because the dynamics couple the two states.

Problem 3 (Hidden Mode): Let \( \mathbf{A}=\operatorname{diag}(-1,-4,2) \) and \( \mathbf{C}=\begin{bmatrix}1&0&0\\0&1&0\end{bmatrix} \). Find the rank of the observability matrix and the unobservable subspace.

Solution: Since \( \mathbf{A} \) is diagonal,

\[ \mathcal{O}_3= \begin{bmatrix} 1&0&0\\ 0&1&0\\ -1&0&0\\ 0&-4&0\\ 1&0&0\\ 0&16&0 \end{bmatrix}. \]

The third column is zero and the first two columns are independent. Therefore \( \operatorname{rank}(\mathcal{O}_3)=2 \). The unobservable subspace is

\[ \mathcal{N}(\mathcal{O}_3)= \operatorname{span}\left\{\begin{bmatrix}0\\0\\1\end{bmatrix}\right\}. \]

The state component associated with the third coordinate is never measured and never dynamically coupled into the measured coordinates.

Problem 4 (Similarity Invariance): Let \( \tilde{\mathbf{A} }=\mathbf{T}^{-1}\mathbf{A}\mathbf{T} \) and \( \tilde{\mathbf{C} }=\mathbf{C}\mathbf{T} \), where \( \mathbf{T} \) is nonsingular. Prove that \( (\mathbf{A},\mathbf{C}) \) is observable iff \( (\tilde{\mathbf{A} },\tilde{\mathbf{C} }) \) is observable.

Solution: The transformed observability matrix is

\[ \tilde{\mathcal{O} }_n= \begin{bmatrix} \tilde{\mathbf{C} }\\ \tilde{\mathbf{C} }\tilde{\mathbf{A} }\\ \vdots\\ \tilde{\mathbf{C} }\tilde{\mathbf{A} }^{n-1} \end{bmatrix} = \mathcal{O}_n\mathbf{T}. \]

Because \( \mathbf{T} \) is nonsingular, \( \operatorname{rank}(\mathcal{O}_n\mathbf{T}) =\operatorname{rank}(\mathcal{O}_n) \). Therefore the full-rank condition is preserved by any nonsingular coordinate change.

Problem 5 (Why Only n Blocks Are Needed): Explain why the observability matrix stops at \( \mathbf{C}\mathbf{A}^{n-1} \), rather than requiring infinitely many rows \( \mathbf{C}\mathbf{A}^k \).

Solution: By the Cayley-Hamilton theorem, \( \mathbf{A} \) satisfies its characteristic equation:

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

Therefore \( \mathbf{A}^n \), and hence every higher power \( \mathbf{A}^k \) for \( k\ge n \), is a linear combination of lower powers. Thus rows beyond \( \mathbf{C}\mathbf{A}^{n-1} \) cannot add new independent information about \( \mathbf{x}_0 \).

The Kalman observability matrix converts the question of state reconstruction into a finite-dimensional rank test. For an \( n \)-state LTI system, the pair \( (\mathbf{A},\mathbf{C}) \) is observable iff \( \mathcal{O}_n \) has rank \( n \). The null space of \( \mathcal{O}_n \) contains exactly the initial state directions that cannot be detected from the output. The test is invariant under similarity transformations, but numerical implementations must inspect singular values and scaling.