Observability asks whether the internal state can be inferred from output measurements. In this chapter we focus on the 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). \]

For observability of the initial state, the input is set to zero because the purpose is to isolate the output caused by \( \mathbf{x}(t_0) \). Thus

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

The observability Gramian is the matrix representation of the mapping from the initial state to the entire output record over an interval \( [t_0,t_f] \), where \( t_f > t_0 \). It turns a signal question into a quadratic-form question.

flowchart TD
  A["Choose observation interval [t0, tf]"] --> B["Set input equal to zero"]
  B --> C["Propagate state with exp(A time)"]
  C --> D["Map state to output with C"]
  D --> E["Square and integrate output magnitude"]
  E --> F["Collect result as the matrix Wo"]
  F --> G["Analyze rank, eigenvalues, and null space"]
        

Let the state dimension be \( n \) and the output dimension be \( p \), with \( \mathbf{A}\in\mathbb{R}^{n\times n} \) and \( \mathbf{C}\in\mathbb{R}^{p\times n} \). For an LTI system observed over an interval of length \( T=t_f-t_0 \), the finite-horizon observability Gramian is

\[ \boxed{\mathbf{W}_o(T)=\int_{0}^{T} e^{\mathbf{A}^{T}\tau}\mathbf{C}^{T}\mathbf{C} e^{\mathbf{A}\tau}\,d\tau}, \qquad T > 0. \]

Equivalently, written directly on \( [t_0,t_f] \),

\[ \boxed{\mathbf{W}_o(t_0,t_f)=\int_{t_0}^{t_f} \boldsymbol{\Phi}^{T}(t,t_0)\mathbf{C}^{T}\mathbf{C} \boldsymbol{\Phi}(t,t_0)\,dt}, \]

where \( \boldsymbol{\Phi}(t,t_0)=e^{\mathbf{A}(t-t_0)} \) for an LTI system. The matrix \( \mathbf{C}^{T}\mathbf{C} \) weights state directions according to how they are seen by the output matrix \( \mathbf{C} \).

Entrywise, if \( \mathbf{c}_i(t) \) denotes the \( i \)-th column of \( \mathbf{C}e^{\mathbf{A}t} \), then

\[ [\mathbf{W}_o(T)]_{ij}=\int_0^T \mathbf{c}_i^{T}(t)\mathbf{c}_j(t)\,dt. \]

Hence the Gramian contains inner products between output signatures generated by basis-state perturbations.

Although Lesson 2 develops output energy more deeply, the defining identity should be established here. For zero input and initial state \( \mathbf{x}_0=\mathbf{x}(0) \),

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

The squared output signal over \( [0,T] \) is

\[ \int_0^T \|\mathbf{y}(t)\|_2^2\,dt =\int_0^T \mathbf{x}_0^{T}e^{\mathbf{A}^{T}t} \mathbf{C}^{T}\mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0\,dt =\mathbf{x}_0^{T}\mathbf{W}_o(T)\mathbf{x}_0. \]

Therefore, \( \mathbf{W}_o(T) \) is not only a rank test object; it also quantifies how visible each initial-state direction is through the measured output.

flowchart TD
  X["Initial state x0"] --> D["Zero-input dynamics"]
  D --> S["State trajectory x(t)"]
  S --> Y["Output trajectory y(t)"]
  Y --> E["Integrate squared output"]
  E --> Q["Quadratic form x0^T Wo x0"]
        

The integrand in the Gramian definition has the form \( \mathbf{M}^{T}(t)\mathbf{M}(t) \), where \( \mathbf{M}(t)=\mathbf{C}e^{\mathbf{A}t} \). Hence the observability Gramian is symmetric:

\[ \mathbf{W}_o^{T}(T) =\left(\int_0^T e^{\mathbf{A}^{T}t}\mathbf{C}^{T}\mathbf{C} e^{\mathbf{A}t}\,dt\right)^{T} =\mathbf{W}_o(T). \]

For any nonzero or zero vector \( \boldsymbol{\xi}\in\mathbb{R}^{n} \),

\[ \boldsymbol{\xi}^{T}\mathbf{W}_o(T)\boldsymbol{\xi} =\int_0^T \|\mathbf{C}e^{\mathbf{A}t}\boldsymbol{\xi}\|_2^2\,dt \ge 0. \]

Thus \( \mathbf{W}_o(T) \) is always positive semidefinite:

\[ \boxed{\mathbf{W}_o(T)=\mathbf{W}_o^{T}(T),\qquad \mathbf{W}_o(T)\succeq \mathbf{0}.} \]

Positive semidefiniteness is guaranteed even if the system is not observable. Strict positive definiteness is stronger and is equivalent to complete observability for the LTI pair \( (\mathbf{A},\mathbf{C}) \).

A direction \( \boldsymbol{\xi} \) lies in the null space of the finite-horizon Gramian precisely when it produces zero output on the whole observation interval:

\[ \boldsymbol{\xi}\in\mathcal{N}(\mathbf{W}_o(T)) \quad\Longleftrightarrow\quad \mathbf{C}e^{\mathbf{A}t}\boldsymbol{\xi}=\mathbf{0} \;\; \text{for all } t\in[0,T]. \]

The implication follows from the nonnegative integral \( \int_0^T \|\mathbf{C}e^{\mathbf{A}t}\boldsymbol{\xi}\|_2^2\,dt \). If the integral is zero, the continuous nonnegative integrand must be zero on the interval.

The Kalman observability matrix from Chapter 14 is

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

For an LTI system, the Gramian test and the Kalman rank test are equivalent:

\[ \boxed{(\mathbf{A},\mathbf{C})\text{ observable} \quad\Longleftrightarrow\quad \operatorname{rank}\mathcal{O}=n \quad\Longleftrightarrow\quad \mathbf{W}_o(T)\succ \mathbf{0}\;\text{ for every }T > 0.} \]

The proof uses analyticity of \( \mathbf{C}e^{\mathbf{A}t}\boldsymbol{\xi} \). If it is zero on a nonzero interval, then all derivatives at \( t=0 \) are zero:

\[ \frac{d^k}{dt^k} \mathbf{C}e^{\mathbf{A}t}\boldsymbol{\xi}\bigg|_{t=0} =\mathbf{C}\mathbf{A}^{k}\boldsymbol{\xi}=\mathbf{0}, \qquad k=0,1,\ldots,n-1. \]

Hence \( \mathcal{O}\boldsymbol{\xi}=\mathbf{0} \). If \( \mathcal{O} \) has full column rank, this forces \( \boldsymbol{\xi}=\mathbf{0} \), so the Gramian has no nonzero null vector and is positive definite.

The finite-horizon Gramian can be computed from the matrix differential equation

\[ \boxed{\dot{\mathbf{W} }(t)=\mathbf{A}^{T}\mathbf{W}(t) +\mathbf{W}(t)\mathbf{A}+\mathbf{C}^{T}\mathbf{C},\qquad \mathbf{W}(0)=\mathbf{0}.} \]

Its solution is exactly \( \mathbf{W}(T)=\mathbf{W}_o(T) \). This follows by differentiating the integral form and by using

\[ \frac{d}{dt}\left(e^{\mathbf{A}^{T}t}\mathbf{C}^{T}\mathbf{C} e^{\mathbf{A}t}\right) =\mathbf{A}^{T}e^{\mathbf{A}^{T}t}\mathbf{C}^{T}\mathbf{C} e^{\mathbf{A}t} +e^{\mathbf{A}^{T}t}\mathbf{C}^{T}\mathbf{C} e^{\mathbf{A}t}\mathbf{A}. \]

If \( \mathbf{A} \) is Hurwitz, meaning every eigenvalue satisfies \( \operatorname{Re}(\lambda_i(\mathbf{A})) < 0 \), then the infinite-horizon Gramian exists:

\[ \boxed{\mathbf{W}_o=\int_{0}^{\infty} e^{\mathbf{A}^{T}t}\mathbf{C}^{T}\mathbf{C} e^{\mathbf{A}t}\,dt}. \]

In that case, it solves the continuous-time Lyapunov equation

\[ \boxed{\mathbf{A}^{T}\mathbf{W}_o+\mathbf{W}_o\mathbf{A} +\mathbf{C}^{T}\mathbf{C}=\mathbf{0}.} \]

This algebraic equation is the standard numerical route for stable continuous-time systems. If \( \mathbf{A} \) is not Hurwitz, the infinite integral may diverge, but the finite-horizon Gramian remains well-defined for every finite \( T > 0 \).

Let \( \mathbf{x}=\mathbf{T}\mathbf{z} \), where \( \mathbf{T} \) is nonsingular. Then

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

If \( \mathbf{W}_{o,x} \) is the Gramian in the \( \mathbf{x} \)-coordinates and \( \mathbf{W}_{o,z} \) is the Gramian in the \( \mathbf{z} \)-coordinates, then

\[ \boxed{\mathbf{W}_{o,z}=\mathbf{T}^{T}\mathbf{W}_{o,x}\mathbf{T}.} \]

This identity is obtained by preserving the same output energy:

\[ \mathbf{x}_0^{T}\mathbf{W}_{o,x}\mathbf{x}_0 =\mathbf{z}_0^{T}\mathbf{T}^{T}\mathbf{W}_{o,x}\mathbf{T}\mathbf{z}_0 =\mathbf{z}_0^{T}\mathbf{W}_{o,z}\mathbf{z}_0. \]

Therefore, rank and positive definiteness are coordinate-invariant, but raw eigenvalues can change under non-orthogonal scaling. In numerical work, state scaling should be treated carefully before interpreting Gramian eigenvalues as physical visibility measures.

Three common computational approaches are used in practice:

1. Quadrature: compute \( e^{\mathbf{A}^{T}t}\mathbf{C}^{T}\mathbf{C} e^{\mathbf{A}t} \) at many time samples and integrate.
2. Finite-horizon Lyapunov ODE: integrate \( \dot{\mathbf{W} }=\mathbf{A}^{T}\mathbf{W} +\mathbf{W}\mathbf{A}+\mathbf{C}^{T}\mathbf{C} \).
3. Infinite-horizon Lyapunov equation: if \( \mathbf{A} \) is Hurwitz, solve \( \mathbf{A}^{T}\mathbf{W} +\mathbf{W}\mathbf{A}+\mathbf{C}^{T}\mathbf{C}=\mathbf{0} \).

Relevant software ecosystems include NumPy/SciPy and python-control in Python, Eigen or Armadillo in C++, EJML or Apache Commons Math in Java, the MATLAB Control System Toolbox and Simulink, and Wolfram Mathematica for symbolic matrix exponentials and Lyapunov equations.

File: Chapter15_Lesson1.py

# Chapter15_Lesson1.py
# Observability Gramian for continuous-time LTI systems.
#
# Required libraries:
#   pip install numpy scipy matplotlib
#
# Optional control library:
#   pip install control
#
# The example uses the zero-input system:
#   x_dot = A x,  y = C x
# and computes:
#   W_o(T) = integral_0^T exp(A^T t) C^T C exp(A t) dt
# and, when A is Hurwitz:
#   A^T W_o + W_o A + C^T C = 0.

import numpy as np
from scipy.linalg import expm, solve_continuous_lyapunov, eigvals
from numpy.linalg import matrix_rank, eigvalsh


def observability_matrix(A: np.ndarray, C: np.ndarray) -> np.ndarray:
    """Return O = [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 finite_observability_gramian(A: np.ndarray, C: np.ndarray, T: float, steps: int = 5000) -> np.ndarray:
    """
    Numerically approximate W_o(T) using the trapezoidal rule.

    W_o(T) = integral_0^T exp(A^T t) C^T C exp(A t) dt.
    """
    if T <= 0:
        raise ValueError("T must be positive.")
    n = A.shape[0]
    Q = C.T @ C
    ts = np.linspace(0.0, T, steps + 1)
    W = np.zeros((n, n), dtype=float)

    for k, t in enumerate(ts):
        E = expm(A * t)
        integrand = E.T @ Q @ E
        weight = 0.5 if k == 0 or k == steps else 1.0
        W += weight * integrand

    W *= T / steps
    return 0.5 * (W + W.T)


def finite_observability_gramian_ode(A: np.ndarray, C: np.ndarray, T: float) -> np.ndarray:
    """
    Compute W_o(T) through the equivalent matrix differential equation:
      dW/dt = A^T W + W A + C^T C,  W(0)=0.

    For compactness, this routine uses scipy.integrate.solve_ivp.
    """
    from scipy.integrate import solve_ivp

    n = A.shape[0]
    Q = C.T @ C

    def rhs(_, w_flat):
        W = w_flat.reshape(n, n)
        Wdot = A.T @ W + W @ A + Q
        return Wdot.reshape(-1)

    sol = solve_ivp(rhs, (0.0, T), np.zeros(n * n), rtol=1e-10, atol=1e-12)
    W = sol.y[:, -1].reshape(n, n)
    return 0.5 * (W + W.T)


def infinite_observability_gramian(A: np.ndarray, C: np.ndarray) -> np.ndarray:
    """
    Compute W_o for Hurwitz A using:
      A^T W_o + W_o A + C^T C = 0.

    SciPy's solve_continuous_lyapunov solves:
      M X + X M^T = Q.
    Choose M=A^T and Q=-(C^T C).
    """
    Q = C.T @ C
    W = solve_continuous_lyapunov(A.T, -Q)
    return 0.5 * (W + W.T)


def output_energy(A: np.ndarray, C: np.ndarray, x0: np.ndarray, T: float) -> float:
    """Compute x0^T W_o(T) x0."""
    W = finite_observability_gramian(A, C, T)
    return float(x0.T @ W @ x0)


def main() -> None:
    # Example 1: observable stable second-order system
    A = np.array([[0.0, 1.0],
                  [-2.0, -3.0]])
    C = np.array([[1.0, 0.0]])

    print("Eigenvalues of A:", eigvals(A))
    O = observability_matrix(A, C)
    print("\nObservability matrix O:")
    print(O)
    print("rank(O) =", matrix_rank(O))

    T = 4.0
    W_T_quad = finite_observability_gramian(A, C, T)
    W_T_ode = finite_observability_gramian_ode(A, C, T)
    W_inf = infinite_observability_gramian(A, C)

    print(f"\nFinite-horizon W_o({T}) by quadrature:")
    print(W_T_quad)
    print("eigenvalues:", eigvalsh(W_T_quad))

    print(f"\nFinite-horizon W_o({T}) by Lyapunov ODE:")
    print(W_T_ode)
    print("difference norm:", np.linalg.norm(W_T_quad - W_T_ode))

    print("\nInfinite-horizon W_o:")
    print(W_inf)
    print("eigenvalues:", eigvalsh(W_inf))

    x0 = np.array([1.0, -0.5])
    print("\nInitial state x0:", x0)
    print("Output energy over [0,T] = x0^T W_o(T) x0 =", output_energy(A, C, x0, T))

    # Example 2: unobservable stable system
    A2 = np.array([[-1.0, 0.0],
                   [0.0, -2.0]])
    C2 = np.array([[1.0, 0.0]])
    W2 = finite_observability_gramian(A2, C2, T)
    print("\nUnobservable example W_o(T):")
    print(W2)
    print("rank(O2) =", matrix_rank(observability_matrix(A2, C2)))
    print("eigenvalues:", eigvalsh(W2))


if __name__ == "__main__":
    main()

      

File: Chapter15_Lesson1.cpp

// Chapter15_Lesson1.cpp
// Finite-horizon observability Gramian by integrating:
//   dW/dt = A^T W + W A + C^T C, W(0)=0.
//
// Library used:
//   Eigen 3
//
// Compile, for example:
//   g++ Chapter15_Lesson1.cpp -std=c++17 -O2 -I /path/to/eigen -o Chapter15_Lesson1

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

using Matrix = Eigen::MatrixXd;

Matrix observabilityMatrix(const Matrix& A, const Matrix& C) {
    const int n = A.rows();
    std::vector<Matrix> blocks;
    Matrix Ak = Matrix::Identity(n, n);

    for (int k = 0; k < n; ++k) {
        blocks.push_back(C * Ak);
        Ak = Ak * A;
    }

    Matrix O(C.rows() * n, n);
    for (int k = 0; k < n; ++k) {
        O.block(k * C.rows(), 0, C.rows(), n) = blocks[k];
    }
    return O;
}

Matrix gramianRhs(const Matrix& A, const Matrix& Q, const Matrix& W) {
    return A.transpose() * W + W * A + Q;
}

Matrix finiteObservabilityGramianRK4(const Matrix& A, const Matrix& C, double T, int steps) {
    const int n = A.rows();
    const double h = T / static_cast<double>(steps);
    const Matrix Q = C.transpose() * C;
    Matrix W = Matrix::Zero(n, n);

    for (int k = 0; k < steps; ++k) {
        Matrix K1 = gramianRhs(A, Q, W);
        Matrix K2 = gramianRhs(A, Q, W + 0.5 * h * K1);
        Matrix K3 = gramianRhs(A, Q, W + 0.5 * h * K2);
        Matrix K4 = gramianRhs(A, Q, W + h * K3);
        W += (h / 6.0) * (K1 + 2.0 * K2 + 2.0 * K3 + K4);
    }

    return 0.5 * (W + W.transpose());
}

int main() {
    Matrix A(2, 2);
    A << 0.0, 1.0,
        -2.0, -3.0;

    Matrix C(1, 2);
    C << 1.0, 0.0;

    Matrix O = observabilityMatrix(A, C);
    std::cout << "Observability matrix O:\n" << O << "\n\n";
    std::cout << "rank(O) = " << O.fullPivLu().rank() << "\n\n";

    const double T = 4.0;
    Matrix W = finiteObservabilityGramianRK4(A, C, T, 20000);

    std::cout << "Finite-horizon observability Gramian W_o(T):\n" << W << "\n\n";

    Eigen::SelfAdjointEigenSolver<Matrix> solver(W);
    std::cout << "Eigenvalues of W_o(T):\n" << solver.eigenvalues() << "\n\n";

    Eigen::Vector2d x0;
    x0 << 1.0, -0.5;
    double energy = (x0.transpose() * W * x0)(0, 0);
    std::cout << "Output energy x0^T W_o(T) x0 = " << energy << "\n";

    return 0;
}

      

File: Chapter15_Lesson1.java

// Chapter15_Lesson1.java
// Finite-horizon observability Gramian with only core Java.
//
// The code integrates:
//   dW/dt = A^T W + W A + C^T C, W(0)=0.
// For production-scale work, use EJML, Apache Commons Math, or ojAlgo.

public class Chapter15_Lesson1 {
    static double[][] transpose(double[][] A) {
        int m = A.length;
        int n = A[0].length;
        double[][] T = new double[n][m];
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                T[j][i] = A[i][j];
            }
        }
        return T;
    }

    static double[][] multiply(double[][] A, double[][] B) {
        int m = A.length;
        int p = A[0].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[][] add(double[][] A, double[][] B) {
        int m = A.length;
        int n = A[0].length;
        double[][] C = new double[m][n];

        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = A[i][j] + B[i][j];
            }
        }
        return C;
    }

    static double[][] scale(double[][] A, double s) {
        int m = A.length;
        int n = A[0].length;
        double[][] C = new double[m][n];

        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = s * A[i][j];
            }
        }
        return C;
    }

    static double[][] rhs(double[][] A, double[][] Q, double[][] W) {
        double[][] AT = transpose(A);
        return add(add(multiply(AT, W), multiply(W, A)), Q);
    }

    static double[][] finiteObservabilityGramianRK4(double[][] A, double[][] C, double T, int steps) {
        int n = A.length;
        double h = T / steps;
        double[][] W = new double[n][n];
        double[][] Q = multiply(transpose(C), C);

        for (int k = 0; k < steps; k++) {
            double[][] K1 = rhs(A, Q, W);
            double[][] K2 = rhs(A, Q, add(W, scale(K1, 0.5 * h)));
            double[][] K3 = rhs(A, Q, add(W, scale(K2, 0.5 * h)));
            double[][] K4 = rhs(A, Q, add(W, scale(K3, h)));

            double[][] increment = scale(add(add(K1, scale(K2, 2.0)), add(scale(K3, 2.0), K4)), h / 6.0);
            W = add(W, increment);
        }

        return W;
    }

    static void printMatrix(String name, double[][] A) {
        System.out.println(name);
        for (double[] row : A) {
            for (double value : row) {
                System.out.printf("% .8f ", value);
            }
            System.out.println();
        }
        System.out.println();
    }

    static double quadraticForm(double[] x, double[][] W) {
        double value = 0.0;
        for (int i = 0; i < x.length; i++) {
            for (int j = 0; j < x.length; j++) {
                value += x[i] * W[i][j] * x[j];
            }
        }
        return value;
    }

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

        double[][] C = {
            {1.0, 0.0}
        };

        double T = 4.0;
        double[][] W = finiteObservabilityGramianRK4(A, C, T, 20000);

        printMatrix("Finite-horizon observability Gramian W_o(T):", W);

        double[] x0 = {1.0, -0.5};
        System.out.printf("Output energy x0^T W_o(T) x0 = %.10f%n", quadraticForm(x0, W));
    }
}

      

File: Chapter15_Lesson1.m

% Chapter15_Lesson1.m
% Observability Gramian for continuous-time LTI systems.
%
% Libraries/toolboxes:
%   Base MATLAB: expm, integral, eig
%   Control System Toolbox: gram, ss, lyap
%   Simulink: optional zero-input output simulation model
%
% System:
%   x_dot = A*x, y = C*x

clear; clc;

A = [0 1; -2 -3];
B = [0; 0];      % not used for observability, included for ss object
C = [1 0];
D = 0;

T = 4.0;
Q = C' * C;

% Finite-horizon Gramian by numerical quadrature
integrand = @(t) expm(A' * t) * Q * expm(A * t);
Wo_T = integral(integrand, 0, T, 'ArrayValued', true);
Wo_T = 0.5 * (Wo_T + Wo_T');

disp('Finite-horizon observability Gramian Wo(T):');
disp(Wo_T);
disp('Eigenvalues of Wo(T):');
disp(eig(Wo_T));

% Observability matrix and rank
O = obsv(A, C);
disp('Observability matrix:');
disp(O);
fprintf('rank(O) = %d\n', rank(O));

% Infinite-horizon Gramian for Hurwitz A:
% lyap(A',Q) solves A'*W + W*A + Q = 0.
Wo_inf = lyap(A', Q);
Wo_inf = 0.5 * (Wo_inf + Wo_inf');

disp('Infinite-horizon observability Gramian:');
disp(Wo_inf);
disp('Eigenvalues of Wo_inf:');
disp(eig(Wo_inf));

% Same infinite-horizon result using the Control System Toolbox gram command.
sys = ss(A, B, C, D);
Wo_control_toolbox = gram(sys, 'o');
disp('Observability Gramian from gram(sys,''o''):');
disp(Wo_control_toolbox);

x0 = [1; -0.5];
energy_T = x0' * Wo_T * x0;
fprintf('Output energy over [0,T] = %.10f\n', energy_T);

% Optional Simulink model: zero-input state-space block with initial condition.
% This block simulates y(t)=C*x(t) for the chosen x0.
createSimulinkModel = false;

if createSimulinkModel
    model = 'Chapter15_Lesson1_ZeroInputOutput';
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    new_system(model);
    open_system(model);

    add_block('simulink/Continuous/State-Space', [model '/ZeroInputStateSpace']);
    set_param([model '/ZeroInputStateSpace'], ...
        'A', mat2str(A), ...
        'B', mat2str(B), ...
        'C', mat2str(C), ...
        'D', mat2str(D), ...
        'X0', mat2str(x0));

    add_block('simulink/Sources/Constant', [model '/ZeroInput']);
    set_param([model '/ZeroInput'], 'Value', '0');

    add_block('simulink/Sinks/Scope', [model '/OutputScope']);

    add_line(model, 'ZeroInput/1', 'ZeroInputStateSpace/1');
    add_line(model, 'ZeroInputStateSpace/1', 'OutputScope/1');

    set_param(model, 'StopTime', num2str(T));
    save_system(model);
end

      

File: Chapter15_Lesson1.nb

Notebook[{
Cell["Chapter15_Lesson1.nb", "Title"],
Cell["Observability Gramian for a continuous-time LTI system", "Text"],
Cell[BoxData[
"ClearAll[A, Cmat, T, Q, t, WoT, WoInf, Omat, x0, energy];\n\
A = { {0, 1}, {-2, -3} };\n\
Cmat = { {1, 0} };\n\
T = 4;\n\
Q = Transpose[Cmat].Cmat;\n\
\n\
WoT = Integrate[Transpose[MatrixExp[A t]].Q.MatrixExp[A t], {t, 0, T}] // FullSimplify;\n\
WoTNumeric = N[WoT];\n\
Print[\"Finite-horizon observability Gramian Wo(T):\"];\n\
Print[MatrixForm[WoTNumeric]];\n\
Print[\"Eigenvalues of Wo(T):\"];\n\
Print[N[Eigenvalues[WoT]]];\n\
\n\
Omat = Join[Cmat, Cmat.A];\n\
Print[\"Observability matrix:\"];\n\
Print[MatrixForm[Omat]];\n\
Print[\"rank(O) = \", MatrixRank[Omat]];\n\
\n\
WoInf = LyapunovSolve[Transpose[A], Q];\n\
Print[\"Infinite-horizon observability Gramian:\"];\n\
Print[MatrixForm[N[WoInf]]];\n\
Print[\"Residual A^T W + W A + C^T C:\"];\n\
Print[MatrixForm[N[Transpose[A].WoInf + WoInf.A + Q]]];\n\
\n\
x0 = {1, -1/2};\n\
energy = x0.WoT.x0;\n\
Print[\"Output energy over [0,T] = \", N[energy]];"
], "Input"]
}]

      

Problem 1 (Double Integrator Gramian): Consider \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \), \( \mathbf{y}=\mathbf{C}\mathbf{x} \), with

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

Compute \( \mathbf{W}_o(T) \) and determine whether the pair is observable for \( T > 0 \).

Solution: Since

\[ e^{\mathbf{A}t}=\begin{bmatrix}1&t\\0&1\end{bmatrix}, \qquad \mathbf{C}e^{\mathbf{A}t}=\begin{bmatrix}1&t\end{bmatrix}, \]

we obtain

\[ \mathbf{W}_o(T)=\int_0^T \begin{bmatrix}1\\t\end{bmatrix} \begin{bmatrix}1&t\end{bmatrix}\,dt = \begin{bmatrix} T & \frac{T^2}{2}\\ \frac{T^2}{2} & \frac{T^3}{3} \end{bmatrix}. \]

Its determinant is

\[ \det\mathbf{W}_o(T)=\frac{T^4}{12}. \]

Thus \( \det\mathbf{W}_o(T) > 0 \) for \( T > 0 \), so \( \mathbf{W}_o(T)\succ\mathbf{0} \) and the pair is observable.

Problem 2 (Unobservable Direction): Let

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

Compute the finite-horizon Gramian and identify the unobservable subspace.

Solution: We have

\[ \mathbf{C}e^{\mathbf{A}t}=\begin{bmatrix}e^{-t}&0\end{bmatrix}. \]

Therefore,

\[ \mathbf{W}_o(T)= \int_0^T \begin{bmatrix}e^{-t}\\0\end{bmatrix} \begin{bmatrix}e^{-t}&0\end{bmatrix}\,dt = \begin{bmatrix} \frac{1-e^{-2T} }{2} & 0\\ 0 & 0 \end{bmatrix}. \]

The null space is \( \operatorname{span}\{[0\;\;1]^T\} \). The second state never appears in the output, so the system is not observable.

Problem 3 (Infinite-Horizon Lyapunov Equation): For

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

compute the infinite-horizon observability Gramian.

Solution: Since \( \mathbf{A} \) is Hurwitz, solve \( \mathbf{A}^{T}\mathbf{W}+\mathbf{W}\mathbf{A} +\mathbf{C}^{T}\mathbf{C}=\mathbf{0} \). Let \( \mathbf{W}=[w_{ij}] \). Since \( \mathbf{C}^{T}\mathbf{C}= \begin{bmatrix}1&1\\1&1\end{bmatrix} \), the scalar equations are

\[ -2w_{11}+1=0,\qquad -3w_{12}+1=0,\qquad -4w_{22}+1=0. \]

Thus

\[ \mathbf{W}_o= \begin{bmatrix} \frac{1}{2} & \frac{1}{3}\\ \frac{1}{3} & \frac{1}{4} \end{bmatrix}. \]

The determinant is \( \frac{1}{72} \), so the Gramian is positive definite and the pair is observable.

Problem 4 (Null-Space Equivalence): Prove that \( \mathcal{N}(\mathbf{W}_o(T)) \) is the unobservable subspace for an LTI system.

Solution: If \( \boldsymbol{\xi}\in\mathcal{N}(\mathbf{W}_o(T)) \), then

\[ 0=\boldsymbol{\xi}^{T}\mathbf{W}_o(T)\boldsymbol{\xi} =\int_0^T \|\mathbf{C}e^{\mathbf{A}t}\boldsymbol{\xi}\|_2^2\,dt. \]

The integrand is continuous and nonnegative, so it is identically zero on the interval. Hence \( \mathbf{C}e^{\mathbf{A}t}\boldsymbol{\xi}=\mathbf{0} \) for all measured times. This means the initial condition \( \boldsymbol{\xi} \) generates no output and cannot be distinguished from the zero initial state. Conversely, any initial-state direction that generates identically zero output makes the above integral zero, so it belongs to \( \mathcal{N}(\mathbf{W}_o(T)) \).

Problem 5 (Coordinate Change): Let \( \mathbf{x}=\mathbf{T}\mathbf{z} \). Prove that the Gramian transforms as \( \mathbf{W}_{o,z}=\mathbf{T}^{T} \mathbf{W}_{o,x}\mathbf{T} \).

Solution: The output generated by \( \mathbf{x}_0=\mathbf{T}\mathbf{z}_0 \) is physically the same signal, regardless of coordinates. Therefore,

\[ \mathbf{x}_0^{T}\mathbf{W}_{o,x}\mathbf{x}_0 =(\mathbf{T}\mathbf{z}_0)^{T}\mathbf{W}_{o,x} (\mathbf{T}\mathbf{z}_0) =\mathbf{z}_0^{T}\mathbf{T}^{T}\mathbf{W}_{o,x} \mathbf{T}\mathbf{z}_0. \]

Since this equality holds for every \( \mathbf{z}_0 \), the transformed Gramian must be \( \mathbf{W}_{o,z}=\mathbf{T}^{T} \mathbf{W}_{o,x}\mathbf{T} \).

The observability Gramian \( \mathbf{W}_o(T)=\int_0^T e^{\mathbf{A}^{T}t}\mathbf{C}^{T}\mathbf{C}e^{\mathbf{A}t}\,dt \) converts the question of state visibility into a symmetric positive-semidefinite matrix. Its quadratic form equals the zero-input output energy generated by an initial state. For LTI systems, complete observability is equivalent to \( \mathbf{W}_o(T)\succ\mathbf{0} \) for every \( T > 0 \), and if \( \mathbf{A} \) is Hurwitz, the infinite-horizon Gramian is obtained from a continuous-time Lyapunov equation.