Consider the continuous-time LTI state equation \( \dot{x}(t)=Ax(t)+Bu(t) \), where \( x(t)\in\mathbb{R}^n \) and \( u(t)\in\mathbb{R}^m \). Given an initial state \( x(0)=x_0 \), a final time \( T>0 \), and a desired terminal state \( x_f \), the state-transition formula gives

\[ x(T)=e^{AT}x_0+\int_0^T e^{A(T-t)}B u(t)\,dt. \]

The terminal constraint can therefore be written as \( \mathcal{L}u=z \), where

\[ z=x_f-e^{AT}x_0, \qquad \mathcal{L}u=\int_0^T e^{A(T-t)}B u(t)\,dt. \]

The minimum-energy steering problem is

\[ \min_{u\in L_2[0,T]} J(u), \qquad J(u)=\int_0^T u(t)^\top u(t)\,dt, \qquad \text{subject to }\quad \mathcal{L}u=z. \]

This is a quadratic optimization problem over a function space. The remarkable result is that its solution is completely determined by the finite-horizon controllability Gramian.

flowchart TD
  A["Start with x0, target xf, horizon T"] --> B["Compute free response exp(A T) x0"]
  B --> C["Required displacement z = xf minus free response"]
  C --> D["Compute finite-horizon Gramian Wc(T)"]
  D --> E["Solve Wc(T) lambda = z"]
  E --> F["Build u_star(t) = B^T exp(A^T (T-t)) lambda"]
  F --> G["Minimum energy = z^T lambda"]
  G --> H["Simulate x_dot = A x + B u_star(t) and verify x(T) = xf"]
        

From Lesson 1, the finite-horizon controllability Gramian is

\[ W_c(T)=\int_0^T e^{A(T-t)}BB^\top e^{A^\top(T-t)}\,dt. \]

By the change of variable \( s=T-t \), this is equivalently

\[ W_c(T)=\int_0^T e^{As}BB^\top e^{A^\top s}\,ds. \]

The matrix \( W_c(T) \) is symmetric positive semidefinite:

\[ q^\top W_c(T)q =\int_0^T \left\|B^\top e^{A^\top(T-t)}q\right\|_2^2\,dt \ge 0. \]

If \( W_c(T) \) is nonsingular, then every \( z\in\mathbb{R}^n \) can be represented as \( \mathcal{L}u=z \), and every state can be reached in time \( T \). If it is singular, then only terminal displacements in \( \operatorname{range}(W_c(T)) \) are reachable.

Assume \( W_c(T) \) is nonsingular. Define \( z=x_f-e^{AT}x_0 \). Then the unique minimum-energy input that steers \( x_0 \) to \( x_f \) over \( [0,T] \) is

\[ u^\star(t)=B^\top e^{A^\top(T-t)}W_c(T)^{-1}z. \]

The corresponding minimum energy is

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

Thus, the inverse Gramian acts as an energy metric on terminal displacements. A small Gramian eigenvalue means that motion in the associated state direction is energetically expensive.

Define the linear endpoint map \( \mathcal{L}:L_2[0,T]\to\mathbb{R}^n \) by

\[ \mathcal{L}u=\int_0^T e^{A(T-t)}B u(t)\,dt. \]

Its adjoint is the map \( \mathcal{L}^\ast:\mathbb{R}^n\to L_2[0,T] \) satisfying \( \int_0^T u(t)^\top(\mathcal{L}^\ast q)(t)\,dt =q^\top\mathcal{L}u \). Direct calculation gives

\[ (\mathcal{L}^\ast q)(t)=B^\top e^{A^\top(T-t)}q. \]

Therefore,

\[ \mathcal{L}\mathcal{L}^\ast q =\int_0^T e^{A(T-t)}BB^\top e^{A^\top(T-t)}q\,dt =W_c(T)q. \]

To solve \( \mathcal{L}u=z \) with minimum \( L_2 \) norm, choose \( u^\star=\mathcal{L}^\ast\lambda \) and impose the endpoint constraint:

\[ \mathcal{L}u^\star =\mathcal{L}\mathcal{L}^\ast\lambda =W_c(T)\lambda=z. \]

Since \( W_c(T) \) is nonsingular, \( \lambda=W_c(T)^{-1}z \), and hence \( u^\star(t)=B^\top e^{A^\top(T-t)}W_c(T)^{-1}z \). The energy becomes

\[ J(u^\star)=\int_0^T (u^\star(t))^\top u^\star(t)\,dt =\lambda^\top W_c(T)\lambda =z^\top W_c(T)^{-1}z. \]

For any other feasible input \( u=u^\star+v \), we have \( \mathcal{L}v=0 \). Then

\[ \int_0^T (u^\star(t))^\top v(t)\,dt =\lambda^\top\mathcal{L}v=0, \]

and therefore

\[ J(u)=J(u^\star)+\int_0^T v(t)^\top v(t)\,dt \ge J(u^\star). \]

This proves both optimality and uniqueness of the minimum-energy input.

With \( x_0=0 \), the set of states reachable with energy at most \( \rho \) is

\[ \mathcal{E}_\rho(T)=\left\{x_f\in\mathbb{R}^n: x_f^\top W_c(T)^{-1}x_f \le \rho\right\}. \]

If \( W_c(T)=V\Lambda V^\top \), where \( \Lambda=\operatorname{diag}(\lambda_1,\dots,\lambda_n) \), and \( y=V^\top x_f \), then

\[ x_f^\top W_c(T)^{-1}x_f =\sum_{i=1}^n \frac{y_i^2}{\lambda_i}. \]

Hence the ellipsoid semi-axis length along eigenvector \( v_i \) is \( \sqrt{\rho\lambda_i} \). Directions with large \( \lambda_i \) are easy to reach; directions with small \( \lambda_i \) are hard to reach.

flowchart TD
  A["Gramian Wc(T)"] --> B["Eigen-decomposition Wc = V Lambda V^T"]
  B --> C["Large eigenvalue direction"]
  B --> D["Small eigenvalue direction"]
  C --> E["Long ellipsoid axis: low required energy"]
  D --> F["Short ellipsoid axis: high required energy"]
  E --> G["Energy formula: sum y_i^2 / lambda_i"]
  F --> G
        

In multi-input systems, different actuators may have different physical costs. Let \( R=R^\top\succ 0 \). The weighted energy criterion is

\[ J_R(u)=\int_0^T u(t)^\top R u(t)\,dt. \]

The weighted Gramian is

\[ W_R(T)=\int_0^T e^{A(T-t)}BR^{-1}B^\top e^{A^\top(T-t)}\,dt. \]

If \( W_R(T) \) is nonsingular, the weighted minimum-energy input and energy are

\[ u_R^\star(t)=R^{-1}B^\top e^{A^\top(T-t)}W_R(T)^{-1}z, \qquad J_{R,\min}=z^\top W_R(T)^{-1}z. \]

This formula is especially useful when actuator channels have different saturation limits, power usage, or reliability. In this lesson, however, we keep the unweighted case as the central object because it directly exposes the role of \( W_c(T) \).

Consider the double integrator

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

Since \( e^{As}=\begin{bmatrix}1&s\\0&1\end{bmatrix} \), the Gramian is

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

Its inverse is

\[ W_c(T)^{-1}= \begin{bmatrix}12/T^3&-6/T^2\\-6/T^2&4/T\end{bmatrix}. \]

For \( x_0=0 \) and \( x_f=\begin{bmatrix}p_f&v_f\end{bmatrix}^\top \), the minimum energy is

\[ J_{\min} =\frac{12p_f^2}{T^3}-\frac{12p_fv_f}{T^2}+\frac{4v_f^2}{T}. \]

The corresponding optimal input is

\[ u^\star(t)= \frac{6(T-2t)}{T^3}p_f +\frac{6t-2T}{T^2}v_f. \]

Notice that requiring a position change becomes much less expensive as \( T \) increases because the leading term scales as \( T^{-3} \). Requiring a final velocity scales as \( T^{-1} \). This difference is not arbitrary; it is a direct consequence of how the input enters the second derivative of position.

In practice, one rarely forms \( W_c(T)^{-1} \) explicitly. Instead, compute \( \lambda \) from the linear system

\[ W_c(T)\lambda=z. \]

Then use \( u^\star(t)=B^\top e^{A^\top(T-t)}\lambda \) and \( J_{\min}=z^\top\lambda \). If \( W_c(T) \) is ill-conditioned, small numerical errors in \( z \) may cause large errors in \( \lambda \) and in the computed energy. A useful diagnostic is

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

Large \( \kappa(W_c(T)) \) means that the system may be mathematically controllable but practically difficult to steer in certain directions. This observation prepares the next lessons, where Gramian rank and Gramian conditioning will be related to controllability tests and easy-versus-hard state directions.

The following implementation computes \( W_c(T) \), solves \( W_c(T)\lambda=z \), constructs the minimum-energy input, and simulates the terminal state. It uses numpy and scipy, which are standard numerical libraries for state-space computation in Python.

Chapter12_Lesson2.py

# Chapter12_Lesson2.py
# Minimum-energy steering for continuous-time LTI systems.
# Requires: numpy, scipy
#
# Model:
#   x_dot(t) = A x(t) + B u(t)
# Goal:
#   steer x(0)=x0 to x(T)=xf while minimizing integral_0^T u(t)^T u(t) dt.

import numpy as np
from scipy.linalg import expm
from scipy.integrate import solve_ivp


def controllability_gramian(A: np.ndarray, B: np.ndarray, T: float, N: int = 4001) -> np.ndarray:
    """Finite-horizon controllability Gramian by trapezoidal quadrature.

    Wc(T) = integral_0^T exp(A (T-t)) B B^T exp(A^T (T-t)) dt.
    """
    n = A.shape[0]
    W = np.zeros((n, n), dtype=float)
    ts = np.linspace(0.0, T, N)
    dt = T / (N - 1)

    for k, t in enumerate(ts):
        E = expm(A * (T - t))
        integrand = E @ B @ B.T @ E.T
        weight = 0.5 if k == 0 or k == N - 1 else 1.0
        W += weight * integrand

    return W * dt


def minimum_energy_control(A: np.ndarray, B: np.ndarray, x0: np.ndarray, xf: np.ndarray, T: float):
    """Return u_star(t), Wc(T), lambda, and minimum energy.

    If Wc(T) is nonsingular:
        u_star(t) = B^T exp(A^T (T-t)) Wc(T)^(-1) (xf - exp(A T) x0)
        E_min     = z^T Wc(T)^(-1) z
    """
    W = controllability_gramian(A, B, T)
    z = xf - expm(A * T) @ x0

    # solve(W, z) is numerically better than inv(W) @ z
    lam = np.linalg.solve(W, z)

    def u_star(t: float) -> np.ndarray:
        return B.T @ expm(A.T * (T - t)) @ lam

    E_min = float(z.T @ lam)
    return u_star, W, lam, E_min


def simulate_closed_loop(A: np.ndarray, B: np.ndarray, u_star, x0: np.ndarray, T: float, N: int = 501):
    """Simulate x_dot = A x + B u_star(t)."""
    def rhs(t, x):
        return A @ x + B @ u_star(t)

    sol = solve_ivp(rhs, (0.0, T), x0, t_eval=np.linspace(0.0, T, N), rtol=1e-9, atol=1e-11)
    return sol.t, sol.y.T


if __name__ == "__main__":
    # Double integrator: x1_dot = x2, x2_dot = u.
    A = np.array([[0.0, 1.0],
                  [0.0, 0.0]])
    B = np.array([[0.0],
                  [1.0]])

    T = 2.0
    x0 = np.array([0.0, 0.0])
    xf = np.array([1.0, 0.0])

    u_star, W, lam, E_min = minimum_energy_control(A, B, x0, xf, T)
    t, x = simulate_closed_loop(A, B, u_star, x0, T)

    print("Wc(T) =")
    print(W)
    print("condition number of Wc(T):", np.linalg.cond(W))
    print("lambda = Wc(T)^(-1) z =", lam)
    print("minimum energy =", E_min)
    print("terminal state reached =", x[-1])
    print("target state =", xf)
    print("terminal error norm =", np.linalg.norm(x[-1] - xf))

    # Sample the optimal input.
    for ti in np.linspace(0.0, T, 5):
        print(f"u_star({ti:.2f}) = {u_star(float(ti))[0]: .6f}")

The C++ version uses Eigen for matrix operations and the unsupported matrix-functions module for the matrix exponential. In large control codes, this computation should be combined with robust linear solvers rather than explicit matrix inversion.

Chapter12_Lesson2.cpp

// Chapter12_Lesson2.cpp
// Minimum-energy steering for continuous-time LTI systems.
// Requires: Eigen 3 with unsupported MatrixFunctions module.
// Example compile command:
//   g++ -std=c++17 Chapter12_Lesson2.cpp -I /path/to/eigen -O2 -o Chapter12_Lesson2

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

using Eigen::MatrixXd;
using Eigen::VectorXd;

MatrixXd controllabilityGramian(const MatrixXd& A, const MatrixXd& B, double T, int N = 4001) {
    const int n = A.rows();
    MatrixXd W = MatrixXd::Zero(n, n);
    const double dt = T / static_cast<double>(N - 1);

    for (int k = 0; k < N; ++k) {
        double t = k * dt;
        MatrixXd E = (A * (T - t)).exp();
        MatrixXd F = E * B * B.transpose() * E.transpose();
        double weight = (k == 0 || k == N - 1) ? 0.5 : 1.0;
        W += weight * F;
    }
    return W * dt;
}

VectorXd minimumEnergyInput(
    double t,
    const MatrixXd& A,
    const MatrixXd& B,
    const VectorXd& lambda,
    double T
) {
    MatrixXd E = (A.transpose() * (T - t)).exp();
    return B.transpose() * E * lambda;
}

VectorXd dynamics(
    double t,
    const VectorXd& x,
    const MatrixXd& A,
    const MatrixXd& B,
    const VectorXd& lambda,
    double T
) {
    VectorXd u = minimumEnergyInput(t, A, B, lambda, T);
    return A * x + B * u;
}

int main() {
    MatrixXd A(2, 2);
    A << 0.0, 1.0,
         0.0, 0.0;

    MatrixXd B(2, 1);
    B << 0.0,
         1.0;

    double T = 2.0;
    VectorXd x0(2), xf(2);
    x0 << 0.0, 0.0;
    xf << 1.0, 0.0;

    MatrixXd W = controllabilityGramian(A, B, T);
    VectorXd z = (xf - (A * T).exp() * x0);
    VectorXd lambda = W.ldlt().solve(z);
    double Emin = z.dot(lambda);

    std::cout << std::setprecision(10);
    std::cout << "Wc(T):\n" << W << "\n\n";
    std::cout << "lambda:\n" << lambda << "\n\n";
    std::cout << "minimum energy: " << Emin << "\n\n";

    // RK4 simulation of x_dot = A x + B u_star(t)
    int N = 1000;
    double h = T / static_cast<double>(N);
    VectorXd x = x0;
    double t = 0.0;

    for (int k = 0; k < N; ++k) {
        VectorXd k1 = dynamics(t, x, A, B, lambda, T);
        VectorXd k2 = dynamics(t + 0.5*h, x + 0.5*h*k1, A, B, lambda, T);
        VectorXd k3 = dynamics(t + 0.5*h, x + 0.5*h*k2, A, B, lambda, T);
        VectorXd k4 = dynamics(t + h, x + h*k3, A, B, lambda, T);
        x += (h / 6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4);
        t += h;
    }

    std::cout << "terminal state reached:\n" << x << "\n\n";
    std::cout << "target state:\n" << xf << "\n\n";
    std::cout << "terminal error norm: " << (x - xf).norm() << "\n";

    return 0;
}

Java has no built-in matrix exponential in the standard library, so this example implements the double-integrator formula directly. This is useful pedagogically because students can compare the closed-form input with the Gramian-derived expression.

Chapter12_Lesson2.java

// Chapter12_Lesson2.java
// Minimum-energy steering for the double-integrator example.
// No external library is required.

public class Chapter12_Lesson2 {
    static double T = 2.0;
    static double pf = 1.0;  // desired final position
    static double vf = 0.0;  // desired final velocity

    // For double integrator from x0 = 0:
    // u*(t) = 6(T - 2t) pf / T^3 + (6t - 2T) vf / T^2
    static double uStar(double t) {
        return 6.0 * (T - 2.0 * t) * pf / Math.pow(T, 3)
             + (6.0 * t - 2.0 * T) * vf / Math.pow(T, 2);
    }

    static double[] f(double t, double[] x) {
        double[] dx = new double[2];
        dx[0] = x[1];
        dx[1] = uStar(t);
        return dx;
    }

    static double[] add(double[] a, double scale, double[] b) {
        return new double[] { a[0] + scale * b[0], a[1] + scale * b[1] };
    }

    static double[] rk4Step(double t, double[] x, double h) {
        double[] k1 = f(t, x);
        double[] k2 = f(t + 0.5 * h, add(x, 0.5 * h, k1));
        double[] k3 = f(t + 0.5 * h, add(x, 0.5 * h, k2));
        double[] k4 = f(t + h, add(x, h, k3));

        return new double[] {
            x[0] + h * (k1[0] + 2.0*k2[0] + 2.0*k3[0] + k4[0]) / 6.0,
            x[1] + h * (k1[1] + 2.0*k2[1] + 2.0*k3[1] + k4[1]) / 6.0
        };
    }

    static double minimumEnergy() {
        // E_min = 12 pf^2/T^3 - 12 pf vf/T^2 + 4 vf^2/T
        return 12.0 * pf * pf / Math.pow(T, 3)
             - 12.0 * pf * vf / Math.pow(T, 2)
             + 4.0 * vf * vf / T;
    }

    public static void main(String[] args) {
        int N = 2000;
        double h = T / N;
        double[] x = {0.0, 0.0};
        double t = 0.0;

        for (int k = 0; k <= 4; i++) {
            double ti = i * T / 4.0;
            System.out.printf("u_star(%.2f) = %.12f%n", ti, uStar(ti));
        }
    }
}

MATLAB provides expm, ode45, and reliable dense linear algebra, making it convenient for finite-horizon Gramian calculations. The same workspace variables can be used in a Simulink State-Space block and a MATLAB Function block.

Chapter12_Lesson2.m

% Chapter12_Lesson2.m
% Minimum-energy steering for continuous-time LTI systems.
% MATLAB implementation plus Simulink-oriented setup variables.

clear; clc;

A = [0 1;
     0 0];
B = [0;
     1];

T  = 2.0;
x0 = [0; 0];
xf = [1; 0];

% Compute Wc(T) from the differential Lyapunov equation:
% dW/dt = A W + W A' + B B', W(0) = 0.
n = size(A, 1);
gramian_ode = @(t, w) reshape(A*reshape(w, n, n) + reshape(w, n, n)*A' + B*B', [], 1);
[~, Wvec] = ode45(gramian_ode, [0 T], zeros(n*n, 1));
W = reshape(Wvec(end, :).', n, n);

z = xf - expm(A*T)*x0;
lambda = W \ z;

u_star = @(t) B' * expm(A'*(T - t)) * lambda;
E_min = z' * lambda;

fprintf('Wc(T) =\n');
disp(W);
fprintf('condition number of Wc(T) = %.6e\n', cond(W));
fprintf('minimum energy = %.12f\n', E_min);

% Simulate x_dot = A x + B u_star(t).
rhs = @(t, x) A*x + B*u_star(t);
[tgrid, xgrid] = ode45(rhs, linspace(0, T, 501), x0);

fprintf('terminal state reached =\n');
disp(xgrid(end, :).');
fprintf('target state =\n');
disp(xf);
fprintf('terminal error norm = %.6e\n', norm(xgrid(end, :).'-xf));

% Plot state and optimal input.
figure;
plot(tgrid, xgrid, 'LineWidth', 1.5);
grid on;
xlabel('time');
ylabel('states');
legend('x_1 position', 'x_2 velocity');

uvals = arrayfun(u_star, tgrid);
figure;
plot(tgrid, uvals, 'LineWidth', 1.5);
grid on;
xlabel('time');
ylabel('u^*(t)');

% Simulink implementation idea:
% 1. Use a State-Space block with matrices A, B, C = eye(2), D = zeros(2,1).
% 2. Use a MATLAB Function block to compute:
%       u = B' * expm(A'*(T - t)) * lambda
%    where T, A, B, and lambda are workspace variables.
% 3. Feed Clock into the MATLAB Function block and feed u into the State-Space block.
% 4. Use Scope blocks to compare x(T) with xf.

Mathematica can compute the double-integrator Gramian symbolically and then verify the resulting input numerically. This is useful for checking closed-form formulas before implementing them in lower-level languages.

Chapter12_Lesson2.nb

(* Chapter12_Lesson2.nb *)
(* Minimum-energy steering for a double-integrator system. *)

Clear["Global`*"];

A = { {0, 1}, {0, 0} };
B = { {0}, {1} };
Tfinal = 2;
x0 = {0, 0};
xf = {1, 0};

Wc = Integrate[
   MatrixExp[A s].B.Transpose[B].Transpose[MatrixExp[A s]],
   {s, 0, Tfinal}
];

z = xf - MatrixExp[A Tfinal].x0;
lambda = LinearSolve[Wc, z];

uStar[t_] := First[Transpose[B].MatrixExp[Transpose[A] (Tfinal - t)].lambda];

Emin = z.lambda;

xSol = NDSolveValue[
   {
    x1'[t] == x2[t],
    x2'[t] == uStar[t],
    x1[0] == x0[[1]],
    x2[0] == x0[[2]]
   },
   {x1, x2},
   {t, 0, Tfinal}
];

terminalState = {xSol[[1]][Tfinal], xSol[[2]][Tfinal]};
terminalError = Norm[terminalState - xf];

Print["Wc(T) = ", MatrixForm[Wc]];
Print["lambda = ", lambda];
Print["minimum energy = ", Emin];
Print["terminal state = ", terminalState];
Print["terminal error norm = ", terminalError];

Plot[uStar[t], {t, 0, Tfinal},
 AxesLabel -> {"t", "u*(t)"},
 PlotLabel -> "Minimum-energy input"]

Problem 1 (Endpoint Constraint): For \( \dot{x}=Ax+Bu \), derive the terminal constraint used in the minimum-energy problem.

Solution: The state-transition formula gives

\[ x(T)=e^{AT}x_0+\int_0^T e^{A(T-t)}B u(t)\,dt. \]

Therefore, reaching \( x_f \) is equivalent to

\[ \int_0^T e^{A(T-t)}B u(t)\,dt =x_f-e^{AT}x_0=z. \]

Problem 2 (Positive Semidefiniteness): Prove that \( W_c(T) \) is positive semidefinite.

Solution: For any \( q\in\mathbb{R}^n \),

\[ q^\top W_c(T)q =\int_0^T q^\top e^{A(T-t)}BB^\top e^{A^\top(T-t)}q\,dt =\int_0^T \left\|B^\top e^{A^\top(T-t)}q\right\|_2^2\,dt \ge 0. \]

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

Problem 3 (Minimum-Energy Formula): Assume \( W_c(T) \) is nonsingular. Show that the energy of \( u^\star(t)=B^\top e^{A^\top(T-t)}W_c(T)^{-1}z \) equals \( z^\top W_c(T)^{-1}z \).

Solution: Let \( \lambda=W_c(T)^{-1}z \). Then

\[ J(u^\star)=\int_0^T \lambda^\top e^{A(T-t)}BB^\top e^{A^\top(T-t)}\lambda\,dt =\lambda^\top W_c(T)\lambda =z^\top W_c(T)^{-1}z. \]

Problem 4 (Double Integrator Gramian): 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^{As}=\begin{bmatrix}1&s\\0&1\end{bmatrix} \),

\[ e^{As}BB^\top e^{A^\top s} = \begin{bmatrix}s^2&s\\s&1\end{bmatrix}. \]

Integrating entrywise gives

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

Problem 5 (Energy Scaling): For the double integrator with \( x_0=0 \) and \( x_f=\begin{bmatrix}p_f&0\end{bmatrix}^\top \), show how the minimum energy scales with \( T \).

Solution: Substituting \( v_f=0 \) in the formula from Section 7 gives

\[ J_{\min}=\frac{12p_f^2}{T^3}. \]

Therefore, doubling the allowed transfer time reduces the position-only steering energy by a factor of \( 8 \).

Problem 6 (Reachability Ellipsoid): Let \( W_c(T)=V\Lambda V^\top \). Explain why small eigenvalues indicate difficult state directions.

Solution: For \( x_f=Vy \),

\[ J_{\min}=x_f^\top W_c(T)^{-1}x_f =\sum_{i=1}^n \frac{y_i^2}{\lambda_i}. \]

If \( \lambda_i \) is small, even a modest component \( y_i \) in that eigenvector direction produces a large energy contribution.

The controllability Gramian provides an exact energy metric for finite horizon state steering. For a desired terminal displacement \( z=x_f-e^{AT}x_0 \), the minimum energy is \( z^\top W_c(T)^{-1}z \), and the optimal input is \( B^\top e^{A^\top(T-t)}W_c(T)^{-1}z \). Gramian eigenvalues define the reachability ellipsoid: large eigenvalues correspond to easy directions and small eigenvalues correspond to difficult directions. Numerical work should solve \( W_c(T)\lambda=z \) rather than forming \( W_c(T)^{-1} \) explicitly.