In the previous controllability lessons, we learned that an LTI system can be checked by the Kalman matrix \( \mathcal{C}=[B\; AB\; A^2B\; \cdots\; A^{n-1}B] \). That rank test answers a qualitative question: can every state direction be reached? The controllability Gramian refines this idea by assigning a matrix measure to the accumulated input influence over a finite time interval.

Consider the continuous-time LTI system \( \dot{x}(t)=Ax(t)+Bu(t) \), where \( x(t)\in\mathbb{R}^n \) and \( u(t)\in\mathbb{R}^m \). For a fixed horizon \( T \), the Gramian summarizes all vectors of the form \( e^{A\tau}Bv \) generated as the input direction \( v\in\mathbb{R}^m \) varies and \( \tau \) sweeps through the time interval.

flowchart TD
  A["Input directions: columns of B"] --> B["Dynamic propagation: exp(A tau) B"]
  B --> C["Accumulate all propagated directions over 0 to T"]
  C --> D["Integral matrix Wc(T)"]
  D --> E["Rank tells reachable dimension"]
  D --> F["Eigenvalues show easy and difficult directions"]
        

The key conceptual shift is from a finite list of algebraic directions \( B,AB,\ldots,A^{n-1}B \) to a time-weighted continuum of directions \( e^{A\tau}B \).

For zero initial condition \( x(0)=0 \), the forced response at time \( T \) is

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

It is often convenient to change variables \( \tau=T-s \). Then \( s=T-\tau \) and the same reachable terminal state can be written as

\[ x(T)=\int_0^T e^{A\tau}B\,u(T-\tau)\,d\tau . \]

Define the finite-horizon reachability operator \( \mathcal{R}_T \) by

\[ \mathcal{R}_T u = \int_0^T e^{A\tau}B\,u(T-\tau)\,d\tau . \]

The range of \( \mathcal{R}_T \) is the reachable subspace at time \( T \) from the origin. The Gramian is the finite-dimensional matrix representation of \( \mathcal{R}_T\mathcal{R}_T^* \), where \( \mathcal{R}_T^* \) is the Hilbert-space adjoint of the reachability operator.

The finite-horizon controllability Gramian of \( (A,B) \) on \( [0,T] \) is the matrix

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

For a general interval \( [t_0,t_1] \) with \( t_1 > t_0 \), the corresponding LTI Gramian is

\[ W_c(t_0,t_1)=\int_{t_0}^{t_1} e^{A(t_1-s)}BB^T e^{A^T(t_1-s)}\,ds . \]

The dimensions are important: \( A\in\mathbb{R}^{n\times n} \), \( B\in\mathbb{R}^{n\times m} \), and therefore \( W_c(T)\in\mathbb{R}^{n\times n} \). The matrix \( BB^T \) describes instantaneous input authority, while multiplication by \( e^{A\tau} \) propagates that authority through the natural dynamics.

The finite-horizon Gramian also satisfies the matrix differential equation

\[ \dot{W}_c(t)=AW_c(t)+W_c(t)A^T+BB^T, \qquad W_c(0)=0 . \]

This differential equation is frequently the most practical implementation route when direct quadrature of the matrix exponential is expensive.

The Gramian is symmetric because the integrand is symmetric:

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

To prove positive semidefiniteness, take an arbitrary vector \( z\in\mathbb{R}^n \). Then

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

Therefore, \( W_c(T)=W_c(T)^T \) and \( W_c(T)\succeq 0 \). If \( W_c(T)\succ 0 \), every nonzero direction \( z \) receives nonzero accumulated input influence over the interval.

The null space has a direct interpretation:

\[ z\in\ker W_c(T) \quad\Longleftrightarrow\quad B^T e^{A^T\tau}z=0 \quad \text{for all }\tau\in[0,T]. \]

Thus a vector in the null space is orthogonal to every dynamically propagated input direction. This links the Gramian to the reachable subspace without introducing any new controllability test beyond the rank ideas already studied.

For continuous-time LTI systems, the following statements are equivalent for any finite horizon \( T > 0 \):

\[ \operatorname{rank}\,[B\;AB\;\cdots\;A^{n-1}B]=n \quad\Longleftrightarrow\quad W_c(T)\succ 0 . \]

Proof sketch. If \( W_c(T) \) is singular, then there exists nonzero \( z \) such that \( z^T W_c(T)z=0 \). From the quadratic form identity, this implies

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

Since \( e^{A^T\tau} \) is analytic in \( \tau \), all derivatives at \( \tau=0 \) must also vanish:

\[ B^T(A^T)^k z=0, \qquad k=0,1,2,\ldots . \]

Equivalently, \( z^T A^kB=0 \) for \( k=0,1,2,\ldots \). By the Cayley-Hamilton theorem, it is sufficient to check \( k=0,1,\ldots,n-1 \). Hence \( z \) is orthogonal to the columns of \( [B\;AB\;\cdots\;A^{n-1}B] \), so the Kalman matrix does not have full row rank. The converse follows by reversing the argument.

flowchart TD
  A["Kalman matrix full rank"] --> B["No nonzero vector orthogonal to B, AB, ..., A^(n-1)B"]
  B --> C["No nonzero vector kills B^T exp(A^T tau) z"]
  C --> D["Quadratic form z^T Wc(T) z positive"]
  D --> E["Wc(T) is positive definite"]
        

Suppose we change state coordinates using \( x=P\bar{x} \), where \( P \) is nonsingular. Then

\[ \dot{\bar{x} } = \bar{A}\bar{x}+\bar{B}u, \qquad \bar{A}=P^{-1}AP, \qquad \bar{B}=P^{-1}B . \]

The transformed Gramian is

\[ \begin{aligned} \bar{W}_c(T) &= \int_0^T e^{\bar{A}\tau}\bar{B}\bar{B}^T e^{\bar{A}^T\tau}\,d\tau \\ &= P^{-1}W_c(T)P^{-T} . \end{aligned} \]

Thus the numerical entries of a Gramian depend on the chosen state coordinates, but the rank property does not:

\[ \operatorname{rank}\,\bar{W}_c(T) = \operatorname{rank}\,W_c(T). \]

This is why controllability is a coordinate-invariant system property, while Gramian eigenvalues and condition numbers must be interpreted with attention to state scaling.

If \( A \) is asymptotically stable, meaning \( \operatorname{Re}(\lambda_i(A)) < 0 \) for all eigenvalues, then the infinite-horizon controllability Gramian is defined by

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

Differentiate the finite-horizon identity and take the limit as \( T\to\infty \). Since \( e^{AT}BB^T e^{A^TT} \to 0 \) for Hurwitz \( A \), the infinite-horizon Gramian solves the continuous-time Lyapunov equation:

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

This equation is central in numerical control because solving a Lyapunov equation is usually more stable and efficient than explicitly computing an improper integral.

There are three common computational routes:

Direct quadrature: evaluate \( e^{A\tau}BB^T e^{A^T\tau} \) over the interval and integrate numerically. This is conceptually closest to the definition but can be expensive for large systems.

Gramian differential equation: integrate \( \dot{W}=AW+WA^T+BB^T \) with \( W(0)=0 \). This avoids repeated quadrature but the matrix state has \( n^2 \) entries.

Lyapunov equation: for stable \( A \), solve \( AW+WA^T+BB^T=0 \). This is the standard infinite-horizon approach.

In exact mathematics, controllability is determined by rank. In computation, small singular values of \( W_c(T) \) indicate directions that are nearly unreachable over the chosen horizon or are poorly scaled. Thus one often inspects \( \sigma_{\min}(W_c(T)) \), \( \sigma_{\max}(W_c(T)) \), and \( \kappa(W_c(T)) \) rather than relying only on a rank command.

Python implementations usually use NumPy for matrix operations and SciPy for matrix exponentials, quadrature, eigenvalue routines, and Lyapunov equations. For control-specific workflows, the python-control package also provides state-space and Gramian-related utilities.

# Chapter12_Lesson1.py
# Definition and numerical computation of the finite-horizon controllability Gramian.
#
# Required packages:
#   pip install numpy scipy
#
# This script demonstrates:
#   1. Wc(T) = integral_0^T exp(A tau) B B^T exp(A^T tau) d tau
#   2. the equivalent differential equation dW/dt = A W + W A^T + B B^T
#   3. symmetry and positive semidefiniteness of Wc(T)
#   4. rank-based Gramian controllability test on a finite horizon
#   5. infinite-horizon Gramian for asymptotically stable A via the Lyapunov equation

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


def controllability_matrix(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """Return C = [B, AB, ..., A^(n-1)B]."""
    n = A.shape[0]
    blocks = []
    Apow = np.eye(n)
    for _ in range(n):
        blocks.append(Apow @ B)
        Apow = Apow @ A
    return np.hstack(blocks)


def finite_horizon_gramian(A: np.ndarray, B: np.ndarray, T: float) -> np.ndarray:
    """
    Compute Wc(T) from the equivalent differential equation:
        dW/dt = A W + W A^T + B B^T,   W(0) = 0.

    This ODE is mathematically equivalent to the integral definition:
        Wc(T) = integral_0^T exp(A tau) B B^T exp(A^T tau) d tau.
    """
    if T <= 0:
        raise ValueError("T must be positive.")

    n = A.shape[0]
    BBt = B @ B.T

    def rhs(_t: float, w_flat: np.ndarray) -> np.ndarray:
        W = w_flat.reshape(n, n)
        dW = A @ W + W @ A.T + BBt
        return dW.reshape(-1)

    sol = solve_ivp(
        rhs,
        (0.0, T),
        np.zeros(n * n),
        method="DOP853",
        rtol=1e-10,
        atol=1e-12,
    )
    if not sol.success:
        raise RuntimeError(sol.message)

    W = sol.y[:, -1].reshape(n, n)
    return 0.5 * (W + W.T)  # remove numerical asymmetry


def infinite_horizon_gramian(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """
    For Hurwitz A, Wc_inf solves:
        A W + W A^T + B B^T = 0.
    SciPy uses A X + X A^T = Q, hence Q = -B B^T.
    """
    W = solve_continuous_lyapunov(A, -(B @ B.T))
    return 0.5 * (W + W.T)


def print_matrix(name: str, M: np.ndarray) -> None:
    print(f"\n{name} =")
    print(np.array2string(M, precision=6, suppress_small=True))


def main() -> None:
    # Example 1: controllable second-order system.
    A1 = np.array([[0.0, 1.0],
                   [-2.0, -3.0]])
    B1 = np.array([[0.0],
                   [1.0]])
    T = 2.0

    C1 = controllability_matrix(A1, B1)
    W1 = finite_horizon_gramian(A1, B1, T)
    Winf1 = infinite_horizon_gramian(A1, B1)

    print("Example 1: controllable system")
    print_matrix("Kalman controllability matrix", C1)
    print("rank(C) =", np.linalg.matrix_rank(C1))
    print_matrix(f"Wc({T})", W1)
    print("eigenvalues of Wc(T) =", eigvalsh(W1))
    print("rank(Wc(T)) =", np.linalg.matrix_rank(W1, tol=1e-9))
    print_matrix("Wc(infinity)", Winf1)

    # Example 2: uncontrollable system. Second state is not affected by input or A-propagation.
    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)
    C2 = controllability_matrix(A2, B2)

    print("\n\nExample 2: uncontrollable system")
    print_matrix("Kalman controllability matrix", C2)
    print("rank(C) =", np.linalg.matrix_rank(C2))
    print_matrix(f"Wc({T})", W2)
    print("eigenvalues of Wc(T) =", eigvalsh(W2))
    print("rank(Wc(T)) =", np.linalg.matrix_rank(W2, tol=1e-9))

    # A simple geometric interpretation:
    # larger eigenvalue means a direction accumulates more control influence.
    vals, vecs = np.linalg.eigh(W1)
    print("\nPrincipal Gramian axes for Example 1:")
    for i, val in enumerate(vals):
        print(f"axis {i+1}: eigenvalue={val:.6f}, direction={vecs[:, i]}")


if __name__ == "__main__":
    main()

      

In C++, large-scale control programs typically use numerical linear algebra libraries such as Eigen, Armadillo, BLAS/LAPACK, or SLICOT. The following file intentionally uses a compact from-scratch 2 by 2 implementation to make the Gramian differential equation transparent.

// Chapter12_Lesson1.cpp
// Finite-horizon controllability Gramian for a 2 by 2 example.
// Compile:
//   g++ -std=c++17 Chapter12_Lesson1.cpp -O2 -o Chapter12_Lesson1
//
// This from-scratch implementation avoids external dependencies and uses RK4
// to integrate the matrix differential equation:
//   dW/dt = A W + W A^T + B B^T,   W(0) = 0.
// For production work, consider Eigen, Armadillo, or the SLICOT routines.

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

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 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 mul(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 gramian_rhs(const Mat2& A, const Vec2& b, const Mat2& W) {
    Mat2 AT = transpose(A);
    Mat2 BBt = outer(b);
    return add(add(mul(A, W), mul(W, AT)), BBt);
}

Mat2 rk4_step(const Mat2& A, const Vec2& b, const Mat2& W, double h) {
    Mat2 k1 = gramian_rhs(A, b, W);
    Mat2 k2 = gramian_rhs(A, b, add(W, scale(k1, 0.5 * h)));
    Mat2 k3 = gramian_rhs(A, b, add(W, scale(k2, 0.5 * h)));
    Mat2 k4 = gramian_rhs(A, b, add(W, scale(k3, h)));

    Mat2 update = add(add(k1, scale(k2, 2.0)), add(scale(k3, 2.0), k4));
    return add(W, scale(update, h / 6.0));
}

Mat2 finite_horizon_gramian(const Mat2& A, const Vec2& b, double T, int steps) {
    Mat2 W = zero2();
    double h = T / static_cast<double>(steps);
    for (int k = 0; k < steps; ++k) {
        W = rk4_step(A, b, W, h);
    }

    // Symmetrize to reduce small numerical skew.
    W[0][1] = W[1][0] = 0.5 * (W[0][1] + W[1][0]);
    return W;
}

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

void print_matrix(const char* name, const Mat2& X) {
    std::cout << "\n" << name << " =\n";
    for (int i = 0; i < 2; ++i) {
        std::cout << "[ ";
        for (int j = 0; j < 2; ++j) {
            std::cout << std::setw(12) << std::setprecision(7) << std::fixed << X[i][j] << " ";
        }
        std::cout << "]\n";
    }
}

int main() {
    Mat2 A = { { {0.0, 1.0}, {-2.0, -3.0} } };
    Vec2 b = {0.0, 1.0};

    double T = 2.0;
    int steps = 20000;

    Mat2 W = finite_horizon_gramian(A, b, T, steps);
    print_matrix("Wc(T)", W);

    std::cout << "\ndet(Wc(T)) = " << std::setprecision(12) << det2(W) << "\n";
    if (det2(W) > 1e-10) {
        std::cout << "The 2 by 2 Gramian is nonsingular: the system is controllable on this horizon.\n";
    } else {
        std::cout << "The Gramian is singular or nearly singular: controllability is lost or ill-conditioned.\n";
    }

    return 0;
}

      

For Java control computation, libraries such as EJML, Apache Commons Math, and JBLAS can be used. The code below again uses a small from-scratch RK4 implementation so that students see the matrix differential equation directly.

// Chapter12_Lesson1.java
// Finite-horizon controllability Gramian for a 2 by 2 example.
// Compile and run:
//   javac Chapter12_Lesson1.java
//   java Chapter12_Lesson1
//
// The code integrates:
//   dW/dt = A W + W A^T + B B^T,   W(0) = 0
// using fourth-order Runge-Kutta.

public class Chapter12_Lesson1 {
    static double[][] zero2() {
        return new double[][]{ {0.0, 0.0}, {0.0, 0.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[][] mul(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[][] rhs(double[][] A, double[] b, double[][] W) {
        double[][] AT = transpose(A);
        return add(add(mul(A, W), mul(W, AT)), outer(b));
    }

    static double[][] rk4Step(double[][] A, double[] b, double[][] W, double h) {
        double[][] k1 = rhs(A, b, W);
        double[][] k2 = rhs(A, b, add(W, scale(k1, 0.5 * h)));
        double[][] k3 = rhs(A, b, add(W, scale(k2, 0.5 * h)));
        double[][] k4 = rhs(A, b, add(W, scale(k3, h)));

        double[][] update = add(add(k1, scale(k2, 2.0)), add(scale(k3, 2.0), k4));
        return add(W, scale(update, h / 6.0));
    }

    static double[][] finiteHorizonGramian(double[][] A, double[] b, double T, int steps) {
        double[][] W = zero2();
        double h = T / steps;
        for (int k = 0; k < steps; k++) {
            W = rk4Step(A, b, W, h);
        }
        double off = 0.5 * (W[0][1] + W[1][0]);
        W[0][1] = off;
        W[1][0] = off;
        return W;
    }

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

    static void printMatrix(String name, double[][] X) {
        System.out.println("\n" + name + " =");
        for (int i = 0; i < 2; i++) {
            System.out.printf("[ %12.7f %12.7f ]%n", X[i][0], X[i][1]);
        }
    }

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

        double T = 2.0;
        int steps = 20000;

        double[][] W = finiteHorizonGramian(A, b, T, steps);
        printMatrix("Wc(T)", W);

        double determinant = det2(W);
        System.out.printf("%ndet(Wc(T)) = %.12f%n", determinant);

        if (determinant > 1e-10) {
            System.out.println("The Gramian is nonsingular: the system is controllable on this horizon.");
        } else {
            System.out.println("The Gramian is singular or nearly singular.");
        }
    }
}

      

MATLAB is especially convenient for this lesson because the Control System Toolbox provides ctrb, ss, gram, and lyap. The script also shows how a simple Simulink state-space model can be created programmatically when Simulink is available.

% Chapter12_Lesson1.m
% Definition and computation of the controllability Gramian.
%
% This script demonstrates:
%   1. finite-horizon Gramian by numerical integration,
%   2. the differential equation dW/dt = A W + W A' + B B',
%   3. infinite-horizon Gramian by the Lyapunov equation,
%   4. optional Control System Toolbox commands,
%   5. optional Simulink model creation for the state equation.

clear; clc;

A = [0 1; -2 -3];
B = [0; 1];
T = 2.0;

fprintf('Kalman controllability matrix:\n');
Cmat = ctrb(A, B);
disp(Cmat);
fprintf('rank(C) = %d\n', rank(Cmat));

% Finite-horizon definition:
% Wc(T) = integral_0^T expm(A*tau)*B*B'*expm(A'*tau) dtau
integrand = @(tau) reshape(expm(A*tau) * B * B' * expm(A'*tau), [], 1);
wvec = integral(integrand, 0, T, 'ArrayValued', true, 'AbsTol', 1e-11, 'RelTol', 1e-11);
WcT = reshape(wvec, size(A));
WcT = 0.5 * (WcT + WcT');

fprintf('\nFinite-horizon Gramian Wc(T):\n');
disp(WcT);
fprintf('eig(Wc(T)):\n');
disp(eig(WcT));
fprintf('rank(Wc(T)) = %d\n', rank(WcT));

% Alternative: integrate the Gramian differential equation.
odefun = @(t, w) reshape(A*reshape(w, 2, 2) + reshape(w, 2, 2)*A' + B*B', [], 1);
[~, Wtraj] = ode45(odefun, [0 T], zeros(4,1));
Wode = reshape(Wtraj(end, :).', 2, 2);
Wode = 0.5 * (Wode + Wode');

fprintf('\nFinite-horizon Gramian from ODE integration:\n');
disp(Wode);

% Infinite-horizon Gramian for Hurwitz A:
% A*W + W*A' + B*B' = 0
Winf = lyap(A, B*B');
fprintf('\nInfinite-horizon Gramian Wc(infinity):\n');
disp(Winf);

% Optional Control System Toolbox command:
sys = ss(A, B, eye(2), zeros(2,1));
try
    Wtoolbox = gram(sys, 'c');
    fprintf('\nControl System Toolbox gram(sys, ''c'') result:\n');
    disp(Wtoolbox);
catch ME
    fprintf('\ngram(sys, ''c'') was not available: %s\n', ME.message);
end

% Optional Simulink construction:
% The Gramian is not a standard Simulink block; however, one may simulate
% xdot = A*x + B*u and estimate reachability effects from impulse-like inputs.
% This block creates a minimal state-space model if Simulink is installed.
try
    model = 'Chapter12_Lesson1_Simulink_StateSpace';
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    new_system(model);
    add_block('simulink/Sources/Step', [model '/Step']);
    add_block('simulink/Continuous/State-Space', [model '/State-Space']);
    add_block('simulink/Sinks/Scope', [model '/Scope']);

    set_param([model '/State-Space'], ...
        'A', mat2str(A), ...
        'B', mat2str(B), ...
        'C', mat2str(eye(2)), ...
        'D', mat2str(zeros(2,1)));

    add_line(model, 'Step/1', 'State-Space/1');
    add_line(model, 'State-Space/1', 'Scope/1');
    save_system(model);
    fprintf('\nCreated optional Simulink model: %s.slx\n', model);
catch ME
    fprintf('\nSimulink model creation skipped: %s\n', ME.message);
end

      

Mathematica is useful for symbolic and high-precision matrix computation. The downloadable notebook contains the same workflow as the plain Wolfram Language code shown here.

(* Chapter12_Lesson1.nb *)
(* Wolfram Language notebook-style code for controllability Gramians. *)

ClearAll["Global`*"];

A = { {0, 1}, {-2, -3} };
B = { {0}, {1} };
Tfinal = 2.0;

controllabilityMatrix[A_, B_] := Module[{n = Length[A]},
  ArrayFlatten[{Table[MatrixPower[A, k].B, {k, 0, n - 1}]}]
];

finiteHorizonGramian[A_, B_, T_] := Module[{n = Length[A], tau},
  NIntegrate[
    MatrixExp[A tau].B.Transpose[B].Transpose[MatrixExp[A tau]],
    {tau, 0, T},
    WorkingPrecision -> 30
  ]
];

Cmat = controllabilityMatrix[A, B];
MatrixRank[Cmat]

WcT = finiteHorizonGramian[A, B, Tfinal] // N;
WcT // MatrixForm
Eigenvalues[WcT] // N
MatrixRank[WcT]

(* Infinite-horizon Gramian by solving A.W + W.Transpose[A] + B.Transpose[B] == 0. *)
n = Length[A];
Wvars = Array[w, {n, n}];
eqns = Flatten[A.Wvars + Wvars.Transpose[A] + B.Transpose[B] == ConstantArray[0, {n, n}]];
sol = Solve[eqns, Flatten[Wvars]][[1]];
Winf = Wvars /. sol;
Winf // MatrixForm
Eigenvalues[Winf] // N

(* Uncontrollable example. *)
A2 = { {0, 0}, {0, -1} };
B2 = { {1}, {0} };
Cmat2 = controllabilityMatrix[A2, B2];
WcT2 = finiteHorizonGramian[A2, B2, Tfinal] // N;

MatrixRank[Cmat2]
WcT2 // MatrixForm
Eigenvalues[WcT2] // N
MatrixRank[WcT2]

      

Problem 1 (Symmetry and Positive Semidefiniteness): Let \( W_c(T)=\int_0^T e^{A\tau}BB^T e^{A^T\tau}d\tau \). Prove that \( W_c(T) \) is symmetric and positive semidefinite.

Solution: Symmetry follows from transposing the integrand. For positive semidefiniteness, take any \( z\in\mathbb{R}^n \):

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

Hence \( W_c(T)=W_c(T)^T \) and \( W_c(T)\succeq 0 \).

Problem 2 (Scalar Gramian): For the scalar system \( \dot{x}=ax+bu \), compute \( W_c(T) \).

Solution: The definition gives

\[ W_c(T)=\int_0^T e^{2a\tau}b^2d\tau . \]

Therefore,

\[ W_c(T)= \begin{cases} \dfrac{b^2}{2a}\left(e^{2aT}-1\right), & a\ne 0,\\ b^2T, & a=0. \end{cases} \]

The scalar system is controllable exactly when \( b\ne 0 \), and in that case \( W_c(T)>0 \) for every \( T>0 \).

Problem 3 (Double Integrator): For \( A=\begin{bmatrix}0&1\\0&0\end{bmatrix} \) and \( B=\begin{bmatrix}0\\1\end{bmatrix} \), compute \( W_c(T) \).

Solution: Since \( e^{A\tau}=I+A\tau =\begin{bmatrix}1&\tau\\0&1\end{bmatrix} \),

\[ e^{A\tau}B= \begin{bmatrix}\tau\\1\end{bmatrix} . \]

Hence

\[ W_c(T)= \int_0^T \begin{bmatrix} \tau^2 & \tau\\ \tau & 1 \end{bmatrix}d\tau = \begin{bmatrix} T^3/3 & T^2/2\\ T^2/2 & T \end{bmatrix} . \]

Its determinant is

\[ \det W_c(T)=\frac{T^4}{12}, \]

which is positive for every \( T>0 \). Therefore the double integrator is controllable on any positive finite horizon.

Problem 4 (Uncontrollable Example): Let \( A=\begin{bmatrix}0&0\\0&-1\end{bmatrix} \) and \( B=\begin{bmatrix}1\\0\end{bmatrix} \). Compute \( W_c(T) \) and determine whether it is nonsingular.

Solution: Since \( e^{A\tau}B=\begin{bmatrix}1\\0\end{bmatrix} \),

\[ W_c(T)=\int_0^T \begin{bmatrix}1&0\\0&0\end{bmatrix}d\tau = \begin{bmatrix}T&0\\0&0\end{bmatrix} . \]

This matrix has rank one, so it is singular. The second state direction is never affected by the input or by dynamic propagation from the first input direction.

Problem 5 (Lyapunov Equation Derivation): Assume \( A \) is Hurwitz and define \( W_c=\int_0^\infty e^{A\tau}BB^T e^{A^T\tau}d\tau \). Prove that \( AW_c+W_cA^T+BB^T=0 \).

Solution: Observe that

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

Integrating this identity from \( 0 \) to \( \infty \) and using \( e^{A\tau}BB^T e^{A^T\tau}\to 0 \) gives

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

which is the desired Lyapunov equation.

The controllability Gramian is the matrix integral that accumulates how input directions are propagated by the state transition matrix. It is symmetric, positive semidefinite, coordinate-dependent in its numerical entries, but coordinate-invariant in rank. For finite horizons, it characterizes the reachable subspace through its range and rank. For stable systems, its infinite-horizon version solves a continuous-time Lyapunov equation. These facts prepare the next lesson, where the Gramian will be used to quantify the control energy required to reach a specified terminal state.