In Lesson 1, controllability was introduced as the ability to steer the state by choosing an input. We now make the idea precise. Consider the continuous-time LTI state equation \( \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t) \), with \( \mathbf{x}(t)\in\mathbb{R}^n \) and \( \mathbf{u}(t)\in\mathbb{R}^m \). If the system starts from the origin, the forced response at time \( T \) is

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

Therefore, an actuator does not directly create arbitrary states. It creates only those states that can be represented by the convolution of \( e^{\mathbf{A}t}\mathbf{B} \) with some admissible input. The reachable subspace is the geometric object that collects all such states.

flowchart TD
  A["Input signal u(t)"] --> B["Actuator matrix B"]
  B --> C["State dynamics shaped by A"]
  C --> D["Integral response over 0 to T"]
  D --> E["Final reachable state x(T)"]
  E --> F["Reachable subspace"]
        

Fix a terminal time \( T > 0 \). The finite-time reachable set from the origin is

\[ \mathcal{R}(T)=\left\{\mathbf{x}_T\in\mathbb{R}^n: \mathbf{x}_T=\int_0^T e^{\mathbf{A}(T-\tau)}\mathbf{B}\mathbf{u}(\tau)\,d\tau \text{ for some }\mathbf{u}\right\}. \]

The notation emphasizes that reachability is defined relative to a time horizon. If a state lies in \( \mathcal{R}(T) \), then there exists an input history that transfers the origin to that state exactly at time \( T \).

For LTI systems with standard unconstrained inputs, the finite-time set \( \mathcal{R}(T) \) is a linear subspace. To see this, let \( \mathbf{x}_1,\mathbf{x}_2\in\mathcal{R}(T) \) be generated by \( \mathbf{u}_1 \) and \( \mathbf{u}_2 \). For scalars \( \alpha,\beta\in\mathbb{R} \),

\[ \alpha\mathbf{x}_1+\beta\mathbf{x}_2 =\int_0^T e^{\mathbf{A}(T-\tau)}\mathbf{B} \left(\alpha\mathbf{u}_1(\tau)+\beta\mathbf{u}_2(\tau)\right)d\tau. \]

Since \( \alpha\mathbf{u}_1+\beta\mathbf{u}_2 \) is again an admissible input, the linear combination is reachable. Hence \( \mathcal{R}(T) \) is closed under addition and scalar multiplication.

The convolution definition is physically meaningful but not always convenient for computation. For an \( n \)-state LTI system, the same reachable directions are captured by the algebraic subspace

\[ \mathscr{R}=\operatorname{span}\{\operatorname{im}\mathbf{B}, \operatorname{im}\mathbf{AB},\dots, \operatorname{im}\mathbf{A}^{n-1}\mathbf{B}\} =\sum_{k=0}^{n-1}\operatorname{im}(\mathbf{A}^k\mathbf{B}). \]

Equivalently, if \( \mathbf{B}=[\mathbf{b}_1\;\dots\;\mathbf{b}_m] \), then \( \mathscr{R} \) is spanned by all vectors \( \mathbf{A}^k\mathbf{b}_j \) for \( k=0,\dots,n-1 \) and \( j=1,\dots,m \). These are the actuator directions and their images after repeated propagation through the internal dynamics.

The matrix collecting these columns is

\[ \mathbf{\mathcal{C} }= \begin{bmatrix} \mathbf{B} & \mathbf{AB} & \mathbf{A}^2\mathbf{B} & \cdots & \mathbf{A}^{n-1}\mathbf{B} \end{bmatrix}. \]

This lesson uses the matrix only as a compact representation of the reachable subspace. The rank-based controllability test will be treated systematically in Chapter 11.

We prove that the convolution-generated reachable directions coincide with \( \mathscr{R} \). First, by the matrix exponential series,

\[ e^{\mathbf{A}s}\mathbf{B} =\sum_{q=0}^\infty \frac{s^q}{q!}\mathbf{A}^q\mathbf{B}. \]

Powers beyond \( n-1 \) introduce no new subspace directions because the characteristic polynomial \( p_\mathbf{A}(\lambda) \) satisfies the Cayley-Hamilton identity:

\[ p_\mathbf{A}(\lambda)=\lambda^n+a_{n-1}\lambda^{n-1}+\cdots+a_1\lambda+a_0, \qquad p_\mathbf{A}(\mathbf{A})=\mathbf{0}. \]

Hence \( \mathbf{A}^n \), and therefore every higher power \( \mathbf{A}^q \), can be expressed as a linear combination of \( \mathbf{I},\mathbf{A},\dots,\mathbf{A}^{n-1} \). Thus every vector \( e^{\mathbf{A}s}\mathbf{B}\mathbf{v} \) lies in \( \mathscr{R} \). Integrals and linear combinations of such vectors also lie in \( \mathscr{R} \), so \( \mathcal{R}(T)\subseteq\mathscr{R} \).

Conversely, suppose a vector \( \mathbf{z} \) is orthogonal to every finite-time response direction:

\[ \mathbf{z}^T e^{\mathbf{A}s}\mathbf{B}=\mathbf{0}^T \quad \text{for all }s\text{ in an interval}. \]

Differentiating at \( s=0 \) gives

\[ \mathbf{z}^T\mathbf{A}^k\mathbf{B}=\mathbf{0}^T, \qquad k=0,1,\dots,n-1. \]

Therefore the orthogonal complement of the convolution-generated reachable directions is the same as the orthogonal complement of \( \mathscr{R} \). The two subspaces must be equal.

If the initial state is not zero, the solution at time \( T \) 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. \]

Thus the input can only modify the free response by an element of the reachable subspace:

\[ \mathbf{x}_T \text{ is reachable from } \mathbf{x}_0 \quad\Longleftrightarrow\quad \mathbf{x}_T-e^{\mathbf{A}T}\mathbf{x}_0\in\mathscr{R}. \]

This formula separates two mechanisms: the natural motion \( e^{\mathbf{A}T}\mathbf{x}_0 \) produced by internal dynamics, and the forced motion produced by the input channel. In continuous-time LTI systems, the state transition matrix is nonsingular, so reachability and steering concepts are closely related. More detailed distinctions will be developed in later lessons.

Let \( \mathbf{Q} \) contain an orthonormal basis for \( \mathscr{R} \). Then any desired target \( \mathbf{x}_d \) decomposes into

\[ \mathbf{x}_d=\underbrace{\mathbf{Q}\mathbf{Q}^T\mathbf{x}_d}_{\text{reachable component} } +\underbrace{\left(\mathbf{I}-\mathbf{Q}\mathbf{Q}^T\right)\mathbf{x}_d}_{\text{unreachable component} }. \]

The first term can be produced by some input. The second term is orthogonal to every reachable state from the origin, so no input acting through \( \mathbf{B} \) can create it. This decomposition is the geometric meaning of limited actuator authority.

flowchart TD
  X["Desired target state"] --> R["Projection onto reachable subspace"]
  X --> U["Projection onto unreachable complement"]
  R --> A["Can be generated by input"]
  U --> B["Cannot be generated by input channel"]
  A --> C["State steering is possible for reachable part"]
  B --> D["State error remains unless actuator structure changes"]
        

Consider the three-state system

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

The actuator affects the first two states through the second-order subsystem, while the third state is dynamically isolated from the input channel. The algebraic reachable-subspace matrix is

\[ \mathbf{\mathcal{C} }= \begin{bmatrix} \mathbf{B} & \mathbf{AB} & \mathbf{A}^2\mathbf{B} \end{bmatrix} = \begin{bmatrix} 0&1&-3\\ 1&-3&7\\ 0&0&0 \end{bmatrix}. \]

Its columns span the plane \( \{[x_1,x_2,0]^T:x_1,x_2\in\mathbb{R}\} \). Therefore every state with nonzero third component contains an unreachable part. For example, for \( \mathbf{x}_d=[1,-0.5,2]^T \), the reachable component is \( [1,-0.5,0]^T \) and the unreachable component is \( [0,0,2]^T \).

For computation, one may approximate the input as piecewise constant on intervals \( [t_j,t_{j+1}] \). The terminal state is then approximated by

\[ \mathbf{x}_N(T)= \sum_{j=0}^{N-1} \left(\int_{t_j}^{t_{j+1} } e^{\mathbf{A}(T-\tau)}\mathbf{B}\,d\tau\right) \mathbf{u}_j = \mathbf{G}_N\bar{\mathbf{u} }. \]

Solving \( \mathbf{G}_N\bar{\mathbf{u} }\approx\mathbf{x}_d \) by least squares gives the closest terminal state within the numerical reachable subspace. If the desired target has an unreachable component, least squares can match only the reachable projection.

This Python program computes the reachable-subspace basis with SVD, projects a target state, and constructs a finite-time piecewise-constant input map.

# Chapter10_Lesson2.py
# Reachable States and Reachable Subspace for continuous-time LTI systems
#
# Model:
#   x_dot(t) = A x(t) + B u(t)
#   x(0) = 0
#   x(T) = integral_0^T exp(A(T-tau)) B u(tau) d tau
#
# This script computes:
#   1) the algebraic reachable subspace span{B, AB, ..., A^(n-1)B}
#   2) a finite-time numerical input-to-state map using piecewise-constant inputs
#   3) the projection of a desired target state onto the reachable subspace

import numpy as np
from scipy.linalg import expm, svd, null_space
from scipy.integrate import quad_vec


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


def orthonormal_column_basis(M: np.ndarray, tol: float = 1e-10) -> np.ndarray:
    """Return an orthonormal basis for the column space of M using SVD."""
    U, s, _ = svd(M, full_matrices=False)
    r = int(np.sum(s > tol))
    return U[:, :r]


def project_onto_subspace(x: np.ndarray, Q: np.ndarray) -> np.ndarray:
    """Project x onto span(Q), where Q has orthonormal columns."""
    return Q @ (Q.T @ x)


def finite_time_input_map(A: np.ndarray, B: np.ndarray, T: float, N: int) -> np.ndarray:
    """
    Construct the matrix G_N mapping piecewise-constant controls to x(T).

    If u(t) = u_j on interval [t_j, t_{j+1}), then:
        x(T) ≈ sum_j integral_{t_j}^{t_{j+1} } exp(A(T-tau))B d tau * u_j

    G_N has shape n x (mN).
    """
    n, m = B.shape
    grid = np.linspace(0.0, T, N + 1)
    blocks = []
    for j in range(N):
        a, b = grid[j], grid[j + 1]

        def integrand(tau):
            return expm(A * (T - tau)) @ B

        block, _ = quad_vec(integrand, a, b)
        blocks.append(block.reshape(n, m))
    return np.hstack(blocks)


def least_squares_piecewise_control(G: np.ndarray, x_target: np.ndarray) -> np.ndarray:
    """
    Minimum-norm vector of piecewise-constant control samples that best reaches x_target.
    """
    u_star, *_ = np.linalg.lstsq(G, x_target, rcond=None)
    return u_star


def demo():
    # Third state is dynamically decoupled from the actuator, so it is not reachable.
    A = np.array([
        [0.0, 1.0, 0.0],
        [-2.0, -3.0, 0.0],
        [0.0, 0.0, -1.0]
    ])
    B = np.array([
        [0.0],
        [1.0],
        [0.0]
    ])

    R = reachability_matrix(A, B)
    Q = orthonormal_column_basis(R)

    print("A =\n", A)
    print("B =\n", B)
    print("\nReachability matrix [B AB A^2B] =\n", R)
    print("Dimension of reachable subspace:", Q.shape[1])
    print("Orthonormal basis for reachable subspace =\n", Q)

    x_desired = np.array([1.0, -0.5, 2.0])
    x_reachable_part = project_onto_subspace(x_desired, Q)
    x_unreachable_part = x_desired - x_reachable_part

    print("\nDesired target:", x_desired)
    print("Reachable component:", x_reachable_part)
    print("Unreachable component:", x_unreachable_part)

    T = 3.0
    N = 60
    G = finite_time_input_map(A, B, T, N)
    u_piecewise = least_squares_piecewise_control(G, x_desired)
    xT = G @ u_piecewise

    print("\nFinite-time map rank:", np.linalg.matrix_rank(G))
    print("Reached final state using least-squares input:", xT)
    print("Final error:", x_desired - xT)

    # Left annihilator of the reachable subspace: any vector here cannot be influenced by u.
    left_annihilator = null_space(R.T)
    print("\nBasis for orthogonal complement of reachable subspace =\n", left_annihilator)


if __name__ == "__main__":
    demo()

The C++ implementation avoids external libraries and computes the reachable-subspace dimension using Gaussian elimination.

// Chapter10_Lesson2.cpp
// Reachable States and Reachable Subspace for continuous-time LTI systems
//
// This from-scratch C++ example builds the algebraic reachable subspace matrix
//   R = [B, AB, A^2B, ..., A^(n-1)B]
// and estimates independent columns using Gaussian elimination.
//
// Compile:
//   g++ -std=c++17 Chapter10_Lesson2.cpp -o Chapter10_Lesson2

#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

using Matrix = std::vector<std::vector<double>>;

Matrix zeros(int r, int c) {
    return Matrix(r, std::vector<double>(c, 0.0));
}

Matrix identity(int n) {
    Matrix I = zeros(n, n);
    for (int i = 0; i < n; ++i) I[i][i] = 1.0;
    return I;
}

Matrix multiply(const Matrix& A, const Matrix& B) {
    int r = static_cast<int>(A.size());
    int k = static_cast<int>(A[0].size());
    int c = static_cast<int>(B[0].size());
    Matrix C = zeros(r, c);
    for (int i = 0; i < r; ++i) {
        for (int j = 0; j < c; ++j) {
            for (int ell = 0; ell < k; ++ell) {
                C[i][j] += A[i][ell] * B[ell][j];
            }
        }
    }
    return C;
}

Matrix hstack(const std::vector<Matrix>& blocks) {
    int rows = static_cast<int>(blocks[0].size());
    int totalCols = 0;
    for (const auto& M : blocks) totalCols += static_cast<int>(M[0].size());

    Matrix H = zeros(rows, totalCols);
    int offset = 0;
    for (const auto& M : blocks) {
        int cols = static_cast<int>(M[0].size());
        for (int i = 0; i < rows; ++i)
            for (int j = 0; j < cols; ++j)
                H[i][offset + j] = M[i][j];
        offset += cols;
    }
    return H;
}

Matrix reachabilityMatrix(const Matrix& A, const Matrix& B) {
    int n = static_cast<int>(A.size());
    Matrix Ak = identity(n);
    std::vector<Matrix> blocks;
    for (int k = 0; k < n; ++k) {
        blocks.push_back(multiply(Ak, B));
        Ak = multiply(Ak, A);
    }
    return hstack(blocks);
}

int rankGaussian(Matrix M, double tol = 1e-10) {
    int rows = static_cast<int>(M.size());
    int cols = static_cast<int>(M[0].size());
    int rank = 0;

    for (int col = 0; col < cols && rank < rows; ++col) {
        int pivot = rank;
        for (int r = rank + 1; r < rows; ++r) {
            if (std::fabs(M[r][col]) > std::fabs(M[pivot][col])) pivot = r;
        }
        if (std::fabs(M[pivot][col]) <= tol) continue;

        std::swap(M[pivot], M[rank]);
        double pivotValue = M[rank][col];
        for (int j = col; j < cols; ++j) M[rank][j] /= pivotValue;

        for (int r = 0; r < rows; ++r) {
            if (r == rank) continue;
            double factor = M[r][col];
            for (int j = col; j < cols; ++j) {
                M[r][j] -= factor * M[rank][j];
            }
        }
        ++rank;
    }
    return rank;
}

void printMatrix(const Matrix& M, const std::string& name) {
    std::cout << name << " =\n";
    for (const auto& row : M) {
        for (double value : row) {
            std::cout << std::setw(12) << std::setprecision(6) << value << " ";
        }
        std::cout << "\n";
    }
}

int main() {
    Matrix A = {
        {0.0, 1.0, 0.0},
        {-2.0, -3.0, 0.0},
        {0.0, 0.0, -1.0}
    };

    Matrix B = {
        {0.0},
        {1.0},
        {0.0}
    };

    Matrix R = reachabilityMatrix(A, B);

    printMatrix(A, "A");
    printMatrix(B, "B");
    printMatrix(R, "R = [B AB A^2B]");

    int reachableDimension = rankGaussian(R);
    std::cout << "\nEstimated dimension of reachable subspace: "
              << reachableDimension << "\n";

    if (reachableDimension < static_cast<int>(A.size())) {
        std::cout << "Some state directions are not reachable from the actuator.\n";
    } else {
        std::cout << "The reachable subspace is the whole state space.\n";
    }

    return 0;
}

The Java version mirrors the C++ implementation and is suitable for teaching matrix operations from scratch.

// Chapter10_Lesson2.java
// Reachable States and Reachable Subspace for continuous-time LTI systems
//
// This from-scratch Java example builds the algebraic reachable subspace matrix
//   R = [B, AB, A^2B, ..., A^(n-1)B]
// and estimates the dimension of its column space using Gaussian elimination.
//
// Compile and run:
//   javac Chapter10_Lesson2.java
//   java Chapter10_Lesson2

public class Chapter10_Lesson2 {
    static double[][] zeros(int r, int c) {
        return new double[r][c];
    }

    static double[][] identity(int n) {
        double[][] I = zeros(n, n);
        for (int i = 0; i < n; i++) I[i][i] = 1.0;
        return I;
    }

    static double[][] multiply(double[][] A, double[][] B) {
        int r = A.length;
        int k = A[0].length;
        int c = B[0].length;
        double[][] C = zeros(r, c);
        for (int i = 0; i < r; i++) {
            for (int j = 0; j < c; j++) {
                for (int ell = 0; ell < k; ell++) {
                    C[i][j] += A[i][ell] * B[ell][j];
                }
            }
        }
        return C;
    }

    static double[][] hstack(double[][][] blocks) {
        int rows = blocks[0].length;
        int totalCols = 0;
        for (double[][] block : blocks) totalCols += block[0].length;

        double[][] H = zeros(rows, totalCols);
        int offset = 0;
        for (double[][] block : blocks) {
            int cols = block[0].length;
            for (int i = 0; i < rows; i++) {
                for (int j = 0; j < cols; j++) {
                    H[i][offset + j] = block[i][j];
                }
            }
            offset += cols;
        }
        return H;
    }

    static double[][] reachabilityMatrix(double[][] A, double[][] B) {
        int n = A.length;
        double[][] Ak = identity(n);
        double[][][] blocks = new double[n][][];
        for (int k = 0; k < n; k++) {
            blocks[k] = multiply(Ak, B);
            Ak = multiply(Ak, A);
        }
        return hstack(blocks);
    }

    static int rankGaussian(double[][] input, double tol) {
        int rows = input.length;
        int cols = input[0].length;
        double[][] M = new double[rows][cols];

        for (int i = 0; i < rows; i++) {
            System.arraycopy(input[i], 0, M[i], 0, cols);
        }

        int rank = 0;
        for (int col = 0; col < cols && rank < rows; col++) {
            int pivot = rank;
            for (int r = rank + 1; r < rows; r++) {
                if (Math.abs(M[r][col]) > Math.abs(M[pivot][col])) pivot = r;
            }
            if (Math.abs(M[pivot][col]) <= tol) continue;

            double[] temp = M[pivot];
            M[pivot] = M[rank];
            M[rank] = temp;

            double pivotValue = M[rank][col];
            for (int j = col; j < cols; j++) M[rank][j] /= pivotValue;

            for (int r = 0; r < rows; r++) {
                if (r == rank) continue;
                double factor = M[r][col];
                for (int j = col; j < cols; j++) {
                    M[r][j] -= factor * M[rank][j];
                }
            }
            rank++;
        }
        return rank;
    }

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

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

        double[][] B = {
            {0.0},
            {1.0},
            {0.0}
        };

        double[][] R = reachabilityMatrix(A, B);
        printMatrix("A", A);
        printMatrix("B", B);
        printMatrix("R = [B AB A^2B]", R);

        int reachableDimension = rankGaussian(R, 1e-10);
        System.out.println("\nEstimated dimension of reachable subspace: " + reachableDimension);

        if (reachableDimension < A.length) {
            System.out.println("Some state directions are not reachable from the actuator.");
        } else {
            System.out.println("The reachable subspace is the whole state space.");
        }
    }
}

The MATLAB script computes the reachable basis and also creates a simple Simulink model containing a State-Space block.

% Chapter10_Lesson2.m
% Reachable States and Reachable Subspace for continuous-time LTI systems
%
% This script computes:
%   1) R = [B AB ... A^(n-1)B]
%   2) an orthonormal basis for the reachable subspace
%   3) projection of a target state onto that subspace
%   4) a simple state-space simulation
%   5) an optional Simulink model built programmatically

clear; clc; close all;

A = [ 0  1  0;
     -2 -3  0;
      0  0 -1 ];

B = [0; 1; 0];
C = eye(3);
D = zeros(3,1);

n = size(A,1);

% Algebraic reachable subspace.
R = [];
Ak = eye(n);
for k = 0:n-1
    R = [R, Ak*B]; %#ok<AGROW>
    Ak = Ak*A;
end

reachableDimension = rank(R);
Q = orth(R);

disp('Reachability matrix R = [B AB A^2B]:');
disp(R);
disp(['Dimension of reachable subspace = ', num2str(reachableDimension)]);
disp('Orthonormal basis Q for reachable subspace:');
disp(Q);

xDesired = [1; -0.5; 2];
xReachable = Q*(Q'*xDesired);
xUnreachable = xDesired - xReachable;

disp('Desired target:');
disp(xDesired);
disp('Reachable component:');
disp(xReachable);
disp('Unreachable component:');
disp(xUnreachable);

% Continuous-time state-space simulation with a smooth input.
sys = ss(A,B,C,D);
t = linspace(0,8,500);
u = sin(2*t) + 0.4*cos(5*t);
[y,t,x] = lsim(sys,u,t,zeros(n,1));

figure;
plot(t,x,'LineWidth',1.5);
grid on;
xlabel('Time (s)');
ylabel('States');
legend('x1','x2','x3');
title('State response: x3 remains unaffected by the actuator from zero initial condition');

% Optional Simulink model generation.
% This creates a simple State-Space block driven by a Sine Wave block.
modelName = 'Chapter10_Lesson2_Simulink';
if bdIsLoaded(modelName)
    close_system(modelName,0);
end
new_system(modelName);
open_system(modelName);

add_block('simulink/Sources/Sine Wave',[modelName '/Input u']);
add_block('simulink/Continuous/State-Space',[modelName '/LTI System']);
add_block('simulink/Sinks/Scope',[modelName '/Scope']);

set_param([modelName '/LTI System'], ...
    'A','[0 1 0; -2 -3 0; 0 0 -1]', ...
    'B','[0; 1; 0]', ...
    'C','eye(3)', ...
    'D','zeros(3,1)');

add_line(modelName,'Input u/1','LTI System/1');
add_line(modelName,'LTI System/1','Scope/1');

set_param(modelName,'StopTime','8');
save_system(modelName);
disp(['Created Simulink model: ', modelName, '.slx']);

The Mathematica notebook-style script symbolically constructs the reachable-subspace matrix and evaluates a finite-time convolution.

(* Chapter10_Lesson2.nb *)
(* Reachable States and Reachable Subspace for continuous-time LTI systems *)

ClearAll["Global`*"];

A = { {0, 1, 0}, {-2, -3, 0}, {0, 0, -1} };
B = { {0}, {1}, {0} };
n = Length[A];

ReachabilityMatrix[A_, B_] := Module[{Ak = IdentityMatrix[Length[A]], blocks = {} },
  Do[
    AppendTo[blocks, Ak.B];
    Ak = Ak.A;
    ,
    {k, 0, Length[A] - 1}
  ];
  ArrayFlatten[{blocks}]
];

R = ReachabilityMatrix[A, B];

reachableDimension = MatrixRank[R];
reachableBasis = Orthogonalize[Transpose[R]];

Print["A = ", MatrixForm[A]];
Print["B = ", MatrixForm[B]];
Print["R = [B AB A^2B] = ", MatrixForm[R]];
Print["Dimension of reachable subspace = ", reachableDimension];
Print["Basis vectors for reachable subspace = ", reachableBasis];

xDesired = {1, -1/2, 2};
q = Orthogonalize[Transpose[R]];
projectionMatrix = Transpose[q].q;
xReachable = projectionMatrix.xDesired;
xUnreachable = xDesired - xReachable;

Print["Desired target = ", xDesired];
Print["Reachable component = ", xReachable];
Print["Unreachable component = ", xUnreachable];

(* Finite-time reachable-state expression from x(0)=0 *)
xT[T_] := Integrate[MatrixExp[A (T - tau)].B*u[tau], {tau, 0, T},
  Assumptions -> T > 0];

Print["Formal reachable state x(T) = "];
Print[MatrixForm[xT[T]]];

(* Numerical illustration with u(t)=Sin[2t] *)
u[t_] := Sin[2 t];
xNumeric[T_?NumericQ] := NIntegrate[
  Flatten[MatrixExp[A (T - tau)].B]*u[tau],
  {tau, 0, T}
];

TableForm[Table[{T, xNumeric[T]}, {T, 0, 5, 1}],
  TableHeadings -> {None, {"T", "x(T)"} }]

Problem 1 (Subspace Property): Let \( \mathcal{R}(T) \) be the finite-time reachable set from the origin for an LTI system with unconstrained inputs. Prove that \( \mathcal{R}(T) \) is a vector subspace.

Solution: The zero state is reachable by choosing \( \mathbf{u}(t)=\mathbf{0} \). If \( \mathbf{x}_1 \) and \( \mathbf{x}_2 \) are reachable by \( \mathbf{u}_1 \) and \( \mathbf{u}_2 \), then for arbitrary scalars \( \alpha \) and \( \beta \), the input \( \alpha\mathbf{u}_1+\beta\mathbf{u}_2 \) reaches \( \alpha\mathbf{x}_1+\beta\mathbf{x}_2 \) by linearity of the integral. Hence the set is closed under linear combinations.

Problem 2 (Computing a Reachable Subspace): For \( \mathbf{A}=\begin{bmatrix}0&1\\0&0\end{bmatrix} \) and \( \mathbf{B}=\begin{bmatrix}0\\1\end{bmatrix} \), find the reachable subspace.

Solution: Compute \( \mathbf{AB}=\begin{bmatrix}1\\0\end{bmatrix} \). Therefore,

\[ \mathscr{R}= \operatorname{span}\left\{ \begin{bmatrix}0\\1\end{bmatrix}, \begin{bmatrix}1\\0\end{bmatrix} \right\}=\mathbb{R}^2. \]

The input directly affects the second state, and the system dynamics transfer that effect into the first state.

Problem 3 (Unreachable Direction): For the numerical example in Section 7, show that the vector \( \mathbf{z}=[0,0,1]^T \) is orthogonal to all reachable states.

Solution: The columns of \( \mathbf{\mathcal{C} } \) all have zero third component. Hence \( \mathbf{z}^T\mathbf{B}=0 \), \( \mathbf{z}^T\mathbf{AB}=0 \), and \( \mathbf{z}^T\mathbf{A}^2\mathbf{B}=0 \). By the algebraic characterization, every reachable state lies in the plane with zero third component, so \( \mathbf{z}^T\mathbf{x}=0 \) for every \( \mathbf{x}\in\mathscr{R} \).

Problem 4 (Steering from a Nonzero Initial State): Let \( \mathbf{x}(0)=\mathbf{x}_0 \). Derive the condition for reaching \( \mathbf{x}_T \) at time \( T \).

Solution: From the state-transition formula,

\[ \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. \]

The right-hand side is an element of the reachable subspace. Therefore, the target is reachable from the initial condition exactly when the displacement after removing the free response lies in \( \mathscr{R} \).

Reachable states are terminal states generated by the forced response of an LTI system. From the origin, the reachable set is a vector subspace. For an \( n \)-state system, this subspace is spanned by the actuator directions and their propagated images \( \mathbf{B},\mathbf{AB},\dots,\mathbf{A}^{n-1}\mathbf{B} \). The reachable-unreachable decomposition explains which state components can be influenced by the input and which components remain beyond actuator authority.