In the previous lessons, controllability was introduced as the ability to move the state of a continuous-time LTI system from one point to another by choosing a suitable input. For the system \( \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t) \), 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 . \]

The input can only create state increments in directions generated by \(\mathbf{B}\) and by the propagation of those directions through \(\mathbf{A}\). For an \(n\)-state system, the standard finite-dimensional reachability directions are collected in

\[ \mathcal{C}= \begin{bmatrix} \mathbf{B} & \mathbf{A}\mathbf{B} & \mathbf{A}^2\mathbf{B} & \cdots & \mathbf{A}^{n-1}\mathbf{B} \end{bmatrix}, \qquad \mathcal{R}=\operatorname{range}(\mathcal{C}). \]

In Chapter 11 this matrix becomes the formal Kalman controllability test. Here, we use it as an example-building tool: if \(\operatorname{rank}(\mathcal{C})=n\), the examples will be called controllable; if the rank is smaller than \(n\), there are unreachable state directions.

flowchart TD
  A["State equation: xdot = A x + B u"] --> B["Direct input directions: columns of B"]
  B --> C["Propagated directions: A B, A^2 B, ..., A^(n-1) B"]
  C --> D["Reachable subspace = span of all generated directions"]
  D --> E{"Does dimension \nequal n?"}
  E -->|"yes"| F["Controllable example"]
  E -->|"no"| G["Uncontrollable example"]
  G --> H["Find unreachable states or modes"]
  

An example is easier to understand when we can exhibit a direction that no input can affect. Suppose there exists a nonzero row vector \(\mathbf{w}^T\) satisfying

\[ \mathbf{w}^T\mathbf{B}=\mathbf{0},\quad \mathbf{w}^T\mathbf{A}\mathbf{B}=\mathbf{0},\quad \ldots,\quad \mathbf{w}^T\mathbf{A}^{n-1}\mathbf{B}=\mathbf{0}. \]

Then \(\mathbf{w}^T\mathcal{C}=\mathbf{0}\), so every reachable increment is orthogonal to \(\mathbf{w}\). The system cannot independently change the component of the state in that missing direction.

A useful equivalent interpretation uses the power-series expansion of the matrix exponential:

\[ \mathbf{w}^T e^{\mathbf{A}t}\mathbf{B} = \sum_{k=0}^\infty \frac{t^k}{k!}\mathbf{w}^T\mathbf{A}^k\mathbf{B}=0 . \]

By the Cayley-Hamilton theorem, powers \(\mathbf{A}^k\) for \(k\ge n\) are linear combinations of lower powers. Therefore, if the first \(n\) propagated input directions miss \(\mathbf{w}\), then all later directions also miss it.

The double integrator is the simplest model of a point mass with input acceleration. Let \(x_1\) be position and \(x_2\) be velocity:

\[ \dot{x}_1=x_2,\qquad \dot{x}_2=u,\qquad \mathbf{A}=\begin{bmatrix}0&1\\0&0\end{bmatrix},\quad \mathbf{B}=\begin{bmatrix}0\\1\end{bmatrix}. \]

The input acts directly on velocity, but velocity integrates into position through the state matrix. Thus

\[ \mathcal{C}= \begin{bmatrix}\mathbf{B}&\mathbf{A}\mathbf{B}\end{bmatrix} = \begin{bmatrix} 0&1\\ 1&0 \end{bmatrix}, \qquad \det(\mathcal{C})=-1\ne 0 . \]

Hence \(\operatorname{rank}(\mathcal{C})=2\). A force input can first change velocity, and the changed velocity changes position. This is a typical example where the actuator is not attached directly to every state, yet the system is still controllable because dynamics propagate authority.

Consider two independent double integrators, but apply force only to the first mass:

\[ \mathbf{x}= \begin{bmatrix}q_1&v_1&q_2&v_2\end{bmatrix}^T,\quad \mathbf{A}= \begin{bmatrix} 0&1&0&0\\ 0&0&0&0\\ 0&0&0&1\\ 0&0&0&0 \end{bmatrix},\quad \mathbf{B}= \begin{bmatrix} 0\\1\\0\\0 \end{bmatrix}. \]

The controllability matrix is

\[ \mathcal{C}= \begin{bmatrix} 0&1&0&0\\ 1&0&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{bmatrix}, \qquad \operatorname{rank}(\mathcal{C})=2<4 . \]

The first mass is controllable, but the second mass is dynamically disconnected from the actuator. The missing directions are the \(q_2\) and \(v_2\) directions. This example emphasizes that more states do not automatically become reachable unless the state matrix provides coupling paths from the actuator to those states.

Now connect the two masses with a spring. The input still acts only on mass 1, but the coupling spring transmits motion to mass 2:

\[ \mathbf{x}= \begin{bmatrix}q_1&v_1&q_2&v_2\end{bmatrix}^T, \quad \mathbf{B}= \begin{bmatrix}0\\1/m_1\\0\\0\end{bmatrix}, \]

\[ \mathbf{A}= \begin{bmatrix} 0&1&0&0\\ -\frac{k_1+k_c}{m_1}&-\frac{c_1}{m_1}&\frac{k_c}{m_1}&0\\ 0&0&0&1\\ \frac{k_c}{m_2}&0&-\frac{k_2+k_c}{m_2}&-\frac{c_2}{m_2} \end{bmatrix}. \]

The input enters \(v_1\) directly. Then \(v_1\) changes \(q_1\). Through the coupling spring term \(k_c(q_1-q_2)\), changes in \(q_1\) affect \(v_2\), and \(v_2\) affects \(q_2\). Therefore the propagated directions \(\mathbf{B},\mathbf{A}\mathbf{B},\mathbf{A}^2\mathbf{B},\mathbf{A}^3\mathbf{B}\) can span the four-dimensional state space for generic nonzero coupling.

flowchart TD
  U["Input force u"] --> V1["Velocity v1"]
  V1 --> Q1["Position q1"]
  Q1 --> KC["Coupling spring kc"]
  KC --> V2["Velocity v2"]
  V2 --> Q2["Position q2"]
  Q2 --> KC
        

If \(k_c=0\), the diagram path from mass 1 to mass 2 disappears and the example reduces to the uncontrollable decoupled case. Thus controllability is a property of the pair \((\mathbf{A},\mathbf{B})\), not of the input matrix alone.

Suppose the state coordinates are already modal coordinates:

\[ \mathbf{A}=\operatorname{diag}(\lambda_1,\lambda_2,\ldots,\lambda_n), \qquad \mathbf{B}= \begin{bmatrix}b_1&b_2&\cdots&b_n\end{bmatrix}^T . \]

Then the controllability matrix becomes

\[ \mathcal{C}= \begin{bmatrix} b_1&\lambda_1 b_1&\lambda_1^2 b_1&\cdots&\lambda_1^{n-1}b_1\\ b_2&\lambda_2 b_2&\lambda_2^2 b_2&\cdots&\lambda_2^{n-1}b_2\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ b_n&\lambda_n b_n&\lambda_n^2 b_n&\cdots&\lambda_n^{n-1}b_n \end{bmatrix}. \]

Its determinant factors as a diagonal input-participation term times a Vandermonde determinant:

\[ \det(\mathcal{C}) = \left(\prod_{i=1}^n b_i\right) \left(\prod_{1\le i<j\le n}(\lambda_j-\lambda_i)\right). \]

Therefore, if the eigenvalues are distinct and every \(b_i\ne 0\), the system is controllable. If \(b_i=0\) for some mode, that mode does not see the input and is unreachable.

For example,

\[ \mathbf{A}= \begin{bmatrix}-1&0&0\\0&-2&0\\0&0&-3\end{bmatrix}, \quad \mathbf{B}= \begin{bmatrix}1\\0\\1\end{bmatrix} \quad\Longrightarrow\quad \operatorname{rank}(\mathcal{C})=2<3 . \]

The second mode is stable, but it is still uncontrollable. Stability and controllability answer different questions: stability asks what the state does naturally, while controllability asks what the input can command.

Consider the diagonal system

\[ \dot{x}_1=-2x_1+u,\qquad \dot{x}_2=\alpha x_2,\qquad \mathbf{A}= \begin{bmatrix}-2&0\\0&\alpha\end{bmatrix},\quad \mathbf{B}= \begin{bmatrix}1\\0\end{bmatrix}. \]

The second state is never affected by the input. The controllability matrix is

\[ \mathcal{C}= \begin{bmatrix} 1&-2\\ 0&0 \end{bmatrix}, \qquad \operatorname{rank}(\mathcal{C})=1<2 . \]

If \(\alpha<0\), the hidden mode decays naturally. If \(\alpha>0\), the hidden mode grows and no state-feedback input through \(u\) can stabilize it. This distinction will become essential later when we study stabilizability and pole placement.

Multi-input systems may be controllable even if no single actuator is sufficient. Consider

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

The first input controls the double-integrator chain \((x_1,x_2)\), while the second input controls \(x_3\). The first input alone cannot move \(x_3\); the second input alone cannot move \((x_1,x_2)\). Together,

\[ \operatorname{rank} \begin{bmatrix} \mathbf{B}&\mathbf{A}\mathbf{B}&\mathbf{A}^2\mathbf{B} \end{bmatrix} =3 . \]

This example distinguishes actuator redundancy from actuator coverage. Adding more inputs is useful only when the new input columns add missing directions to the reachable subspace.

The following implementations compute \(\mathcal{C}=[\mathbf{B},\mathbf{A}\mathbf{B},\ldots,\mathbf{A}^{n-1}\mathbf{B}]\) and compare its rank with the state dimension. Python uses NumPy and optionally the python-control package. MATLAB uses the Control System Toolbox function ctrb. C++ and Java include from-scratch Gaussian-elimination rank routines.

Chapter10_Lesson5.py

# Chapter10_Lesson5.py
"""
Modern Control - Chapter 10, Lesson 5
Examples of controllable and uncontrollable systems.

Libraries:
    numpy: matrix construction and rank computation.
    scipy.linalg: optional matrix utilities.
    python-control: optional; if installed, control.ctrb(A, B) gives the
    same controllability matrix used here.

Run:
    python Chapter10_Lesson5.py
"""

import numpy as np
from numpy.linalg import matrix_rank


def controllability_matrix(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """Return C = [B, AB, A^2 B, ..., 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 analyze_system(name: str, A: np.ndarray, B: np.ndarray, tol: float = 1e-10) -> None:
    C = controllability_matrix(A, B)
    r = matrix_rank(C, tol=tol)
    n = A.shape[0]
    print(f"\n{name}")
    print("-" * len(name))
    print("A =\n", A)
    print("B =\n", B)
    print("Controllability matrix C =\n", C)
    print(f"rank(C) = {r} out of n = {n}")
    print("Conclusion:", "controllable" if r == n else "uncontrollable")


def main() -> None:
    # Example 1: double integrator, input is acceleration.
    A1 = np.array([[0.0, 1.0],
                   [0.0, 0.0]])
    B1 = np.array([[0.0],
                   [1.0]])
    analyze_system("Example 1: double integrator", A1, B1)

    # Example 2: two decoupled integrator chains, only the first chain is actuated.
    A2 = np.array([[0.0, 1.0, 0.0, 0.0],
                   [0.0, 0.0, 0.0, 0.0],
                   [0.0, 0.0, 0.0, 1.0],
                   [0.0, 0.0, 0.0, 0.0]])
    B2 = np.array([[0.0],
                   [1.0],
                   [0.0],
                   [0.0]])
    analyze_system("Example 2: two decoupled masses, one actuator", A2, B2)

    # Example 3: coupled two-mass system. Input acts on mass 1.
    # State x = [q1, v1, q2, v2]^T.
    m1, m2 = 1.0, 1.0
    k1, k2, kc = 1.0, 1.2, 0.8
    c1, c2 = 0.1, 0.2
    A3 = np.array([[0.0, 1.0, 0.0, 0.0],
                   [-(k1 + kc) / m1, -c1 / m1, kc / m1, 0.0],
                   [0.0, 0.0, 0.0, 1.0],
                   [kc / m2, 0.0, -(k2 + kc) / m2, -c2 / m2]])
    B3 = np.array([[0.0],
                   [1.0 / m1],
                   [0.0],
                   [0.0]])
    analyze_system("Example 3: coupled two-mass oscillator", A3, B3)

    # Example 4: diagonal modal system.
    # A mode is unreachable if its corresponding entry in B is zero.
    A4 = np.diag([-1.0, -2.0, -3.0])
    B4a = np.array([[1.0],
                    [1.0],
                    [1.0]])
    B4b = np.array([[1.0],
                    [0.0],
                    [1.0]])
    analyze_system("Example 4a: diagonal system, all modes actuated", A4, B4a)
    analyze_system("Example 4b: diagonal system, middle mode not actuated", A4, B4b)

    # Example 5: stable uncontrollable hidden mode.
    A5 = np.array([[-2.0, 0.0],
                   [0.0, -0.5]])
    B5 = np.array([[1.0],
                   [0.0]])
    analyze_system("Example 5: stable but uncontrollable hidden mode", A5, B5)

    # Optional: compare with python-control if installed.
    try:
        import control
        C_control = control.ctrb(A3, B3)
        print("\npython-control check for Example 3:")
        print("rank(control.ctrb(A3, B3)) =", matrix_rank(C_control))
    except ImportError:
        print("\nOptional package 'control' is not installed; custom implementation was used.")


if __name__ == "__main__":
    main()

Chapter10_Lesson5.cpp

// Chapter10_Lesson5.cpp
// Modern Control - Chapter 10, Lesson 5
// Controllability examples using a small from-scratch matrix implementation.
// Compile:
//     g++ -std=c++17 Chapter10_Lesson5.cpp -o Chapter10_Lesson5
// Run:
//     ./Chapter10_Lesson5

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

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

Matrix zeros(int rows, int cols) {
    return Matrix(rows, std::vector<double>(cols, 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 m = 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 k = 0; k < m; ++k)
            for (int j = 0; j < c; ++j)
                C[i][j] += A[i][k] * B[k][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;
}

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[rank], M[pivot]);
        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;
}

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

void printMatrix(const Matrix& M) {
    for (const auto& row : M) {
        for (double value : row) {
            std::cout << std::setw(10) << std::setprecision(4) << value << " ";
        }
        std::cout << "\n";
    }
}

void analyze(const std::string& name, const Matrix& A, const Matrix& B) {
    Matrix C = controllabilityMatrix(A, B);
    int r = rankGaussian(C);
    int n = static_cast<int>(A.size());

    std::cout << "\n" << name << "\n";
    std::cout << std::string(name.size(), '-') << "\n";
    std::cout << "Controllability matrix:\n";
    printMatrix(C);
    std::cout << "rank(C) = " << r << " out of n = " << n << "\n";
    std::cout << "Conclusion: " << (r == n ? "controllable" : "uncontrollable") << "\n";
}

int main() {
    Matrix A1 = { {0, 1},
                 {0, 0} };
    Matrix B1 = { {0},
                 {1} };
    analyze("Example 1: double integrator", A1, B1);

    Matrix A2 = { {0, 1, 0, 0},
                 {0, 0, 0, 0},
                 {0, 0, 0, 1},
                 {0, 0, 0, 0} };
    Matrix B2 = { {0},
                 {1},
                 {0},
                 {0} };
    analyze("Example 2: two decoupled masses, one actuator", A2, B2);

    double m1 = 1.0, m2 = 1.0, k1 = 1.0, k2 = 1.2, kc = 0.8, c1 = 0.1, c2 = 0.2;
    Matrix A3 = { {0, 1, 0, 0},
                 {-(k1 + kc) / m1, -c1 / m1, kc / m1, 0},
                 {0, 0, 0, 1},
                 {kc / m2, 0, -(k2 + kc) / m2, -c2 / m2} };
    Matrix B3 = { {0},
                 {1.0 / m1},
                 {0},
                 {0} };
    analyze("Example 3: coupled two-mass oscillator", A3, B3);

    Matrix A4 = { {-1, 0, 0},
                 {0, -2, 0},
                 {0, 0, -3} };
    Matrix B4 = { {1},
                 {0},
                 {1} };
    analyze("Example 4: diagonal system with one unactuated mode", A4, B4);

    return 0;
}

Chapter10_Lesson5.java

// Chapter10_Lesson5.java
// Modern Control - Chapter 10, Lesson 5
// Controllability examples using from-scratch matrix operations.
// Compile:
//     javac Chapter10_Lesson5.java
// Run:
//     java Chapter10_Lesson5

public class Chapter10_Lesson5 {
    static double[][] zeros(int rows, int cols) {
        return new double[rows][cols];
    }

    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 m = A[0].length;
        int c = B[0].length;
        double[][] C = zeros(r, c);
        for (int i = 0; i < r; i++) {
            for (int k = 0; k < m; k++) {
                for (int j = 0; j < c; j++) {
                    C[i][j] += A[i][k] * B[k][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 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[rank];
            M[rank] = M[pivot];
            M[pivot] = 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 double[][] controllabilityMatrix(double[][] A, double[][] B) {
        int n = A.length;
        double[][][] blocks = new double[n][][];
        double[][] Ak = identity(n);
        for (int k = 0; k < n; k++) {
            blocks[k] = multiply(Ak, B);
            Ak = multiply(Ak, A);
        }
        return hstack(blocks);
    }

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

    static void analyze(String name, double[][] A, double[][] B) {
        double[][] C = controllabilityMatrix(A, B);
        int r = rankGaussian(C, 1e-10);
        int n = A.length;

        System.out.println("\n" + name);
        System.out.println("-".repeat(name.length()));
        System.out.println("Controllability matrix:");
        printMatrix(C);
        System.out.println("rank(C) = " + r + " out of n = " + n);
        System.out.println("Conclusion: " + (r == n ? "controllable" : "uncontrollable"));
    }

    public static void main(String[] args) {
        double[][] A1 = { {0, 1},
                         {0, 0} };
        double[][] B1 = { {0},
                         {1} };
        analyze("Example 1: double integrator", A1, B1);

        double[][] A2 = { {0, 1, 0, 0},
                         {0, 0, 0, 0},
                         {0, 0, 0, 1},
                         {0, 0, 0, 0} };
        double[][] B2 = { {0},
                         {1},
                         {0},
                         {0} };
        analyze("Example 2: two decoupled masses, one actuator", A2, B2);

        double m1 = 1.0, m2 = 1.0, k1 = 1.0, k2 = 1.2, kc = 0.8, c1 = 0.1, c2 = 0.2;
        double[][] A3 = { {0, 1, 0, 0},
                         {-(k1 + kc) / m1, -c1 / m1, kc / m1, 0},
                         {0, 0, 0, 1},
                         {kc / m2, 0, -(k2 + kc) / m2, -c2 / m2} };
        double[][] B3 = { {0},
                         {1.0 / m1},
                         {0},
                         {0} };
        analyze("Example 3: coupled two-mass oscillator", A3, B3);

        double[][] A4 = { {-1, 0, 0},
                         {0, -2, 0},
                         {0, 0, -3} };
        double[][] B4 = { {1},
                         {0},
                         {1} };
        analyze("Example 4: diagonal system with one unactuated mode", A4, B4);
    }
}

Chapter10_Lesson5.m

% Chapter10_Lesson5.m
% Modern Control - Chapter 10, Lesson 5
% Examples of controllable and uncontrollable systems.
%
% MATLAB libraries:
%   Control System Toolbox: ctrb, ss, rank, eig.
%   Simulink connection: the ss(A,B,C,D) objects below can be used with
%   State-Space blocks or lsim for time-domain checks.

clear; clc;

function analyzeSystem(name, A, B)
    n = size(A, 1);
    Ctr = ctrb(A, B);
    r = rank(Ctr);
    fprintf('\n%s\n', name);
    fprintf('%s\n', repmat('-', 1, strlength(name)));
    disp('Controllability matrix:');
    disp(Ctr);
    fprintf('rank(C) = %d out of n = %d\n', r, n);
    if r == n
        disp('Conclusion: controllable');
    else
        disp('Conclusion: uncontrollable');
    end
end

% Example 1: double integrator.
A1 = [0 1;
      0 0];
B1 = [0; 1];
analyzeSystem("Example 1: double integrator", A1, B1);

% Example 2: two decoupled masses, one actuator.
A2 = [0 1 0 0;
      0 0 0 0;
      0 0 0 1;
      0 0 0 0];
B2 = [0; 1; 0; 0];
analyzeSystem("Example 2: two decoupled masses, one actuator", A2, B2);

% Example 3: coupled two-mass oscillator.
m1 = 1.0; m2 = 1.0;
k1 = 1.0; k2 = 1.2; kc = 0.8;
c1 = 0.1; c2 = 0.2;

A3 = [0, 1, 0, 0;
     -(k1+kc)/m1, -c1/m1, kc/m1, 0;
      0, 0, 0, 1;
      kc/m2, 0, -(k2+kc)/m2, -c2/m2];
B3 = [0; 1/m1; 0; 0];
analyzeSystem("Example 3: coupled two-mass oscillator", A3, B3);

% Example 4: diagonal modal systems.
A4 = diag([-1, -2, -3]);
B4a = [1; 1; 1];
B4b = [1; 0; 1];

analyzeSystem("Example 4a: diagonal system, all modes actuated", A4, B4a);
analyzeSystem("Example 4b: diagonal system, middle mode not actuated", A4, B4b);

% State-space model object for simulation or Simulink State-Space block data.
C3 = eye(4);
D3 = zeros(4, 1);
sys3 = ss(A3, B3, C3, D3);

% Simple time response for Example 3 with a sinusoidal input.
t = linspace(0, 15, 1000);
u = sin(t);
[y, tOut, x] = lsim(sys3, u, t);

figure;
plot(tOut, x);
grid on;
xlabel('time (s)');
ylabel('states');
title('Chapter10 Lesson5: State response of coupled two-mass system');
legend('q1', 'v1', 'q2', 'v2');

Chapter10_Lesson5.nb

(* Chapter10_Lesson5.nb *)
(* Modern Control - Chapter 10, Lesson 5 *)
(* Examples of controllable and uncontrollable systems in Wolfram Mathematica. *)

ClearAll[controllabilityMatrix, analyzeSystem];

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

analyzeSystem[name_String, A_, B_] := Module[
  {Ctr, r, n},
  Ctr = controllabilityMatrix[A, B];
  r = MatrixRank[Ctr];
  n = Length[A];
  Print["\n", name];
  Print[StringRepeat["-", StringLength[name]]];
  Print["Controllability matrix:"];
  Print[MatrixForm[Ctr]];
  Print["rank(C) = ", r, " out of n = ", n];
  If[r == n, Print["Conclusion: controllable"], Print["Conclusion: uncontrollable"]];
];

(* Example 1: double integrator. *)
A1 = { {0, 1}, {0, 0} };
B1 = { {0}, {1} };
analyzeSystem["Example 1: double integrator", A1, B1];

(* Example 2: two decoupled masses, one actuator. *)
A2 = { {0, 1, 0, 0},
      {0, 0, 0, 0},
      {0, 0, 0, 1},
      {0, 0, 0, 0} };
B2 = { {0}, {1}, {0}, {0} };
analyzeSystem["Example 2: two decoupled masses, one actuator", A2, B2];

(* Example 3: coupled two-mass oscillator. *)
m1 = 1.0; m2 = 1.0;
k1 = 1.0; k2 = 1.2; kc = 0.8;
c1 = 0.1; c2 = 0.2;

A3 = { {0, 1, 0, 0},
      {-(k1 + kc)/m1, -c1/m1, kc/m1, 0},
      {0, 0, 0, 1},
      {kc/m2, 0, -(k2 + kc)/m2, -c2/m2} };
B3 = { {0}, {1/m1}, {0}, {0} };
analyzeSystem["Example 3: coupled two-mass oscillator", A3, B3];

(* Example 4: diagonal modal systems. *)
A4 = DiagonalMatrix[{-1, -2, -3}];
B4a = { {1}, {1}, {1} };
B4b = { {1}, {0}, {1} };

analyzeSystem["Example 4a: diagonal system, all modes actuated", A4, B4a];
analyzeSystem["Example 4b: diagonal system, middle mode not actuated", A4, B4b];

(* Symbolic diagonal modal determinant: if lambdas are distinct and all b_i are nonzero,
   the controllability matrix has full rank. *)
A3sym = DiagonalMatrix[{\[Lambda]1, \[Lambda]2, \[Lambda]3}];
B3sym = { {b1}, {b2}, {b3} };
Ctr3sym = controllabilityMatrix[A3sym, B3sym];
FullSimplify[Det[Ctr3sym]]

Most elementary examples fall into one of the following patterns:

1. Integrator-chain propagation: the input enters a derivative state, and integration propagates authority to lower-order states.
2. Disconnected subsystem: a block of states has no path from the input, so that block is unreachable.
3. Dynamic coupling: springs, inertias, electrical interconnections, or fluid coupling transmit input authority to states that are not directly actuated.
4. Modal participation: in modal coordinates, a mode is controllable only if the input has nonzero projection onto that mode.
5. Multi-input coverage: several weak actuators may collectively span all state directions.

These patterns provide the engineering intuition behind the formal tests developed in the next chapter.

Problem 1 (Double Integrator): For \(\dot{x}_1=x_2,\;\dot{x}_2=u\), compute the controllability matrix and decide whether the system is controllable.

Solution:

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

\[ \mathcal{C}=\begin{bmatrix}0&1\\1&0\end{bmatrix}, \qquad \operatorname{rank}(\mathcal{C})=2 . \]

Since the rank equals the state dimension, the system is controllable.

Problem 2 (Disconnected Fourth-Order System): Consider

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

Find the reachable subspace.

Solution:

\[ \mathbf{B}= \begin{bmatrix}0\\1\\0\\0\end{bmatrix}, \quad \mathbf{A}\mathbf{B}= \begin{bmatrix}1\\0\\0\\0\end{bmatrix}, \quad \mathbf{A}^2\mathbf{B}=\mathbf{0}, \quad \mathbf{A}^3\mathbf{B}=\mathbf{0}. \]

\[ \mathcal{R}= \operatorname{span} \left\{ \begin{bmatrix}0\\1\\0\\0\end{bmatrix}, \begin{bmatrix}1\\0\\0\\0\end{bmatrix} \right\} = \left\{\begin{bmatrix}a\\b\\0\\0\end{bmatrix}:a,b\in\mathbb{R}\right\}. \]

The coordinates \(q_2\) and \(v_2\) are unreachable.

Problem 3 (Diagonal Modal System): Let \(\mathbf{A}=\operatorname{diag}(-1,-2,-3)\) and \(\mathbf{B}=\begin{bmatrix}1&1&1\end{bmatrix}^T\). Is the system controllable?

Solution:

\[ \mathcal{C}= \begin{bmatrix} 1&-1&1\\ 1&-2&4\\ 1&-3&9 \end{bmatrix}, \qquad \det(\mathcal{C})=-2\ne 0 . \]

The rank is \(3\), so all three modes are controllable.

Problem 4 (Missing Modal Input): Repeat Problem 3 with \(\mathbf{B}=\begin{bmatrix}1&0&1\end{bmatrix}^T\).

Solution:

\[ \mathcal{C}= \begin{bmatrix} 1&-1&1\\ 0&0&0\\ 1&-3&9 \end{bmatrix}, \qquad \operatorname{rank}(\mathcal{C})=2<3 . \]

The second modal state is unreachable because the second entry of \(\mathbf{B}\) is zero.

Problem 5 (Left-Null Certificate): For the system in Problem 4, find a nonzero vector \(\mathbf{w}\) such that \(\mathbf{w}^T\mathcal{C}=\mathbf{0}\).

Solution:

\[ \mathbf{w}= \begin{bmatrix}0\\1\\0\end{bmatrix} \quad\Longrightarrow\quad \mathbf{w}^T\mathcal{C} = \begin{bmatrix}0&1&0\end{bmatrix} \begin{bmatrix} 1&-1&1\\ 0&0&0\\ 1&-3&9 \end{bmatrix} = \begin{bmatrix}0&0&0\end{bmatrix}. \]

This vector selects the second modal coordinate, proving that the input cannot affect that coordinate.

Controllability examples reveal how input authority moves through a system. A double integrator is controllable because velocity propagates into position. Decoupled unactuated subsystems are uncontrollable because no input path reaches them. Coupling can restore controllability by transmitting actuator authority across states. In modal coordinates, a mode is reachable only when the input has nonzero participation in that mode. These examples prepare the formal rank, PBH, canonical-form, and graph-based tests in the next chapter.