Consider the continuous-time LTI state equation already introduced in earlier lessons: \( \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t) \), where \( \mathbf{x}(t)\in\mathbb{R}^n \) and \( \mathbf{u}(t)\in\mathbb{R}^m \). The finite-time reachability question is not only whether a state can be reached, but whether it can be reached at a prescribed terminal time \( T \).

The steering problem is:

\[ \text{Given } \mathbf{x}_0,\,\mathbf{x}_T,\,T>0, \text{ find } \mathbf{u}:[0,T]\mapsto\mathbb{R}^m \text{ such that } \mathbf{x}(0)=\mathbf{x}_0, \;\mathbf{x}(T)=\mathbf{x}_T. \]

From the variation-of-constants formula, the terminal state is

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

Therefore the input must generate the displacement \( \mathbf{d}_T=\mathbf{x}_T-e^{\mathbf{A}T}\mathbf{x}_0 \) through the finite-time input-to-state operator \( \mathcal{L}_T \):

\[ \mathcal{L}_T\mathbf{u} = \int_0^T e^{\mathbf{A}(T-\tau)}\mathbf{B}\mathbf{u}(\tau)\,d\tau, \qquad \mathbf{d}_T\in\operatorname{Range}(\mathcal{L}_T). \]

flowchart TD
  A["Initial state x0"] --> B["Free motion: exp(A T) x0"]
  C["Input signal u(t), 0 <= t <= T"] --> D["Forced displacement via finite-time endpoint map"]
  B --> E["Terminal state x(T)"]
  D --> E
  E --> F["Reach target xT?"]
  F --> G["Need displacement dT = xT - exp(A T) x0 in reachable range"]
        

For a fixed duration \( T>0 \), define the finite-time reachable displacement set from the origin by

\[ \mathcal{R}(T) = \left\{ \int_0^T e^{\mathbf{A}(T-\tau)}\mathbf{B}\mathbf{u}(\tau)\,d\tau : \mathbf{u}\in L_2([0,T],\mathbb{R}^m) \right\}. \]

The full reachable set from \( \mathbf{x}_0 \) at exactly time \( T \) is the affine set

\[ \mathcal{X}(T;\mathbf{x}_0) = e^{\mathbf{A}T}\mathbf{x}_0+\mathcal{R}(T). \]

Thus finite-time reachability is naturally an affine geometry problem. The input determines only the displacement subspace, while the initial condition contributes the free-motion offset \( e^{\mathbf{A}T}\mathbf{x}_0 \).

A useful finite-horizon matrix, later studied in more detail as a controllability Gramian, is

\[ \mathbf{W}(T)= \int_0^T e^{\mathbf{A}s}\mathbf{B}\mathbf{B}^\top e^{\mathbf{A}^\top s}\,ds. \]

The range of \( \mathbf{W}(T) \) equals the reachable displacement subspace \( \mathcal{R}(T) \). Consequently, \( \mathbf{x}_T \) is reachable from \( \mathbf{x}_0 \) at time \( T \) if and only if

\[ \mathbf{d}_T = \mathbf{x}_T-e^{\mathbf{A}T}\mathbf{x}_0 \in \operatorname{Range}(\mathbf{W}(T)). \]

The key identity is obtained from orthogonality. For any vector \( \mathbf{v}\in\mathbb{R}^n \),

\[ \mathbf{v}^\top \mathbf{W}(T)\mathbf{v} = \int_0^T \mathbf{v}^\top e^{\mathbf{A}s}\mathbf{B}\mathbf{B}^\top e^{\mathbf{A}^\top s}\mathbf{v}\,ds = \int_0^T \left\|\mathbf{B}^\top e^{\mathbf{A}^\top s}\mathbf{v}\right\|_2^2\,ds. \]

Hence \( \mathbf{v}\in\operatorname{Null}(\mathbf{W}(T)) \) if and only if \( \mathbf{B}^\top e^{\mathbf{A}^\top s}\mathbf{v}=\mathbf{0} \) for every \( s\in[0,T] \). By differentiating at \( s=0 \), this implies

\[ \mathbf{B}^\top(\mathbf{A}^\top)^k\mathbf{v}=\mathbf{0}, \qquad k=0,1,2,\ldots. \]

Equivalently,

\[ \mathbf{v}^\top \mathbf{A}^k\mathbf{B}=\mathbf{0}, \qquad k=0,1,2,\ldots. \]

Because higher powers of \( \mathbf{A} \) are linearly dependent on the first \( n \) powers by the Cayley-Hamilton theorem, the finite set \( \mathbf{B},\mathbf{A}\mathbf{B},\ldots, \mathbf{A}^{n-1}\mathbf{B} \) determines the same orthogonal complement. Therefore

\[ \operatorname{Range}(\mathbf{W}(T)) = \operatorname{span}\left\{ \mathbf{B},\mathbf{A}\mathbf{B},\ldots, \mathbf{A}^{n-1}\mathbf{B} \right\}, \qquad T>0. \]

This proof explains an important fact: for continuous-time LTI systems, the reachable subspace itself does not change with the positive horizon length. What changes strongly with \( T \) is the required input magnitude or energy.

A state can be reachable for every positive terminal time but still require impractically large input when the available time is short. Among all inputs satisfying the terminal constraint, the minimum-energy input minimizes

\[ J(\mathbf{u})=\int_0^T \mathbf{u}^\top(t)\mathbf{u}(t)\,dt. \]

When \( \mathbf{d}_T\in\operatorname{Range}(\mathbf{W}(T)) \), one minimum-norm steering input is

\[ \mathbf{u}^\star(t) = \mathbf{B}^\top e^{\mathbf{A}^\top(T-t)} \mathbf{W}^\dagger(T)\mathbf{d}_T, \qquad 0\le t\le T, \]

where \( \mathbf{W}^\dagger(T) \) is the Moore-Penrose inverse. If \( \mathbf{W}(T) \) is nonsingular, this becomes the ordinary inverse. The minimum energy is

\[ J^\star = \mathbf{d}_T^\top\mathbf{W}^\dagger(T)\mathbf{d}_T. \]

Shorter horizons usually reduce the entries of \( \mathbf{W}(T) \) and amplify \( \mathbf{W}^\dagger(T) \). Therefore finite-time reachability must be separated from practical reachability: the latter also depends on actuator saturation, bandwidth, and available energy.

flowchart TD
  A["Choose terminal time T"] --> B["Compute displacement dT"]
  B --> C["Check dT in Range W(T)"]
  C -->|not reachable| D["No admissible input exists"]
  C -->|reachable| E["Compute u*(t) using W(T)^dagger"]
  E --> F["Energy J* = dT^T W(T)^dagger dT"]
  F --> G["If J* too large, increase T or redesign actuator placement"]
        

The double integrator, \( \dot{x}_1=x_2,\;\dot{x}_2=u \), has

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

The finite-time matrix is

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

Its determinant is

\[ \det \mathbf{W}(T)=\frac{T^4}{12},\qquad T>0. \]

Thus the double integrator can reach any state in any positive time, but the energy can grow rapidly as \( T \) becomes small. For example, steering from \( [0,0]^\top \) to \( [1,0]^\top \) gives

\[ J^\star = \begin{bmatrix}1&0\end{bmatrix} \mathbf{W}^{-1}(T) \begin{bmatrix}1\\0\end{bmatrix} = \frac{12}{T^3}. \]

Halving the available duration multiplies this particular transfer energy by \( 2^3=8 \). This is the basic mathematical reason why aggressive finite-time maneuvers require large actuator authority.

The statement \( \mathbf{x}_T\in\mathcal{X}(T;\mathbf{x}_0) \) answers an existence question. Engineering design requires additional duration-dependent checks:

\[ \|\mathbf{u}(t)\|_2\le u_{\max},\qquad \int_0^T \|\mathbf{u}(t)\|_2^2dt\le E_{\max},\qquad T_{\min}\le T\le T_{\max}. \]

These inequalities are not part of the basic linear reachability definition, but they are essential in implementation. If a target is mathematically reachable but violates input limits, a longer duration, a different actuator matrix \( \mathbf{B} \), or a different trajectory strategy is needed.

For a reachable target, increasing \( T \) usually enlarges the finite-horizon matrix in a positive-semidefinite sense:

\[ \mathbf{W}(T_2)-\mathbf{W}(T_1) = \int_{T_1}^{T_2} e^{\mathbf{A}s}\mathbf{B}\mathbf{B}^\top e^{\mathbf{A}^\top s}ds \succeq \mathbf{0}, \qquad 0<T_1<T_2. \]

This does not mean every transfer becomes easy, especially for unstable \( \mathbf{A} \) or poorly aligned inputs, but it gives a rigorous first explanation of why duration helps.

This Python implementation uses numpy and scipy to compute the finite-time matrix, check reachability of the requested terminal state, and sample the minimum-energy input.

# Chapter10_Lesson3.py
# Finite-time reachability and minimum-energy steering for continuous-time LTI systems.
# Requires: numpy, scipy

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


def controllability_matrix(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """Return [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_time_gramian(A: np.ndarray, B: np.ndarray, T: float) -> np.ndarray:
    """Compute W(T)=int_0^T exp(A s) B B' exp(A' s) ds."""
    if T <= 0:
        raise ValueError("T must be positive.")

    def integrand(s):
        E = expm(A * s)
        return E @ B @ B.T @ E.T

    W, _ = quad_vec(integrand, 0.0, T)
    return 0.5 * (W + W.T)


def min_energy_control(A: np.ndarray, B: np.ndarray, x0: np.ndarray, xT: np.ndarray, T: float):
    """Return a callable u(t), the Gramian W(T), displacement d, and minimum energy."""
    PhiT = expm(A * T)
    W = finite_time_gramian(A, B, T)
    d = xT - PhiT @ x0
    Winv = np.linalg.pinv(W)

    if np.linalg.norm(W @ Winv @ d - d) > 1e-8:
        raise ValueError("The requested terminal state is not reachable over this horizon.")

    def u(t):
        return B.T @ expm(A.T * (T - t)) @ Winv @ d

    energy = float(d.T @ Winv @ d)
    return u, W, d, energy


def simulate_closed_form(A, B, x0, u, T, steps=400):
    """Approximate x(T)=Phi(T)x0+int_0^T Phi(T-t)B u(t) dt by trapezoidal integration."""
    ts = np.linspace(0.0, T, steps + 1)
    total = np.zeros_like(x0, dtype=float)
    for k, t in enumerate(ts):
        weight = 0.5 if k == 0 or k == len(ts) - 1 else 1.0
        total += weight * (expm(A * (T - t)) @ B @ np.atleast_1d(u(t)))
    total *= T / steps
    return expm(A * T) @ x0 + total


def demo():
    # 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]])
    x0 = np.array([0.0, 0.0])
    xT = np.array([1.0, 0.0])

    C = controllability_matrix(A, B)
    print("Controllability matrix:\n", C)
    print("rank =", np.linalg.matrix_rank(C))

    for T in [0.5, 1.0, 2.0, 4.0]:
        u, W, d, energy = min_energy_control(A, B, x0, xT, T)
        xT_num = simulate_closed_form(A, B, x0, u, T)
        print(f"\nT = {T:.2f}")
        print("W(T) =\n", W)
        print("det(W) =", np.linalg.det(W))
        print("minimum energy =", energy)
        print("reconstructed x(T) =", xT_num)
        print("sample u(0), u(T/2), u(T) =", u(0.0), u(T / 2.0), u(T))


if __name__ == "__main__":
    demo()

The C++ version uses the Eigen library. The finite-time matrix is computed by trapezoidal quadrature, while the matrix exponential is provided by Eigen's unsupported matrix-functions module.

// Chapter10_Lesson3.cpp
// Finite-time reachability and minimum-energy steering using Eigen.
// Compile example:
// g++ -std=c++17 Chapter10_Lesson3.cpp -I /path/to/eigen -O2 -o Chapter10_Lesson3

#include <iostream>
#include <vector>
#include <cmath>
#include <Eigen/Dense>
#include <unsupported/Eigen/MatrixFunctions>

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

Matrix controllabilityMatrix(const Matrix& A, const Matrix& B) {
    const int n = A.rows();
    Matrix C(n, n * B.cols());
    Matrix Apow = Matrix::Identity(n, n);
    for (int k = 0; k < n; ++k) {
        C.block(0, k * B.cols(), n, B.cols()) = Apow * B;
        Apow = Apow * A;
    }
    return C;
}

Matrix finiteTimeGramian(const Matrix& A, const Matrix& B, double T, int steps = 5000) {
    const int n = A.rows();
    Matrix W = Matrix::Zero(n, n);
    const double h = T / steps;
    for (int k = 0; k <= steps; ++k) {
        double s = k * h;
        double weight = (k == 0 || k == steps) ? 0.5 : 1.0;
        Matrix E = (A * s).exp();
        W += weight * (E * B * B.transpose() * E.transpose());
    }
    W *= h;
    return 0.5 * (W + W.transpose());
}

Matrix pseudoInverse(const Matrix& M, double tol = 1e-10) {
    Eigen::JacobiSVD<Matrix> svd(M, Eigen::ComputeThinU | Eigen::ComputeThinV);
    Vector s = svd.singularValues();
    Matrix Sinv = Matrix::Zero(svd.matrixV().cols(), svd.matrixU().cols());
    for (int i = 0; i < s.size(); ++i) {
        if (s(i) > tol) {
            Sinv(i, i) = 1.0 / s(i);
        }
    }
    return svd.matrixV() * Sinv * svd.matrixU().transpose();
}

Vector minimumEnergyInput(const Matrix& A, const Matrix& B, const Matrix& Winv,
                          const Vector& d, double T, double t) {
    Matrix E = (A.transpose() * (T - t)).exp();
    return B.transpose() * E * Winv * d;
}

int main() {
    Matrix A(2, 2);
    A << 0.0, 1.0,
         0.0, 0.0;
    Matrix B(2, 1);
    B << 0.0,
         1.0;
    Vector x0(2), xT(2);
    x0 << 0.0, 0.0;
    xT << 1.0, 0.0;

    Matrix C = controllabilityMatrix(A, B);
    Eigen::FullPivLU<Matrix> lu(C);
    std::cout << "Controllability matrix:\n" << C << "\n";
    std::cout << "rank = " << lu.rank() << "\n\n";

    for (double T : {0.5, 1.0, 2.0, 4.0}) {
        Matrix PhiT = (A * T).exp();
        Matrix W = finiteTimeGramian(A, B, T);
        Matrix Winv = pseudoInverse(W);
        Vector d = xT - PhiT * x0;
        double energy = d.transpose() * Winv * d;

        std::cout << "T = " << T << "\n";
        std::cout << "W(T):\n" << W << "\n";
        std::cout << "det(W) = " << W.determinant() << "\n";
        std::cout << "minimum energy = " << energy << "\n";
        std::cout << "u(0) = " << minimumEnergyInput(A, B, Winv, d, T, 0.0).transpose() << "\n";
        std::cout << "u(T/2) = " << minimumEnergyInput(A, B, Winv, d, T, T / 2.0).transpose() << "\n";
        std::cout << "u(T) = " << minimumEnergyInput(A, B, Winv, d, T, T).transpose() << "\n\n";
    }
    return 0;
}

The Java version uses the analytic finite-time matrix for the double integrator, avoiding external numerical libraries while preserving the main finite-time reachability computation.

// Chapter10_Lesson3.java
// Analytic finite-time reachability example for the double integrator.
// Compile and run:
// javac Chapter10_Lesson3.java && java Chapter10_Lesson3

public class Chapter10_Lesson3 {
    static double[][] gramianDoubleIntegrator(double T) {
        return new double[][] {
            {Math.pow(T, 3.0) / 3.0, Math.pow(T, 2.0) / 2.0},
            {Math.pow(T, 2.0) / 2.0, T}
        };
    }

    static double[][] inverse2x2(double[][] M) {
        double det = M[0][0] * M[1][1] - M[0][1] * M[1][0];
        if (Math.abs(det) < 1e-12) {
            throw new IllegalArgumentException("Matrix is numerically singular.");
        }
        return new double[][] {
            { M[1][1] / det, -M[0][1] / det},
            {-M[1][0] / det,  M[0][0] / det}
        };
    }

    static double[] matVec(double[][] M, double[] v) {
        return new double[] {
            M[0][0] * v[0] + M[0][1] * v[1],
            M[1][0] * v[0] + M[1][1] * v[1]
        };
    }

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

    static double uMinEnergy(double T, double t, double[] lambda) {
        // For A=[[0,1],[0,0]], B=[[0],[1]], B' exp(A'(T-t)) = [T-t, 1].
        return (T - t) * lambda[0] + lambda[1];
    }

    public static void main(String[] args) {
        // Steering x0=[0,0] to xT=[1,0] for the double integrator.
        double[] d = {1.0, 0.0};

        for (double T : new double[] {0.5, 1.0, 2.0, 4.0}) {
            double[][] W = gramianDoubleIntegrator(T);
            double[][] Winv = inverse2x2(W);
            double[] lambda = matVec(Winv, d);
            double energy = dot(d, lambda);

            System.out.println("T = " + T);
            System.out.printf("W(T) = [[%.6f, %.6f], [%.6f, %.6f]]%n",
                    W[0][0], W[0][1], W[1][0], W[1][1]);
            System.out.println("minimum energy = " + energy);
            System.out.println("u(0) = " + uMinEnergy(T, 0.0, lambda));
            System.out.println("u(T/2) = " + uMinEnergy(T, T / 2.0, lambda));
            System.out.println("u(T) = " + uMinEnergy(T, T, lambda));
            System.out.println();
        }
    }
}

The MATLAB script uses expm, numerical integration, and ctrb. The comments also describe a direct Simulink block-level implementation.

% Chapter10_Lesson3.m
% Finite-time reachability and minimum-energy steering in MATLAB/Simulink.
% Uses Control System Toolbox for ctrb; the core Gramian computation is explicit.

clear; clc;

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

Cmat = ctrb(A, B);
disp("Controllability matrix:");
disp(Cmat);
disp("rank = " + rank(Cmat));

for T = [0.5 1.0 2.0 4.0]
    PhiT = expm(A*T);
    W = integral(@(s) gramianIntegrand(A, B, s), 0, T, ...
                 "ArrayValued", true, "RelTol", 1e-10, "AbsTol", 1e-12);
    W = 0.5*(W + W');
    d = xT - PhiT*x0;

    if norm(W*pinv(W)*d - d) > 1e-8
        error("Requested terminal state is not reachable over this horizon.");
    end

    Winv = pinv(W);
    energy = d' * Winv * d;

    u = @(t) B' * expm(A'*(T-t)) * Winv * d;

    fprintf("\nT = %.2f\n", T);
    disp("W(T) = "); disp(W);
    fprintf("det(W) = %.12g\n", det(W));
    fprintf("minimum energy = %.12g\n", energy);
    fprintf("u(0) = %.12g, u(T/2) = %.12g, u(T) = %.12g\n", ...
            u(0), u(T/2), u(T));
end

% Simulink note:
% Build a State-Space block with A=[0 1;0 0], B=[0;1], C=eye(2), D=zeros(2,1).
% Drive it with a MATLAB Function block implementing u(t)=B'*expm(A'*(T-t))*pinv(W)*d.
% Use a Clock block as the function input t and stop the simulation at t=T.

function G = gramianIntegrand(A, B, s)
    E = expm(A*s);
    G = E*B*B'*E';
end

The Mathematica notebook expression below computes the same quantities using symbolic-style matrix operations, numerical integration, and PseudoInverse.

Notebook[{
Cell["Chapter10_Lesson3.nb", "Title"],
Cell["Finite-Time Reachability and Minimum-Energy Control Duration", "Section"],
Cell[BoxData[
 RowBox[{
  RowBox[{"A", "=", RowBox[{"{ {", RowBox[{"0", ",", "1"}], "}", ",", RowBox[{"{", RowBox[{"0", ",", "0"}], "} }"}]}], ";"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"B", "=", RowBox[{"{ {", "0", "}", ",", RowBox[{"{", "1", "}"}]}]}], ";"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"x0", "=", RowBox[{"{", RowBox[{"0", ",", "0"}], "}"}]}], ";",
  RowBox[{"xT", "=", RowBox[{"{", RowBox[{"1", ",", "0"}], "}"}]}], ";"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"controllabilityMatrix", "[", RowBox[{"A_", ",", "B_"}], "]"}], ":=",
  RowBox[{"ArrayFlatten", "[", RowBox[{"{", RowBox[{"Table", "[", RowBox[{
  RowBox[{"MatrixPower", "[", RowBox[{"A", ",", "k"}], "]"}], ".", "B"}], ",",
  RowBox[{"{", RowBox[{"k", ",", "0", ",", RowBox[{
  RowBox[{"Length", "[", "A", "]"}], "-", "1"}]}], "}"}]}], "]"}], "}"}], "]"}]}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"finiteTimeGramian", "[", RowBox[{"A_", ",", "B_", ",", "T_"}], "]"}], ":=",
  RowBox[{"NIntegrate", "[", RowBox[{
  RowBox[{"MatrixExp", "[", RowBox[{"A", " ", "s"}], "]"}], ".", "B", ".",
  RowBox[{"Transpose", "[", "B", "]"}], ".", RowBox[{"Transpose", "[",
  RowBox[{"MatrixExp", "[", RowBox[{"A", " ", "s"}], "]"}], "]"}]}], ",",
  RowBox[{"{", RowBox[{"s", ",", "0", ",", "T"}], "}"}]}], "]"}]}]], "Input"],
Cell[BoxData[
 RowBox[{"MatrixRank", "[", RowBox[{"controllabilityMatrix", "[", RowBox[{"A", ",", "B"}], "]"}], "]"}]], "Input"],
Cell[BoxData[
 RowBox[{"Table", "[", RowBox[{
  RowBox[{"Module", "[", RowBox[{
  RowBox[{"{", RowBox[{"W", ",", "d", ",", "energy", ",", "u"}], "}"}], ",",
  RowBox[{
  RowBox[{"W", "=", RowBox[{"finiteTimeGramian", "[", RowBox[{"A", ",", "B", ",", "T"}], "]"}]}], ";",
  RowBox[{"d", "=", RowBox[{"xT", "-", RowBox[{
  RowBox[{"MatrixExp", "[", RowBox[{"A", " ", "T"}], "]"}], ".", "x0"}]}]}], ";",
  RowBox[{"energy", "=", RowBox[{"d", ".", RowBox[{"PseudoInverse", "[", "W", "]"}], ".", "d"}]}], ";",
  RowBox[{"u", "=", RowBox[{"Function", "[", RowBox[{"t", ",",
  RowBox[{
  RowBox[{"Transpose", "[", "B", "]"}], ".", RowBox[{"Transpose", "[",
  RowBox[{"MatrixExp", "[", RowBox[{"A", " ", RowBox[{"(", RowBox[{"T", "-", "t"}], ")"}]}], "]"}], "]"}], ".",
  RowBox[{"PseudoInverse", "[", "W", "]"}], ".", "d"}]}], "]"}]}], ";",
  RowBox[{"{", RowBox[{"T", ",", "W", ",", "energy", ",",
  RowBox[{"u", "[", "0", "]"}], ",", RowBox[{"u", "[", RowBox[{"T", "/", "2"}], "]"}], ",", RowBox[{"u", "[", "T", "]"}]}], "}"}]}]}],
  "]"}], ",", RowBox[{"{", RowBox[{"T", ",", RowBox[{"{", RowBox[{"0.5", ",", "1.0", ",", "2.0", ",", "4.0"}], "}"}]}], "}"}]}], "]"}]], "Input"]
}]

Problem 1 (Scalar System): Consider \( \dot{x}=ax+bu \) with \( b\ne0 \). Derive the finite-time condition for reaching \( x_T \) from \( x_0 \) at time \( T>0 \).

Solution: The terminal state is \( x_T=e^{aT}x_0+\int_0^T e^{a(T-\tau)}b u(\tau)d\tau \). Since \( b\ne0 \), the one-dimensional displacement space is all of \( \mathbb{R} \). Therefore every scalar terminal state is reachable for every \( T>0 \). The finite-time matrix is

\[ W(T)=\int_0^T e^{2as}b^2ds = \begin{cases} b^2\dfrac{e^{2aT}-1}{2a}, & a\ne0,\\ b^2T, & a=0. \end{cases} \]

Because \( W(T)>0 \) for \( T>0 \), all scalar targets are reachable.

Problem 2 (Double-Integrator Energy Scaling): For the double integrator, show that steering \( [0,0]^\top \) to \( [p,0]^\top \) in time \( T \) requires minimum energy \( 12p^2/T^3 \).

Solution: From Section 5, \( \mathbf{W}(T)=\begin{bmatrix}T^3/3&T^2/2\\T^2/2&T\end{bmatrix} \). Its inverse is

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

With \( \mathbf{d}=[p,0]^\top \),

\[ J^\star=\mathbf{d}^\top\mathbf{W}^{-1}(T)\mathbf{d} = \begin{bmatrix}p&0\end{bmatrix} \begin{bmatrix} 12/T^3 & -6/T^2\\ -6/T^2 & 4/T \end{bmatrix} \begin{bmatrix}p\\0\end{bmatrix} = \frac{12p^2}{T^3}. \]

Problem 3 (Unreachable Component): Let \( \mathbf{A}=\mathbf{0}_{2\times2} \) and \( \mathbf{B}=\begin{bmatrix}1\\0\end{bmatrix} \). Characterize all states reachable from the origin at time \( T>0 \).

Solution: Since \( e^{\mathbf{A}s}=\mathbf{I} \),

\[ \mathbf{W}(T) = \int_0^T \begin{bmatrix}1\\0\end{bmatrix} \begin{bmatrix}1&0\end{bmatrix}ds = \begin{bmatrix}T&0\\0&0\end{bmatrix}. \]

Therefore \( \operatorname{Range}(\mathbf{W}(T)) \) is the first coordinate axis. From the origin, reachable terminal states are exactly \( [\alpha,0]^\top \), where \( \alpha\in\mathbb{R} \). The second coordinate is unreachable because no input enters that state equation.

Problem 4 (Monotonicity of the Finite-Time Matrix): Prove that \( \mathbf{W}(T_2)-\mathbf{W}(T_1)\succeq\mathbf{0} \) whenever \( 0<T_1<T_2 \).

Solution: Direct subtraction gives

\[ \mathbf{W}(T_2)-\mathbf{W}(T_1) = \int_{T_1}^{T_2} e^{\mathbf{A}s}\mathbf{B}\mathbf{B}^\top e^{\mathbf{A}^\top s}ds. \]

For any \( \mathbf{v} \),

\[ \mathbf{v}^\top(\mathbf{W}(T_2)-\mathbf{W}(T_1))\mathbf{v} = \int_{T_1}^{T_2} \left\|\mathbf{B}^\top e^{\mathbf{A}^\top s}\mathbf{v}\right\|_2^2ds \ge0. \]

Hence the difference is positive semidefinite.

Problem 5 (Transfer Feasibility Test): Suppose \( \mathbf{W}(T) \) is singular. Give a computational test for whether a requested terminal state is reachable.

Solution: Compute \( \mathbf{d}_T=\mathbf{x}_T-e^{\mathbf{A}T}\mathbf{x}_0 \) and check whether projection onto the range of \( \mathbf{W}(T) \) leaves the displacement unchanged:

\[ \mathbf{d}_T\in\operatorname{Range}(\mathbf{W}(T)) \quad\Longleftrightarrow\quad \mathbf{W}(T)\mathbf{W}^\dagger(T)\mathbf{d}_T = \mathbf{d}_T. \]

Numerically one uses \( \|\mathbf{W}(T)\mathbf{W}^\dagger(T)\mathbf{d}_T-\mathbf{d}_T\|_2 \le \varepsilon \) for a small tolerance \( \varepsilon \).

Finite-time reachability asks whether an input can drive a system from \( \mathbf{x}_0 \) to \( \mathbf{x}_T \) exactly at a prescribed time. For continuous-time LTI systems, the reachable subspace for any positive horizon is the same span generated by \( \mathbf{B},\mathbf{A}\mathbf{B},\ldots, \mathbf{A}^{n-1}\mathbf{B} \), but the required energy depends strongly on the chosen duration. The minimum-energy formula based on \( \mathbf{W}^\dagger(T) \) connects the abstract reachability question to actuator sizing and feasible maneuver timing.