Consider the continuous-time LTI system \( \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t) \) on the finite interval \( [t_0,t_f] \), where \( t_f > t_0 \). From earlier lessons, the forced solution is

\[ \mathbf{x}(t_f)=e^{\mathbf{A}(t_f-t_0)}\mathbf{x}(t_0)+ \int_{t_0}^{t_f} e^{\mathbf{A}(t_f-\tau)}\mathbf{B}\mathbf{u}(\tau)\,d\tau. \]

Define the required displacement from the uncontrolled terminal state as \( \boldsymbol{\delta} = \mathbf{x}_f - e^{\mathbf{A}(t_f-t_0)}\mathbf{x}_0 \). The controllability Gramian is

\[ \mathbf{W}_c(t_0,t_f)=\int_{t_0}^{t_f} e^{\mathbf{A}(t_f-\tau)}\mathbf{B}\mathbf{B}^T e^{\mathbf{A}^T(t_f-\tau)}\,d\tau. \]

If \( \mathbf{W}_c \) is nonsingular, every \( \boldsymbol{\delta}\in\mathbb{R}^n \) can be produced. The stronger question in this lesson is not simply whether the state is reachable, but how much input energy it requires.

flowchart TD
  A["Target state xf"] --> B["Required displacement delta"]
  B --> C["Compute finite-horizon Gramian Wc"]
  C --> D["Eigen-decompose Wc"]
  D --> E["Large eigenvalue direction: low energy"]
  D --> F["Small eigenvalue direction: high energy"]
  E --> G["Easy states to reach"]
  F --> H["Difficult states to reach"]
        

Use the quadratic input-energy functional \( J(\mathbf{u})=\int_{t_0}^{t_f}\mathbf{u}^T(t)\mathbf{u}(t)\,dt \). Among all controls that steer \( \mathbf{x}_0 \) to \( \mathbf{x}_f \), the unique minimum-energy control is

\[ \boxed{\mathbf{u}^*(t)=\mathbf{B}^T e^{\mathbf{A}^T(t_f-t)} \mathbf{W}_c^{-1}(t_0,t_f)\boldsymbol{\delta} }. \]

The corresponding minimum energy is the Gramian inverse quadratic form

\[ \boxed{J_{\min}(\boldsymbol{\delta})= \boldsymbol{\delta}^T\mathbf{W}_c^{-1}(t_0,t_f)\boldsymbol{\delta} }. \]

Proof. Define the linear input-to-displacement operator \( \mathcal{L}\mathbf{u}=\int_{t_0}^{t_f} e^{\mathbf{A}(t_f-\tau)}\mathbf{B}\mathbf{u}(\tau)d\tau \). Its adjoint is \( (\mathcal{L}^*\mathbf{p})(t)=\mathbf{B}^T e^{\mathbf{A}^T(t_f-t)}\mathbf{p} \). Since \( \mathcal{L}\mathcal{L}^*=\mathbf{W}_c \), the minimum-norm solution of \( \mathcal{L}\mathbf{u}=\boldsymbol{\delta} \) is \( \mathbf{u}^*=\mathcal{L}^*\mathbf{W}_c^{-1}\boldsymbol{\delta} \). Substituting this expression into \( J \) gives

\[ \begin{aligned} J(\mathbf{u}^*) &=\int_{t_0}^{t_f} \boldsymbol{\delta}^T\mathbf{W}_c^{-1} e^{\mathbf{A}(t_f-t)}\mathbf{B}\mathbf{B}^T e^{\mathbf{A}^T(t_f-t)}\mathbf{W}_c^{-1}\boldsymbol{\delta}\,dt \\ &=\boldsymbol{\delta}^T\mathbf{W}_c^{-1}\mathbf{W}_c \mathbf{W}_c^{-1}\boldsymbol{\delta} =\boldsymbol{\delta}^T\mathbf{W}_c^{-1}\boldsymbol{\delta}. \end{aligned} \]

Since \( \mathbf{W}_c \) is symmetric positive definite for a completely controllable finite-horizon system, it admits the orthonormal decomposition

\[ \mathbf{W}_c=\mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^T, \quad \boldsymbol{\Lambda}=\operatorname{diag}(\lambda_1,\ldots,\lambda_n), \quad 0<\lambda_1\le\cdots\le\lambda_n. \]

Write the displacement in Gramian principal coordinates as \( \boldsymbol{\delta}=\sum_{i=1}^n \alpha_i\mathbf{q}_i \). Then

\[ J_{\min}(\boldsymbol{\delta})= \sum_{i=1}^n \frac{\alpha_i^2}{\lambda_i}. \]

This formula is the precise mathematical statement behind the words easy and difficult. A unit displacement along \( \mathbf{q}_i \) costs \( 1/\lambda_i \). Hence large \( \lambda_i \) means a direction is reached cheaply, while small \( \lambda_i \) means that even a modest terminal displacement may require large input energy.

flowchart TD
  A["Input energy budget rho squared"] --> B["Reachable ellipsoid"]
  B --> C["Long axes: large Gramian eigenvalues"]
  B --> D["Short axes: small Gramian eigenvalues"]
  C --> E["Low energy directions"]
  D --> F["High energy directions"]
  F --> G["Weak actuator authority, fast decay, or short horizon"]
        

For a fixed energy budget \( \rho^2 \), the set of reachable displacements is the ellipsoid

\[ \begin{aligned} \mathcal{E}_\rho &= \left\{\boldsymbol{\delta}: \right. \\ &\quad \left. \boldsymbol{\delta}^{T}\mathbf{W}_c^{-1}\boldsymbol{\delta} \le \rho^{2} \right\}. \end{aligned} \]

In principal coordinates \( \boldsymbol{\delta}=\mathbf{Q}\mathbf{z} \), this becomes

\[ \sum_{i=1}^n \frac{z_i^2}{\rho^2\lambda_i}\le 1. \]

Therefore the semi-axis length in direction \( \mathbf{q}_i \) is \( \rho\sqrt{\lambda_i} \). Long axes correspond to states that the actuators can create with small energy; short axes correspond to states that are barely affected by the available input over the interval.

A useful numerical summary is the Gramian condition number \( \kappa(\mathbf{W}_c)=\lambda_{\max}/\lambda_{\min} \). A large condition number does not mean the system is mathematically uncontrollable, but it means controllability is energetically ill conditioned: one direction may require orders of magnitude more energy than another.

Suppose \( \mathbf{A} \) is diagonalizable and use modal coordinates \( \mathbf{x}=\mathbf{V}\mathbf{z} \), so \( \dot{\mathbf{z} }=\boldsymbol{\Lambda}_A\mathbf{z}+\tilde{\mathbf{B} }\mathbf{u} \), where \( \tilde{\mathbf{B} }=\mathbf{V}^{-1}\mathbf{B} \). The modal Gramian is

\[ \tilde{\mathbf{W} }_c=\int_0^T e^{\boldsymbol{\Lambda}_A s}\tilde{\mathbf{B} } \tilde{\mathbf{B} }^T e^{\boldsymbol{\Lambda}_A^T s}ds. \]

For a diagonal stable example with \( \mathbf{A}=\operatorname{diag}(-a_1,\ldots,-a_n) \), \( a_i>0 \), and scalar input distribution \( \mathbf{b}=[b_1,\ldots,b_n]^T \), the entries are

\[ (\mathbf{W}_c)_{ij}=b_i b_j \frac{1-e^{-(a_i+a_j)T} }{a_i+a_j}. \]

This expression reveals three common sources of difficult states: weak modal input coefficient \( |b_i| \), fast modal decay \( a_i \), and short control horizon \( T \). A state direction may be controllable in the Kalman rank sense and still be extremely expensive if the corresponding Gramian eigenvalue is very small.

The Gramian depends on the state coordinates. If \( \mathbf{x}=\mathbf{T}\bar{\mathbf{x} } \), then the transformed Gramian is

\[ \bar{\mathbf{W} }_c=\mathbf{T}^{-1}\mathbf{W}_c\mathbf{T}^{-T}. \]

The energy for the same physical terminal state is invariant because \( \bar{\boldsymbol{\delta} }=\mathbf{T}^{-1}\boldsymbol{\delta} \) gives

\[ \bar{\boldsymbol{\delta} }^T\bar{\mathbf{W} }_c^{-1} \bar{\boldsymbol{\delta} }=\boldsymbol{\delta}^T\mathbf{W}_c^{-1} \boldsymbol{\delta}. \]

However, the visual impression of eigenvectors and ellipsoid axes can be distorted by arbitrary unit choices. For engineering interpretation, states should be scaled to comparable physical units or assessed using a chosen state-weighting matrix.

\[ J_{\min,Q}(\boldsymbol{\delta})= \boldsymbol{\delta}^T\mathbf{Q}_s^{1/2} \left(\mathbf{Q}_s^{1/2}\mathbf{W}_c\mathbf{Q}_s^{1/2}\right)^{-1} \mathbf{Q}_s^{1/2}\boldsymbol{\delta}. \]

Here \( \mathbf{Q}_s \) encodes the relative importance or normalization of state variables. Without such scaling, comparing a position measured in meters with an angular velocity measured in radians per second can be misleading.

Consider \( \mathbf{A}=\operatorname{diag}(-1,-4) \) and \( \mathbf{B}=[1,\varepsilon]^T \), with \( 0<\varepsilon\ll 1 \). For horizon \( T \), the Gramian is

\[ \mathbf{W}_c(T)= \begin{bmatrix} \dfrac{1-e^{-2T} }{2} & \varepsilon\dfrac{1-e^{-5T} }{5} \\ \varepsilon\dfrac{1-e^{-5T} }{5} & \varepsilon^2\dfrac{1-e^{-8T} }{8} \end{bmatrix}. \]

The second coordinate is difficult for two reasons: it has small direct input coefficient \( \varepsilon \), and its natural dynamics decay quickly. As \( \varepsilon\to0 \), the smallest Gramian eigenvalue goes to zero approximately on the order of \( \varepsilon^2 \), so the energy for unit displacement in the difficult direction grows on the order of \( 1/\varepsilon^2 \).

\[ J_{\min}(\mathbf{e}_2)=\mathbf{e}_2^T\mathbf{W}_c^{-1} \mathbf{e}_2=\mathcal{O}\left(\frac{1}{\varepsilon^2}\right). \]

In exact mathematics, controllability is binary: the Gramian is either nonsingular or singular. In computation and engineering design, it is more useful to inspect the spectrum \( \lambda_1,\ldots,\lambda_n \). A very small \( \lambda_{\min} \) means that finite actuator limits, saturation, sensor noise, or modeling error may make the state practically unreachable even when rank tests pass.

\[ \frac{J_{\min}(\mathbf{q}_{\min})}{J_{\min}(\mathbf{q}_{\max})} =\frac{\lambda_{\max} }{\lambda_{\min} }=\kappa(\mathbf{W}_c). \]

Thus the Gramian condition number measures the spread between the most difficult and easiest normalized principal directions. In actuator placement problems, one seeks input matrices \( \mathbf{B} \) that increase small Gramian eigenvalues, not merely those that make the Kalman matrix full rank.

Chapter12_Lesson4.py

# Chapter12_Lesson4.py
# Modern Control — Chapter 12, Lesson 4
# Insights into "Easy" vs "Difficult" States to Reach using controllability Gramians
# Libraries: NumPy, SciPy, Matplotlib

import numpy as np
from scipy.linalg import expm, eigh
import matplotlib.pyplot as plt


def finite_horizon_gramian(A, B, T, steps=4000):
    """Compute Wc(T)=int_0^T exp(A s) B B^T exp(A^T s) ds by trapezoidal quadrature."""
    n = A.shape[0]
    W = np.zeros((n, n), dtype=float)
    ds = T / steps
    for k in range(steps + 1):
        s = k * ds
        Phi = expm(A * s)
        integrand = Phi @ B @ B.T @ Phi.T
        weight = 0.5 if k in (0, steps) else 1.0
        W += weight * integrand * ds
    return 0.5 * (W + W.T)


def minimum_energy(W, delta):
    """Return E_min = delta^T W^{-1} delta."""
    return float(delta.T @ np.linalg.solve(W, delta))


def optimal_input(A, B, W, delta, T, t_grid):
    """u*(t)=B^T exp(A^T(T-t)) W^{-1} delta for steering x(0)=0 to delta."""
    alpha = np.linalg.solve(W, delta)
    values = []
    for t in t_grid:
        values.append((B.T @ expm(A.T * (T - t)) @ alpha).reshape(-1))
    return np.vstack(values)


def plot_energy_ellipse(W, rho=1.0, filename="Chapter12_Lesson4_energy_ellipse.png"):
    """Plot the reachable energy ellipsoid delta^T W^{-1} delta <= rho^2 for n=2."""
    eigvals, eigvecs = eigh(W)
    theta = np.linspace(0, 2 * np.pi, 400)
    circle = np.vstack([np.cos(theta), np.sin(theta)])
    ellipse = eigvecs @ np.diag(rho * np.sqrt(eigvals)) @ circle

    plt.figure(figsize=(6, 5))
    plt.plot(ellipse[0, :], ellipse[1, :], label="energy boundary")
    plt.axhline(0, linewidth=0.8)
    plt.axvline(0, linewidth=0.8)
    for i in range(2):
        axis = eigvecs[:, i] * rho * np.sqrt(eigvals[i])
        plt.plot([0, axis[0]], [0, axis[1]], linewidth=2, label=f"axis {i+1}")
    plt.gca().set_aspect("equal", adjustable="box")
    plt.xlabel("state component x1")
    plt.ylabel("state component x2")
    plt.title("Finite-horizon controllability energy ellipsoid")
    plt.legend()
    plt.tight_layout()
    plt.savefig(filename, dpi=160)
    print(f"Saved {filename}")


if __name__ == "__main__":
    # Example: the second state is directly actuated only weakly.
    A = np.array([[-1.0, 0.0],
                  [ 0.0, -4.0]])
    B = np.array([[1.0],
                  [0.08]])
    T = 2.0

    W = finite_horizon_gramian(A, B, T)
    eigvals, eigvecs = eigh(W)
    cond_W = eigvals[-1] / eigvals[0]

    print("Wc(T) =")
    print(W)
    print("Eigenvalues of Wc(T):", eigvals)
    print("Condition number:", cond_W)
    print("Eigenvectors (columns):")
    print(eigvecs)

    # Energy required to reach unit displacement along each Gramian eigenvector.
    for i in range(2):
        q = eigvecs[:, i]
        print(f"Energy to reach q_{i+1}: {minimum_energy(W, q):.6g}")
        print(f"Theoretical value 1/lambda_{i+1}: {1.0 / eigvals[i]:.6g}")

    # Energy for physical coordinate targets.
    e1 = np.array([1.0, 0.0])
    e2 = np.array([0.0, 1.0])
    print("Energy to reach e1:", minimum_energy(W, e1))
    print("Energy to reach e2:", minimum_energy(W, e2))

    # Verify optimal input energy by numerical integration.
    delta = e2
    t_grid = np.linspace(0.0, T, 2001)
    u = optimal_input(A, B, W, delta, T, t_grid)
    numerical_energy = np.trapezoid(np.sum(u * u, axis=1), t_grid)
    print("Numerical energy for target e2:", numerical_energy)
    print("Analytic minimum energy for target e2:", minimum_energy(W, delta))

    plot_energy_ellipse(W, rho=1.0)

      

Chapter12_Lesson4.cpp

// Chapter12_Lesson4.cpp
// Modern Control — Chapter 12, Lesson 4
// Finite-horizon controllability Gramian and energy directions
// Requires Eigen, including unsupported MatrixFunctions module.
// Example compile command:
// g++ -std=c++17 Chapter12_Lesson4.cpp -I /path/to/eigen -O2 -o Chapter12_Lesson4

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

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

Matrix finiteHorizonGramian(const Matrix& A, const Matrix& B, double T, int steps) {
    const int n = A.rows();
    Matrix W = Matrix::Zero(n, n);
    const double ds = T / static_cast<double>(steps);

    for (int k = 0; k <= steps; ++k) {
        const double s = k * ds;
        Matrix Phi = (A * s).exp();
        Matrix integrand = Phi * B * B.transpose() * Phi.transpose();
        double weight = (k == 0 || k == steps) ? 0.5 : 1.0;
        W += weight * integrand * ds;
    }
    return 0.5 * (W + W.transpose());
}

double minimumEnergy(const Matrix& W, const Vector& delta) {
    Vector alpha = W.ldlt().solve(delta);
    return delta.dot(alpha);
}

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

    Matrix B(2, 1);
    B << 1.0,
         0.08;

    const double T = 2.0;
    Matrix W = finiteHorizonGramian(A, B, T, 4000);

    Eigen::SelfAdjointEigenSolver<Matrix> solver(W);
    Vector lambda = solver.eigenvalues();
    Matrix Q = solver.eigenvectors();

    std::cout << std::setprecision(10);
    std::cout << "Wc(T) =\n" << W << "\n\n";
    std::cout << "Eigenvalues = " << lambda.transpose() << "\n";
    std::cout << "Condition number = " << lambda(lambda.size() - 1) / lambda(0) << "\n\n";
    std::cout << "Eigenvectors (columns) =\n" << Q << "\n\n";

    for (int i = 0; i < 2; ++i) {
        Vector q = Q.col(i);
        std::cout << "Energy to reach q_" << (i + 1) << " = "
                  << minimumEnergy(W, q) << " ; expected = " << 1.0 / lambda(i) << "\n";
    }

    Vector e1(2), e2(2);
    e1 << 1.0, 0.0;
    e2 << 0.0, 1.0;
    std::cout << "\nEnergy to reach e1 = " << minimumEnergy(W, e1) << "\n";
    std::cout << "Energy to reach e2 = " << minimumEnergy(W, e2) << "\n";
    return 0;
}

      

Chapter12_Lesson4.java

// Chapter12_Lesson4.java
// Modern Control — Chapter 12, Lesson 4
// From-scratch finite-horizon controllability Gramian for a 2-state SISO example.
// No external Java library is required.

public class Chapter12_Lesson4 {
    static double[][] A = { {-1.0, 0.0}, {0.0, -4.0} };
    static double[][] B = { {1.0}, {0.08} };

    static double[][] add(double[][] X, double[][] Y) {
        int r = X.length, c = X[0].length;
        double[][] Z = new double[r][c];
        for (int i = 0; i < r; i++) {
            for (int j = 0; j < c; j++) Z[i][j] = X[i][j] + Y[i][j];
        }
        return Z;
    }

    static double[][] scale(double[][] X, double a) {
        int r = X.length, c = X[0].length;
        double[][] Z = new double[r][c];
        for (int i = 0; i < r; i++) {
            for (int j = 0; j < c; j++) Z[i][j] = a * X[i][j];
        }
        return Z;
    }

    static double[][] mul(double[][] X, double[][] Y) {
        int r = X.length, c = Y[0].length, inner = Y.length;
        double[][] Z = new double[r][c];
        for (int i = 0; i < r; i++) {
            for (int j = 0; j < c; j++) {
                for (int k = 0; k < inner; k++) Z[i][j] += X[i][k] * Y[k][j];
            }
        }
        return Z;
    }

    static double[][] transpose(double[][] X) {
        int r = X.length, c = X[0].length;
        double[][] Z = new double[c][r];
        for (int i = 0; i < r; i++) {
            for (int j = 0; j < c; j++) Z[j][i] = X[i][j];
        }
        return Z;
    }

    static double[][] gramianDerivative(double[][] W) {
        // dW/ds = A W + W A^T + B B^T, W(0)=0.
        return add(add(mul(A, W), mul(W, transpose(A))), mul(B, transpose(B)));
    }

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

    static double[][] finiteHorizonGramian(double T, int steps) {
        double h = T / steps;
        double[][] W = { {0.0, 0.0}, {0.0, 0.0} };
        for (int i = 0; i < steps; i++) W = rk4Step(W, h);
        // Symmetrize to remove tiny numerical asymmetry.
        double avg = 0.5 * (W[0][1] + W[1][0]);
        W[0][1] = avg;
        W[1][0] = avg;
        return W;
    }

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

    static double minimumEnergy(double[][] W, double[] delta) {
        double det = W[0][0] * W[1][1] - W[0][1] * W[1][0];
        double inv00 = W[1][1] / det;
        double inv01 = -W[0][1] / det;
        double inv10 = -W[1][0] / det;
        double inv11 = W[0][0] / det;
        return delta[0] * (inv00 * delta[0] + inv01 * delta[1])
             + delta[1] * (inv10 * delta[0] + inv11 * delta[1]);
    }

    static void printMatrix(String name, double[][] X) {
        System.out.println(name + " =");
        for (double[] row : X) {
            for (double v : row) System.out.printf("%14.8f ", v);
            System.out.println();
        }
    }

    public static void main(String[] args) {
        double T = 2.0;
        double[][] W = finiteHorizonGramian(T, 4000);
        printMatrix("Wc(T)", W);
        double[] lambda = eigenvalues2x2Symmetric(W);
        System.out.printf("Eigenvalues: %.10f %.10f%n", lambda[0], lambda[1]);
        System.out.printf("Condition number: %.6f%n", lambda[1] / lambda[0]);
        System.out.printf("Energy to reach e1: %.8f%n", minimumEnergy(W, new double[]{1.0, 0.0}));
        System.out.printf("Energy to reach e2: %.8f%n", minimumEnergy(W, new double[]{0.0, 1.0}));
        System.out.printf("Best-axis energy: %.8f%n", 1.0 / lambda[1]);
        System.out.printf("Worst-axis energy: %.8f%n", 1.0 / lambda[0]);
    }
}

      

Chapter12_Lesson4.m

% Chapter12_Lesson4.m
% Modern Control — Chapter 12, Lesson 4
% Finite-horizon controllability Gramian, energy ellipsoid, and optional Simulink model.
% Requires base MATLAB. Control System Toolbox and Simulink sections are optional.

clear; clc; close all;

A = [-1 0; 0 -4];
B = [1; 0.08];
T = 2.0;
N = 4000;
dt = T/N;

% Wc(T)=int_0^T exp(A*s) B B' exp(A'*s) ds by trapezoidal quadrature.
W = zeros(2,2);
for k = 0:N
    s = k*dt;
    Phi = expm(A*s);
    G = Phi*B*B'*Phi';
    w = 1.0;
    if k == 0 || k == N
        w = 0.5;
    end
    W = W + w*G*dt;
end
W = 0.5*(W + W');

[V,D] = eig(W);
lambda = diag(D);
[lambda,idx] = sort(lambda,'ascend');
V = V(:,idx);

fprintf('Wc(T)=\n'); disp(W);
fprintf('Eigenvalues:\n'); disp(lambda.');
fprintf('Condition number = %.6g\n', lambda(end)/lambda(1));

E = @(delta) delta'*(W\delta);
e1 = [1;0];
e2 = [0;1];
fprintf('Energy to reach e1 = %.8f\n', E(e1));
fprintf('Energy to reach e2 = %.8f\n', E(e2));
fprintf('Best-axis energy = %.8f\n', 1/lambda(end));
fprintf('Worst-axis energy = %.8f\n', 1/lambda(1));

% Plot reachable energy ellipsoid delta' inv(W) delta <= rho^2.
rho = 1;
theta = linspace(0,2*pi,400);
circle = [cos(theta); sin(theta)];
ellipse = V*diag(rho*sqrt(lambda))*circle;
figure;
plot(ellipse(1,:), ellipse(2,:), 'LineWidth', 1.5); grid on; axis equal;
xlabel('state component x_1'); ylabel('state component x_2');
title('Finite-horizon controllability energy ellipsoid');
hold on;
for i = 1:2
    axisVector = V(:,i)*rho*sqrt(lambda(i));
    plot([0 axisVector(1)], [0 axisVector(2)], 'LineWidth', 2);
end

% Optional Control System Toolbox comparison for infinite-horizon stable systems.
if exist('ss','file') == 2 && exist('gram','file') == 2
    sys = ss(A,B,eye(2),zeros(2,1));
    Winf = gram(sys,'c');
    fprintf('Infinite-horizon Wc from Control System Toolbox:\n'); disp(Winf);
end

% Optional Simulink skeleton: create a state-space plant block for simulation.
% Feed a designed input signal u(t) to the input port and inspect state outputs.
if exist('new_system','file') == 2
    model = 'Chapter12_Lesson4_Simulink';
    if bdIsLoaded(model)
        close_system(model,0);
    end
    new_system(model);
    add_block('simulink/Sources/In1', [model '/u']);
    add_block('simulink/Continuous/State-Space', [model '/Plant']);
    add_block('simulink/Sinks/Out1', [model '/x']);
    set_param([model '/Plant'], 'A', mat2str(A), 'B', mat2str(B), ...
        'C', mat2str(eye(2)), 'D', mat2str(zeros(2,1)));
    add_line(model, 'u/1', 'Plant/1');
    add_line(model, 'Plant/1', 'x/1');
    save_system(model);
    fprintf('Created optional Simulink model: %s.slx\n', model);
end

      

Chapter12_Lesson4.nb

(* Chapter12_Lesson4.nb *)
(* Modern Control — Chapter 12, Lesson 4 *)
(* Wolfram Mathematica implementation of finite-horizon controllability energy directions. *)

ClearAll["Global`*"];
A = { {-1, 0}, {0, -4} };
B = { {1}, {0.08} };
T = 2;

W = NIntegrate[
   MatrixExp[A s].B.Transpose[B].Transpose[MatrixExp[A s]],
   {s, 0, T}
];
W = (W + Transpose[W])/2;

{lambda, Q} = Eigensystem[W];
ord = Ordering[lambda];
lambda = lambda[[ord]];
Q = Transpose[Q[[ord]]];

energy[delta_] := delta.LinearSolve[W, delta];

e1 = {1, 0};
e2 = {0, 1};

Print["Wc(T) = ", MatrixForm[W]];
Print["Eigenvalues = ", lambda];
Print["Condition number = ", Max[lambda]/Min[lambda]];
Print["Energy to reach e1 = ", energy[e1]];
Print["Energy to reach e2 = ", energy[e2]];
Print["Best-axis energy = ", 1/Max[lambda]];
Print["Worst-axis energy = ", 1/Min[lambda]];

rho = 1;
ellipse[theta_] := Q.DiagonalMatrix[rho Sqrt[lambda]].{Cos[theta], Sin[theta]};
ParametricPlot[ellipse[theta], {theta, 0, 2 Pi},
 AxesLabel -> {"x1", "x2"},
 PlotLabel -> "Finite-horizon controllability energy ellipsoid",
 AspectRatio -> Automatic,
 GridLines -> Automatic
]

      

Problem 1 (Principal-Axis Energy): Let \( \mathbf{W}_c=\mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^T \) with positive eigenvalues. Show that a unit displacement in direction \( \mathbf{q}_i \) requires energy \( 1/\lambda_i \).

Solution: Set \( \boldsymbol{\delta}=\mathbf{q}_i \). Then

\[ J_{\min}=\mathbf{q}_i^T\mathbf{W}_c^{-1}\mathbf{q}_i =\mathbf{q}_i^T\mathbf{Q}\boldsymbol{\Lambda}^{-1} \mathbf{Q}^T\mathbf{q}_i=\frac{1}{\lambda_i}. \]

Therefore the largest Gramian eigenvalue gives the easiest unit principal direction, and the smallest eigenvalue gives the most difficult unit principal direction.

Problem 2 (Energy Ellipsoid Axes): Prove that \( \boldsymbol{\delta}^T\mathbf{W}_c^{-1}\boldsymbol{\delta}\le\rho^2 \) has semi-axis lengths \( \rho\sqrt{\lambda_i} \).

Solution: Substitute \( \boldsymbol{\delta}=\mathbf{Q}\mathbf{z} \):

\[ \mathbf{z}^T\boldsymbol{\Lambda}^{-1}\mathbf{z}\le\rho^2 \quad\Longleftrightarrow\quad \sum_{i=1}^n \frac{z_i^2}{\rho^2\lambda_i}\le 1. \]

This is the standard ellipsoid equation; the axis associated with \( \mathbf{q}_i \) has length \( \rho\sqrt{\lambda_i} \).

Problem 3 (Weak Actuation): For \( \mathbf{A}=\operatorname{diag}(-1,-4) \) and \( \mathbf{B}=[1,\varepsilon]^T \), derive the finite-horizon Gramian and explain why the second state is difficult when \( \varepsilon \) is small.

Solution: Since \( e^{\mathbf{A}s}=\operatorname{diag}(e^{-s},e^{-4s}) \),

\[ \mathbf{W}_c(T)=\int_0^T \begin{bmatrix}e^{-s}\\0\end{bmatrix} \begin{bmatrix}e^{-s}&0\end{bmatrix}ds +\text{cross and second-mode terms induced by }\varepsilon, \]

\[ \mathbf{W}_c(T)= \begin{bmatrix} \dfrac{1-e^{-2T} }{2} & \varepsilon\dfrac{1-e^{-5T} }{5} \\ \varepsilon\dfrac{1-e^{-5T} }{5} & \varepsilon^2\dfrac{1-e^{-8T} }{8} \end{bmatrix}. \]

The second diagonal term is proportional to \( \varepsilon^2 \), so a unit target with strong second-state component generally requires energy proportional to \( 1/\varepsilon^2 \).

Problem 4 (Coordinate Transformation): If \( \mathbf{x}=\mathbf{T}\bar{\mathbf{x} } \), show that the minimum energy for a fixed physical terminal state is invariant.

Solution: The transformed Gramian is \( \bar{\mathbf{W} }_c=\mathbf{T}^{-1}\mathbf{W}_c\mathbf{T}^{-T} \) and \( \bar{\boldsymbol{\delta} }=\mathbf{T}^{-1}\boldsymbol{\delta} \). Thus

\[ \bar{\boldsymbol{\delta} }^T\bar{\mathbf{W} }_c^{-1} \bar{\boldsymbol{\delta} } =\boldsymbol{\delta}^T\mathbf{T}^{-T} \left(\mathbf{T}^T\mathbf{W}_c^{-1}\mathbf{T}\right) \mathbf{T}^{-1}\boldsymbol{\delta} =\boldsymbol{\delta}^T\mathbf{W}_c^{-1}\boldsymbol{\delta}. \]

Problem 5 (Finite Horizon Effect): For the scalar system \( \dot{x}=-a x+b u \), \( a>0 \), compute \( W_c(T) \) and discuss what happens as \( T\to0 \) and \( T\to\infty \).

Solution:

\[ W_c(T)=\int_0^T b^2 e^{-2as}ds= b^2\frac{1-e^{-2aT} }{2a}. \]

For short horizons, \( W_c(T)\approx b^2T \), so \( J_{\min}(\delta)\approx\delta^2/(b^2T) \), which becomes very large as \( T\to0 \). For long horizons, \( W_c(T)\to b^2/(2a) \), so the energy approaches a finite limit.

The controllability Gramian refines the binary rank-based concept of controllability into a quantitative energy geometry. The formula \( J_{\min}=\boldsymbol{\delta}^T\mathbf{W}_c^{-1}\boldsymbol{\delta} \) shows exactly why Gramian eigenvectors with large eigenvalues are easy to reach and eigenvectors with small eigenvalues are difficult. The reachable set under a finite energy budget is an ellipsoid whose axes reveal actuator authority, modal accessibility, and finite-horizon limitations.