In Lesson 1, observability was tested using the Kalman observability matrix. For an LTI 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), \]

observability depends only on the pair \( (\mathbf{A},\mathbf{C}) \). The Kalman matrix is

\[ \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}, \qquad (\mathbf{A},\mathbf{C})\text{ observable} \Longleftrightarrow \operatorname{rank}(\mathcal{O}_n)=n. \]

The PBH test gives an equivalent condition using the eigenstructure of \( \mathbf{A} \). Instead of checking all powers of \( \mathbf{A} \), it asks whether every eigenmode of the state matrix is visible through \( \mathbf{C} \).

flowchart TD
  A["Start with pair (A, C)"] --> B["Compute eigenvalues of A"]
  B --> C["For each eigenvalue lambda"]
  C --> D["Build matrix [lambda I - A; C]"]
  D --> E{"Rank equals n?"}
  E -->|"yes for this lambda"| F["Mode is visible"]
  E -->|"no"| G["Unobservable mode exists"]
  F --> H{"All eigenvalues passed?"}
  G --> I["Pair (A, C) is not observable"]
  H -->|"yes"| J["Pair (A, C) is observable"]
  H -->|"no"| I
        

Let \( \mathbf{A}\in\mathbb{R}^{n\times n} \) and \( \mathbf{C}\in\mathbb{R}^{p\times n} \). The PBH test states that \( (\mathbf{A},\mathbf{C}) \) is observable if and only if, for every eigenvalue \( \lambda\in\sigma(\mathbf{A}) \),

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

Since \( \lambda\mathbf{I}-\mathbf{A} \) is singular precisely at eigenvalues, it is unnecessary to test values of \( \lambda \) that are not eigenvalues. If \( \lambda\notin\sigma(\mathbf{A}) \), then \( \lambda\mathbf{I}-\mathbf{A} \) already has full rank \( n \), so adding \( \mathbf{C} \) cannot reduce the rank.

The equivalent left-eigenvector form is often the most interpretable:

\[ \boxed{(\mathbf{A},\mathbf{C})\text{ observable} \Longleftrightarrow \nexists\,\mathbf{q}\ne\mathbf{0},\;\lambda\in\mathbb{C} \text{ such that } \mathbf{q}^*\mathbf{A}=\lambda\mathbf{q}^*,\quad \mathbf{q}^*\mathbf{C}^*=\mathbf{0}^*.} \]

Equivalently, if \( \mathbf{v} \) is a right eigenvector of \( \mathbf{A}^* \), then \( \mathbf{C}\mathbf{v}\ne\mathbf{0} \) must hold for every eigenvector direction. A zero value \( \mathbf{C}\mathbf{v}=\mathbf{0} \) means that this modal state direction produces no output signature.

The unobservable subspace is the null space of the observability matrix:

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

By the Cayley-Hamilton theorem, if the output derivatives vanish for powers \( k=0,1,\dots,n-1 \), then they vanish for all higher powers. Therefore \( \mathbf{x}_0\in\mathcal{N}_o \) produces zero output under the zero-input response:

\[ \mathbf{y}(t)=\mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0 =\mathbf{0} \quad \text{for all admissible times}. \]

Now suppose \( \mathbf{q}^*\mathbf{A}=\lambda\mathbf{q}^* \) and \( \mathbf{q}^*\mathbf{C}^*=\mathbf{0}^* \). Taking conjugate transpose gives \( \mathbf{A}^*\mathbf{q}=\bar{\lambda}\mathbf{q} \) and \( \mathbf{C}\mathbf{q}=\mathbf{0} \). This modal direction is invisible at the output. In rank form, invisibility means that there exists a nonzero vector \( \mathbf{z} \) satisfying

\[ (\lambda\mathbf{I}-\mathbf{A})\mathbf{z}=\mathbf{0}, \qquad \mathbf{C}\mathbf{z}=\mathbf{0}. \]

These two equations can be stacked as

\[ \begin{bmatrix} \lambda\mathbf{I}-\mathbf{A} \\ \mathbf{C} \end{bmatrix}\mathbf{z}=\mathbf{0}. \]

Hence the stacked matrix fails to have full column rank if and only if an eigenvector direction of \( \mathbf{A} \) is also in the null space of \( \mathbf{C} \). Therefore,

\[ \ker\begin{bmatrix}\lambda\mathbf{I}-\mathbf{A}\\ \mathbf{C}\end{bmatrix}=\left\{\mathbf{0}\right\} \quad \Longleftrightarrow \quad \operatorname{rank}\begin{bmatrix}\lambda\mathbf{I}-\mathbf{A}\\ \mathbf{C}\end{bmatrix}=n. \]

This proves the PBH condition. The proof also explains why the test is a modal visibility test rather than merely a mechanical rank calculation.

If \( \mathbf{A}\mathbf{v}=\lambda\mathbf{v} \) and the initial condition is \( \mathbf{x}_0=\alpha\mathbf{v} \), then the zero-input state trajectory is

\[ \mathbf{x}(t)=\alpha e^{\lambda t}\mathbf{v},\qquad \mathbf{y}(t)=\alpha e^{\lambda t}\mathbf{C}\mathbf{v}. \]

Thus, if \( \mathbf{C}\mathbf{v}=\mathbf{0} \), that eigenmode never appears in the output, regardless of its amplitude \( \alpha \). In a diagonalizable system, this gives a particularly direct interpretation: every modal coordinate must be seen by at least one output channel.

flowchart LR
  M1["Eigenmode lambda1"] --> C1["Output matrix C"]
  M2["Eigenmode lambda2"] --> C1
  M3["Eigenmode lambda3"] --> C1
  C1 --> Y["Measured output y"]
  M1 -. "C times v1 nonzero" .-> Y
  M2 -. "C times v2 equals zero" .-> U["Invisible mode"]
  M3 -. "C times v3 nonzero" .-> Y
        

For repeated eigenvalues, it is not enough to check one eigenvector arbitrarily. The PBH rank condition automatically handles eigenspaces of dimension greater than one and defective cases. For each eigenvalue \( \lambda \), the test requires the stacked matrix to have full column rank, so every direction in the eigenspace associated with \( \lambda \) must be output-visible.

Detectability was introduced in Chapter 13 as the weaker property that unobservable modes are allowed only if they are internally stable. The PBH test gives a compact detectability test. For a continuous-time system, \( (\mathbf{A},\mathbf{C}) \) is detectable if all PBH rank failures occur only at eigenvalues whose real parts are negative:

\[ (\mathbf{A},\mathbf{C})\text{ detectable} \Longleftrightarrow \operatorname{rank}\begin{bmatrix}\lambda\mathbf{I}-\mathbf{A}\\ \mathbf{C}\end{bmatrix}=n \quad \text{for every}\quad \lambda\in\sigma(\mathbf{A})\text{ with }\operatorname{Re}(\lambda)\ge 0. \]

The observability test checks all modes; the detectability test checks only modes that can persist or grow. In this lesson, the primary focus remains observability, but the detectability version is important because it shows how PBH separates modal visibility from modal stability.

Example 1: observable second-order system. Let

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

The characteristic polynomial is \( \lambda^2+3\lambda+2=(\lambda+1)(\lambda+2) \), so \( \lambda_1=-1 \) and \( \lambda_2=-2 \). For \( \lambda=-1 \),

\[ \begin{bmatrix}\lambda\mathbf{I}-\mathbf{A}\\\mathbf{C}\end{bmatrix} =\begin{bmatrix}-1 & -1\\ 2 & 2\\ 1 & 0\end{bmatrix}, \qquad \operatorname{rank}=2. \]

For \( \lambda=-2 \),

\[ \begin{bmatrix}\lambda\mathbf{I}-\mathbf{A}\\\mathbf{C}\end{bmatrix} =\begin{bmatrix}-2 & -1\\ 2 & 1\\ 1 & 0\end{bmatrix}, \qquad \operatorname{rank}=2. \]

Since the state dimension is \( n=2 \) and both modal PBH matrices have rank 2, the pair is observable.

Example 2: unobservable diagonal system. Let

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

For \( \lambda=-2 \),

\[ \begin{bmatrix}\lambda\mathbf{I}-\mathbf{A}\\\mathbf{C}\end{bmatrix} =\begin{bmatrix}-1 & 0\\ 0 & 0\\ 1 & 0\end{bmatrix}, \qquad \operatorname{rank}=1<2. \]

The second state direction is an eigenmode of \( \mathbf{A} \) but is not measured by \( \mathbf{C} \). Therefore the pair is not observable. However, because the unobservable eigenvalue is \( -2 \), this specific pair is detectable in continuous time.

The PBH test is algebraically exact, but numerical implementation needs a rank tolerance. A practical algorithm is:

1. Compute eigenvalues \( \lambda_i \) of \( \mathbf{A} \).
2. For each eigenvalue, form \( \mathbf{M}_i=\begin{bmatrix}\lambda_i\mathbf{I}-\mathbf{A}\\\mathbf{C}\end{bmatrix} \).
3. Compute singular values \( \sigma_1\ge\cdots\ge\sigma_n \) of \( \mathbf{M}_i \).
4. Declare full column rank if the smallest singular value is above a tolerance.

\[ \mathbf{M}_i\text{ has full column rank} \Longleftrightarrow \sigma_n(\mathbf{M}_i)>0. \]

In computation, one usually replaces the exact condition by

\[ \sigma_n(\mathbf{M}_i)>\varepsilon\,\max(m,n)\,\sigma_1(\mathbf{M}_i), \]

where \( m=n+p \) is the number of rows of the stacked PBH matrix and \( \varepsilon \) is machine precision or a user-selected engineering tolerance. If the smallest singular value is very small but nonzero, the system is mathematically observable but weakly observable; state reconstruction may be highly sensitive to noise.

Python libraries commonly used for this lesson are NumPy for matrix operations, SciPy for eigenvalue routines, and python-control for state-space models and observability matrices. The implementation below keeps the PBH rank test explicit.

File: Chapter14_Lesson2.py

"""
Chapter14_Lesson2.py
PBH Test for Observability for continuous-time LTI systems.

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

Observability depends only on (A, C). The PBH test says that (A, C) is
observable iff rank([lambda*I - A; C]) = n for every eigenvalue lambda of A.
"""

from __future__ import annotations

import numpy as np
from numpy.linalg import matrix_rank

try:
    from scipy.linalg import eigvals
except Exception:  # pragma: no cover
    eigvals = None


def observability_matrix(A: np.ndarray, C: np.ndarray) -> np.ndarray:
    """Return O = [C; C A; ...; C A^(n-1)]."""
    A = np.asarray(A, dtype=float)
    C = np.atleast_2d(np.asarray(C, dtype=float))
    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 pbh_observability_test(A: np.ndarray, C: np.ndarray, tol: float = 1e-9) -> dict:
    """Run the PBH observability test and return a diagnostic dictionary."""
    A = np.asarray(A, dtype=float)
    C = np.atleast_2d(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 = A.shape[0]
    if eigvals is not None:
        lambdas = eigvals(A)
    else:
        lambdas = np.linalg.eigvals(A)

    rank_O = matrix_rank(observability_matrix(A, C), tol=tol)
    mode_reports = []
    observable = True

    for lam in lambdas:
        pbh_matrix = np.vstack((lam * np.eye(n, dtype=complex) - A, C.astype(complex)))
        r = matrix_rank(pbh_matrix, tol=tol)
        passed = (r == n)
        observable = observable and passed
        mode_reports.append(
            {
                "lambda": lam,
                "rank": int(r),
                "passed": bool(passed),
            }
        )

    return {
        "observable_by_pbh": bool(observable),
        "observability_matrix_rank": int(rank_O),
        "state_dimension": int(n),
        "mode_reports": mode_reports,
    }


def left_eigenvector_check(A: np.ndarray, C: np.ndarray, tol: float = 1e-8) -> list[dict]:
    """
    Check the equivalent condition: no nonzero left eigenvector q* A = lambda q*
    satisfies q* C? For observability, use A.T right eigenvectors v with C v != 0.
    """
    A = np.asarray(A, dtype=float)
    C = np.atleast_2d(np.asarray(C, dtype=float))
    vals, vecs = np.linalg.eig(A.T)
    reports = []
    for k, lam in enumerate(vals):
        v = vecs[:, k]
        cv = C.astype(complex) @ v.astype(complex)
        reports.append(
            {
                "lambda": lam,
                "norm_Cv": float(np.linalg.norm(cv)),
                "visible": bool(np.linalg.norm(cv) > tol),
            }
        )
    return reports


if __name__ == "__main__":
    A1 = np.array([[0.0, 1.0], [-2.0, -3.0]])
    C1 = np.array([[1.0, 0.0]])

    A2 = np.array([[-1.0, 0.0], [0.0, -2.0]])
    C2 = np.array([[1.0, 0.0]])

    for name, A, C in [("Example 1", A1, C1), ("Example 2", A2, C2)]:
        print("=" * 70)
        print(name)
        result = pbh_observability_test(A, C)
        print("PBH observable:", result["observable_by_pbh"])
        print("Rank of observability matrix:", result["observability_matrix_rank"])
        for report in result["mode_reports"]:
            print(
                f"lambda={report['lambda']}, "
                f"rank={report['rank']}, passed={report['passed']}"
            )
        print("Left-eigenvector visibility:")
        for report in left_eigenvector_check(A, C):
            print(
                f"lambda={report['lambda']}, "
                f"||C v||={report['norm_Cv']:.3e}, visible={report['visible']}"
            )

For C++, the most convenient modern numerical library is Eigen. It provides dense matrices, eigensolvers, and rank estimation through decompositions such as full-pivot LU.

File: Chapter14_Lesson2.cpp

/*
Chapter14_Lesson2.cpp
PBH Test for Observability using Eigen.

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

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

using Complex = std::complex<double>;
using MatrixXd = Eigen::MatrixXd;
using MatrixXcd = Eigen::MatrixXcd;

int complexRank(const MatrixXcd& M, double tol = 1e-9) {
    Eigen::FullPivLU<MatrixXcd> lu(M);
    lu.setThreshold(tol);
    return lu.rank();
}

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;
}

bool pbhObservable(const MatrixXd& A, const MatrixXd& C, double tol = 1e-9) {
    const int n = static_cast<int>(A.rows());
    const int p = static_cast<int>(C.rows());

    Eigen::EigenSolver<MatrixXd> solver(A);
    const auto eigenvalues = solver.eigenvalues();

    bool observable = true;
    for (int k = 0; k < eigenvalues.size(); ++k) {
        Complex lambda = eigenvalues(k);
        MatrixXcd M(n + p, n);
        M.topRows(n) = lambda * MatrixXcd::Identity(n, n) - A.cast<Complex>();
        M.bottomRows(p) = C.cast<Complex>();

        int r = complexRank(M, tol);
        bool passed = (r == n);
        observable = observable && passed;

        std::cout << "lambda = " << lambda
                  << ", rank = " << r
                  << ", PBH mode passed = " << std::boolalpha << passed
                  << "\n";
    }
    return observable;
}

int main() {
    MatrixXd A1(2, 2);
    A1 << 0.0, 1.0,
         -2.0, -3.0;
    MatrixXd C1(1, 2);
    C1 << 1.0, 0.0;

    MatrixXd A2(2, 2);
    A2 << -1.0, 0.0,
           0.0, -2.0;
    MatrixXd C2(1, 2);
    C2 << 1.0, 0.0;

    std::vector<std::pair<MatrixXd, MatrixXd>> examples = { {A1, C1}, {A2, C2} };

    for (std::size_t i = 0; i < examples.size(); ++i) {
        std::cout << "================ Example " << i + 1 << " ================\n";
        const MatrixXd& A = examples[i].first;
        const MatrixXd& C = examples[i].second;
        MatrixXd O = observabilityMatrix(A, C);
        Eigen::FullPivLU<MatrixXd> lu(O);
        std::cout << "Rank of observability matrix = " << lu.rank() << "\n";
        bool ok = pbhObservable(A, C);
        std::cout << "Observable by PBH = " << std::boolalpha << ok << "\n";
    }

    return 0;
}

Java does not include a built-in numerical linear algebra package. This example uses Apache Commons Math for eigenvalues and a small from-scratch Gaussian-elimination rank routine. Complex PBH matrices are handled through realification.

File: Chapter14_Lesson2.java

/*
Chapter14_Lesson2.java
PBH Test for Observability using Apache Commons Math for eigenvalues
and a from-scratch Gaussian-elimination rank routine.

Compile example:
    javac -cp commons-math3-3.6.1.jar Chapter14_Lesson2.java
Run example:
    java -cp .:commons-math3-3.6.1.jar Chapter14_Lesson2
*/

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.linear.RealMatrix;

public class Chapter14_Lesson2 {
    static double[][] vstack(double[][] top, double[][] bottom) {
        int rows = top.length + bottom.length;
        int cols = top[0].length;
        double[][] out = new double[rows][cols];
        for (int i = 0; i < top.length; i++) {
            System.arraycopy(top[i], 0, out[i], 0, cols);
        }
        for (int i = 0; i < bottom.length; i++) {
            System.arraycopy(bottom[i], 0, out[top.length + i], 0, cols);
        }
        return out;
    }

    static int rankReal(double[][] input, double tol) {
        int m = input.length;
        int n = input[0].length;
        double[][] a = new double[m][n];
        for (int i = 0; i < m; i++) {
            System.arraycopy(input[i], 0, a[i], 0, n);
        }

        int rank = 0;
        int row = 0;
        for (int col = 0; col < n && row < m; col++) {
            int pivot = row;
            for (int i = row + 1; i < m; i++) {
                if (Math.abs(a[i][col]) > Math.abs(a[pivot][col])) pivot = i;
            }
            if (Math.abs(a[pivot][col]) <= tol) continue;

            double[] temp = a[row];
            a[row] = a[pivot];
            a[pivot] = temp;

            double pivotValue = a[row][col];
            for (int j = col; j < n; j++) a[row][j] /= pivotValue;

            for (int i = 0; i < m; i++) {
                if (i == row) continue;
                double factor = a[i][col];
                for (int j = col; j < n; j++) a[i][j] -= factor * a[row][j];
            }
            row++;
            rank++;
        }
        return rank;
    }

    static int rankComplexByRealification(double[][] real, double[][] imag, double tol) {
        int m = real.length;
        int n = real[0].length;
        double[][] R = new double[2 * m][2 * n];
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                R[i][j] = real[i][j];
                R[i][j + n] = -imag[i][j];
                R[i + m][j] = imag[i][j];
                R[i + m][j + n] = real[i][j];
            }
        }
        return rankReal(R, tol);
    }

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

    static boolean pbhObservable(double[][] A, double[][] C, double tol) {
        int n = A.length;
        int p = C.length;
        RealMatrix Am = new Array2DRowRealMatrix(A);
        EigenDecomposition eig = new EigenDecomposition(Am);
        boolean observable = true;

        for (int k = 0; k < n; k++) {
            double lr = eig.getRealEigenvalue(k);
            double li = eig.getImagEigenvalue(k);

            double[][] real = new double[n + p][n];
            double[][] imag = new double[n + p][n];

            for (int i = 0; i < n; i++) {
                for (int j = 0; j < n; j++) {
                    real[i][j] = -A[i][j];
                    imag[i][j] = 0.0;
                }
                real[i][i] += lr;
                imag[i][i] += li;
            }
            for (int i = 0; i < p; i++) {
                for (int j = 0; j < n; j++) {
                    real[n + i][j] = C[i][j];
                    imag[n + i][j] = 0.0;
                }
            }

            int realifiedRank = rankComplexByRealification(real, imag, tol);
            boolean passed = (realifiedRank == 2 * n);
            observable = observable && passed;
            System.out.printf("lambda = %.6f%+.6fi, realified rank = %d, PBH mode passed = %b%n",
                    lr, li, realifiedRank, passed);
        }
        return observable;
    }

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

    static double[][] multiply(double[][] A, double[][] B) {
        int m = A.length;
        int n = B[0].length;
        int r = 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 < r; k++) sum += A[i][k] * B[k][j];
                C[i][j] = sum;
            }
        }
        return C;
    }

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

        double[][] A2 = { {-1.0, 0.0}, {0.0, -2.0} };
        double[][] C2 = { {1.0, 0.0} };

        double[][][] As = {A1, A2};
        double[][][] Cs = {C1, C2};

        for (int i = 0; i < As.length; i++) {
            System.out.println("================ Example " + (i + 1) + " ================");
            int rankO = rankReal(observabilityMatrix(As[i], Cs[i]), 1e-9);
            System.out.println("Rank of observability matrix = " + rankO);
            boolean ok = pbhObservable(As[i], Cs[i], 1e-9);
            System.out.println("Observable by PBH = " + ok);
        }
    }
}

MATLAB directly supports eig, rank, and, when the Control System Toolbox is available, obsv. The script below implements the PBH rank test explicitly and also creates a simple Simulink model containing a State-Space block when Simulink is available.

File: Chapter14_Lesson2.m

% Chapter14_Lesson2.m
% PBH Test for Observability and optional Simulink State-Space model creation.
%
% System:
%   x_dot = A x + B u
%   y     = C x + D u
%
% Observability depends only on (A, C).

clear; clc;

A1 = [0 1; -2 -3];
C1 = [1 0];

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

tol = 1e-9;

fprintf('================ Example 1 ================\n');
pbh_observability_report(A1, C1, tol);

fprintf('================ Example 2 ================\n');
pbh_observability_report(A2, C2, tol);

% Optional: create a minimal Simulink model with a State-Space block.
% This requires Simulink. The model is for simulation; the PBH test itself is
% algebraic and is implemented above.
if exist('simulink', 'file') == 4
    model = 'Chapter14_Lesson2_PBH_Observability_Model';
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    new_system(model);
    open_system(model);

    B = [0; 1];
    D = 0;
    assignin('base', 'A_sim', A1);
    assignin('base', 'B_sim', B);
    assignin('base', 'C_sim', C1);
    assignin('base', 'D_sim', D);

    add_block('simulink/Sources/Step', [model '/Step Input'], ...
        'Position', [80 90 120 120]);
    add_block('simulink/Continuous/State-Space', [model '/State-Space Plant'], ...
        'A', 'A_sim', 'B', 'B_sim', 'C', 'C_sim', 'D', 'D_sim', ...
        'Position', [180 75 330 135]);
    add_block('simulink/Sinks/Scope', [model '/Output Scope'], ...
        'Position', [400 85 440 125]);
    add_line(model, 'Step Input/1', 'State-Space Plant/1');
    add_line(model, 'State-Space Plant/1', 'Output Scope/1');
    save_system(model);
    fprintf('Created Simulink model: %s.slx\n', model);
else
    fprintf('Simulink not available; algebraic PBH test completed.\n');
end

function pbh_observability_report(A, C, tol)
    n = size(A, 1);
    O = observability_matrix(A, C);
    rankO = rank(O, tol);
    fprintf('Rank of observability matrix = %d of %d\n', rankO, n);

    lambdas = eig(A);
    observable = true;
    for k = 1:length(lambdas)
        lambda = lambdas(k);
        M = [lambda*eye(n) - A; C];
        r = rank(M, tol);
        passed = (r == n);
        observable = observable && passed;
        fprintf('lambda = %.6f%+.6fi, PBH rank = %d, passed = %d\n', ...
            real(lambda), imag(lambda), r, passed);
    end
    fprintf('Observable by PBH = %d\n', observable);
end

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

Wolfram Mathematica is suitable for symbolic PBH calculations because it can compute eigenvalues, ranks, and exact matrix expressions when the entries are rational or symbolic.

File: Chapter14_Lesson2.nb

Notebook[{
Cell["Chapter14_Lesson2.nb", "Title"],
Cell["PBH Test for Observability", "Section"],
Cell[BoxData[
 RowBox[{"ClearAll", "[", RowBox[{"pbhObservable", ",", "observabilityMatrix"}], "]"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"observabilityMatrix", "[", RowBox[{"A_", ",", "C_"}], "]"}], ":=",
  RowBox[{"Join", "@@", RowBox[{"Table", "[",
    RowBox[{RowBox[{"C", ".", RowBox[{"MatrixPower", "[", RowBox[{"A", ",", "k"}], "]"}]}], ",",
    RowBox[{"{", RowBox[{"k", ",", "0", ",", RowBox[{RowBox[{"Length", "[", "A", "]"}], "-", "1"}]}], "}"}]}], "]"}]}]}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"pbhObservable", "[", RowBox[{"A_", ",", "C_"}], "]"}], ":=", RowBox[{"Module", "[",
   RowBox[{RowBox[{"{", RowBox[{"n", ",", "lambdas", ",", "tests"}], "}"}], ",",
    RowBox[{
     RowBox[{"n", "=", RowBox[{"Length", "[", "A", "]"}]}], ";",
     RowBox[{"lambdas", "=", RowBox[{"Eigenvalues", "[", "A", "]"}]}], ";",
     RowBox[{"tests", "=", RowBox[{"Table", "[",
       RowBox[{RowBox[{"{", RowBox[{"lambda", ",",
          RowBox[{"MatrixRank", "[", RowBox[{"Join", "[",
            RowBox[{RowBox[{"lambda", " ", RowBox[{"IdentityMatrix", "[", "n", "]"}], "-", "A"}], ",", "C"}], "]"}], "]"}]}], "}"}], ",",
        RowBox[{"{", RowBox[{"lambda", ",", "lambdas"}], "}"}]}], "]"}]}], ";",
     RowBox[{"{", RowBox[{RowBox[{""ObservableByPBH"", "->", RowBox[{"AllTrue", "[", RowBox[{"tests", ",", RowBox[{RowBox[{RowBox[{"#[[2]]", "==", "n"}], "&"}]}], "]"}]}], ",",
       RowBox[{""ModeRanks"", "->", "tests"}], ",",
       RowBox[{""ObservabilityMatrixRank"", "->", RowBox[{"MatrixRank", "[", RowBox[{"observabilityMatrix", "[", RowBox[{"A", ",", "C"}], "]"}], "]"}]}]}], "}"}]}]}], "]"}]}]], "Input"],
Cell[BoxData[
 RowBox[{RowBox[{"A1", "=", RowBox[{"{ {", RowBox[{RowBox[{"0", ",", "1"}], "},", RowBox[{"{", RowBox[{RowBox[{"-", "2"}], ",", RowBox[{"-", "3"}]}], "}"}]}], "}"}]}], ";",
  RowBox[{"C1", "=", RowBox[{"{ {", RowBox[{"1", ",", "0"}], "} }"}]}], ";",
  RowBox[{"pbhObservable", "[", RowBox[{"A1", ",", "C1"}], "]"}]}]], "Input"],
Cell[BoxData[
 RowBox[{RowBox[{"A2", "=", RowBox[{"{ {", RowBox[{RowBox[{"-", "1"}], ",", "0"}], "},", RowBox[{"{", RowBox[{"0", ",", RowBox[{"-", "2"}]}], "}"}], "}"}]}], ";",
  RowBox[{"C2", "=", RowBox[{"{ {", RowBox[{"1", ",", "0"}], "} }"}]}], ";",
  RowBox[{"pbhObservable", "[", RowBox[{"A2", ",", "C2"}], "]"}]}]], "Input"]
}]

Problem 1 (PBH rank test): Consider \( \mathbf{A}=\begin{bmatrix}0&1\\-6&-5\end{bmatrix} \) and \( \mathbf{C}=\begin{bmatrix}1&0\end{bmatrix} \). Use the PBH test to determine whether the pair is observable.

Solution: The characteristic polynomial is

\[ \det(\lambda\mathbf{I}-\mathbf{A})=\lambda^2+5\lambda+6 = (\lambda+2)(\lambda+3). \]

For \( \lambda=-2 \),

\[ \begin{bmatrix}-2\mathbf{I}-\mathbf{A}\\\mathbf{C}\end{bmatrix} =\begin{bmatrix}-2&-1\\6&3\\1&0\end{bmatrix}, \qquad \operatorname{rank}=2. \]

For \( \lambda=-3 \),

\[ \begin{bmatrix}-3\mathbf{I}-\mathbf{A}\\\mathbf{C}\end{bmatrix} =\begin{bmatrix}-3&-1\\6&2\\1&0\end{bmatrix}, \qquad \operatorname{rank}=2. \]

Both ranks equal \( n=2 \), hence the pair is observable.

Problem 2 (unobservable eigenmode): Let \( \mathbf{A}=\operatorname{diag}(-1,-4,2) \) and \( \mathbf{C}=\begin{bmatrix}1&0&0\\0&1&0\end{bmatrix} \). Identify the unobservable mode and decide detectability.

Solution: The eigenvalues are \( -1,-4,2 \). The third eigenvector is \( \mathbf{e}_3 \), and

\[ \mathbf{C}\mathbf{e}_3=\begin{bmatrix}1&0&0\\0&1&0\end{bmatrix} \begin{bmatrix}0\\0\\1\end{bmatrix}=\begin{bmatrix}0\\0\end{bmatrix}. \]

Therefore the mode associated with \( \lambda=2 \) is unobservable. Since \( \operatorname{Re}(2)>0 \), the pair is not detectable.

Problem 3 (equivalence with Kalman matrix): Show that if \( \operatorname{rank}(\mathcal{O}_n)<n \), then there exists at least one eigenvalue \( \lambda \) for which the PBH matrix loses rank.

Solution: If the observability matrix loses rank, then \( \mathcal{N}_o=\ker(\mathcal{O}_n) \) contains a nonzero vector. The subspace \( \mathcal{N}_o \) is invariant under \( \mathbf{A} \) because \( \mathbf{C}\mathbf{A}^k\mathbf{x}=\mathbf{0} \) for the required powers implies the same condition for \( \mathbf{A}\mathbf{x} \), using Cayley-Hamilton for the last power. Any nonzero finite-dimensional invariant subspace contains an eigenvector of the restricted operator. Thus there is a nonzero \( \mathbf{z}\in\mathcal{N}_o \) and an eigenvalue \( \lambda \) such that \( \mathbf{A}\mathbf{z}=\lambda\mathbf{z} \). Since \( \mathbf{C}\mathbf{z}=\mathbf{0} \),

\[ \begin{bmatrix}\lambda\mathbf{I}-\mathbf{A}\\\mathbf{C}\end{bmatrix} \mathbf{z}=\mathbf{0},\qquad \mathbf{z}\ne\mathbf{0}. \]

Hence the stacked PBH matrix does not have full column rank.

Problem 4 (sensor selection): Let \( \mathbf{A}=\operatorname{diag}(-1,-2,-3) \). Compare \( \mathbf{C}_1=\begin{bmatrix}1&1&0\end{bmatrix} \) and \( \mathbf{C}_2=\begin{bmatrix}1&1&1\end{bmatrix} \). Which output matrix gives an observable pair?

Solution: For a diagonal matrix with distinct eigenvalues, the eigenvectors are the coordinate vectors \( \mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3 \). The first output matrix satisfies

\[ \mathbf{C}_1\mathbf{e}_3=0, \]

so the \( -3 \) mode is not visible and \( (\mathbf{A},\mathbf{C}_1) \) is not observable. For \( \mathbf{C}_2 \),

\[ \mathbf{C}_2\mathbf{e}_1=1, \qquad \mathbf{C}_2\mathbf{e}_2=1, \qquad \mathbf{C}_2\mathbf{e}_3=1. \]

Every eigenmode is visible, so \( (\mathbf{A},\mathbf{C}_2) \) is observable.

Problem 5 (repeated eigenvalue): Let \( \mathbf{A}=\begin{bmatrix}-1&0\\0&-1\end{bmatrix} \) and \( \mathbf{C}=\begin{bmatrix}1&0\end{bmatrix} \). Use PBH to test observability.

Solution: There is only one eigenvalue, \( \lambda=-1 \), but its eigenspace is two-dimensional. The PBH matrix is

\[ \begin{bmatrix}-1\mathbf{I}-\mathbf{A}\\\mathbf{C}\end{bmatrix} =\begin{bmatrix}0&0\\0&0\\1&0\end{bmatrix}, \qquad \operatorname{rank}=1<2. \]

Therefore the pair is not observable. This example shows why repeated eigenvalues require a full rank test rather than a single scalar check.

The PBH test states that \( (\mathbf{A},\mathbf{C}) \) is observable exactly when \( \operatorname{rank}\begin{bmatrix}\lambda\mathbf{I}-\mathbf{A}\\\mathbf{C}\end{bmatrix}=n \) for every eigenvalue of \( \mathbf{A} \). The test is equivalent to the Kalman observability rank condition but gives a more modal interpretation: no eigenmode of the state matrix may lie entirely in the null space of the output map. Numerically, the PBH test is best implemented through singular-value-based rank decisions with explicit tolerances.