Consider the continuous-time LTI system already used in the previous lesson:

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

In this lesson we study only the output produced by the initial state, so the input is set to zero: \( \mathbf{u}(t)=\mathbf{0} \). Then

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

The natural measurement of how strongly the initial state appears in the output over the interval \( [0,T] \) is the squared signal norm:

\[ E_o(\mathbf{x}_0,T)=\int_0^T \mathbf{y}(t)^T\mathbf{y}(t)\,dt,\qquad T>0. \]

Large output energy means the initial condition leaves a strong trace in the measured signal. Small output energy means the initial condition is weakly visible; zero output energy means it is completely invisible over that interval.

flowchart TD
  A["Initial state x0"] --> B["Zero-input state motion"]
  B --> C["Output signal y(t)"]
  C --> D["Integrate squared output over 0 to T"]
  D --> E["Output energy Eo(x0,T)"]
  E --> F["Large: strongly visible \nstate direction"]
  E --> G["Small: weakly visible \nstate direction"]
  E --> H["Zero: invisible \nstate direction"]
        

For a finite interval \( [0,T] \), the observability Gramian is

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

This matrix is symmetric and positive semidefinite. Symmetry follows because the integrand is symmetric:

\[ \left(e^{\mathbf{A}^T \tau}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}\tau}\right)^T =e^{\mathbf{A}^T \tau}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}\tau}. \]

Positive semidefiniteness follows because, for any vector \( \boldsymbol{\xi} \),

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

Therefore \( \mathbf{W}_o(0,T) \) defines an energy geometry on the state space. It is not an arbitrary matrix: it is the accumulated sensitivity of the output to the initial condition.

Theorem. For the zero-input LTI system \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \), \( \mathbf{y}=\mathbf{C}\mathbf{x} \), the output energy generated by \( \mathbf{x}_0 \) on \( [0,T] \) satisfies

\[ E_o(\mathbf{x}_0,T)=\mathbf{x}_0^T \mathbf{W}_o(0,T)\mathbf{x}_0. \]

Proof. Since the input is zero,

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

Hence

\[ \begin{aligned} E_o(\mathbf{x}_0,T) &=\int_0^T \mathbf{y}(t)^T\mathbf{y}(t)\,dt \\ &=\int_0^T \left(\mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0\right)^T \left(\mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0\right)\,dt \\ &=\int_0^T \mathbf{x}_0^T e^{\mathbf{A}^Tt}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}t}\mathbf{x}_0\,dt \\ &=\mathbf{x}_0^T \left(\int_0^T e^{\mathbf{A}^Tt}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}t}\,dt\right)\mathbf{x}_0 \\ &=\mathbf{x}_0^T\mathbf{W}_o(0,T)\mathbf{x}_0. \end{aligned} \]

This completes the proof. The theorem is important because it turns a signal-space integral into a finite-dimensional quadratic form in the initial state.

Because the integrand is nonnegative, zero energy has a strong implication:

\[ E_o(\mathbf{x}_0,T)=0 \quad \Longleftrightarrow \quad \mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0=\mathbf{0} \text{ for all } t\in[0,T]. \]

In other words, an initial state produces zero output energy exactly when its entire zero-input output trajectory is identically zero. Using the Gramian:

\[ E_o(\mathbf{x}_0,T)=0 \quad \Longleftrightarrow \quad \mathbf{x}_0\in\mathcal{N}\bigl(\mathbf{W}_o(0,T)\bigr). \]

If \( \mathbf{W}_o(0,T) \) is positive definite, then the only zero-energy initial state is the zero state:

\[ \mathbf{W}_o(0,T)>0 \quad \Longrightarrow \quad E_o(\mathbf{x}_0,T)>0 \text{ for every } \mathbf{x}_0\ne\mathbf{0}. \]

This is the energy version of observability. If every nonzero initial state produces nonzero output energy on the observation interval, then no nonzero state is hidden from the sensors.

The finite-horizon Gramian and the finite-dimensional observability matrix are consistent. For LTI systems, when \( T>0 \),

\[ \operatorname{rank}\mathbf{W}_o(0,T)=n \quad \Longleftrightarrow \quad \operatorname{rank}\begin{bmatrix} \mathbf{C} \\ \mathbf{C}\mathbf{A} \\ \mathbf{C}\mathbf{A}^2 \\ \vdots \\ \mathbf{C}\mathbf{A}^{n-1} \end{bmatrix}=n. \]

Suppose \( \mathbf{W}_o(0,T) \) has eigenvalue decomposition

\[ \mathbf{W}_o(0,T)=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^T,\qquad \boldsymbol{\Lambda}=\operatorname{diag}(\lambda_1,\ldots,\lambda_n), \qquad \lambda_i\ge 0. \]

If \( \mathbf{x}_0=\sum_i \alpha_i\mathbf{v}_i \), then

\[ E_o(\mathbf{x}_0,T)=\sum_{i=1}^n \lambda_i\alpha_i^2. \]

Therefore an eigenvector with large \( \lambda_i \) is a state direction that creates large output energy. An eigenvector with very small \( \lambda_i \) is weakly observed. If \( \lambda_i=0 \), that direction is unobservable on the interval.

flowchart TD
  A["Gramian eigenpairs"] --> B["Large eigenvalue direction"]
  A --> C["Small eigenvalue direction"]
  A --> D["Zero eigenvalue direction"]
  B --> E["Large output energy for \nsame state amplitude"]
  C --> F["Weak output trace and \nsensitive reconstruction"]
  D --> G["No output trace: \nunobservable direction"]
  E --> H["Energy level set: \nx0 prime W x0 equals rho squared"]
  F --> H
  G --> H
        

The energy-bounded set

\[ \mathcal{E}_\rho= \left\{\mathbf{x}_0: \mathbf{x}_0^T\mathbf{W}_o(0,T)\mathbf{x}_0\le\rho^2 \right\} \]

is an ellipsoid when the Gramian is positive definite. Directions with large Gramian eigenvalues have short ellipsoid axes because a small initial displacement already produces noticeable output energy.

If \( \mathbf{A} \) is asymptotically stable, the infinite-horizon output energy is finite:

\[ E_o(\mathbf{x}_0,\infty)= \int_0^\infty \left\|\mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0\right\|_2^2\,dt. \]

The infinite-horizon observability Gramian is

\[ \mathbf{W}_o= \int_0^\infty e^{\mathbf{A}^T \tau} \mathbf{C}^T\mathbf{C}e^{\mathbf{A}\tau}\,d\tau. \]

It satisfies the continuous-time Lyapunov equation

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

Proof. Differentiate the integrand:

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

Integrating from \( 0 \) to \( \infty \) and using stability of \( \mathbf{A} \),

\[ \begin{aligned} \mathbf{A}^T\mathbf{W}_o+\mathbf{W}_o\mathbf{A} &=\int_0^\infty \frac{d}{d\tau} \left(e^{\mathbf{A}^T \tau} \mathbf{C}^T\mathbf{C}e^{\mathbf{A}\tau}\right)d\tau \\ &=\lim_{\tau\to\infty} e^{\mathbf{A}^T \tau}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}\tau}-\mathbf{C}^T\mathbf{C} \\ &=-\mathbf{C}^T\mathbf{C}. \end{aligned} \]

Hence for stable systems, \( E_o(\mathbf{x}_0,\infty)= \mathbf{x}_0^T\mathbf{W}_o\mathbf{x}_0 \).

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

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

If \( \bar{\mathbf{W} }_o \) denotes the observability Gramian in the \( \mathbf{z} \) coordinates, then

\[ \bar{\mathbf{W} }_o=\mathbf{T}^T\mathbf{W}_o\mathbf{T}. \]

The numerical matrix changes, but the physical output energy does not:

\[ \mathbf{z}_0^T\bar{\mathbf{W} }_o\mathbf{z}_0 =\mathbf{z}_0^T\mathbf{T}^T\mathbf{W}_o\mathbf{T}\mathbf{z}_0 =\mathbf{x}_0^T\mathbf{W}_o\mathbf{x}_0. \]

The practical computation depends on the horizon and on whether \( \mathbf{A} \) is stable.

\[ \begin{array}{ll} \text{finite horizon:} & \mathbf{W}_o(0,T)=\int_0^T e^{\mathbf{A}^T \tau}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}\tau}d\tau, \\[4pt] \text{stable infinite horizon:} & \mathbf{A}^T\mathbf{W}_o+\mathbf{W}_o\mathbf{A} =-\mathbf{C}^T\mathbf{C}. \end{array} \]

Useful software support includes Python scipy.linalg.solve_continuous_lyapunov and scipy.linalg.expm; C++ libraries such as Eigen, Armadillo, and SLICOT/LAPACK wrappers; Java libraries such as EJML and Apache Commons Math; MATLAB/Simulink tools such as lyap, gram, ss, and initial; and Wolfram Mathematica functions such as MatrixExp, NIntegrate, and KroneckerProduct.

Chapter15_Lesson2.py

# Chapter15_Lesson2.py
# Output Energy Generated by Initial States
# Requires: numpy, scipy

import numpy as np
from scipy.linalg import expm, solve_continuous_lyapunov, eigvalsh
from scipy.integrate import quad


def finite_observability_gramian(A: np.ndarray, C: np.ndarray, T: float, n_grid: int = 2000) -> np.ndarray:
    """Approximate W_o(0,T)=int_0^T exp(A^T t) C^T C exp(A t) dt."""
    if T <= 0:
        raise ValueError("T must be positive.")
    Q = C.T @ C
    grid = np.linspace(0.0, T, n_grid)
    dt = grid[1] - grid[0]
    W = np.zeros((A.shape[0], A.shape[0]))

    for k, t in enumerate(grid):
        Phi = expm(A * t)
        weight = 0.5 if k == 0 or k == n_grid - 1 else 1.0
        W += weight * (Phi.T @ Q @ Phi)

    return dt * W


def output_energy_by_quadrature(A: np.ndarray, C: np.ndarray, x0: np.ndarray, T: float) -> float:
    """Compute int_0^T ||C exp(A t) x0||_2^2 dt."""
    def integrand(t: float) -> float:
        y = C @ expm(A * t) @ x0
        return float(y.T @ y)

    value, _ = quad(integrand, 0.0, T, epsabs=1e-10, epsrel=1e-10, limit=200)
    return value


def main() -> None:
    A = np.array([
        [-1.0,  0.2,  0.0],
        [ 0.0, -0.6,  1.0],
        [ 0.0, -1.2, -0.8]
    ], dtype=float)

    C = np.array([
        [1.0, 0.0, 0.0],
        [0.0, 0.5, 1.0]
    ], dtype=float)

    x0 = np.array([1.0, -0.4, 0.8], dtype=float)
    T = 5.0

    W_T = finite_observability_gramian(A, C, T)
    E_T_quadratic = float(x0.T @ W_T @ x0)
    E_T_integral = output_energy_by_quadrature(A, C, x0, T)

    # Infinite-horizon Gramian for stable A:
    # A^T W + W A = -C^T C
    W_inf = solve_continuous_lyapunov(A.T, -(C.T @ C))
    E_inf = float(x0.T @ W_inf @ x0)

    print("Finite-horizon observability Gramian W_o(0,T):")
    print(W_T)
    print("\nEigenvalues of W_o(0,T):", eigvalsh(W_T))
    print("\nEnergy from x0^T W_o x0:", E_T_quadratic)
    print("Energy from direct output integration:", E_T_integral)
    print("Absolute difference:", abs(E_T_quadratic - E_T_integral))

    print("\nInfinite-horizon observability Gramian W_o:")
    print(W_inf)
    print("Infinite-horizon output energy:", E_inf)


if __name__ == "__main__":
    main()
      

Chapter15_Lesson2.cpp

// Chapter15_Lesson2.cpp
// Output Energy Generated by Initial States
// Requires Eigen: https://eigen.tuxfamily.org
// Compile: g++ -std=c++17 Chapter15_Lesson2.cpp -O2 -I /path/to/eigen -o Chapter15_Lesson2

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

using Matrix = Eigen::MatrixXd;
using Vector = Eigen::VectorXd;

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

Matrix rk4_step_gramian(const Matrix& A, const Matrix& Q, const Matrix& W, double h) {
    Matrix k1 = gramian_rhs(A, Q, W);
    Matrix k2 = gramian_rhs(A, Q, W + 0.5 * h * k1);
    Matrix k3 = gramian_rhs(A, Q, W + 0.5 * h * k2);
    Matrix k4 = gramian_rhs(A, Q, W + h * k3);
    return W + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
}

Vector rk4_step_state(const Matrix& A, const Vector& x, double h) {
    Vector k1 = A * x;
    Vector k2 = A * (x + 0.5 * h * k1);
    Vector k3 = A * (x + 0.5 * h * k2);
    Vector k4 = A * (x + h * k3);
    return x + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
}

double output_density(const Matrix& C, const Vector& x) {
    Vector y = C * x;
    return y.squaredNorm();
}

int main() {
    Matrix A(3, 3);
    A << -1.0,  0.2,  0.0,
          0.0, -0.6,  1.0,
          0.0, -1.2, -0.8;

    Matrix C(2, 3);
    C << 1.0, 0.0, 0.0,
         0.0, 0.5, 1.0;

    Vector x0(3);
    x0 << 1.0, -0.4, 0.8;

    const double T = 5.0;
    const int steps = 10000;
    const double h = T / static_cast<double>(steps);

    Matrix Q = C.transpose() * C;
    Matrix W = Matrix::Zero(3, 3);
    Vector x = x0;

    double E_direct = 0.0;
    double previous = output_density(C, x);

    for (int k = 0; k < steps; ++k) {
        W = rk4_step_gramian(A, Q, W, h);
        Vector x_next = rk4_step_state(A, x, h);
        double next = output_density(C, x_next);
        E_direct += 0.5 * h * (previous + next);
        x = x_next;
        previous = next;
    }

    double E_quad = x0.transpose() * W * x0;

    std::cout << std::fixed << std::setprecision(10);
    std::cout << "Finite-horizon observability Gramian W_o(0,T):\n" << W << "\n\n";
    std::cout << "Energy from x0^T W_o x0: " << E_quad << "\n";
    std::cout << "Energy from direct output integration: " << E_direct << "\n";
    std::cout << "Absolute difference: " << std::abs(E_quad - E_direct) << "\n";

    return 0;
}
      

Chapter15_Lesson2.java

// Chapter15_Lesson2.java
// Output Energy Generated by Initial States
// From-scratch RK4 integration; no external Java libraries required.
// Compile: javac Chapter15_Lesson2.java
// Run:     java Chapter15_Lesson2

public class Chapter15_Lesson2 {
    static final int N = 3;
    static final int P = 2;

    static double[][] zeros() {
        return new double[N][N];
    }

    static double[][] transpose(double[][] A) {
        double[][] T = zeros();
        for (int i = 0; i < N; i++)
            for (int j = 0; j < N; j++)
                T[i][j] = A[j][i];
        return T;
    }

    static double[][] add(double[][] A, double[][] B) {
        double[][] R = zeros();
        for (int i = 0; i < N; i++)
            for (int j = 0; j < N; j++)
                R[i][j] = A[i][j] + B[i][j];
        return R;
    }

    static double[][] scale(double[][] A, double s) {
        double[][] R = zeros();
        for (int i = 0; i < N; i++)
            for (int j = 0; j < N; j++)
                R[i][j] = s * A[i][j];
        return R;
    }

    static double[][] multiply(double[][] A, double[][] B) {
        double[][] R = zeros();
        for (int i = 0; i < N; i++)
            for (int j = 0; j < N; j++)
                for (int k = 0; k < N; k++)
                    R[i][j] += A[i][k] * B[k][j];
        return R;
    }

    static double[] matvec(double[][] A, double[] x) {
        double[] y = new double[N];
        for (int i = 0; i < N; i++)
            for (int j = 0; j < N; j++)
                y[i] += A[i][j] * x[j];
        return y;
    }

    static double[][] ctC(double[][] C) {
        double[][] Q = zeros();
        for (int i = 0; i < N; i++)
            for (int j = 0; j < N; j++)
                for (int r = 0; r < P; r++)
                    Q[i][j] += C[r][i] * C[r][j];
        return Q;
    }

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

    static double[][] rk4W(double[][] A, double[][] Q, double[][] W, double h) {
        double[][] k1 = rhsW(A, Q, W);
        double[][] k2 = rhsW(A, Q, add(W, scale(k1, 0.5 * h)));
        double[][] k3 = rhsW(A, Q, add(W, scale(k2, 0.5 * h)));
        double[][] k4 = rhsW(A, Q, add(W, scale(k3, h)));
        return add(W, scale(add(add(k1, scale(k2, 2.0)), add(scale(k3, 2.0), k4)), h / 6.0));
    }

    static double[] rk4X(double[][] A, double[] x, double h) {
        double[] k1 = matvec(A, x);
        double[] x2 = new double[N];
        for (int i = 0; i < N; i++) x2[i] = x[i] + 0.5 * h * k1[i];

        double[] k2 = matvec(A, x2);
        double[] x3 = new double[N];
        for (int i = 0; i < N; i++) x3[i] = x[i] + 0.5 * h * k2[i];

        double[] k3 = matvec(A, x3);
        double[] x4 = new double[N];
        for (int i = 0; i < N; i++) x4[i] = x[i] + h * k3[i];

        double[] k4 = matvec(A, x4);
        double[] xn = new double[N];
        for (int i = 0; i < N; i++)
            xn[i] = x[i] + h * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
        return xn;
    }

    static double density(double[][] C, double[] x) {
        double sum = 0.0;
        for (int r = 0; r < P; r++) {
            double yr = 0.0;
            for (int j = 0; j < N; j++) yr += C[r][j] * x[j];
            sum += yr * yr;
        }
        return sum;
    }

    static double quadratic(double[][] W, double[] x) {
        double[] Wx = matvec(W, x);
        double s = 0.0;
        for (int i = 0; i < N; i++) s += x[i] * Wx[i];
        return s;
    }

    static void printMatrix(double[][] W) {
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++)
                System.out.printf("%14.8f ", W[i][j]);
            System.out.println();
        }
    }

    public static void main(String[] args) {
        double[][] A = {
            {-1.0,  0.2,  0.0},
            { 0.0, -0.6,  1.0},
            { 0.0, -1.2, -0.8}
        };

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

        double[] x0 = {1.0, -0.4, 0.8};
        double T = 5.0;
        int steps = 10000;
        double h = T / steps;

        double[][] Q = ctC(C);
        double[][] W = zeros();
        double[] x = x0.clone();

        double EDirect = 0.0;
        double previous = density(C, x);

        for (int k = 0; k < steps; k++) {
            W = rk4W(A, Q, W, h);
            double[] xNext = rk4X(A, x, h);
            double next = density(C, xNext);
            EDirect += 0.5 * h * (previous + next);
            x = xNext;
            previous = next;
        }

        double EQuad = quadratic(W, x0);
        System.out.println("Finite-horizon observability Gramian W_o(0,T):");
        printMatrix(W);
        System.out.printf("%nEnergy from x0^T W_o x0: %.10f%n", EQuad);
        System.out.printf("Energy from direct output integration: %.10f%n", EDirect);
        System.out.printf("Absolute difference: %.10e%n", Math.abs(EQuad - EDirect));
    }
}
      

The script computes the Gramian through both a Lyapunov equation and a finite-horizon matrix differential equation. The initial command is the MATLAB/Simulink-style zero-input response calculation: it is the same response that a Simulink State-Space block would produce when initialized at \( \mathbf{x}_0 \) with zero input.

Chapter15_Lesson2.m

% Chapter15_Lesson2.m
% Output Energy Generated by Initial States
% Requires: Control System Toolbox for lyap, ss, and initial.

clear; clc; close all;

A = [-1.0  0.2  0.0;
      0.0 -0.6  1.0;
      0.0 -1.2 -0.8];

C = [1.0 0.0 0.0;
     0.0 0.5 1.0];

x0 = [1.0; -0.4; 0.8];
T  = 5.0;
Q  = C' * C;
n  = size(A, 1);

% Infinite-horizon observability Gramian:
% A' * W + W * A + C' * C = 0
Wo_inf = lyap(A', Q);
E_inf = x0' * Wo_inf * x0;

% Finite-horizon Gramian from dW/dt = A'W + WA + C'C, W(0)=0.
odefun = @(t,wvec) reshape(A' * reshape(wvec,n,n) + reshape(wvec,n,n) * A + Q, [], 1);
[tW, wsol] = ode45(odefun, [0 T], zeros(n*n, 1));
Wo_T = reshape(wsol(end,:).', n, n);
Wo_T = 0.5 * (Wo_T + Wo_T');

E_T_quad = x0' * Wo_T * x0;

% Simulink-equivalent zero-input initial response.
% A Simulink State-Space block with this A,C and initial state x0 gives the same y(t).
Bdummy = zeros(n, 1);
Ddummy = zeros(size(C, 1), 1);
sys = ss(A, Bdummy, C, Ddummy);

time = linspace(0, T, 1000);
[y, t, x] = initial(sys, x0, time);
E_T_initial = trapz(t, sum(y.^2, 2));

disp('Finite-horizon observability Gramian W_o(0,T):');
disp(Wo_T);
disp('Eigenvalues of W_o(0,T):');
disp(eig(Wo_T).');
fprintf('Energy from x0'' W_o x0: %.12f\n', E_T_quad);
fprintf('Energy from initial-response integration: %.12f\n', E_T_initial);
fprintf('Absolute difference: %.3e\n', abs(E_T_quad - E_T_initial));

disp('Infinite-horizon observability Gramian W_o:');
disp(Wo_inf);
fprintf('Infinite-horizon output energy: %.12f\n', E_inf);

figure;
plot(t, y, 'LineWidth', 1.5);
grid on;
xlabel('Time');
ylabel('Output components');
title('Zero-input output response y(t)=C exp(A t) x_0');
legend('y_1', 'y_2');
      

Chapter15_Lesson2.nb

(* Chapter15_Lesson2.nb *)
(* Output Energy Generated by Initial States *)

ClearAll["Global`*"];

A = { {-1.0, 0.2, 0.0}, {0.0, -0.6, 1.0}, {0.0, -1.2, -0.8} };
Cmat = { {1.0, 0.0, 0.0}, {0.0, 0.5, 1.0} };
x0 = {1.0, -0.4, 0.8};
Tfinal = 5.0;
n = Length[A];
Q = Transpose[Cmat].Cmat;

vec[M_] := Flatten[Transpose[M]];
unvec[v_] := Transpose[Partition[v, n]];

(* Infinite-horizon observability Gramian:
   Transpose[A].W + W.A == -Transpose[Cmat].Cmat *)
K = KroneckerProduct[IdentityMatrix[n], Transpose[A]] +
    KroneckerProduct[Transpose[A], IdentityMatrix[n]];
WoInf = unvec[LinearSolve[K, -vec[Q]]];

Phi[t_] := MatrixExp[A t];

WoFinite = Table[
   NIntegrate[(Transpose[Phi[s]].Q.Phi[s])[[i, j]], {s, 0, Tfinal},
    AccuracyGoal -> 10, PrecisionGoal -> 10],
   {i, n}, {j, n}
];

EnergyFromGramian = x0.WoFinite.x0;

y[t_] := Cmat.Phi[t].x0;
EnergyFromOutput = NIntegrate[y[s].y[s], {s, 0, Tfinal},
   AccuracyGoal -> 10, PrecisionGoal -> 10];

{MatrixForm[WoFinite],
 Eigenvalues[WoFinite],
 EnergyFromGramian,
 EnergyFromOutput,
 Abs[EnergyFromGramian - EnergyFromOutput],
 MatrixForm[WoInf],
 x0.WoInf.x0}
      

Problem 1 (Energy Equality): For \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \), \( \mathbf{y}=\mathbf{C}\mathbf{x} \), prove that \( E_o(\mathbf{x}_0,T)=\mathbf{x}_0^T\mathbf{W}_o(0,T)\mathbf{x}_0 \).

Solution: The zero-input solution is \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0 \). Thus

\[ \begin{aligned} E_o(\mathbf{x}_0,T) &=\int_0^T \left(\mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0\right)^T \left(\mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0\right)dt \\ &=\mathbf{x}_0^T \left(\int_0^T e^{\mathbf{A}^Tt}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}t}dt\right)\mathbf{x}_0 =\mathbf{x}_0^T\mathbf{W}_o(0,T)\mathbf{x}_0. \end{aligned} \]

Problem 2 (Unobservable Direction): Let \( \mathbf{v}\ne\mathbf{0} \) satisfy \( \mathbf{W}_o(0,T)\mathbf{v}=\mathbf{0} \). Show that \( \mathbf{v} \) produces zero output on \( [0,T] \).

Solution: Since the Gramian is positive semidefinite,

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

A continuous nonnegative function with zero integral over a finite interval must be zero everywhere on that interval. Hence \( \mathbf{C}e^{\mathbf{A}t}\mathbf{v}=\mathbf{0} \) for all \( t\in[0,T] \).

Problem 3 (Scalar System): Consider \( \dot{x}=ax \), \( y=cx \), where \( a\ne0 \). Find \( W_o(0,T) \) and \( E_o(x_0,T) \).

Solution:

\[ W_o(0,T)=\int_0^T c^2e^{2a\tau}d\tau =\frac{c^2}{2a}\left(e^{2aT}-1\right). \]

\[ E_o(x_0,T)=x_0^2W_o(0,T) =\frac{c^2x_0^2}{2a}\left(e^{2aT}-1\right). \]

If \( a<0 \), the infinite-horizon value is

\[ W_o=\int_0^\infty c^2e^{2a\tau}d\tau =-\frac{c^2}{2a}. \]

Problem 4 (Lyapunov Equation): Assume \( \mathbf{A} \) is asymptotically stable. Prove that the infinite-horizon observability Gramian satisfies \( \mathbf{A}^T\mathbf{W}_o+\mathbf{W}_o\mathbf{A} =-\mathbf{C}^T\mathbf{C} \).

Solution: Starting from

\[ \mathbf{W}_o=\int_0^\infty e^{\mathbf{A}^T\tau}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}\tau}d\tau, \]

multiply by \( \mathbf{A}^T \) on the left and by \( \mathbf{A} \) on the right:

\[ \mathbf{A}^T\mathbf{W}_o+\mathbf{W}_o\mathbf{A} = \int_0^\infty \frac{d}{d\tau} \left(e^{\mathbf{A}^T\tau}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}\tau}\right)d\tau. \]

Since \( \mathbf{A} \) is stable, the upper boundary term vanishes, while the lower boundary term is \( \mathbf{C}^T\mathbf{C} \). Therefore

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

Problem 5 (Eigenvalue Meaning): Let \( \mathbf{W}_o=\operatorname{diag}(100,1,0.01) \). Compare the output energy for unit initial states along the three coordinate axes.

Solution: For \( \mathbf{e}_i \), \( E_o(\mathbf{e}_i)=\mathbf{e}_i^T\mathbf{W}_o\mathbf{e}_i \). Thus

\[ E_o(\mathbf{e}_1)=100,\qquad E_o(\mathbf{e}_2)=1,\qquad E_o(\mathbf{e}_3)=0.01. \]

The first direction is strongly visible, the second is moderately visible, and the third is weakly visible. Reconstructing the third component from noisy output data would be much more sensitive.

The output energy generated by an initial condition is exactly the quadratic form induced by the observability Gramian: \( E_o(\mathbf{x}_0,T)=\mathbf{x}_0^T\mathbf{W}_o(0,T)\mathbf{x}_0 \). The Gramian is symmetric positive semidefinite, its null space contains initial states invisible to the output, and its eigenvalues quantify how strongly each state direction appears in the measured signal. For stable systems, the infinite-horizon Gramian solves a Lyapunov equation, which is the basis of efficient computation in control software.