In earlier chapters, the internal state vector \( \mathbf{x}(t) \in \mathbb{R}^n \) was used to describe the system's internal condition. In practice, however, sensors usually measure only a few combinations of the state. For a continuous-time LTI model,

\[ \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) \]

the sensor matrix \( \mathbf{C} \) decides which state combinations are directly visible. Observability asks a precise inverse question: from the observed signal \( \mathbf{y}(t) \), and assuming the input \( \mathbf{u}(t) \) and matrices \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) are known, can we recover the unknown initial state \( \mathbf{x}(t_0) \)?

flowchart TD
  X0["Unknown initial state x(t0)"] --> DYN["Internal dynamics: xdot = A x + B u"]
  U["Known input u(t)"] --> DYN
  DYN --> SENS["Sensors: y = C x + D u"]
  SENS --> Y["Measured output history y(t)"]
  Y --> Q["Can x(t0) be uniquely inferred?"]
  Q --> YES["Observable: no hidden ambiguity"]
  Q --> NO["Unobservable: at least one hidden direction"]
        

This is the measurement counterpart of state steering: controllability studies actuator authority, while observability studies sensor information. This lesson avoids formal duality until Chapter 14, but the intuitive symmetry is useful.

By the solution formula from previous chapters, for \( t \in [t_0,t_f] \),

\[ \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}(t_0) +\int_{t_0}^t e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds . \]

Substitution into the output equation gives

\[ \mathbf{y}(t)= \mathbf{C}e^{\mathbf{A}(t-t_0)}\mathbf{x}(t_0) +\mathbf{C}\int_{t_0}^t e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds +\mathbf{D}\mathbf{u}(t). \]

The second and third terms are known when the input is known. Therefore the part of the output that contains information about the unknown initial state is

\[ \mathbf{z}(t)=\mathbf{y}(t) -\mathbf{C}\int_{t_0}^t e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds -\mathbf{D}\mathbf{u}(t) = \mathbf{C}e^{\mathbf{A}(t-t_0)}\mathbf{x}(t_0). \]

Thus observability is fundamentally a question about whether the map \( \mathbf{x}(t_0) \mapsto \mathbf{z}(t) \) is one-to-one over the observation interval.

Take two possible initial states \( \mathbf{x}_a(t_0) \) and \( \mathbf{x}_b(t_0) \). With the same known input, the difference between their input-corrected outputs is

\[ \Delta \mathbf{z}(t) = \mathbf{C}e^{\mathbf{A}(t-t_0)} \left(\mathbf{x}_a(t_0)-\mathbf{x}_b(t_0)\right). \]

Let \( \boldsymbol{\eta}=\mathbf{x}_a(t_0)-\mathbf{x}_b(t_0) \). If

\[ \mathbf{C}e^{\mathbf{A}(t-t_0)}\boldsymbol{\eta}=\mathbf{0} \quad \text{for every } t\in [t_0,t_f], \]

then the two initial states generate exactly the same measured output history. The vector \( \boldsymbol{\eta} \neq \mathbf{0} \) is an unobservable direction: movement along that direction is hidden from the sensors for all measured times.

\[ \boxed{\text{Observable means: the only output-invisible initial-state difference is } \boldsymbol{\eta}=\mathbf{0}.} \]

To see why powers of \( \mathbf{A} \) appear in observability, consider the zero-input case with \( \mathbf{D}=\mathbf{0} \):

\[ \mathbf{y}(t)=\mathbf{C}e^{\mathbf{A}(t-t_0)}\mathbf{x}(t_0). \]

Differentiating at \( t=t_0 \) gives

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

Stacking the first \( n \) derivative equations gives the measurement algebraic system

\[ \underbrace{\begin{bmatrix} \mathbf{y}(t_0)\\ \dot{\mathbf{y} }(t_0)\\ \vdots\\ \mathbf{y}^{(n-1)}(t_0) \end{bmatrix} }_{\mathbf{s} } = \underbrace{\begin{bmatrix} \mathbf{C}\\ \mathbf{C}\mathbf{A}\\ \vdots\\ \mathbf{C}\mathbf{A}^{n-1} \end{bmatrix} }_{\mathbf{O}_n} \mathbf{x}(t_0). \]

The matrix \( \mathbf{O}_n \) is the finite derivative version of the output-history map. If its columns are linearly independent, the stacked measurements contain enough independent equations to determine the initial state.

\[ \boxed{\operatorname{rank}\mathbf{O}_n=n \quad \Longrightarrow \quad \mathbf{x}(t_0) \text{ is uniquely determined by these derivatives.} } \]

The formal rank criterion is developed in Chapter 14. Here, the key intuition is that the sensor first sees \( \mathbf{C}\mathbf{x} \), then sees how dynamics propagate hidden components into measured coordinates through \( \mathbf{C}\mathbf{A}\mathbf{x} \), \( \mathbf{C}\mathbf{A}^2\mathbf{x} \), and so on.

Each measurement derivative removes possible ambiguity. The first measurement requires an invisible difference \( \boldsymbol{\eta} \) to satisfy \( \mathbf{C}\boldsymbol{\eta}=\mathbf{0} \). The first derivative additionally requires \( \mathbf{C}\mathbf{A}\boldsymbol{\eta}=\mathbf{0} \). Continuing gives the common hidden subspace

\[ \mathcal{N}_o = \ker(\mathbf{C}) \cap \ker(\mathbf{C}\mathbf{A}) \cap \cdots \cap \ker(\mathbf{C}\mathbf{A}^{n-1}). \]

If \( \mathcal{N}_o=\{\mathbf{0}\} \), there is no nonzero hidden initial-state direction. If \( \mathcal{N}_o \) contains a nonzero vector, then initial states separated by that vector cannot be distinguished from the measured output history.

flowchart TD
  ETA["Candidate ambiguity eta"] --> C0["Must satisfy C eta = 0"]
  C0 --> C1["Must satisfy C A eta = 0"]
  C1 --> C2["Must satisfy C A^2 eta = 0"]
  C2 --> CN["Continue up to C A^(n-1) eta = 0"]
  CN --> Z["Only eta = 0?"]
  Z -->|yes| OBS["No hidden direction: observable intuition"]
  Z -->|no| UNOBS["Nonzero hidden direction: unobservable intuition"]
        

Consider a second-order system with state \( \mathbf{x}=\begin{bmatrix}x_1 & x_2\end{bmatrix}^T \), where \( x_1 \) is position and \( x_2 \) is velocity:

\[ \dot{\mathbf{x} } = \begin{bmatrix}0 & 1\\ -2 & -3\end{bmatrix}\mathbf{x}, \qquad y= \begin{bmatrix}1 & 0\end{bmatrix}\mathbf{x}=x_1. \]

The sensor directly measures only position. But since \( \dot{y}=\dot{x}_1=x_2 \), the velocity is visible through the time variation of the measured position. The first two derivative equations are

\[ \begin{bmatrix}y(0)\\ \dot{y}(0)\end{bmatrix} = \begin{bmatrix}1 & 0\\ 0 & 1\end{bmatrix} \begin{bmatrix}x_1(0)\\x_2(0)\end{bmatrix}. \]

Hence \( x_1(0)=y(0) \) and \( x_2(0)=\dot{y}(0) \). Even one physical sensor can reveal multiple states if the dynamics couple those states into the measured coordinate.

Now consider

\[ \dot{\mathbf{x} } = \begin{bmatrix}-1 & 0\\0 & -2\end{bmatrix}\mathbf{x}, \qquad y= \begin{bmatrix}1 & 0\end{bmatrix}\mathbf{x}=x_1. \]

The output is

\[ y(t)=e^{-t}x_1(0). \]

The second initial coordinate \( x_2(0) \) never appears in the output. The derivative stack is

\[ \begin{bmatrix}y(0)\\\dot{y}(0)\end{bmatrix} = \begin{bmatrix}1 & 0\\-1 & 0\end{bmatrix} \begin{bmatrix}x_1(0)\\x_2(0)\end{bmatrix}. \]

The columns are dependent, and all states of the form \( \begin{bmatrix}x_1(0) & \alpha\end{bmatrix}^T \) produce the same measured signal for any value of \( \alpha \). The second mode is internally present but externally invisible to this sensor.

Notice the important separation: this hidden mode is stable because it decays as \( e^{-2t} \), but it is still not observable. Stability and observability are different properties.

Observability is normally defined with the same known input applied to all candidate initial states. The input does not hide the initial state because its contribution to the output is computable from the model:

\[ \mathbf{y}(t) = \underbrace{\mathbf{C}e^{\mathbf{A}(t-t_0)}\mathbf{x}(t_0)}_{\text{depends on initial state} } + \underbrace{\mathbf{C}\int_{t_0}^t e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds +\mathbf{D}\mathbf{u}(t)}_{\text{known input contribution} }. \]

In this lesson, we focus on reconstructing the initial state from output measurements under a known input. Problems involving unknown disturbances or unknown inputs require additional ideas and are outside the scope of this first observability lesson.

Computationally, the simplest observability check builds the stacked matrix \( \mathbf{O}_n \) and computes its rank. In floating-point arithmetic, the practical question is not only whether the rank is full, but also whether the matrix is well-conditioned. Nearly dependent columns mean that state reconstruction is sensitive to measurement noise.

\[ \widehat{\mathbf{x} }(t_0) = \mathbf{O}_n^\dagger \mathbf{s}, \qquad \mathbf{s}= \begin{bmatrix} \mathbf{y}(t_0) & \dot{\mathbf{y} }(t_0) & \cdots & \mathbf{y}^{(n-1)}(t_0) \end{bmatrix}^T. \]

Here \( \mathbf{O}_n^\dagger \) denotes the Moore-Penrose pseudoinverse. If the system is observable and the derivative data are exact, this reconstruction is exact. With noisy data, later lessons will motivate observer design rather than direct numerical differentiation.

Chapter13_Lesson1.py

# Chapter13_Lesson1.py
# Intuitive observability from output measurements for continuous-time LTI systems.
# Requirements: numpy, scipy. Optional: python-control for comparison.

import numpy as np
from numpy.linalg import matrix_rank, pinv, norm
from scipy.linalg import expm


def observability_matrix(A: np.ndarray, C: np.ndarray) -> np.ndarray:
    """Build O_n = [C; C A; ...; C A^(n-1)]."""
    A = np.asarray(A, dtype=float)
    C = 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 output_from_initial_state(A, C, x0, times):
    """Zero-input output y(t)=C exp(A t) x0 for t0=0."""
    return np.array([(C @ expm(A * t) @ x0).ravel() for t in times])


def derivative_stack_measurement(A, C, x0):
    """Stack y(0), y_dot(0), ..., y^(n-1)(0) = O_n x0."""
    O = observability_matrix(A, C)
    return O @ x0


def reconstruct_x0_from_derivatives(A, C, derivative_vector):
    """Least-squares reconstruction of x0 from stacked derivatives."""
    O = observability_matrix(A, C)
    return pinv(O) @ derivative_vector


def report_system(name, A, C, x_true):
    O = observability_matrix(A, C)
    print(f"\n{name}")
    print("A =\n", A)
    print("C =\n", C)
    print("O_n =\n", O)
    print("rank(O_n) =", matrix_rank(O), "out of n =", A.shape[0])

    z = derivative_stack_measurement(A, C, x_true)
    x_hat = reconstruct_x0_from_derivatives(A, C, z)
    print("true x0 =", x_true.ravel())
    print("stacked derivatives =", z.ravel())
    print("least-squares reconstructed x0 =", x_hat.ravel())
    print("reconstruction error norm =", norm(x_hat - x_true))


if __name__ == "__main__":
    # Example 1: measuring position reveals velocity through the derivative of position.
    A1 = np.array([[0.0, 1.0],
                   [-2.0, -3.0]])
    C1 = np.array([[1.0, 0.0]])
    x01 = np.array([[1.0], [2.0]])
    report_system("Example 1: position sensor, observable second-order system", A1, C1, x01)

    # Show output traces for two different initial states.
    times = np.linspace(0.0, 5.0, 101)
    y1 = output_from_initial_state(A1, C1, np.array([[1.0], [2.0]]), times)
    y2 = output_from_initial_state(A1, C1, np.array([[1.0], [-1.0]]), times)
    print("\nFirst five zero-input outputs for x0=[1,2]^T:", y1[:5, 0])
    print("First five zero-input outputs for x0=[1,-1]^T:", y2[:5, 0])

    # Example 2: diagonal dynamics with a sensor that sees only the first coordinate.
    # The second state never affects the output, so it is an unobservable direction.
    A2 = np.array([[-1.0, 0.0],
                   [0.0, -2.0]])
    C2 = np.array([[1.0, 0.0]])
    x02 = np.array([[3.0], [4.0]])
    report_system("Example 2: second mode invisible, unobservable", A2, C2, x02)

    # If desired and installed:
    # import control
    # print(control.obsv(A1, C1))

      

Python libraries commonly used for this topic include numpy for matrices, scipy.linalg.expm for matrix exponentials, scipy.signal.StateSpace for system objects, and the optional python-control package for functions such as obsv.

Chapter13_Lesson1.cpp

// Chapter13_Lesson1.cpp
// From-scratch observability matrix and rank test for small LTI systems.
// Compile: g++ -std=c++17 Chapter13_Lesson1.cpp -o Chapter13_Lesson1

#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>

using Matrix = std::vector<std::vector<double>>;

Matrix identity(int n) {
    Matrix I(n, std::vector<double>(n, 0.0));
    for (int i = 0; i < n; ++i) I[i][i] = 1.0;
    return I;
}

Matrix multiply(const Matrix& A, const Matrix& B) {
    int m = static_cast<int>(A.size());
    int p = static_cast<int>(B.size());
    int n = static_cast<int>(B[0].size());
    Matrix C(m, std::vector<double>(n, 0.0));
    for (int i = 0; i < m; ++i)
        for (int k = 0; k < p; ++k)
            for (int j = 0; j < n; ++j)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}

Matrix vstack(const std::vector<Matrix>& blocks) {
    int cols = static_cast<int>(blocks[0][0].size());
    Matrix out;
    for (const auto& B : blocks) {
        for (const auto& row : B) {
            if (static_cast<int>(row.size()) != cols) {
                throw std::runtime_error("Column mismatch in vstack.");
            }
            out.push_back(row);
        }
    }
    return out;
}

Matrix observabilityMatrix(const Matrix& A, const Matrix& C) {
    int n = static_cast<int>(A.size());
    Matrix Ak = identity(n);
    std::vector<Matrix> blocks;
    for (int k = 0; k < n; ++k) {
        blocks.push_back(multiply(C, Ak));
        Ak = multiply(Ak, A);
    }
    return vstack(blocks);
}

int rankGaussian(Matrix M, double tol = 1e-10) {
    int rows = static_cast<int>(M.size());
    int cols = static_cast<int>(M[0].size());
    int r = 0;
    for (int c = 0; c < cols && r < rows; ++c) {
        int pivot = r;
        for (int i = r + 1; i < rows; ++i)
            if (std::fabs(M[i][c]) > std::fabs(M[pivot][c])) pivot = i;

        if (std::fabs(M[pivot][c]) < tol) continue;

        std::swap(M[r], M[pivot]);
        double div = M[r][c];
        for (int j = c; j < cols; ++j) M[r][j] /= div;

        for (int i = 0; i < rows; ++i) {
            if (i == r) continue;
            double factor = M[i][c];
            for (int j = c; j < cols; ++j) M[i][j] -= factor * M[r][j];
        }
        ++r;
    }
    return r;
}

void printMatrix(const Matrix& M) {
    for (const auto& row : M) {
        for (double x : row) std::cout << std::setw(10) << std::setprecision(5) << x << " ";
        std::cout << "\n";
    }
}

void report(const std::string& name, const Matrix& A, const Matrix& C) {
    Matrix O = observabilityMatrix(A, C);
    std::cout << "\n" << name << "\nO_n =\n";
    printMatrix(O);
    std::cout << "rank(O_n) = " << rankGaussian(O)
              << " out of n = " << A.size() << "\n";
}

int main() {
    Matrix A1 = { {0.0, 1.0},
                 {-2.0, -3.0} };
    Matrix C1 = { {1.0, 0.0} };
    report("Example 1: position sensor, observable", A1, C1);

    Matrix A2 = { {-1.0, 0.0},
                 {0.0, -2.0} };
    Matrix C2 = { {1.0, 0.0} };
    report("Example 2: second state invisible, unobservable", A2, C2);

    return 0;
}

      

For production numerical control software, C++ implementations usually rely on linear algebra libraries such as Eigen, Armadillo, or BLAS/LAPACK wrappers. The example above is intentionally from scratch so that the rank calculation and matrix stacking are transparent.

Chapter13_Lesson1.java

// Chapter13_Lesson1.java
// From-scratch observability matrix and rank test for small LTI systems.
// Compile: javac Chapter13_Lesson1.java
// Run:     java Chapter13_Lesson1

public class Chapter13_Lesson1 {
    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 p = B.length;
        int n = B[0].length;
        double[][] C = new double[m][n];
        for (int i = 0; i < m; i++)
            for (int k = 0; k < p; k++)
                for (int j = 0; j < n; j++)
                    C[i][j] += A[i][k] * B[k][j];
        return C;
    }

    static double[][] observabilityMatrix(double[][] A, double[][] C) {
        int n = A.length;
        int q = C.length;
        double[][] O = new double[q * n][n];
        double[][] Ak = identity(n);

        for (int block = 0; block < n; block++) {
            double[][] CAk = multiply(C, Ak);
            for (int i = 0; i < q; i++)
                for (int j = 0; j < n; j++)
                    O[block * q + i][j] = CAk[i][j];
            Ak = multiply(Ak, A);
        }
        return O;
    }

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

        int r = 0;
        for (int c = 0; c < cols && r < rows; c++) {
            int pivot = r;
            for (int i = r + 1; i < rows; i++)
                if (Math.abs(M[i][c]) > Math.abs(M[pivot][c])) pivot = i;

            if (Math.abs(M[pivot][c]) < tol) continue;

            double[] tmp = M[r];
            M[r] = M[pivot];
            M[pivot] = tmp;

            double div = M[r][c];
            for (int j = c; j < cols; j++) M[r][j] /= div;

            for (int i = 0; i < rows; i++) {
                if (i == r) continue;
                double factor = M[i][c];
                for (int j = c; j < cols; j++) M[i][j] -= factor * M[r][j];
            }
            r++;
        }
        return r;
    }

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

    static void report(String name, double[][] A, double[][] C) {
        double[][] O = observabilityMatrix(A, C);
        System.out.println("\n" + name);
        System.out.println("O_n =");
        printMatrix(O);
        System.out.println("rank(O_n) = " + rank(O, 1e-10) + " out of n = " + A.length);
    }

    public static void main(String[] args) {
        double[][] A1 = { {0.0, 1.0}, {-2.0, -3.0} };
        double[][] C1 = { {1.0, 0.0} };
        report("Example 1: position sensor, observable", A1, C1);

        double[][] A2 = { {-1.0, 0.0}, {0.0, -2.0} };
        double[][] C2 = { {1.0, 0.0} };
        report("Example 2: second state invisible, unobservable", A2, C2);
    }
}

      

Java projects often use libraries such as Apache Commons Math, EJML, or ojAlgo for linear algebra. The code above avoids external dependencies and implements only the operations needed for this lesson.

Chapter13_Lesson1.m

% Chapter13_Lesson1.m
% Intuitive observability from output measurements for continuous-time LTI systems.
% MATLAB Control System Toolbox functions: ss, initial, obsv.
% The first part also works with basic MATLAB except initial/ss/obsv.

clear; clc; close all;

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

A2 = [-1 0; 0 -2];
C2 = [1 0];
x02 = [3; 4];

fprintf('\nExample 1: position sensor, observable\n');
O1 = local_observability_matrix(A1, C1);
disp(O1);
fprintf('rank(O1) = %d out of n = %d\n', rank(O1), size(A1,1));

fprintf('\nExample 2: second state invisible, unobservable\n');
O2 = local_observability_matrix(A2, C2);
disp(O2);
fprintf('rank(O2) = %d out of n = %d\n', rank(O2), size(A2,1));

% Derivative-stack reconstruction: z = O_n x0.
z1 = O1*x01;
x01_hat = pinv(O1)*z1;
fprintf('\nTrue x01 = [%g %g]^T\n', x01(1), x01(2));
fprintf('Recovered x01 from y(0), y_dot(0) = [%g %g]^T\n', x01_hat(1), x01_hat(2));

% Compare output traces for two different initial conditions.
t = linspace(0, 5, 300);
Y1 = zeros(size(t));
Y2 = zeros(size(t));
x0a = [1; 2];
x0b = [1; -1];
for k = 1:numel(t)
    Y1(k) = C1*expm(A1*t(k))*x0a;
    Y2(k) = C1*expm(A1*t(k))*x0b;
end

figure;
plot(t, Y1, 'LineWidth', 1.5); hold on;
plot(t, Y2, '--', 'LineWidth', 1.5);
grid on;
xlabel('time');
ylabel('zero-input output y(t)');
legend('x0 = [1, 2]^T', 'x0 = [1, -1]^T');
title('Different initial states produce distinguishable outputs');

% Optional Control System Toolbox check.
if exist('obsv', 'file') == 2
    fprintf('\nControl System Toolbox obsv(A1,C1):\n');
    disp(obsv(A1,C1));
end

% Optional Simulink model creation if Simulink is installed.
if exist('simulink', 'file') == 4
    modelName = 'Chapter13_Lesson1_Simulink';
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    new_system(modelName);
    open_system(modelName);

    add_block('simulink/Continuous/State-Space', [modelName '/State-Space']);
    set_param([modelName '/State-Space'], 'A', 'A1', 'B', '[0;0]', 'C', 'C1', 'D', '0', 'X0', 'x01');

    add_block('simulink/Sinks/Scope', [modelName '/Scope']);
    add_line(modelName, 'State-Space/1', 'Scope/1');

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

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

      

MATLAB's Control System Toolbox provides ss, initial, and obsv. The script also includes an optional programmatic Simulink model using a State-Space block when Simulink is installed.

Chapter13_Lesson1.nb

(* Chapter13_Lesson1.nb *)
(* Wolfram Mathematica / Wolfram Language code for intuitive observability. *)

ClearAll[observabilityMatrix, reportSystem];

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

reportSystem[name_, A_, C_, x0_] := Module[
  {O, z, xhat},
  O = observabilityMatrix[A, C];
  z = O . x0;
  xhat = PseudoInverse[O] . z;
  Print["\n", name];
  Print["A = ", MatrixForm[A]];
  Print["C = ", MatrixForm[C]];
  Print["O_n = ", MatrixForm[O]];
  Print["Rank[O_n] = ", MatrixRank[O], " out of n = ", Length[A]];
  Print["stacked derivatives = ", z];
  Print["least-squares reconstructed x0 = ", xhat];
];

A1 = { {0, 1}, {-2, -3} };
C1 = { {1, 0} };
x01 = {1, 2};
reportSystem["Example 1: position sensor, observable", A1, C1, x01];

A2 = { {-1, 0}, {0, -2} };
C2 = { {1, 0} };
x02 = {3, 4};
reportSystem["Example 2: second state invisible, unobservable", A2, C2, x02];

(* Zero-input output y(t)=C Exp[A t] x0. *)
y1[t_] := C1 . MatrixExp[A1 t] . {1, 2};
y2[t_] := C1 . MatrixExp[A1 t] . {1, -1};

Plot[
  {First[y1[t]], First[y2[t]]},
  {t, 0, 5},
  PlotLegends -> {"x0 = {1,2}", "x0 = {1,-1}"},
  AxesLabel -> {"time", "y(t)"},
  PlotLabel -> "Distinguishable zero-input output traces"
]

      

Mathematica is especially convenient for symbolic matrix powers, matrix exponentials, exact rank calculations, and symbolic output formulas.

Problem 1 (Distinguishable Initial States): For \( \dot{\mathbf{x} }=\begin{bmatrix}0 & 1\\-2 & -3\end{bmatrix}\mathbf{x} \) and \( y=\begin{bmatrix}1 & 0\end{bmatrix}\mathbf{x} \), determine whether the two initial states \( \mathbf{x}_a(0)=\begin{bmatrix}1\\2\end{bmatrix} \) and \( \mathbf{x}_b(0)=\begin{bmatrix}1\\-1\end{bmatrix} \) are distinguishable from output measurements.

Solution: Their difference is \( \boldsymbol{\eta}=\begin{bmatrix}0\\3\end{bmatrix} \). The first output difference is \( \mathbf{C}\boldsymbol{\eta}=0 \), so the outputs agree at \( t=0 \). However,

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

Therefore \( \dot{y}_a(0)\neq \dot{y}_b(0) \), so the two states are distinguishable.

Problem 2 (Invisible Mode): For \( \mathbf{A}=\operatorname{diag}(-1,-2) \) and \( \mathbf{C}=\begin{bmatrix}1 & 0\end{bmatrix} \), find all initial states that produce the same output as \( \mathbf{x}_0=\begin{bmatrix}5\\0\end{bmatrix} \).

Solution: The output is \( y(t)=5e^{-t} \). Since the second coordinate never appears in the output, every initial state \( \begin{bmatrix}5\\\alpha\end{bmatrix} \), with arbitrary scalar \( \alpha \), produces the same output. The hidden direction is \( \operatorname{span}\{\begin{bmatrix}0\\1\end{bmatrix}\} \).

Problem 3 (Derivative Stack Reconstruction): Suppose \( \mathbf{O}_2=\begin{bmatrix}1 & 0\\0 & 1\end{bmatrix} \) and the measured derivative stack is \( \mathbf{s}=\begin{bmatrix}4\\-2\end{bmatrix} \). Recover \( \mathbf{x}(0) \).

Solution: Since \( \mathbf{s}=\mathbf{O}_2\mathbf{x}(0) \) and \( \mathbf{O}_2=\mathbf{I} \),

\[ \mathbf{x}(0)=\mathbf{O}_2^{-1}\mathbf{s} = \begin{bmatrix}4\\-2\end{bmatrix}. \]

Problem 4 (Known Input Correction): Explain why a known input does not change the intuitive observability question for the initial state.

Solution: The measured output can be decomposed into an initial-state term plus a known input-generated term:

\[ \mathbf{y}(t)= \mathbf{C}e^{\mathbf{A}(t-t_0)}\mathbf{x}(t_0) + \mathbf{g}_u(t), \]

where \( \mathbf{g}_u(t)= \mathbf{C}\int_{t_0}^t e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds +\mathbf{D}\mathbf{u}(t) \) is computable from the known input and model. Subtracting \( \mathbf{g}_u(t) \) leaves exactly the same initial-state information map \( \mathbf{C}e^{\mathbf{A}(t-t_0)}\mathbf{x}(t_0) \).

Problem 5 (Why Only the First n Derivatives Matter): Give an intuitive argument for why the first \( n \) derivative equations are enough for an \( n \)-state LTI system.

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

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

Multiplying by \( \mathbf{C} \) on the left shows that \( \mathbf{C}\mathbf{A}^n \) is a linear combination of \( \mathbf{C},\mathbf{C}\mathbf{A},\dots, \mathbf{C}\mathbf{A}^{n-1} \). The same logic applies to all higher powers. Therefore higher derivatives cannot introduce new independent row information beyond the first \( n \) derivative blocks.

Observability means that output measurements contain enough information to distinguish the initial state. For LTI systems, the output history depends on the initial state through \( \mathbf{C}e^{\mathbf{A}(t-t_0)}\mathbf{x}(t_0) \). Hidden state directions are vectors that lie simultaneously in the null spaces of \( \mathbf{C},\mathbf{C}\mathbf{A}, \dots,\mathbf{C}\mathbf{A}^{n-1} \). If no nonzero hidden direction exists, the state is observable from the measurements.