In previous lessons, the observability Gramian \( \mathbf{W}_o \) was introduced as a matrix that measures how strongly an initial state \( \mathbf{x}_0 \) appears in the output signal \( \mathbf{y}(t) \). Earlier, the controllability Gramian \( \mathbf{W}_c \) measured how easily an input can move the state. The central fact of this lesson is that these two objects are not merely analogous: they are dual mathematical objects.

For the continuous-time 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 of \( (\mathbf{A},\mathbf{C}) \) is equivalent to controllability of the dual pair \( (\mathbf{A}^T,\mathbf{C}^T) \). Therefore, many results about actuation directions can be reinterpreted as results about sensing directions.

flowchart TD
  A["Original System"] --> B["State matrix A"]
  A --> C["Output matrix C"]
  B --> D["Observability pair: (A, C)"]
  C --> D
  D --> E["Observability Gramian Wo"]
  A --> F["Dual System"]
  C --> F
  F --> G["Dual pair: (A^T, C^T)"]
  G --> H["Controllability Gramian Wc_dual"]
  E --> I["Wo = Wc_dual"]
  H --> I
        

Consider the autonomous output response generated only by the initial state:

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

The state trajectory and output are

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

The output energy over the finite interval \( [0,T] \), with \( T>0 \), is

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

Substituting the output expression gives

\[ \begin{aligned} E_y(0,T;\mathbf{x}_0) &= \int_0^T \mathbf{x}_0^T e^{\mathbf{A}^Tt}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}t}\mathbf{x}_0\,dt \\[2mm] &= \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. \end{aligned} \]

Therefore, the finite-horizon observability Gramian is

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

The quadratic form \( \mathbf{x}_0^T\mathbf{W}_o(0,T)\mathbf{x}_0 \) equals the output energy produced by the initial state.

The dual system associated with the sensing pair \( (\mathbf{A},\mathbf{C}) \) is constructed by transposing the state matrix and treating \( \mathbf{C}^T \) as an input matrix:

\[ \dot{\mathbf{z} }(t)=\mathbf{A}^T\mathbf{z}(t)+\mathbf{C}^T\mathbf{v}(t). \]

For a generic system \( \dot{\mathbf{x} }=\mathbf{F}\mathbf{x}+\mathbf{G}\mathbf{u} \), the finite-horizon controllability Gramian is

\[ \mathbf{W}_c^{(\mathbf{F},\mathbf{G})}(0,T)= \int_0^T e^{\mathbf{F}t}\mathbf{G}\mathbf{G}^T e^{\mathbf{F}^Tt}\,dt. \]

Setting \( \mathbf{F}=\mathbf{A}^T \) and \( \mathbf{G}=\mathbf{C}^T \) gives

\[ \begin{aligned} \mathbf{W}_c^{dual}(0,T) &= \int_0^T e^{\mathbf{A}^Tt}\mathbf{C}^T \left(\mathbf{C}^T\right)^T e^{(\mathbf{A}^T)^Tt}\,dt \\[2mm] &= \int_0^T e^{\mathbf{A}^Tt}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}t}\,dt. \end{aligned} \]

Hence,

\[ \boxed{ \mathbf{W}_o^{(\mathbf{A},\mathbf{C})}(0,T)= \mathbf{W}_c^{(\mathbf{A}^T,\mathbf{C}^T)}(0,T) }. \]

This identity is exact. It is not an approximation and does not require diagonalizability of \( \mathbf{A} \).

From the Gramian rank test, the pair \( (\mathbf{A},\mathbf{C}) \) is observable on \( [0,T] \) if and only if

\[ \operatorname{rank}\mathbf{W}_o(0,T)=n. \]

Since \( \mathbf{W}_o(0,T) \) is exactly the controllability Gramian of \( (\mathbf{A}^T,\mathbf{C}^T) \), we also have

\[ \boxed{ (\mathbf{A},\mathbf{C}) \text{ observable} \iff (\mathbf{A}^T,\mathbf{C}^T) \text{ controllable}. } \]

In terms of Kalman matrices,

\[ \mathcal{O}= \begin{bmatrix} \mathbf{C}\\ \mathbf{C}\mathbf{A}\\ \vdots\\ \mathbf{C}\mathbf{A}^{n-1} \end{bmatrix}, \qquad \mathcal{C}_{dual}= \begin{bmatrix} \mathbf{C}^T & \mathbf{A}^T\mathbf{C}^T & \cdots & (\mathbf{A}^T)^{n-1}\mathbf{C}^T \end{bmatrix}. \]

Since \( \mathcal{C}_{dual}=\mathcal{O}^T \), their ranks are equal:

\[ \operatorname{rank}\mathcal{O} = \operatorname{rank}\mathcal{O}^T = \operatorname{rank}\mathcal{C}_{dual}. \]

Thus, Kalman rank duality and Gramian duality are two equivalent views of the same underlying linear-algebraic symmetry.

Because \( \mathbf{W}_o(0,T) \) is symmetric positive semidefinite, loss of observability is equivalent to existence of a nonzero vector in its null space. Let \( \boldsymbol{\xi}\neq \mathbf{0} \). Then

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

Since the integrand is nonnegative,

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

Expanding the matrix exponential,

\[ \mathbf{C}e^{\mathbf{A}t}\boldsymbol{\xi} = \mathbf{C} \left( \mathbf{I}+\mathbf{A}t+ \frac{\mathbf{A}^2t^2}{2!}+\cdots \right) \boldsymbol{\xi}. \]

If this expression vanishes for every \( t \) on an interval, then all Taylor coefficients vanish:

\[ \mathbf{C}\boldsymbol{\xi}=\mathbf{0},\quad \mathbf{C}\mathbf{A}\boldsymbol{\xi}=\mathbf{0},\quad \ldots,\quad \mathbf{C}\mathbf{A}^{n-1}\boldsymbol{\xi}=\mathbf{0}. \]

Therefore,

\[ \boldsymbol{\xi}\in\ker\mathbf{W}_o(0,T) \iff \boldsymbol{\xi}\in\ker\mathcal{O}. \]

Consequently,

\[ \boxed{ \mathbf{W}_o(0,T) \text{ positive definite} \iff \mathcal{O} \text{ has full column rank} }. \]

If \( \mathbf{A} \) is Hurwitz, meaning all eigenvalues have negative real parts, the infinite-horizon observability Gramian is finite:

\[ \mathbf{W}_o= \int_0^\infty e^{\mathbf{A}^Tt}\mathbf{C}^T\mathbf{C}e^{\mathbf{A}t}\,dt. \]

This matrix satisfies the continuous-time Lyapunov equation

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

The controllability Gramian of the dual system \( \dot{\mathbf{z} }=\mathbf{A}^T\mathbf{z}+\mathbf{C}^T\mathbf{v} \) satisfies

\[ \mathbf{A}^T\mathbf{W}_c^{dual} + \mathbf{W}_c^{dual}\mathbf{A} + \mathbf{C}^T\mathbf{C} = \mathbf{0}. \]

The Lyapunov equation is the same. Under the Hurwitz assumption, its symmetric solution is unique. Hence

\[ \boxed{ \mathbf{W}_o=\mathbf{W}_c^{dual} }. \]

This infinite-horizon identity is especially important for numerical computation because standard controllability-Gramian solvers can compute observability Gramians by applying them to the dual pair.

The controllability Gramian measures input energy required to reach a state, while the observability Gramian measures output energy produced by an initial state. Duality says these two energy pictures are transposed versions of each other.

\[ \underbrace{ \mathbf{x}_0^T\mathbf{W}_o\mathbf{x}_0 }_{\text{output energy from initial state} } = \underbrace{ \mathbf{x}_0^T\mathbf{W}_c^{dual}\mathbf{x}_0 }_{\text{dual reachability strength} } . \]

A direction \( \mathbf{v} \) with large \( \mathbf{v}^T\mathbf{W}_o\mathbf{v} \) is strongly visible in the output. The same direction is easily reachable in the dual system. A direction with small \( \mathbf{v}^T\mathbf{W}_o\mathbf{v} \) is weakly visible and corresponds to a difficult-to-reach direction in the dual system.

flowchart TD
  A["Original initial state direction v"] --> B["Output energy: vT Wo v"]
  B --> C["Large value: strongly observable"]
  B --> D["Small value: weakly observable"]
  C --> E["Dual system: easy reachable direction"]
  D --> F["Dual system: difficult reachable direction"]
        

Therefore, eigenvectors of \( \mathbf{W}_o \) identify strongly and weakly sensed state directions:

\[ \mathbf{W}_o\mathbf{q}_i=\lambda_i\mathbf{q}_i, \qquad \lambda_i\ge 0. \]

If \( \lambda_i \) is small, the initial condition in direction \( \mathbf{q}_i \) generates little output energy and is numerically difficult to reconstruct from measurements.

The identity \( \mathbf{W}_o^{(\mathbf{A},\mathbf{C})} = \mathbf{W}_c^{(\mathbf{A}^T,\mathbf{C}^T)} \) gives a practical recipe:

\[ \text{Compute observability Gramian} = \text{compute controllability Gramian of the dual pair}. \]

Common libraries and tools include:

• Python: numpy, scipy.linalg, and python-control. The function scipy.linalg.solve_continuous_lyapunov is useful for infinite-horizon Gramians.
• C++: Eigen, Armadillo, Boost.Odeint, and SLICOT-based wrappers. For teaching, a from-scratch matrix-exponential quadrature is also useful.
• Java: EJML and Apache Commons Math can be used for matrix operations. The example below implements finite-horizon integration directly.
• MATLAB/Simulink: Control System Toolbox provides gram, lyap, and ss. Simulink can represent the dual system using a State-Space block.
• Wolfram Mathematica: MatrixExp, NIntegrate, and LyapunovSolve allow symbolic or numerical verification.

Chapter15_Lesson5.py

# Chapter15_Lesson5.py
"""
Duality between the observability Gramian of (A, C)
and the controllability Gramian of the dual pair (A.T, C.T).

Requires:
    numpy
    scipy

Run:
    python Chapter15_Lesson5.py
"""

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


def finite_horizon_observability_gramian(A: np.ndarray, C: np.ndarray, T: float, steps: int = 4000) -> np.ndarray:
    """Compute W_o(0,T) = integral_0^T exp(A.T t) C.T C exp(A t) dt by trapezoidal quadrature."""
    W = np.zeros((A.shape[0], A.shape[0]))
    dt = T / steps
    Q = C.T @ C
    for k in range(steps + 1):
        t = k * dt
        E = expm(A * t)
        integrand = E.T @ Q @ E
        weight = 0.5 if (k == 0 or k == steps) else 1.0
        W += weight * integrand
    return W * dt


def finite_horizon_controllability_gramian(A: np.ndarray, B: np.ndarray, T: float, steps: int = 4000) -> np.ndarray:
    """Compute W_c(0,T) = integral_0^T exp(A t) B B.T exp(A.T t) dt by trapezoidal quadrature."""
    W = np.zeros((A.shape[0], A.shape[0]))
    dt = T / steps
    R = B @ B.T
    for k in range(steps + 1):
        t = k * dt
        E = expm(A * t)
        integrand = E @ R @ E.T
        weight = 0.5 if (k == 0 or k == steps) else 1.0
        W += weight * integrand
    return W * dt


def infinite_horizon_observability_gramian(A: np.ndarray, C: np.ndarray) -> np.ndarray:
    """
    For Hurwitz A, solve A.T W + W A + C.T C = 0.
    SciPy uses solve_continuous_lyapunov(a, q): a X + X a.T = q.
    Hence use a=A.T and q=-(C.T C).
    """
    return solve_continuous_lyapunov(A.T, -(C.T @ C))


def infinite_horizon_controllability_gramian(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """For Hurwitz A, solve A W + W A.T + B B.T = 0."""
    return solve_continuous_lyapunov(A, -(B @ B.T))


def main() -> None:
    A = np.array([[-1.0, 2.0],
                  [-3.0, -4.0]])
    C = np.array([[1.0, 0.5]])
    T = 3.0

    print("Eigenvalues of A:", eigvals(A))

    Wo_T = finite_horizon_observability_gramian(A, C, T)
    Wc_dual_T = finite_horizon_controllability_gramian(A.T, C.T, T)

    print("\nFinite-horizon W_o(A,C):")
    print(Wo_T)
    print("\nFinite-horizon W_c(A.T,C.T):")
    print(Wc_dual_T)
    print("\nFinite-horizon Frobenius difference:", np.linalg.norm(Wo_T - Wc_dual_T, ord="fro"))

    Wo_inf = infinite_horizon_observability_gramian(A, C)
    Wc_dual_inf = infinite_horizon_controllability_gramian(A.T, C.T)

    print("\nInfinite-horizon W_o(A,C):")
    print(Wo_inf)
    print("\nInfinite-horizon W_c(A.T,C.T):")
    print(Wc_dual_inf)
    print("\nInfinite-horizon Frobenius difference:", np.linalg.norm(Wo_inf - Wc_dual_inf, ord="fro"))

    eig_Wo = np.linalg.eigvalsh(Wo_inf)
    print("\nEigenvalues of W_o:", eig_Wo)
    print("Observable by Gramian test?", np.min(eig_Wo) > 1e-9)


if __name__ == "__main__":
    main()
      

Chapter15_Lesson5.cpp

// Chapter15_Lesson5.cpp
// From-scratch finite-horizon Gramian integration for a 2x2 example.
// Compile:
//     g++ -std=c++17 Chapter15_Lesson5.cpp -O2 -o Chapter15_Lesson5
// Run:
//     ./Chapter15_Lesson5

#include <array>
#include <cmath>
#include <iostream>
#include <iomanip>

using Mat2 = std::array<std::array<double, 2>, 2>;
using Vec2 = std::array<double, 2>;

Mat2 zeros() {
    return { { {0.0, 0.0}, {0.0, 0.0} } };
}

Mat2 add(const Mat2& A, const Mat2& B) {
    Mat2 C = zeros();
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            C[i][j] = A[i][j] + B[i][j];
    return C;
}

Mat2 scale(const Mat2& A, double s) {
    Mat2 C = zeros();
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            C[i][j] = s * A[i][j];
    return C;
}

Mat2 mul(const Mat2& A, const Mat2& B) {
    Mat2 C = zeros();
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            for (int k = 0; k < 2; ++k)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}

Mat2 trans(const Mat2& A) {
    return { { {A[0][0], A[1][0]}, {A[0][1], A[1][1]} } };
}

double frobenius_norm(const Mat2& A) {
    double s = 0.0;
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            s += A[i][j] * A[i][j];
    return std::sqrt(s);
}

Mat2 mat_exp_2x2(const Mat2& A, double t) {
    Mat2 At = scale(A, t);
    int squarings = 6;
    double factor = std::pow(2.0, squarings);
    Mat2 X = scale(At, 1.0 / factor);

    Mat2 E = { { {1.0, 0.0}, {0.0, 1.0} } };
    Mat2 term = E;
    for (int k = 1; k <= 24; ++k) {
        term = scale(mul(term, X), 1.0 / k);
        E = add(E, term);
    }

    for (int i = 0; i < squarings; ++i) {
        E = mul(E, E);
    }
    return E;
}

Mat2 observability_integrand(const Mat2& A, const Vec2& C, double t) {
    Mat2 E = mat_exp_2x2(A, t);
    Mat2 Q = { { {C[0] * C[0], C[0] * C[1]},
               {C[1] * C[0], C[1] * C[1]} } };
    return mul(mul(trans(E), Q), E);
}

Mat2 controllability_integrand(const Mat2& A, const Vec2& B, double t) {
    Mat2 E = mat_exp_2x2(A, t);
    Mat2 R = { { {B[0] * B[0], B[0] * B[1]},
               {B[1] * B[0], B[1] * B[1]} } };
    return mul(mul(E, R), trans(E));
}

Mat2 integrate_observability(const Mat2& A, const Vec2& C, double T, int steps) {
    Mat2 W = zeros();
    double dt = T / steps;
    for (int k = 0; k <= steps; ++k) {
        double weight = (k == 0 || k == steps) ? 0.5 : 1.0;
        W = add(W, scale(observability_integrand(A, C, k * dt), weight));
    }
    return scale(W, dt);
}

Mat2 integrate_controllability(const Mat2& A, const Vec2& B, double T, int steps) {
    Mat2 W = zeros();
    double dt = T / steps;
    for (int k = 0; k <= steps; ++k) {
        double weight = (k == 0 || k == steps) ? 0.5 : 1.0;
        W = add(W, scale(controllability_integrand(A, B, k * dt), weight));
    }
    return scale(W, dt);
}

void print_matrix(const Mat2& A) {
    std::cout << std::fixed << std::setprecision(8);
    for (const auto& row : A) {
        std::cout << "[ " << row[0] << "  " << row[1] << " ]\n";
    }
}

int main() {
    Mat2 A = { { {-1.0, 2.0}, {-3.0, -4.0} } };
    Mat2 AT = trans(A);
    Vec2 C = {1.0, 0.5};
    Vec2 Bdual = {1.0, 0.5};

    double T = 3.0;
    int steps = 4000;

    Mat2 Wo = integrate_observability(A, C, T, steps);
    Mat2 WcDual = integrate_controllability(AT, Bdual, T, steps);

    std::cout << "Finite-horizon observability Gramian W_o(A,C):\n";
    print_matrix(Wo);

    std::cout << "\nFinite-horizon controllability Gramian W_c(A^T,C^T):\n";
    print_matrix(WcDual);

    Mat2 D = add(Wo, scale(WcDual, -1.0));
    std::cout << "\nFrobenius difference: " << frobenius_norm(D) << "\n";

    return 0;
}
      

Chapter15_Lesson5.java

// Chapter15_Lesson5.java
// From-scratch finite-horizon Gramian integration for a 2x2 example.
// Compile:
//     javac Chapter15_Lesson5.java
// Run:
//     java Chapter15_Lesson5

public class Chapter15_Lesson5 {
    static double[][] zeros() {
        return new double[][] { {0.0, 0.0}, {0.0, 0.0} };
    }

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

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

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

    static double[][] trans(double[][] A) {
        return new double[][] { {A[0][0], A[1][0]}, {A[0][1], A[1][1]} };
    }

    static double frobeniusNorm(double[][] A) {
        double s = 0.0;
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                s += A[i][j] * A[i][j];
        return Math.sqrt(s);
    }

    static double[][] expm2x2(double[][] A, double t) {
        int squarings = 6;
        double factor = Math.pow(2.0, squarings);
        double[][] X = scale(A, t / factor);

        double[][] E = { {1.0, 0.0}, {0.0, 1.0} };
        double[][] term = { {1.0, 0.0}, {0.0, 1.0} };

        for (int k = 1; k <= 24; k++) {
            term = scale(mul(term, X), 1.0 / k);
            E = add(E, term);
        }

        for (int i = 0; i < squarings; i++) {
            E = mul(E, E);
        }
        return E;
    }

    static double[][] observabilityIntegrand(double[][] A, double[] C, double t) {
        double[][] E = expm2x2(A, t);
        double[][] Q = { {C[0] * C[0], C[0] * C[1]},
                        {C[1] * C[0], C[1] * C[1]} };
        return mul(mul(trans(E), Q), E);
    }

    static double[][] controllabilityIntegrand(double[][] A, double[] B, double t) {
        double[][] E = expm2x2(A, t);
        double[][] R = { {B[0] * B[0], B[0] * B[1]},
                        {B[1] * B[0], B[1] * B[1]} };
        return mul(mul(E, R), trans(E));
    }

    static double[][] integrateObservability(double[][] A, double[] C, double T, int steps) {
        double[][] W = zeros();
        double dt = T / steps;
        for (int k = 0; k <= steps; k++) {
            double weight = (k == 0 || k == steps) ? 0.5 : 1.0;
            W = add(W, scale(observabilityIntegrand(A, C, k * dt), weight));
        }
        return scale(W, dt);
    }

    static double[][] integrateControllability(double[][] A, double[] B, double T, int steps) {
        double[][] W = zeros();
        double dt = T / steps;
        for (int k = 0; k <= steps; k++) {
            double weight = (k == 0 || k == steps) ? 0.5 : 1.0;
            W = add(W, scale(controllabilityIntegrand(A, B, k * dt), weight));
        }
        return scale(W, dt);
    }

    static void printMatrix(double[][] A) {
        for (int i = 0; i < 2; i++) {
            System.out.printf("[ %.8f  %.8f ]%n", A[i][0], A[i][1]);
        }
    }

    public static void main(String[] args) {
        double[][] A = { {-1.0, 2.0}, {-3.0, -4.0} };
        double[][] AT = trans(A);
        double[] C = {1.0, 0.5};
        double[] Bdual = {1.0, 0.5};

        double T = 3.0;
        int steps = 4000;

        double[][] Wo = integrateObservability(A, C, T, steps);
        double[][] WcDual = integrateControllability(AT, Bdual, T, steps);

        System.out.println("Finite-horizon observability Gramian W_o(A,C):");
        printMatrix(Wo);

        System.out.println("\nFinite-horizon controllability Gramian W_c(A^T,C^T):");
        printMatrix(WcDual);

        double[][] D = add(Wo, scale(WcDual, -1.0));
        System.out.println("\nFrobenius difference: " + frobeniusNorm(D));
    }
}
      

Chapter15_Lesson5.m

% Chapter15_Lesson5.m
% Duality between observability and controllability Gramians.

clear; clc;

A = [-1  2;
     -3 -4];
C = [1 0.5];
T = 3.0;

N = 4000;
dt = T / N;
Wo = zeros(2,2);
WcDual = zeros(2,2);

for k = 0:N
    t = k * dt;
    E = expm(A * t);
    integrand_o = E' * C' * C * E;

    Edual = expm(A' * t);
    Bdual = C';
    integrand_c_dual = Edual * Bdual * Bdual' * Edual';

    weight = 1.0;
    if k == 0 || k == N
        weight = 0.5;
    end

    Wo = Wo + weight * integrand_o;
    WcDual = WcDual + weight * integrand_c_dual;
end

Wo = dt * Wo;
WcDual = dt * WcDual;

disp('Finite-horizon W_o(A,C):');
disp(Wo);

disp('Finite-horizon W_c(A'',C''):');
disp(WcDual);

disp('Finite-horizon Frobenius difference:');
disp(norm(Wo - WcDual, 'fro'));

WoInf = lyap(A', C' * C);
WcDualInf = lyap(A', C' * C);

disp('Infinite-horizon W_o(A,C):');
disp(WoInf);

disp('Infinite-horizon W_c(A'',C''):');
disp(WcDualInf);

disp('Infinite-horizon Frobenius difference:');
disp(norm(WoInf - WcDualInf, 'fro'));

try
    sys = ss(A, [0;0], C, 0);
    sysDual = ss(A', C', eye(2), zeros(2,1));
    WoToolbox = gram(sys, 'o');
    WcDualToolbox = gram(sysDual, 'c');
    disp('Toolbox W_o(A,C):');
    disp(WoToolbox);
    disp('Toolbox W_c(A'',C''):');
    disp(WcDualToolbox);
catch ME
    disp('Control System Toolbox verification skipped:');
    disp(ME.message);
end

try
    modelName = 'Chapter15_Lesson5_Duality_Model';
    new_system(modelName);
    open_system(modelName);
    add_block('simulink/Continuous/State-Space', [modelName '/Dual_State_Space']);
    set_param([modelName '/Dual_State_Space'], ...
        'A', 'A''', ...
        'B', 'C''', ...
        'C', 'eye(2)', ...
        'D', 'zeros(2,1)');
    save_system(modelName);
    disp(['Created Simulink model: ' modelName '.slx']);
catch ME
    disp('Simulink model creation skipped:');
    disp(ME.message);
end
      

Chapter15_Lesson5.nb

(* Chapter15_Lesson5.nb *)

A = { {-1, 2}, {-3, -4} };
Cmat = { {1, 0.5} };
T = 3;

Wo = NIntegrate[
   Transpose[MatrixExp[A t]].Transpose[Cmat].Cmat.MatrixExp[A t],
   {t, 0, T}
];

WcDual = NIntegrate[
   MatrixExp[Transpose[A] t].Transpose[Cmat].Cmat.
    Transpose[MatrixExp[Transpose[A] t]],
   {t, 0, T}
];

MatrixForm[Wo]
MatrixForm[WcDual]
Norm[Wo - WcDual]

WoInf = LyapunovSolve[Transpose[A], -Transpose[Cmat].Cmat];
WcDualInf = LyapunovSolve[Transpose[A], -Transpose[Cmat].Cmat];

MatrixForm[WoInf]
Norm[WoInf - WcDualInf]
      

Problem 1: Let \( \mathbf{W}_o(0,T) \) be the finite-horizon observability Gramian of \( (\mathbf{A},\mathbf{C}) \). Prove that it is symmetric positive semidefinite.

Solution:

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

The integrand is symmetric because

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

For any vector \( \boldsymbol{\xi} \),

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

Therefore, \( \mathbf{W}_o(0,T) \) is symmetric positive semidefinite.

Problem 2: Prove that \( \mathbf{W}_o^{(\mathbf{A},\mathbf{C})}(0,T) = \mathbf{W}_c^{(\mathbf{A}^T,\mathbf{C}^T)}(0,T) \).

Solution:

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

Hence the observability Gramian of the original pair equals the controllability Gramian of the dual pair.

Problem 3: Show that \( (\mathbf{A},\mathbf{C}) \) is observable if and only if \( (\mathbf{A}^T,\mathbf{C}^T) \) is controllable.

Solution:

The observability matrix of \( (\mathbf{A},\mathbf{C}) \) is

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

The controllability matrix of the dual pair is

\[ \mathcal{C}_{dual}= \begin{bmatrix} \mathbf{C}^T & \mathbf{A}^T\mathbf{C}^T & \cdots & (\mathbf{A}^T)^{n-1}\mathbf{C}^T \end{bmatrix} = \mathcal{O}^T. \]

Since transposition preserves rank,

\[ \operatorname{rank}\mathcal{C}_{dual} = \operatorname{rank}\mathcal{O}. \]

Therefore, \( \mathcal{O} \) has full column rank if and only if \( \mathcal{C}_{dual} \) has full row rank, proving the claim.

Problem 4: Consider

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

Determine whether the system is observable using the observability matrix and interpret the result using the Gramian.

Solution:

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

The rank is \( 1 \), not \( 2 \), so the system is not observable. The second state direction \( [0\;1]^T \) never appears in the output because the sensor measures only the first state. Therefore,

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

The zero eigenvalue of \( \mathbf{W}_o \) confirms the unobservable state direction.

Problem 5: Suppose \( \mathbf{A} \) is Hurwitz. Derive the Lyapunov equation for the infinite-horizon observability Gramian.

Solution:

Define

\[ \mathbf{W}_o= \int_0^\infty e^{\mathbf{A}^Tt}\mathbf{C}^T\mathbf{C}e^{\mathbf{A}t}\,dt. \]

Then

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

Since \( \mathbf{A} \) is Hurwitz, the limit term is zero. Thus,

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

or equivalently,

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

The observability Gramian of \( (\mathbf{A},\mathbf{C}) \) is exactly the controllability Gramian of the dual pair \( (\mathbf{A}^T,\mathbf{C}^T) \). This identity connects output-energy analysis to input-energy analysis, transfers controllability rank tests to observability rank tests, and explains why observability Lyapunov equations are transposed versions of controllability Lyapunov equations. Numerically, the duality is useful because one can reuse controllability-Gramian algorithms to compute observability Gramians.