In the previous lessons of this chapter, state-feedback limitations were connected to uncontrollable modes, transmission zeros, effort/sensitivity trade-offs, and model uncertainty. Here the limitation is structural. A full-state feedback law \( u=-Kx+r \) assumes that every actuator channel can use every state variable and that all intended actuator channels exist. In many physical systems this assumption is false. The plant may have only a small number of force/torque inputs, and the controller may have a prescribed communication graph.

\[ \dot{x}=Ax+Bu,\qquad u=-Kx+r,\qquad A_{cl}=A-BK. \]

If \( B \) has too few effective columns, some directions of the state space may be difficult or impossible to move. If \( K \) is sparse, then even a controllable pair \( (A,B) \) may not admit the desired closed-loop eigenvalue assignment using the allowed feedback links.

flowchart TD
  P["Physical plant: x_dot = A x + B u"] --> A1["Available actuator locations"]
  A1 --> B1["Selected input matrix B_S"]
  B1 --> C1["Controllability and actuator authority"]
  C1 --> D1["Allowed feedback links M"]
  D1 --> E1["Sparse gain K with zero entries fixed"]
  E1 --> F1["Closed-loop matrix A - B_S K"]
  F1 --> G1["Check poles, effort, robustness, and feasibility"]
        

Suppose the designer has a library of candidate actuator directions collected in \( B_f\in\mathbb{R}^{n\times q} \). A chosen actuator subset is represented by a selection matrix \( S_a\in\mathbb{R}^{q\times m} \), so that

\[ B_S=B_fS_a,\qquad \dot{x}=Ax+B_Su. \]

The selected actuator set is sufficient for arbitrary pole assignment only if the selected pair \( (A,B_S) \) is controllable. Using the Kalman controllability matrix,

\[ \mathcal{C}(A,B_S)=\begin{bmatrix}B_S & AB_S & A^2B_S & \cdots & A^{n-1}B_S\end{bmatrix},\qquad \operatorname{rank}\mathcal{C}(A,B_S)=n. \]

Equivalently, the PBH condition requires \( \operatorname{rank}\begin{bmatrix}\lambda I-A & B_S\end{bmatrix}=n \) for every eigenvalue \( \lambda\in\sigma(A) \). If this condition fails for an unstable mode, no state-feedback gain can move that unstable mode to the stable half-plane.

\[ \exists q^*\neq 0:\quad q^*A=\lambda q^*,\quad q^*B_S=0 \quad\Longrightarrow\quad q^*(A-B_SK)=\lambda q^*. \]

The formula above is one of the most important structural obstructions: if a left eigenvector is orthogonal to every selected actuator direction, then that mode is invisible to the input channel. Changing \( K \) cannot affect it because \( q^*B_SK=0 \).

Limited actuators restrict \( B \). Sparse feedback restricts \( K \). Let \( M\in\{0,1\}^{m\times n} \) be a structural mask. The admissible feedback set is

\[ \mathcal{S}_M=\left\{K\in\mathbb{R}^{m\times n}: K_{ij}=0\;\text{whenever}\;M_{ij}=0\right\}. \]

Thus the design problem is no longer unconstrained pole placement. It becomes

\[ \text{find } K\in\mathcal{S}_M \quad\text{such that}\quad \sigma(A-BK)=\Lambda_d, \]

where \( \Lambda_d \) denotes the desired multiset of closed-loop eigenvalues. Even when \( (A,B) \) is controllable, this sparse assignment problem can fail because the map \( K\mapsto \det(sI-A+BK) \) is restricted to a lower-dimensional family.

\[ d=\|M\|_0,\qquad F_M:\mathbb{R}^d\longrightarrow \mathbb{R}^n,\qquad \theta\longmapsto \text{coefficients of }\det(sI-A+BK(\theta)). \]

If \( d < n \), then arbitrary assignment of all \( n \) polynomial coefficients is generically impossible. If \( d\ge n \), success is still not guaranteed because the map is nonlinear and may be singular at the available structure.

Sparse feedback also restricts which eigenvalues can be moved strongly. Let \( \lambda \) be a simple eigenvalue of \( A_{cl}=A-BK \) with right and left eigenvectors \( v \) and \( w \), normalized so that \( w^*v=1 \). A small admissible change \( \Delta K\in\mathcal{S}_M \) produces

\[ \Delta\lambda =w^*\Delta A_{cl}v+O(\|\Delta A_{cl}\|^2) =-w^*B\Delta K v+O(\|B\Delta K\|^2). \]

Therefore, if the allowed nonzero entries of \( \Delta K \) make the bilinear term \( w^*B\Delta K v \) small or zero, then the corresponding closed-loop pole has poor structural mobility. This is a precise version of the engineering statement that an actuator cannot stabilize or shape a mode that it cannot effectively influence through the permitted feedback paths.

\[ \mathcal{D}_\lambda(M)=\left\{-w^*B\Delta K v: \Delta K\in\mathcal{S}_M,\;\|\Delta K\|_F=1\right\}. \]

The set \( \mathcal{D}_\lambda(M) \) describes the first-order directions in which the eigenvalue can move under the sparse structure.

Controllability rank is a yes/no test. It does not say whether a state is easy or difficult to reach. For a stable open-loop matrix \( A \), the infinite-horizon controllability Gramian associated with the selected actuator matrix \( B_S \) satisfies

\[ AW_c+W_cA^T+B_SB_S^T=0. \]

The minimum input energy required to reach \( x_f \) from the origin in the corresponding infinite-horizon stable setting is proportional to

\[ E_{\min}(x_f)=x_f^TW_c^{-1}x_f,\qquad W_c\succ 0. \]

Small eigenvalues of \( W_c \) identify state directions that require large input energy. This gives a quantitative method for comparing actuator locations. A common selection objective is

\[ \max_{S_a}\;\log\det(W_c(S_a)+\varepsilon I) \quad\text{subject to}\quad |S_a|=m, \]

where \( \varepsilon I \) regularizes the determinant when the Gramian is nearly singular. This criterion prefers actuator sets that spread authority broadly across state space.

flowchart TD
  S["Candidate actuator columns"] --> R["Rank test: controllable or not"]
  R --> G["Authority test: Gramian size and conditioning"]
  G --> M["Sparse feedback mask"]
  M --> O["Constrained gain search"]
  O --> E["Closed-loop eigenvalue and simulation check"]
  E --> V["Accept design or change actuators / feedback links"]
        

A structurally constrained feedback design is often written as a constrained optimization problem. One pole-oriented formulation is

\[ \min_{K\in\mathcal{S}_M}\; \sum_{i=1}^n |\lambda_i(A-BK)-\lambda_i^d|^2 +\rho\|K\|_F^2. \]

This objective penalizes deviation from the desired poles and also discourages excessive gain. However, eigenvalues are nonlinear functions of \( K \), so the optimization is generally nonconvex. A polynomial-coefficient formulation avoids eigenvalue ordering:

\[ \min_{K\in\mathcal{S}_M}\; \left\|\operatorname{coef}\left(\det(sI-A+BK)\right) -\operatorname{coef}(p_d(s))\right\|_2^2 +\rho\|K\|_F^2. \]

If actuator selection is also part of the design, the problem becomes combinatorial:

\[ \min_{S_a,K}\;J(S_a,K) \quad\text{subject to}\quad B_S=B_fS_a,\quad K\in\mathcal{S}_M,\quad |S_a|=m. \]

In practice, engineers often combine discrete search over actuator subsets with continuous optimization over the allowed entries of \( K \). For small systems, exhaustive subset search is possible. For large systems, greedy Gramian-based selection or graph-structured heuristics are commonly used before numerical gain tuning.

Proposition 1 (Uncontrollable mode invariance): If \( q^*A=\lambda q^* \) and \( q^*B=0 \), then \( \lambda \) is an eigenvalue of \( A-BK \) for every compatible feedback matrix \( K \).

Proof:

\[ q^*(A-BK)=q^*A-q^*BK=\lambda q^*-0\cdot K=\lambda q^*. \]

Hence the same left eigenvector remains a left eigenvector of the closed-loop matrix. Therefore the mode is unassignable by state feedback.

Proposition 2 (Sparse feedback is an affine structural subspace): The admissible set \( \mathcal{S}_M \) is a linear subspace of \( \mathbb{R}^{m\times n} \) with dimension \( \|M\|_0 \).

Proof: If \( K_1,K_2\in\mathcal{S}_M \) and \( a,b\in\mathbb{R} \), then every forbidden entry satisfies \( (aK_1+bK_2)_{ij}=a(K_1)_{ij}+b(K_2)_{ij}=0 \). Therefore \( aK_1+bK_2\in\mathcal{S}_M \). The free coordinates are exactly the entries where \( M_{ij}=1 \), so the dimension is \( \|M\|_0 \).

Proposition 3 (Parameter-count obstruction): If the sparse mask has \( d=\|M\|_0 < n \), then a generic \( n \)-coefficient closed-loop characteristic polynomial cannot be assigned arbitrarily.

Proof idea: The characteristic polynomial coefficients are generated by the smooth map \( F_M:\mathbb{R}^d\longrightarrow\mathbb{R}^n \). The image of a smooth map from a lower-dimensional parameter space cannot contain an open subset of \( \mathbb{R}^n \) under regularity assumptions. Thus arbitrary assignment is generically impossible when the number of free feedback parameters is smaller than the number of polynomial coefficients.

The Python implementation uses numpy and scipy. The optional python-control package can also be used for ctrb, place, and state-space simulation, but the code below implements the controllability matrix explicitly.

Chapter25_Lesson5.py

# Chapter25_Lesson5.py
"""
Structural Constraints: Limited Actuators and Sparse Feedback

This script demonstrates three ideas from Chapter 25, Lesson 5:
1. actuator selection changes the controllability matrix;
2. a sparse feedback mask restricts the admissible K matrix;
3. a constrained feedback gain may stabilize the system, but it generally
   cannot reproduce an arbitrary dense pole-placement design.

Required libraries:
    pip install numpy scipy

Optional modern-control libraries to explore further:
    pip install control slycot
"""

import itertools
import numpy as np
from scipy.linalg import solve_continuous_are
from scipy.integrate import solve_ivp


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


def matrix_rank(M: np.ndarray, tol: float = 1e-9) -> int:
    """Numerical rank through singular values."""
    s = np.linalg.svd(M, compute_uv=False)
    return int(np.sum(s > tol))


def gramian_energy_score(A: np.ndarray, B: np.ndarray) -> float:
    """
    Use a finite-horizon Gramian-like proxy by numerical quadrature.
    This works even when open-loop A is not asymptotically stable.
    Larger trace means more directions are easier to actuate.
    """
    n = A.shape[0]
    T = 5.0
    steps = 400
    dt = T / steps
    W = np.zeros((n, n))
    Phi = np.eye(n)
    # crude Euler approximation of exp(A t) for educational transparency
    for _ in range(steps):
        W += Phi @ B @ B.T @ Phi.T * dt
        Phi = Phi + dt * A @ Phi
    return float(np.trace(W))


def lqr_gain(A: np.ndarray, B: np.ndarray, Q: np.ndarray, R: np.ndarray) -> np.ndarray:
    """Continuous-time LQR gain K = R^(-1) B' P."""
    P = solve_continuous_are(A, B, Q, R)
    return np.linalg.solve(R, B.T @ P)


def apply_feedback_mask(K_dense: np.ndarray, mask: np.ndarray) -> np.ndarray:
    """Project a dense gain to a sparse structural pattern."""
    return K_dense * mask


def simulate_closed_loop(A: np.ndarray, B: np.ndarray, K: np.ndarray, x0: np.ndarray):
    """Simulate x_dot = (A - B K)x."""
    Acl = A - B @ K

    def rhs(t, x):
        return Acl @ x

    return solve_ivp(rhs, (0.0, 10.0), x0, dense_output=True, max_step=0.02)


def main() -> None:
    np.set_printoptions(precision=4, suppress=True)

    # Four-state coupled oscillator-like model.
    A = np.array([
        [0.0,  1.0,  0.0,  0.0],
        [-2.0, -0.25, 0.7,  0.0],
        [0.0,  0.0,  0.0,  1.0],
        [0.6,  0.0, -1.5, -0.20],
    ])

    # Candidate physical actuator columns.
    B_candidates = np.array([
        [0.0, 0.0, 0.0],
        [1.0, 0.2, 0.0],
        [0.0, 0.0, 0.0],
        [0.0, 0.4, 1.0],
    ])

    n = A.shape[0]
    q = B_candidates.shape[1]

    print("Actuator subset analysis")
    print("------------------------")
    for r in range(1, q + 1):
        for subset in itertools.combinations(range(q), r):
            B = B_candidates[:, subset]
            C = controllability_matrix(A, B)
            rank_C = matrix_rank(C)
            score = gramian_energy_score(A, B)
            print(f"subset={subset}, rank(C)={rank_C}/{n}, trace(W_T)≈{score:.3f}")

    # Select two actuators.
    selected = (0, 2)
    B = B_candidates[:, selected]

    # Dense LQR design.
    Q = np.diag([20.0, 1.0, 20.0, 1.0])
    R = np.diag([0.5, 0.5])
    K_dense = lqr_gain(A, B, Q, R)

    # Sparse communication/measurement pattern:
    # actuator 1 sees states x1 and x2; actuator 2 sees states x3 and x4.
    mask = np.array([
        [1.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 1.0, 1.0],
    ])
    K_sparse = apply_feedback_mask(K_dense, mask)

    print("\nDense K:")
    print(K_dense)
    print("Sparse projected K:")
    print(K_sparse)

    eig_dense = np.linalg.eigvals(A - B @ K_dense)
    eig_sparse = np.linalg.eigvals(A - B @ K_sparse)
    eig_open = np.linalg.eigvals(A)

    print("\nOpen-loop eigenvalues:", eig_open)
    print("Dense closed-loop eigenvalues:", eig_dense)
    print("Sparse closed-loop eigenvalues:", eig_sparse)

    x0 = np.array([1.0, 0.0, -0.7, 0.2])
    sol_dense = simulate_closed_loop(A, B, K_dense, x0)
    sol_sparse = simulate_closed_loop(A, B, K_sparse, x0)

    print("\nFinal dense state:", sol_dense.y[:, -1])
    print("Final sparse state:", sol_sparse.y[:, -1])
    print("\nInterpretation:")
    print("A sparse mask preserves only allowed feedback links. Even when the pair")
    print("(A,B) is controllable, the masked gain may have weaker pole motion,")
    print("higher residual oscillation, or even fail to stabilize in harder examples.")


if __name__ == "__main__":
    main()

The C++ implementation uses the Eigen library, which is a common choice for control-oriented matrix computation in C++ projects.

Chapter25_Lesson5.cpp

// Chapter25_Lesson5.cpp
/*
Structural Constraints: Limited Actuators and Sparse Feedback

C++ implementation using Eigen:
    g++ -std=c++17 Chapter25_Lesson5.cpp -I /path/to/eigen -O2 -o Chapter25_Lesson5

Eigen is widely used in modern-control prototyping for matrix operations.
This example computes controllability rank and compares dense/sparse gains.
*/

#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <iostream>
#include <vector>

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

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

int numericalRank(const MatrixXd& M, double tol = 1e-9) {
    Eigen::JacobiSVD<MatrixXd> svd(M);
    const auto& s = svd.singularValues();
    int r = 0;
    for (int i = 0; i < s.size(); ++i) {
        if (s(i) > tol) {
            ++r;
        }
    }
    return r;
}

MatrixXd applyMask(const MatrixXd& K, const MatrixXd& mask) {
    return K.cwiseProduct(mask);
}

VectorXd eulerSimulate(const MatrixXd& Acl, const VectorXd& x0, double dt, int steps) {
    VectorXd x = x0;
    for (int k = 0; k < steps; ++k) {
        x = x + dt * (Acl * x);
    }
    return x;
}

void printEigenvalues(const MatrixXd& Acl, const std::string& label) {
    Eigen::EigenSolver<MatrixXd> solver(Acl);
    std::cout << label << "\n";
    for (int i = 0; i < solver.eigenvalues().size(); ++i) {
        std::cout << "  " << solver.eigenvalues()(i) << "\n";
    }
}

int main() {
    MatrixXd A(4, 4);
    A << 0.0,  1.0,  0.0,  0.0,
        -2.0, -0.25, 0.7,  0.0,
         0.0,  0.0,  0.0,  1.0,
         0.6,  0.0, -1.5, -0.20;

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

    MatrixXd C = controllabilityMatrix(A, B);
    std::cout << "rank(C) = " << numericalRank(C) << " / " << A.rows() << "\n";

    // A representative dense stabilizing gain from an external design step.
    MatrixXd Kdense(2, 4);
    Kdense << 3.0, 2.0, 0.9, 0.4,
              0.6, 0.3, 2.6, 1.8;

    MatrixXd mask(2, 4);
    mask << 1.0, 1.0, 0.0, 0.0,
            0.0, 0.0, 1.0, 1.0;

    MatrixXd Ksparse = applyMask(Kdense, mask);

    std::cout << "\nDense K:\n" << Kdense << "\n";
    std::cout << "\nSparse K:\n" << Ksparse << "\n";

    MatrixXd AclDense = A - B * Kdense;
    MatrixXd AclSparse = A - B * Ksparse;

    printEigenvalues(A, "\nOpen-loop eigenvalues:");
    printEigenvalues(AclDense, "\nDense closed-loop eigenvalues:");
    printEigenvalues(AclSparse, "\nSparse closed-loop eigenvalues:");

    VectorXd x0(4);
    x0 << 1.0, 0.0, -0.7, 0.2;
    VectorXd xfDense = eulerSimulate(AclDense, x0, 0.001, 10000);
    VectorXd xfSparse = eulerSimulate(AclSparse, x0, 0.001, 10000);

    std::cout << "\nFinal dense state:\n" << xfDense.transpose() << "\n";
    std::cout << "Final sparse state:\n" << xfSparse.transpose() << "\n";

    return 0;
}

This Java file implements the matrix operations from scratch. For larger projects, libraries such as EJML or Apache Commons Math are more appropriate for eigenvalue, SVD, and Riccati computations.

Chapter25_Lesson5.java

// Chapter25_Lesson5.java
/*
Structural Constraints: Limited Actuators and Sparse Feedback

Compile and run:
    javac Chapter25_Lesson5.java
    java Chapter25_Lesson5

This pure Java version implements basic matrix multiplication, controllability
matrix construction, numerical rank by Gaussian elimination, and Euler
simulation of the closed loop.
*/

public class Chapter25_Lesson5 {
    static double[][] multiply(double[][] A, double[][] B) {
        int n = A.length;
        int p = B[0].length;
        int m = B.length;
        double[][] C = new double[n][p];
        for (int i = 0; i < n; i++) {
            for (int k = 0; k < m; k++) {
                for (int j = 0; j < p; j++) {
                    C[i][j] += A[i][k] * B[k][j];
                }
            }
        }
        return C;
    }

    static double[][] subtract(double[][] A, double[][] B) {
        int n = A.length;
        int m = A[0].length;
        double[][] C = new double[n][m];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                C[i][j] = A[i][j] - B[i][j];
            }
        }
        return C;
    }

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

    static double[][] controllabilityMatrix(double[][] A, double[][] B) {
        int n = A.length;
        int m = B[0].length;
        double[][] C = new double[n][n * m];
        double[][] Ap = identity(n);
        for (int block = 0; block < n; block++) {
            double[][] ApB = multiply(Ap, B);
            for (int i = 0; i < n; i++) {
                for (int j = 0; j < m; j++) {
                    C[i][block * m + j] = ApB[i][j];
                }
            }
            Ap = multiply(Ap, A);
        }
        return C;
    }

    static int rank(double[][] input, double tol) {
        int rows = input.length;
        int cols = input[0].length;
        double[][] A = new double[rows][cols];
        for (int i = 0; i < rows; i++) {
            System.arraycopy(input[i], 0, A[i], 0, cols);
        }

        int r = 0;
        for (int c = 0; c < cols && r < rows; c++) {
            int pivot = r;
            for (int i = r + 1; i < rows; i++) {
                if (Math.abs(A[i][c]) > Math.abs(A[pivot][c])) {
                    pivot = i;
                }
            }
            if (Math.abs(A[pivot][c]) <= tol) {
                continue;
            }
            double[] temp = A[r];
            A[r] = A[pivot];
            A[pivot] = temp;

            double div = A[r][c];
            for (int j = c; j < cols; j++) {
                A[r][j] /= div;
            }
            for (int i = 0; i < rows; i++) {
                if (i != r) {
                    double factor = A[i][c];
                    for (int j = c; j < cols; j++) {
                        A[i][j] -= factor * A[r][j];
                    }
                }
            }
            r++;
        }
        return r;
    }

    static double[][] applyMask(double[][] K, double[][] mask) {
        int rows = K.length;
        int cols = K[0].length;
        double[][] S = new double[rows][cols];
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                S[i][j] = K[i][j] * mask[i][j];
            }
        }
        return S;
    }

    static double[] matVec(double[][] A, double[] x) {
        int n = A.length;
        int m = A[0].length;
        double[] y = new double[n];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                y[i] += A[i][j] * x[j];
            }
        }
        return y;
    }

    static double[] simulate(double[][] Acl, double[] x0, double dt, int steps) {
        double[] x = x0.clone();
        for (int k = 0; k < steps; k++) {
            double[] dx = matVec(Acl, x);
            for (int i = 0; i < x.length; i++) {
                x[i] += dt * dx[i];
            }
        }
        return x;
    }

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

    static void printVector(String label, double[] v) {
        System.out.print(label);
        for (double value : v) {
            System.out.printf("%10.4f ", value);
        }
        System.out.println();
    }

    public static void main(String[] args) {
        double[][] A = {
            {0.0,  1.0,  0.0,  0.0},
            {-2.0, -0.25, 0.7, 0.0},
            {0.0,  0.0,  0.0,  1.0},
            {0.6,  0.0, -1.5, -0.20}
        };

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

        double[][] Kdense = {
            {3.0, 2.0, 0.9, 0.4},
            {0.6, 0.3, 2.6, 1.8}
        };

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

        double[][] C = controllabilityMatrix(A, B);
        System.out.println("rank(C) = " + rank(C, 1.0e-9) + " / " + A.length);

        double[][] Ksparse = applyMask(Kdense, mask);
        printMatrix("\nDense K:", Kdense);
        printMatrix("\nSparse K:", Ksparse);

        double[][] AclDense = subtract(A, multiply(B, Kdense));
        double[][] AclSparse = subtract(A, multiply(B, Ksparse));

        double[] x0 = {1.0, 0.0, -0.7, 0.2};
        double[] xfDense = simulate(AclDense, x0, 0.001, 10000);
        double[] xfSparse = simulate(AclSparse, x0, 0.001, 10000);

        printVector("\nFinal dense state:  ", xfDense);
        printVector("Final sparse state: ", xfSparse);

        System.out.println("\nUse EJML or Apache Commons Math if eigenvalues and Riccati solvers");
        System.out.println("are needed in a larger Java modern-control project.");
    }
}

The MATLAB script uses Control System Toolbox functions and also creates a simple Simulink model for the optimized sparse closed-loop system.

Chapter25_Lesson5.m

% Chapter25_Lesson5.m
% Structural Constraints: Limited Actuators and Sparse Feedback
%
% MATLAB / Simulink implementation.
% Recommended toolboxes:
%   Control System Toolbox: ctrb, place, lqr, ss, initial
%   Simulink: automatic block-diagram construction at the end of this file

clear; clc; close all;

A = [ 0.0   1.0   0.0   0.0;
     -2.0  -0.25  0.7   0.0;
      0.0   0.0   0.0   1.0;
      0.6   0.0  -1.5  -0.20 ];

B_candidates = [0.0  0.0  0.0;
                1.0  0.2  0.0;
                0.0  0.0  0.0;
                0.0  0.4  1.0];

n = size(A,1);
q = size(B_candidates,2);

fprintf("Actuator subset analysis\n");
fprintf("------------------------\n");
for maskInt = 1:(2^q - 1)
    idx = find(bitget(maskInt,1:q));
    Bsel = B_candidates(:,idx);
    Cmat = ctrb(A,Bsel);
    fprintf("subset = ");
    fprintf("%d ", idx);
    fprintf(", rank(C) = %d / %d\n", rank(Cmat), n);
end

selected = [1 3];
B = B_candidates(:,selected);

Q = diag([20 1 20 1]);
R = diag([0.5 0.5]);

% Dense LQR design.
Kdense = lqr(A,B,Q,R);

% Sparse feedback pattern:
% actuator 1 receives x1,x2; actuator 2 receives x3,x4.
Mask = [1 1 0 0;
        0 0 1 1];

KsparseProjection = Kdense .* Mask;

fprintf("\nDense LQR gain:\n");
disp(Kdense);
fprintf("Sparse projected gain:\n");
disp(KsparseProjection);

fprintf("Dense closed-loop eigenvalues:\n");
disp(eig(A - B*Kdense));
fprintf("Sparse projected closed-loop eigenvalues:\n");
disp(eig(A - B*KsparseProjection));

% Optional nonlinear optimization over only the allowed entries of K.
% This searches for a sparse K whose closed-loop polynomial approximates
% the desired pole polynomial.
desiredPoles = [-1.2+1.0i, -1.2-1.0i, -1.8+0.8i, -1.8-0.8i];
targetPoly = poly(desiredPoles);

allowed = find(Mask(:) == 1);
theta0 = Kdense(allowed);

objective = @(theta) sparsePoleObjective(theta, A, B, Mask, allowed, targetPoly);
thetaStar = fminsearch(objective, theta0);
Kopt = zeros(size(Mask));
Kopt(allowed) = thetaStar;

fprintf("Optimized sparse gain:\n");
disp(Kopt);
fprintf("Optimized sparse closed-loop eigenvalues:\n");
disp(eig(A - B*Kopt));

sysDense = ss(A - B*Kdense, eye(n), eye(n), zeros(n));
sysSparse = ss(A - B*Kopt, eye(n), eye(n), zeros(n));
x0 = [1; 0; -0.7; 0.2];

figure;
initial(sysDense, sysSparse, x0, 10);
grid on;
legend("dense LQR", "optimized sparse");

% Simulink model generation: x_dot = A x + B u, u = -Kopt x.
modelName = "Chapter25_Lesson5_Simulink";
if bdIsLoaded(modelName)
    close_system(modelName,0);
end
new_system(modelName);
open_system(modelName);

add_block("simulink/Sources/Constant", modelName + "/zero_input", ...
    "Value", "0", "Position", [40 80 90 110]);

add_block("simulink/Continuous/State-Space", modelName + "/closed_loop_ss", ...
    "A", "A - B*Kopt", ...
    "B", "zeros(4,1)", ...
    "C", "eye(4)", ...
    "D", "zeros(4,1)", ...
    "X0", "x0", ...
    "Position", [150 55 300 135]);

add_block("simulink/Sinks/Scope", modelName + "/states_scope", ...
    "Position", [370 65 430 125]);

add_line(modelName, "zero_input/1", "closed_loop_ss/1");
add_line(modelName, "closed_loop_ss/1", "states_scope/1");

set_param(modelName, "StopTime", "10");
save_system(modelName);

function J = sparsePoleObjective(theta, A, B, Mask, allowed, targetPoly)
    K = zeros(size(Mask));
    K(allowed) = theta;
    p = poly(eig(A - B*K));
    poleError = norm(real(p - targetPoly),2)^2;
    gainPenalty = 1.0e-3 * norm(K,"fro")^2;
    J = poleError + gainPenalty;
end

The Mathematica notebook-style file uses symbolic and numeric matrix operations to compute controllability ranks, sparse gains, and eigenvalues.

Chapter25_Lesson5.nb

(* Chapter25_Lesson5.nb *)
Notebook[{
Cell["Chapter25_Lesson5.nb", "Title"],
Cell["Structural Constraints: Limited Actuators and Sparse Feedback", "Subtitle"],
Cell["This Wolfram Mathematica notebook-style script studies actuator selection, controllability rank, sparse feedback masks, and closed-loop eigenvalues.", "Text"],
Cell[BoxData[
 ToBoxes[
  ClearAll["Global`*"];
  A = { {0, 1, 0, 0}, {-2, -0.25, 0.7, 0}, {0, 0, 0, 1}, {0.6, 0, -1.5, -0.20} };
  Bcandidates = { {0, 0, 0}, {1, 0.2, 0}, {0, 0, 0}, {0, 0.4, 1} };
  n = Length[A];
  controllabilityMatrix[A_, B_] := ArrayFlatten[{Table[MatrixPower[A, k].B, {k, 0, n - 1}]}];
  subsets = Rest[Subsets[Range[3]]];
  actuatorReport = Table[
    Bsel = Bcandidates[[All, s]];
    {s, MatrixRank[controllabilityMatrix[A, Bsel]]},
    {s, subsets}
  ];
  actuatorReport
 ]
], "Input"],
Cell["A selected two-actuator structure is used below. The dense feedback is compared with a sparse communication mask.", "Text"],
Cell[BoxData[
 ToBoxes[
  B = Bcandidates[[All, {1, 3}]];
  Kdense = { {3.0, 2.0, 0.9, 0.4}, {0.6, 0.3, 2.6, 1.8} };
  mask = { {1, 1, 0, 0}, {0, 0, 1, 1} };
  Ksparse = Kdense*mask;
  {
    "OpenLoopEigenvalues" -> Eigenvalues[A],
    "DenseClosedLoopEigenvalues" -> Eigenvalues[A - B.Kdense],
    "SparseClosedLoopEigenvalues" -> Eigenvalues[A - B.Ksparse]
  }
 ]
], "Input"],
Cell["Finite-horizon controllability Gramian proxy for actuator-authority comparison.", "Text"],
Cell[BoxData[
 ToBoxes[
  gramianApprox[A_, B_, T_, h_] := Module[{W, Phi, steps},
    W = ConstantArray[0, {Length[A], Length[A]}];
    Phi = IdentityMatrix[Length[A]];
    steps = Round[T/h];
    Do[
      W = W + Phi.B.Transpose[B].Transpose[Phi] h;
      Phi = Phi + h A.Phi,
      {steps}
    ];
    W
  ];
  Table[
    Bsel = Bcandidates[[All, s]];
    {s, Tr[gramianApprox[A, Bsel, 5, 0.01]]},
    {s, subsets}
  ]
 ]
], "Input"]
}]

Problem 1 (Actuator selection and controllability): Consider \( A=\begin{bmatrix}0&1&0\\0&0&1\\0&0&0\end{bmatrix} \) and three candidate actuator columns \( b_1=\begin{bmatrix}1\\0\\0\end{bmatrix} \), \( b_2=\begin{bmatrix}0\\1\\0\end{bmatrix} \), and \( b_3=\begin{bmatrix}0\\0\\1\end{bmatrix} \). Which single actuator makes the pair controllable?

Solution: For a single column \( b \), compute \( \mathcal{C}=[b\;Ab\;A^2b] \). If \( b=b_1 \), then \( Ab_1=0 \) and the rank is 1. If \( b=b_2 \), then \( Ab_2=b_1 \) and \( A^2b_2=0 \), so the rank is 2. If \( b=b_3 \), then \( Ab_3=b_2 \) and \( A^2b_3=b_1 \), so

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

Thus only the third actuator location gives controllability with one actuator.

Problem 2 (Unassignable pole proof): Let \( q^*A=\lambda q^* \) and \( q^*B=0 \). Show that no feedback gain \( K \) can move \( \lambda \).

Solution: Multiplying the closed-loop matrix from the left gives

\[ q^*(A-BK)=q^*A-q^*BK=\lambda q^*. \]

Therefore \( \lambda \) remains a closed-loop eigenvalue. If \( \operatorname{Re}(\lambda)>0 \), the system cannot be stabilized through that input matrix.

Problem 3 (Sparse parameter count): A fourth-order SISO system uses the sparse feedback law \( u=-(k_1x_1+k_3x_3) \). Explain why arbitrary fourth-order pole placement is generically impossible.

Solution: The closed-loop characteristic polynomial has four independent coefficients, but the feedback has only two free parameters, \( k_1 \) and \( k_3 \). The coefficient map is \( F:\mathbb{R}^2\longrightarrow\mathbb{R}^4 \). A two-dimensional parameter set cannot generically cover an open subset of the four-dimensional coefficient space, so only a restricted family of pole sets can be achieved.

Problem 4 (Sparse projection is not pole placement): Suppose a dense gain \( K_d \) is designed by LQR or pole placement. A sparse mask is then applied: \( K_s=K_d\circ M \). Is \( K_s \) guaranteed to preserve closed-loop stability?

Solution: No. Projection removes feedback links and changes the closed-loop matrix from \( A-BK_d \) to \( A-B(K_d\circ M) \). The difference is

\[ \Delta A_{cl}=B\left(K_d-K_d\circ M\right). \]

If this perturbation is large in directions associated with lightly damped modes, eigenvalues may cross into the unstable half-plane. A sparse gain should therefore be designed directly under the mask, not merely obtained by deleting entries after a dense design.

Problem 5 (Gramian interpretation): Let \( W_c \) have eigenvalues \( 10,1,0.01 \). Which state direction is hardest to reach and why?

Solution: The direction associated with eigenvalue \( 0.01 \) is hardest to reach. Since \( E_{\min}(x_f)=x_f^TW_c^{-1}x_f \), the energy along a Gramian eigenvector is inversely proportional to the corresponding Gramian eigenvalue. Thus the energy scaling in that direction is approximately \( 1/0.01=100 \).

Structural constraints limit what state feedback can accomplish even before numerical conditioning, uncertainty, or actuator saturation are considered. Limited actuators restrict the input matrix and can create uncontrollable modes. Sparse feedback restricts the gain matrix and can prevent arbitrary pole placement even when the selected pair is controllable. Practical design therefore requires actuator selection, controllability tests, Gramian-based authority checks, sparse gain design, and closed-loop validation.