In the previous lessons, controllability was first expressed through reachable directions and then through the energy needed to steer a state. The controllability Gramian converts those ideas into a matrix test. For the continuous-time LTI system \( \dot x(t)=Ax(t)+Bu(t) \), with \( x(t)\in\mathbb{R}^n \) and \( u(t)\in\mathbb{R}^m \), the finite-horizon controllability Gramian over \( [0,T] \) is

\[ W_c(T)=\int_0^T e^{A\tau}BB^T e^{A^T\tau}\,d\tau,\qquad T > 0. \]

The central test is simple: the pair \( (A,B) \) is controllable on \( [0,T] \) exactly when \( W_c(T) \) is nonsingular, equivalently positive definite. This is a stronger numerical object than the rank of the Kalman matrix because it also reveals which state directions require high or low input energy.

flowchart TD
  A["Given xdot = A x + B u"] --> B["Choose horizon T > 0"]
  B --> C["Compute Wc(T) by integral or ODE"]
  C --> D["Symmetrize Wc(T)"]
  D --> E["Check eigenvalues of Wc(T)"]
  E -->|all positive| F["Controllable on the horizon"]
  E -->|some zero or near zero| G["Unreachable or weakly reachable directions"]
  F --> H["Use Wc inverse for minimum energy"]
  G --> I["Inspect null space and actuator directions"]
        

The matrix \( W_c(T) \) is symmetric and positive semidefinite. For any vector \( v\in\mathbb{R}^n \),

\[ v^T W_c(T)v =\int_0^T v^T e^{A\tau}BB^T e^{A^T\tau}v\,d\tau =\int_0^T \left\|B^T e^{A^T\tau}v\right\|_2^2\,d\tau \ge 0. \]

Therefore, the Gramian fails to be positive definite precisely when there exists a nonzero direction \( v \) such that \( B^T e^{A^T\tau}v=0 \) for every \( \tau\in[0,T] \). Such a vector is orthogonal to every direction that the input can generate through the dynamics.

\[ W_c(T)>0 \quad \Longleftrightarrow \quad v^T W_c(T)v>0 \; \text{for every}\; v\neq 0. \]

Numerically, one usually computes the eigenvalues \( \lambda_1\le\lambda_2\le\cdots\le\lambda_n \) of the symmetric matrix \( W_c(T) \). A practical test is

\[ \lambda_{\min}(W_c(T)) > \varepsilon, \]

where \( \varepsilon \) is a tolerance chosen relative to problem scaling and floating-point precision. Exact mathematics uses \( \varepsilon=0 \); computation cannot.

Define the Kalman controllability matrix \( \mathcal{C}=[B\;AB\;A^2B\;\cdots\;A^{n-1}B] \). For an LTI system, the finite-horizon Gramian test and the Kalman rank test are equivalent.

\[ W_c(T)>0 \quad \Longleftrightarrow \quad \operatorname{rank}\mathcal{C}=n, \qquad T>0. \]

Proof. First suppose that \( W_c(T) \) is not positive definite. Then there is a nonzero vector \( v \) satisfying

\[ 0=v^T W_c(T)v =\int_0^T \left\|B^T e^{A^T\tau}v\right\|_2^2\,d\tau. \]

Since the integrand is continuous and nonnegative, it must vanish for all \( \tau\in[0,T] \):

\[ B^T e^{A^T\tau}v=0 \qquad \text{for all}\quad \tau\in[0,T]. \]

Expanding the matrix exponential in a power series gives

\[ B^T e^{A^T\tau}v =\sum_{k=0}^\infty \frac{\tau^k}{k!}B^T(A^T)^k v=0 \qquad \text{for all}\quad \tau\in[0,T]. \]

A real analytic function that is identically zero on an interval has all Taylor coefficients equal to zero. Thus,

\[ B^T(A^T)^k v=0 \quad \text{for every}\quad k=0,1,2,\ldots. \]

Equivalently, \( v^T A^k B=0 \) for all \( k \). In particular, \( v^T[B\;AB\;\cdots\;A^{n-1}B]=0 \), so the rows of \( \mathcal{C} \) cannot span \( \mathbb{R}^n \); hence \( \operatorname{rank}\mathcal{C}<n \).

Conversely, suppose \( \operatorname{rank}\mathcal{C}<n \). Then there is a nonzero vector \( v \) such that \( v^T A^kB=0 \) for \( k=0,1,\ldots,n-1 \). By the Cayley-Hamilton theorem, every power \( A^k \) with \( k\ge n \) is a linear combination of \( I,A,\ldots,A^{n-1} \). Therefore \( v^T A^kB=0 \) for all \( k\ge0 \), which implies \( B^T e^{A^T\tau}v=0 \) for all \( \tau \). Hence

\[ v^T W_c(T)v =\int_0^T \left\|B^T e^{A^T\tau}v\right\|_2^2\,d\tau=0, \]

so \( W_c(T) \) is not positive definite. This proves the equivalence.

For \( x(0)=0 \), the final state at time \( T \) is

\[ x(T)=\int_0^T e^{A(T-s)}Bu(s)\,ds. \]

Define the linear operator \( \mathcal{R}_T:L_2[0,T]\to\mathbb{R}^n \) by \( \mathcal{R}_Tu=x(T) \). Its adjoint is

\[ (\mathcal{R}_T^*z)(s)=B^T e^{A^T(T-s)}z. \]

Then the Gramian is exactly the operator product \( \mathcal{R}_T\mathcal{R}_T^* \):

\[ \mathcal{R}_T\mathcal{R}_T^*z =\int_0^T e^{A(T-s)}BB^T e^{A^T(T-s)}z\,ds =W_c(T)z. \]

Thus, \( W_c(T) \) being nonsingular means the reachability operator has full range. This operator viewpoint explains why the Gramian is not merely a test but also the natural matrix for minimum-energy steering.

When \( W_c(T)>0 \), the minimum-energy input that transfers \( x(0)=0 \) to \( x(T)=x_f \) is

\[ u^*(t)=B^T e^{A^T(T-t)}W_c(T)^{-1}x_f. \]

The minimum required energy is

\[ J_{\min}(x_f,T)=\int_0^T \|u^*(t)\|_2^2\,dt =x_f^T W_c(T)^{-1}x_f. \]

Therefore, a state direction associated with a small Gramian eigenvalue is controllable in exact arithmetic but difficult to reach in practice. If \( W_c(T)q_i=\lambda_i q_i \) and \( \|q_i\|_2=1 \), then

\[ J_{\min}(q_i,T)=\frac{1}{\lambda_i}. \]

flowchart LR
  A["Large eigenvalue of Wc"] --> B["Low energy direction"]
  C["Small eigenvalue of Wc"] --> D["High energy direction"]
  E["Zero eigenvalue of Wc"] --> F["Unreachable direction"]
  B --> G["Numerically strong controllability"]
  D --> H["Numerically weak controllability"]
  F --> I["Rank test fails"]
        

If \( A \) is Hurwitz, meaning every eigenvalue of \( A \) has negative real part, the infinite-horizon Gramian exists:

\[ W_c=\int_0^\infty e^{At}BB^T e^{A^Tt}\,dt. \]

Differentiating \( e^{At}BB^T e^{A^Tt} \) and integrating from \( 0 \) to \( \infty \) yields the continuous-time Lyapunov equation:

\[ AW_c+W_cA^T+BB^T=0. \]

For Hurwitz \( A \), this Lyapunov equation has a unique symmetric solution. The system is controllable if and only if this solution is positive definite:

\[ (A,B)\;\text{controllable} \quad \Longleftrightarrow \quad W_c>0. \]

The finite-horizon test does not require stability. The infinite-horizon test is computationally convenient but only valid as an integral test when \( A \) is Hurwitz.

For the discrete-time LTI system \( x_{k+1}=Ax_k+Bu_k \), the finite-horizon controllability Gramian over \( N \) steps is

\[ W_c[N]=\sum_{i=0}^{N-1}A^iBB^T(A^T)^i. \]

The corresponding test is

\[ W_c[N]>0 \quad \Longleftrightarrow \quad \operatorname{rank}[B\;AB\;\cdots\;A^{n-1}B]=n, \qquad N\ge n. \]

If \( A \) is Schur stable, the infinite-horizon discrete Gramian satisfies the Stein equation

\[ W_c=AW_cA^T+BB^T. \]

This discrete-time statement is the direct counterpart of the continuous-time Lyapunov equation.

Exact controllability is a rank property, but engineering computation is finite precision. A controllable system may have a Gramian whose smallest eigenvalue is extremely small. Then the input energy needed to reach the corresponding direction is extremely large, and the system behaves as if it were nearly uncontrollable.

\[ \kappa(W_c)=\frac{\lambda_{\max}(W_c)}{\lambda_{\min}(W_c)}. \]

A large condition number indicates anisotropic reachability: some state directions are easy to reach, while others are difficult. A careful Gramian-based workflow should report \( \operatorname{rank}(W_c) \), \( \lambda_{\min}(W_c) \), \( \det(W_c) \), and \( \kappa(W_c) \), rather than only a binary answer.

\[ \text{robust numerical pass} \quad \Longleftrightarrow \quad \lambda_{\min}(W_c)\;\text{is safely above the selected tolerance}. \]

Consider the second-order system

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

The Kalman matrix is

\[ \mathcal{C}=[B\;AB] =\begin{bmatrix}0&1\\1&-3\end{bmatrix},\qquad \det(\mathcal{C})=-1\neq0. \]

Hence the system is controllable. The finite-horizon Gramian over \( T=2 \) must be positive definite:

\[ W_c(2)=\int_0^2 e^{A\tau}BB^T e^{A^T\tau}\,d\tau>0. \]

By contrast, the system

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

has \( \mathcal{C}=[B\;AB]=\begin{bmatrix}1&0\\0&0\end{bmatrix} \), so the second state direction is not reachable. Its Gramian is singular for every finite horizon \( T>0 \).

The following file computes finite-horizon and infinite-horizon Gramians, tests eigenvalues and rank, and samples the minimum-energy control input.

Chapter12_Lesson3.py

# Chapter12_Lesson3.py
# Gramian-based controllability tests for continuous-time LTI systems.
# Requires: numpy, scipy

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


def finite_horizon_gramian(A: np.ndarray, B: np.ndarray, T: float) -> np.ndarray:
    """Compute Wc(T)=int_0^T exp(A tau) B B^T exp(A^T tau) d tau."""
    if T <= 0:
        raise ValueError("T must be positive.")

    A = np.asarray(A, dtype=float)
    B = np.asarray(B, dtype=float)
    BBt = B @ B.T

    def integrand(tau: float) -> np.ndarray:
        E = expm(A * tau)
        return E @ BBt @ E.T

    W, _ = quad_vec(integrand, 0.0, T, epsabs=1e-10, epsrel=1e-10)
    return 0.5 * (W + W.T)


def stable_infinite_horizon_gramian(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """For Hurwitz A, solve A W + W A^T + B B^T = 0."""
    A = np.asarray(A, dtype=float)
    B = np.asarray(B, dtype=float)
    # SciPy solves A X + X A^H = Q, so use Q = -B B^T for
    # A W + W A^T + B B^T = 0.
    W = solve_continuous_lyapunov(A, -(B @ B.T))
    return 0.5 * (W + W.T)


def gramian_test(W: np.ndarray, tol: float = 1e-9) -> dict:
    """Return eigenvalue, rank, determinant, and condition-number diagnostics."""
    lam = eigvalsh(0.5 * (W + W.T))
    rank = int(np.sum(lam > tol))
    det = float(np.linalg.det(W))
    cond = float(lam[-1] / lam[0]) if lam[0] > tol else np.inf
    return {
        "eigenvalues": lam,
        "rank_by_eigenvalues": rank,
        "rank_by_matrix_rank": int(matrix_rank(W, tol=tol)),
        "determinant": det,
        "condition_number": cond,
        "positive_definite": bool(lam[0] > tol),
    }


def minimum_energy_input(A: np.ndarray, B: np.ndarray, xf: np.ndarray, T: float, samples: int = 9):
    """
    Return the finite-horizon minimum-energy control:
    u*(t)=B^T exp(A^T (T-t)) Wc(T)^(-1) xf, for x(0)=0.
    """
    W = finite_horizon_gramian(A, B, T)
    eta = np.linalg.solve(W, xf)
    ts = np.linspace(0.0, T, samples)
    us = []
    for t in ts:
        u = B.T @ expm(A.T * (T - t)) @ eta
        us.append(u)
    return ts, np.array(us), W


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

    T = 2.0
    W1 = finite_horizon_gramian(A1, B1, T)
    print("Controllable example: finite-horizon Wc(T)")
    print(W1)
    print(gramian_test(W1))

    Winf = stable_infinite_horizon_gramian(A1, B1)
    print("\nControllable example: infinite-horizon Wc")
    print(Winf)
    print(gramian_test(Winf))

    xf = np.array([1.0, 0.0])
    ts, us, W_for_energy = minimum_energy_input(A1, B1, xf, T)
    Emin = xf.T @ np.linalg.solve(W_for_energy, xf)
    print("\nMinimum energy to reach xf =", xf, "over T =", T, "is", Emin)
    print("Sampled u*(t):")
    for t, u in zip(ts, us):
        print(f"t={t:5.2f}, u={u.ravel()}")

    # Uncontrollable example: second state is never actuated or dynamically reached.
    A2 = np.array([[0.0, 0.0],
                   [0.0, -1.0]])
    B2 = np.array([[1.0],
                   [0.0]])
    W2 = finite_horizon_gramian(A2, B2, T)
    print("\nUncontrollable example: finite-horizon Wc(T)")
    print(W2)
    print(gramian_test(W2))
      

The C++ implementation below uses a Taylor approximation for the matrix exponential and Simpson integration for the Gramian. It is intentionally written from scratch to expose the numerical structure.

Chapter12_Lesson3.cpp

// Chapter12_Lesson3.cpp
// Finite-horizon Gramian test for a 2-state continuous-time LTI system.
// Standard C++17 only: no external linear algebra library.

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

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

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

Mat2 identity2() {
    return { { {1.0, 0.0}, {0.0, 1.0} } };
}

Mat2 add(const Mat2& X, const Mat2& Y) {
    Mat2 Z = zero2();
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            Z[i][j] = X[i][j] + Y[i][j];
    return Z;
}

Mat2 scale(const Mat2& X, double a) {
    Mat2 Z = zero2();
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            Z[i][j] = a * X[i][j];
    return Z;
}

Mat2 multiply(const Mat2& X, const Mat2& Y) {
    Mat2 Z = zero2();
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            for (int k = 0; k < 2; ++k)
                Z[i][j] += X[i][k] * Y[k][j];
    return Z;
}

Mat2 transpose(const Mat2& X) {
    return { { {X[0][0], X[1][0]}, {X[0][1], X[1][1]} } };
}

Mat2 outer(const Vec2& b) {
    return { { {b[0] * b[0], b[0] * b[1]},
             {b[1] * b[0], b[1] * b[1]} } };
}

Mat2 expm_taylor(const Mat2& A, double t, int terms = 40) {
    Mat2 At = scale(A, t);
    Mat2 E = identity2();
    Mat2 P = identity2();
    double factorial = 1.0;

    for (int k = 1; k <= terms; ++k) {
        P = multiply(P, At);
        factorial *= static_cast<double>(k);
        E = add(E, scale(P, 1.0 / factorial));
    }
    return E;
}

Mat2 integrand(const Mat2& A, const Vec2& b, double tau) {
    Mat2 E = expm_taylor(A, tau);
    Mat2 BBt = outer(b);
    return multiply(multiply(E, BBt), transpose(E));
}

Mat2 finite_horizon_gramian(const Mat2& A, const Vec2& b, double T, int N = 2000) {
    if (N % 2 == 1) ++N; // Simpson rule requires even number of subintervals.
    double h = T / static_cast<double>(N);
    Mat2 W = zero2();

    for (int k = 0; k <= N; ++k) {
        double tau = k * h;
        double coeff = (k == 0 || k == N) ? 1.0 : ((k % 2 == 0) ? 2.0 : 4.0);
        W = add(W, scale(integrand(A, b, tau), coeff));
    }

    W = scale(W, h / 3.0);
    W[0][1] = W[1][0] = 0.5 * (W[0][1] + W[1][0]);
    return W;
}

double det2(const Mat2& W) {
    return W[0][0] * W[1][1] - W[0][1] * W[1][0];
}

void symmetric_eigenvalues(const Mat2& W, double& lambda_min, double& lambda_max) {
    double a = W[0][0];
    double d = W[1][1];
    double c = 0.5 * (W[0][1] + W[1][0]);
    double tr = a + d;
    double disc = std::sqrt((a - d) * (a - d) + 4.0 * c * c);
    lambda_min = 0.5 * (tr - disc);
    lambda_max = 0.5 * (tr + disc);
}

void print_matrix(const Mat2& W) {
    std::cout << std::fixed << std::setprecision(10);
    std::cout << "[" << W[0][0] << ", " << W[0][1] << "]\n";
    std::cout << "[" << W[1][0] << ", " << W[1][1] << "]\n";
}

void report(const Mat2& W, double tol = 1e-8) {
    double lmin, lmax;
    symmetric_eigenvalues(W, lmin, lmax);
    std::cout << "det(W) = " << det2(W) << "\n";
    std::cout << "lambda_min = " << lmin << "\n";
    std::cout << "lambda_max = " << lmax << "\n";
    std::cout << "positive definite? " << (lmin > tol ? "yes" : "no") << "\n";
    if (lmin > tol)
        std::cout << "condition number = " << lmax / lmin << "\n";
    else
        std::cout << "condition number = infinity\n";
}

int main() {
    Mat2 A1 = { { {0.0, 1.0}, {-2.0, -3.0} } };
    Vec2 b1 = {0.0, 1.0};
    double T = 2.0;

    std::cout << "Controllable example Wc(T):\n";
    Mat2 W1 = finite_horizon_gramian(A1, b1, T);
    print_matrix(W1);
    report(W1);

    Mat2 A2 = { { {0.0, 0.0}, {0.0, -1.0} } };
    Vec2 b2 = {1.0, 0.0};

    std::cout << "\nUncontrollable example Wc(T):\n";
    Mat2 W2 = finite_horizon_gramian(A2, b2, T);
    print_matrix(W2);
    report(W2);

    return 0;
}
      

The Java version mirrors the C++ implementation and is useful for students who want a fully transparent implementation without external numerical libraries.

Chapter12_Lesson3.java

// Chapter12_Lesson3.java
// Finite-horizon Gramian test for 2-state continuous-time LTI systems.
// Uses only the Java standard library.

public class Chapter12_Lesson3 {
    static double[][] zero2() {
        return new double[][]{ {0.0, 0.0}, {0.0, 0.0} };
    }

    static double[][] identity2() {
        return new double[][]{ {1.0, 0.0}, {0.0, 1.0} };
    }

    static double[][] add(double[][] X, double[][] Y) {
        double[][] Z = zero2();
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                Z[i][j] = X[i][j] + Y[i][j];
        return Z;
    }

    static double[][] scale(double[][] X, double a) {
        double[][] Z = zero2();
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                Z[i][j] = a * X[i][j];
        return Z;
    }

    static double[][] multiply(double[][] X, double[][] Y) {
        double[][] Z = zero2();
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                for (int k = 0; k < 2; k++)
                    Z[i][j] += X[i][k] * Y[k][j];
        return Z;
    }

    static double[][] transpose(double[][] X) {
        return new double[][]{ {X[0][0], X[1][0]}, {X[0][1], X[1][1]} };
    }

    static double[][] outer(double[] b) {
        return new double[][]{
            {b[0] * b[0], b[0] * b[1]},
            {b[1] * b[0], b[1] * b[1]}
        };
    }

    static double[][] expmTaylor(double[][] A, double t, int terms) {
        double[][] At = scale(A, t);
        double[][] E = identity2();
        double[][] P = identity2();
        double factorial = 1.0;

        for (int k = 1; k <= terms; k++) {
            P = multiply(P, At);
            factorial *= k;
            E = add(E, scale(P, 1.0 / factorial));
        }
        return E;
    }

    static double[][] integrand(double[][] A, double[] b, double tau) {
        double[][] E = expmTaylor(A, tau, 40);
        double[][] BBt = outer(b);
        return multiply(multiply(E, BBt), transpose(E));
    }

    static double[][] finiteHorizonGramian(double[][] A, double[] b, double T, int N) {
        if (N % 2 == 1) N++;
        double h = T / N;
        double[][] W = zero2();

        for (int k = 0; k <= N; k++) {
            double tau = k * h;
            double coeff = (k == 0 || k == N) ? 1.0 : ((k % 2 == 0) ? 2.0 : 4.0);
            W = add(W, scale(integrand(A, b, tau), coeff));
        }

        W = scale(W, h / 3.0);
        double off = 0.5 * (W[0][1] + W[1][0]);
        W[0][1] = off;
        W[1][0] = off;
        return W;
    }

    static double det2(double[][] W) {
        return W[0][0] * W[1][1] - W[0][1] * W[1][0];
    }

    static double[] symmetricEigenvalues(double[][] W) {
        double a = W[0][0];
        double d = W[1][1];
        double c = 0.5 * (W[0][1] + W[1][0]);
        double tr = a + d;
        double disc = Math.sqrt((a - d) * (a - d) + 4.0 * c * c);
        return new double[]{0.5 * (tr - disc), 0.5 * (tr + disc)};
    }

    static void printMatrix(double[][] W) {
        System.out.printf("[%.10f, %.10f]%n", W[0][0], W[0][1]);
        System.out.printf("[%.10f, %.10f]%n", W[1][0], W[1][1]);
    }

    static void report(double[][] W) {
        double[] ev = symmetricEigenvalues(W);
        double tol = 1e-8;
        System.out.println("det(W) = " + det2(W));
        System.out.println("lambda_min = " + ev[0]);
        System.out.println("lambda_max = " + ev[1]);
        System.out.println("positive definite? " + (ev[0] > tol ? "yes" : "no"));
        if (ev[0] > tol)
            System.out.println("condition number = " + (ev[1] / ev[0]));
        else
            System.out.println("condition number = infinity");
    }

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

        System.out.println("Controllable example Wc(T):");
        double[][] W1 = finiteHorizonGramian(A1, b1, T, 2000);
        printMatrix(W1);
        report(W1);

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

        System.out.println("\nUncontrollable example Wc(T):");
        double[][] W2 = finiteHorizonGramian(A2, b2, T, 2000);
        printMatrix(W2);
        report(W2);
    }
}
      

MATLAB provides direct control-analysis functions such as ctrb, rank, lyap, and gram. The script also includes a finite-horizon numerical integration method and a Simulink verification note.

Chapter12_Lesson3.m

% Chapter12_Lesson3.m
% Gramian-based controllability tests in MATLAB/Simulink style.
% Requires Control System Toolbox for ctrb and gram; numerical integration works without gram.

clear; clc;

A = [0 1; -2 -3];
B = [0; 1];
C = eye(2);
D = zeros(2,1);
T = 2.0;

% Classical Kalman rank test for comparison
Co = ctrb(A,B);
fprintf('Rank of Kalman controllability matrix = %d\n', rank(Co));

% Finite-horizon controllability Gramian:
% Wc(T) = int_0^T expm(A*tau)*B*B'*expm(A'*tau) d tau
N = 2000;
tau = linspace(0,T,N+1);
h = T/N;
W = zeros(size(A));
for k = 1:(N+1)
    E = expm(A*tau(k));
    F = E*B*B'*E';
    if k == 1 || k == N+1
        coeff = 1;
    elseif mod(k-1,2) == 0
        coeff = 2;
    else
        coeff = 4;
    end
    W = W + coeff*F;
end
W = (h/3)*W;
W = 0.5*(W + W');

disp('Finite-horizon Wc(T):');
disp(W);
disp('Eigenvalues of Wc(T):');
disp(eig(W));
fprintf('det(Wc(T)) = %.12g\n', det(W));
fprintf('rank(Wc(T)) = %d\n', rank(W,1e-9));

% Infinite-horizon Gramian exists if A is Hurwitz.
if all(real(eig(A)) < 0)
    Winf = lyap(A, B*B');   % solves A*W + W*A' + B*B' = 0
    disp('Infinite-horizon Wc:');
    disp(Winf);
    disp('Eigenvalues of infinite-horizon Wc:');
    disp(eig(Winf));
end

% Minimum-energy steering from x(0)=0 to xf over [0,T].
xf = [1; 0];
eta = W\xf;
Emin = xf'*(W\xf);
fprintf('Minimum energy to reach xf = [1;0] over T=%.2f is %.12g\n', T, Emin);

sample_t = linspace(0,T,9);
fprintf('Sampled minimum-energy u*(t):\n');
for i = 1:length(sample_t)
    t = sample_t(i);
    u = B'*expm(A'*(T-t))*eta;
    fprintf('t = %.3f, u = %.12g\n', t, u);
end

% Simulink note:
% 1) Place a State-Space block with A, B, C = eye(2), D = zeros(2,1).
% 2) Feed the computed u*(t) through a From Workspace block.
% 3) Verify that the final simulated state approaches xf at t = T.
      

Mathematica can compute the Gramian symbolically or numerically. The notebook expression below evaluates finite-horizon and infinite-horizon tests and computes the minimum-energy input samples.

Chapter12_Lesson3.nb

Notebook[{
 Cell["Chapter12_Lesson3.nb", "Title"],
 Cell["Gramian-Based Controllability Tests", "Section"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"A", "=", "{ {0,1},{-2,-3} }"}], ";",
   RowBox[{"B", "=", "{ {0},{1} }"}], ";",
   RowBox[{"T", "=", "2"}], ";"}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"Wc", "=",
    RowBox[{"NIntegrate", "[",
     RowBox[{
      RowBox[{
       RowBox[{"MatrixExp", "[", RowBox[{"A", " ", "tau"}], "]"}], ".",
       "B", ".", RowBox[{"Transpose", "[", "B", "]"}], ".",
       RowBox[{"Transpose", "[",
        RowBox[{"MatrixExp", "[", RowBox[{"A", " ", "tau"}], "]"}], "]"}]}],
      ",", RowBox[{"{", RowBox[{"tau", ",", "0", ",", "T"}], "}"}]}],
     "]"}]}], ";",
   RowBox[{"Wc", "=", RowBox[{"N", "[", RowBox[{
     RowBox[{"(", RowBox[{"Wc", "+", RowBox[{"Transpose", "[", "Wc", "]"}]}],
     ")"}], "/", "2"}], "]"}]}]}]], "Input"],
 Cell[BoxData[RowBox[{"Eigenvalues", "[", "Wc", "]"}]], "Input"],
 Cell[BoxData[RowBox[{"MatrixRank", "[", "Wc", "]"}]], "Input"],
 Cell[BoxData[RowBox[{"Det", "[", "Wc", "]"}]], "Input"],
 Cell["Infinite-horizon Gramian for Hurwitz A", "Section"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"vars", "=", RowBox[{"{", RowBox[{"w11", ",", "w12", ",", "w22"}], "}"}]}], ";",
   RowBox[{"W", "=", RowBox[{"{ {", RowBox[{"w11", ",", "w12"}], "},{", RowBox[{"w12", ",", "w22"}], "} }"}]}], ";",
   RowBox[{"eqs", "=", RowBox[{"Thread", "[",
    RowBox[{"Flatten", "[", RowBox[{
     RowBox[{"A", ".", "W"}], "+", RowBox[{"W", ".", RowBox[{"Transpose", "[", "A", "]"}]}], "+",
     RowBox[{"B", ".", RowBox[{"Transpose", "[", "B", "]"}]}]}], "]"}], "==", "0"}], "]"}]}], ";",
   RowBox[{"Winf", "=", RowBox[{"W", "/.", RowBox[{"First", "[", RowBox[{"Solve", "[", RowBox[{"eqs", ",", "vars"}], "]"}], "]"}]}]}]}]], "Input"],
 Cell[BoxData[RowBox[{"Eigenvalues", "[", "Winf", "]"}]], "Input"],
 Cell["Minimum-energy steering", "Section"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"xf", "=", RowBox[{"{", RowBox[{"1", ",", "0"}], "}"}]}], ";",
   RowBox[{"eta", "=", RowBox[{"LinearSolve", "[", RowBox[{"Wc", ",", "xf"}], "]"}]}], ";",
   RowBox[{"Emin", "=", RowBox[{"xf", ".", "eta"}]}]}]], "Input"],
 Cell[BoxData[
  RowBox[{"Table", "[",
   RowBox[{
    RowBox[{
     RowBox[{"Transpose", "[", "B", "]"}], ".",
     RowBox[{"MatrixExp", "[", RowBox[{
      RowBox[{"Transpose", "[", "A", "]"}], " ", RowBox[{"(", RowBox[{"T", "-", "t"}], ")"}]}], "]"}], ".",
     "eta"}],
    ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "T", ",", "T/8"}], "}"}]}],
   "]"}]], "Input"]
}]
      

Problem 1 (Positive Semidefiniteness): Prove that \( W_c(T) \) is symmetric positive semidefinite for every \( T>0 \).

Solution: Symmetry follows from \( (e^{A\tau}BB^T e^{A^T\tau})^T=e^{A\tau}BB^T e^{A^T\tau} \). For any \( v \),

\[ v^T W_c(T)v =\int_0^T \left\|B^T e^{A^T\tau}v\right\|_2^2\,d\tau\ge0. \]

Hence \( W_c(T) \) is positive semidefinite.

Problem 2 (Singular Gramian and Unreachable Direction): Suppose \( W_c(T) \) is singular. Show that there is a nonzero vector \( v \) orthogonal to every reachable state from the origin.

Solution: Since \( W_c(T) \) is symmetric positive semidefinite and singular, there exists \( v\neq0 \) such that \( W_c(T)v=0 \). Then

\[ 0=v^T W_c(T)v =\int_0^T \left\|B^T e^{A^T\tau}v\right\|_2^2\,d\tau. \]

Thus \( B^T e^{A^T\tau}v=0 \) for all \( \tau\in[0,T] \). For any input \( u \), the final state from the origin satisfies

\[ v^Tx(T)=\int_0^T v^Te^{A(T-s)}Bu(s)\,ds=0. \]

Therefore every reachable state lies in the hyperplane orthogonal to \( v \).

Problem 3 (Minimum Energy Along an Eigenvector): Let \( W_c(T)q_i=\lambda_i q_i \) with \( \|q_i\|_2=1 \) and \( \lambda_i>0 \). Compute the minimum energy needed to reach \( q_i \) from the origin.

Solution: By the minimum-energy formula,

\[ J_{\min}(q_i,T)=q_i^T W_c(T)^{-1}q_i. \]

Since \( W_c(T)^{-1}q_i=\lambda_i^{-1}q_i \),

\[ J_{\min}(q_i,T)=q_i^T\lambda_i^{-1}q_i =\frac{1}{\lambda_i}. \]

Problem 4 (Infinite-Horizon Lyapunov Equation): Assume \( A \) is Hurwitz. Derive \( AW_c+W_cA^T+BB^T=0 \) from the integral definition of the infinite-horizon controllability Gramian.

Solution: Start from \( W_c=\int_0^\infty e^{At}BB^T e^{A^Tt}dt \). Then

\[ AW_c+W_cA^T =\int_0^\infty \frac{d}{dt} \left(e^{At}BB^T e^{A^Tt}\right)dt. \]

Because \( A \) is Hurwitz, \( e^{At}BB^T e^{A^Tt} \) converges to zero as \( t\to\infty \). Therefore,

\[ AW_c+W_cA^T =0-BB^T=-BB^T, \]

which gives \( AW_c+W_cA^T+BB^T=0 \).

Problem 5 (Discrete-Time Test): For \( x_{k+1}=Ax_k+Bu_k \), prove that if \( W_c[n] \) is positive definite, then the Kalman controllability matrix has rank \( n \).

Solution: Suppose the Kalman matrix does not have full rank. Then there exists \( v\neq0 \) such that \( v^T A^iB=0 \) for \( i=0,1,\ldots,n-1 \). Hence

\[ v^T W_c[n]v =\sum_{i=0}^{n-1}v^TA^iBB^T(A^T)^iv =\sum_{i=0}^{n-1}\left\|B^T(A^T)^iv\right\|_2^2=0. \]

This contradicts positive definiteness of \( W_c[n] \). Therefore the Kalman matrix must have full rank.

Gramian-based controllability tests provide an energy-aware alternative to the Kalman rank test. The finite-horizon Gramian is positive definite exactly when the LTI system is controllable, and its eigenstructure quantifies the energetic difficulty of reaching state-space directions. For Hurwitz continuous-time systems, the infinite-horizon Gramian is obtained from a Lyapunov equation; for Schur stable discrete-time systems, it is obtained from a Stein equation. In computation, the smallest eigenvalue and condition number of the Gramian are often more informative than a binary rank declaration.