In the state-space model \( \dot{\mathbf{x} } = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u} \), the matrix \( \mathbf{A} \) describes the internal dynamics, while \( \mathbf{B} \) describes how actuators enter the system. Therefore, actuator placement is represented mathematically by the columns of \( \mathbf{B} \). If \( \mathbf{u}\in\mathbb{R}^m \), then

\[ \mathbf{B} = \begin{bmatrix} \mathbf{b}_1 & \mathbf{b}_2 & \cdots & \mathbf{b}_m \end{bmatrix}, \qquad \mathbf{B}\mathbf{u} = \sum_{j=1}^m \mathbf{b}_j u_j. \]

The vector \( \mathbf{b}_j \) is the immediate state-space direction produced by actuator \( j \). For example, a force actuator usually appears directly in acceleration states rather than position states. A torque actuator appears in angular acceleration equations. A voltage actuator may first affect electrical current, and only later mechanical motion through electromechanical coupling.

The key physical point is that direct actuator direction and dynamic propagation together create possible motion directions. A system may have few actuators, yet be steerable because the natural coupling in \( \mathbf{A} \) spreads those actuator effects into other state directions.

Suppose a short pulse is applied through actuator \( j \) at time \( \tau \). Its immediate effect is along \( \mathbf{b}_j \), but after time \( t-\tau \) this effect becomes \( e^{\mathbf{A}(t-\tau)}\mathbf{b}_j \). Thus the forced part of the state trajectory is

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

This formula shows why actuator placement cannot be evaluated from \( \mathbf{B} \) alone. The columns of \( \mathbf{B} \) are propagated by \( e^{\mathbf{A}s} \). Therefore, the relevant family of input-generated directions includes \( \mathbf{B}, \mathbf{A}\mathbf{B}, \mathbf{A}^2\mathbf{B}, \ldots \). In finite-dimensional systems, only the first \( n \) powers are needed to span all possible directions generated by repeated dynamic propagation.

\[ \mathcal{R} = \operatorname{span}\left\{ \mathbf{B},\mathbf{A}\mathbf{B}, \mathbf{A}^2\mathbf{B},\ldots,\mathbf{A}^{n-1}\mathbf{B} \right\}. \]

Here \( \mathcal{R} \) is the reachable subspace already introduced in earlier lessons. Its physical interpretation is: the actuator does not need to push directly in every state direction; it must push in directions that the system dynamics can spread into the desired state directions.

flowchart TD
  A["Actuator location"] --> B["Input direction b_j"]
  B --> C["Immediate motion direction"]
  C --> D["Dynamic propagation by A"]
  D --> E["Directions: b_j, A b_j, A^2 b_j, ..."]
  E --> F["Reachable subspace"]
  F --> G["State-steering authority"]
        

Consider a linearized second-order mechanical model

\[ \mathbf{M}\ddot{\mathbf{q} } + \mathbf{D}\dot{\mathbf{q} } + \mathbf{K}\mathbf{q} = \mathbf{G}\mathbf{u}, \]

where \( \mathbf{q} \) is the generalized coordinate vector, \( \mathbf{M} \) is the mass or inertia matrix, \( \mathbf{D} \) is damping, \( \mathbf{K} \) is stiffness, and \( \mathbf{G} \) maps actuator forces or torques into generalized coordinates. With \( \mathbf{x}=\begin{bmatrix}\mathbf{q}^T & \dot{\mathbf{q} }^T\end{bmatrix}^T \), we obtain

\[ \dot{\mathbf{x} } = \underbrace{\begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{D} \end{bmatrix} }_{\mathbf{A} } \mathbf{x} + \underbrace{\begin{bmatrix} \mathbf{0} \\ \mathbf{M}^{-1}\mathbf{G} \end{bmatrix} }_{\mathbf{B} } \mathbf{u}. \]

This equation is a direct mathematical expression of actuator placement. If a force is applied to the first mass, the corresponding column of \( \mathbf{G} \) has a nonzero entry in the first generalized-force coordinate. If an actuator is placed at a joint, the column represents the joint torque direction. If a thruster is placed away from the center of mass, its column may contain both force and moment components.

The matrix \( \mathbf{M}^{-1}\mathbf{G} \) also shows why authority depends on inertia. The same actuator force may create strong acceleration in a light direction but weak acceleration in a direction with large effective inertia. Hence actuator placement is not only geometrical; it is also inertial and dynamic.

A useful physical question is not only whether the system can be moved, but whether each natural mode can be influenced. Suppose \( \mathbf{A} \) is diagonalizable:

\[ \mathbf{A} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}, \qquad \mathbf{z}=\mathbf{V}^{-1}\mathbf{x}. \]

In modal coordinates,

\[ \dot{\mathbf{z} } = \boldsymbol{\Lambda}\mathbf{z} + \underbrace{\mathbf{V}^{-1}\mathbf{B} }_{\mathbf{B}_m} \mathbf{u}. \]

The \( i \)-th row of \( \mathbf{B}_m \) tells how strongly the input excites mode \( i \). If \( \mathbf{w}_i^T \) is the left eigenvector of \( \mathbf{A} \), then the scalar or row vector \( \mathbf{w}_i^T\mathbf{B} \) is the modal input gain. A common quantitative modal-authority score is

\[ \alpha_i = \left\|\mathbf{w}_i^T\mathbf{B}\right\|_2. \]

If \( \alpha_i=0 \), input cannot excite mode \( i \) directly in modal coordinates. If a natural mode has little or no modal authority, it may be practically difficult or impossible to move, even when other modes are easily moved.

Proof of the modal statement. For one modal coordinate,

\[ \dot{z}_i = \lambda_i z_i + \mathbf{w}_i^T\mathbf{B}\mathbf{u}, \qquad z_i(t)=e^{\lambda_i t}z_i(0)+ \int_0^t e^{\lambda_i(t-\tau)}\mathbf{w}_i^T\mathbf{B} \mathbf{u}(\tau)\,d\tau. \]

If \( \mathbf{w}_i^T\mathbf{B}=\mathbf{0} \), the integral term is always zero. Thus the input cannot change that modal coordinate beyond its natural homogeneous evolution. This is the physical meaning of an actuator being placed at a modal node or in a direction orthogonal to that mode.

Physical authority is not only binary. Two actuator placements may both reach the full state space, but one may require much more input energy. For a finite duration \( T \), define the finite-horizon controllability Gramian

\[ \mathbf{W}_c(T)= \int_0^T e^{\mathbf{A}\tau}\mathbf{B}\mathbf{B}^T e^{\mathbf{A}^T\tau}\,d\tau. \]

For a unit target direction \( \mathbf{v} \), the scalar \( \mathbf{v}^T\mathbf{W}_c(T)\mathbf{v} \) measures how much input-generated motion appears along \( \mathbf{v} \) during the time interval. Large values indicate high authority; small values indicate weak authority.

\[ a_T(\mathbf{v})=\mathbf{v}^T\mathbf{W}_c(T)\mathbf{v}, \qquad \|\mathbf{v}\|_2=1. \]

When \( \mathbf{W}_c(T) \) is nonsingular, the minimum input energy needed to move from the origin to a target \( \mathbf{x}_f \) at time \( T \) is

\[ E_{\min}(\mathbf{x}_f,T)= \mathbf{x}_f^T\mathbf{W}_c(T)^{-1}\mathbf{x}_f. \]

Therefore, actuator authority has an energy interpretation: directions corresponding to small Gramian eigenvalues are expensive directions. Directions corresponding to large Gramian eigenvalues are easy directions. This explains why practical actuator placement often seeks not merely reachability, but well-conditioned reachability.

Sketch of proof. For zero initial condition, the final state satisfies

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

Applying a Hilbert-space minimum-norm argument gives the energy-optimal input

\[ \mathbf{u}^*(\tau)= \mathbf{B}^T e^{\mathbf{A}^T(T-\tau)} \mathbf{W}_c(T)^{-1}\mathbf{x}_f. \]

Substituting \( \mathbf{u}^* \) into \( \int_0^T \|\mathbf{u}(\tau)\|_2^2d\tau \) yields the stated expression for \( E_{\min} \).

With multiple actuators, \( \mathbf{B} \) may have redundant columns. Redundancy is not useless: it can reduce effort, improve fault tolerance, and allow input allocation under constraints. However, if two actuators produce nearly the same column direction, they may increase force capacity without increasing directional diversity.

If a desired instantaneous generalized input \( \mathbf{v} \) lies in \( \operatorname{range}(\mathbf{B}) \), one unconstrained allocation is

\[ \mathbf{u} = \mathbf{B}^\dagger \mathbf{v} + \left(\mathbf{I}-\mathbf{B}^\dagger\mathbf{B}\right) \boldsymbol{\eta}, \]

where \( \mathbf{B}^\dagger \) is the Moore-Penrose pseudoinverse and \( \boldsymbol{\eta} \) is an arbitrary vector in the redundant-actuation null-space. In real systems, constraints such as \( |u_j|\leq u_{j,\max} \), rate limits, dead zones, and actuator bandwidth reduce practical authority even when the mathematical reachable subspace is large.

Hence the phrase actuator authority has three layers:

1. Geometric authority: which state directions can be generated by \( \mathbf{B} \) and its propagation through \( \mathbf{A} \).
2. Energetic authority: how much input energy is needed to reach specific directions.
3. Constrained authority: what can be achieved under amplitude, rate, thermal, saturation, or safety limits.

In engineering practice, actuator placement is evaluated before final feedback design. The designer first writes the physical model, derives \( \mathbf{B} \) for each candidate placement, and then compares geometric and energetic authority. This procedure is especially important for flexible structures, robotic manipulators, aircraft, spacecraft, marine vehicles, and coupled electromechanical systems.

flowchart TD
  P["Physical model"] --> G["Candidate actuator geometry G"]
  G --> B["Build input matrix B"]
  B --> R["Check generated directions: B, AB, A^2B, ..."]
  R --> M["Compute modal authority scores"]
  M --> W["Compute finite-time energy indicators"]
  W --> C["Check constraints and saturation"]
  C --> D["Choose placement with adequate authority"]
        

This workflow also prevents a common design mistake: selecting an actuator because it has high local force, while ignoring whether that force actually couples into the important modes or state directions.

The following script compares several actuator placements in a two-degree-of-freedom mass-spring-damper model. It uses NumPy and SciPy to compute the reachability matrix, a finite-horizon Gramian, directional authority, and modal authority scores.

Chapter10_Lesson4.py

"""
Chapter10_Lesson4.py

Physical Interpretation: Actuator Placement and Authority
Modern Control - Chapter 10, Lesson 4

This script compares actuator placements for a two-degree-of-freedom
mass-spring-damper system. It computes:
1. The finite-dimensional reachability matrix [B, AB, ..., A^(n-1)B].
2. The finite-horizon controllability Gramian Wc(T).
3. Directional authority along target state directions.
4. A modal projection score based on left eigenvectors of A.

Dependencies:
    pip install numpy scipy
"""

import numpy as np
from scipy.linalg import expm


def controllability_matrix(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """Return C = [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 numerical_gramian(A: np.ndarray, B: np.ndarray, T: float = 5.0, N: int = 1000) -> np.ndarray:
    """
    Approximate Wc(T) = int_0^T exp(A tau) B B^T exp(A^T tau) d tau
    by the trapezoidal rule.
    """
    n = A.shape[0]
    W = np.zeros((n, n), dtype=float)
    dt = T / N

    for k in range(N + 1):
        tau = k * dt
        Phi = expm(A * tau)
        integrand = Phi @ B @ B.T @ Phi.T
        weight = 0.5 if k == 0 or k == N else 1.0
        W += weight * integrand

    return W * dt


def modal_authority(A: np.ndarray, B: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """
    Compute ||w_i^T B||_2 for each left eigenvector w_i^T of A.
    If this value is near zero, mode i has weak or zero direct modal authority.
    """
    eigvals, V = np.linalg.eig(A)
    Vinv = np.linalg.inv(V)
    scores = np.linalg.norm(Vinv @ B, axis=1)
    return eigvals, scores


def summarize_placement(name: str, A: np.ndarray, B: np.ndarray, target_directions: list[tuple[str, np.ndarray]]) -> None:
    print("\n" + "=" * 72)
    print(f"Actuator placement: {name}")
    print("=" * 72)

    Ctrb = controllability_matrix(A, B)
    rank_C = np.linalg.matrix_rank(Ctrb, tol=1e-9)
    print("Rank of [B, AB, ..., A^(n-1)B]:", rank_C, "out of", A.shape[0])

    W = numerical_gramian(A, B, T=5.0, N=800)
    eigW = np.linalg.eigvalsh(W)
    print("Eigenvalues of finite-horizon Gramian Wc(5):")
    print(np.round(eigW, 8))

    print("Directional authority v^T Wc v for selected target directions:")
    for label, v in target_directions:
        v = v.reshape(-1, 1)
        score = (v.T @ W @ v).item()
        print(f"  {label:18s}: {score:.8f}")

    eigvals, scores = modal_authority(A, B)
    print("Modal authority scores ||w_i^T B||_2:")
    for lam, score in zip(eigvals, scores):
        print(f"  lambda={lam.real:+.4f}{lam.imag:+.4f}j, score={score:.8f}")


def main() -> None:
    # Two coupled unit masses with spring coupling and small damping:
    # qdd = -K q - D qdot + G u
    K = np.array([[2.0, -1.0],
                  [-1.0, 2.0]])
    D = np.array([[0.08, 0.00],
                  [0.00, 0.08]])
    M = np.eye(2)

    # State x = [q1, q2, q1dot, q2dot]^T
    Z = np.zeros((2, 2))
    I = np.eye(2)
    A = np.block([[Z, I],
                  [-np.linalg.solve(M, K), -np.linalg.solve(M, D)]])

    # Candidate actuator force-placement matrices G in qdd = ... + M^{-1}G u
    candidates = {
        "force on mass 1 only": np.array([[1.0], [0.0]]),
        "force on mass 2 only": np.array([[0.0], [1.0]]),
        "same force on both masses": np.array([[1.0], [1.0]]),
        "independent forces on both masses": np.eye(2),
    }

    target_directions = [
        ("displace mass 1", np.array([1.0, 0.0, 0.0, 0.0])),
        ("displace mass 2", np.array([0.0, 1.0, 0.0, 0.0])),
        ("relative motion", np.array([1.0, -1.0, 0.0, 0.0]) / np.sqrt(2)),
        ("common motion", np.array([1.0, 1.0, 0.0, 0.0]) / np.sqrt(2)),
    ]

    for name, G in candidates.items():
        B = np.vstack([np.zeros((2, G.shape[1])), np.linalg.solve(M, G)])
        summarize_placement(name, A, B, target_directions)


if __name__ == "__main__":
    main()

The C++ implementation avoids external numerical libraries and builds the reachability matrix from scratch. This is useful pedagogically because students can see how repeated multiplication by \( \mathbf{A} \) produces dynamically propagated actuator directions.

Chapter10_Lesson4.cpp

/*
Chapter10_Lesson4.cpp

Physical Interpretation: Actuator Placement and Authority
Modern Control - Chapter 10, Lesson 4

This from-scratch C++ example computes the reachability matrix
[B, AB, A^2B, ..., A^(n-1)B] and its numerical rank for several actuator
placements in a two-degree-of-freedom mechanical system.

Build:
    g++ -std=c++17 Chapter10_Lesson4.cpp -O2 -o Chapter10_Lesson4
*/

#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>
#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 kdim = static_cast<int>(B.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 < kdim; ++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& B : blocks) totalCols += static_cast<int>(B[0].size());

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

Matrix controllabilityMatrix(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 numericalRank(Matrix M, double tol = 1e-9) {
    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 i = rank + 1; i < rows; ++i) {
            if (std::fabs(M[i][col]) > std::fabs(M[pivot][col])) pivot = i;
        }

        if (std::fabs(M[pivot][col]) <= tol) continue;
        std::swap(M[pivot], M[rank]);

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

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

double columnEnergyScore(const Matrix& C) {
    double sum = 0.0;
    for (const auto& row : C) {
        for (double value : row) sum += value * value;
    }
    return std::sqrt(sum);
}

void analyzePlacement(const std::string& name, const Matrix& A, const Matrix& B) {
    Matrix C = controllabilityMatrix(A, B);
    int rank = numericalRank(C);
    double score = columnEnergyScore(C);

    std::cout << "\nActuator placement: " << name << "\n";
    std::cout << "Rank of [B, AB, ..., A^(n-1)B]: " << rank << " out of "
              << A.size() << "\n";
    std::cout << "Frobenius norm authority score of reachability matrix: "
              << std::fixed << std::setprecision(6) << score << "\n";
}

int main() {
    // State x = [q1, q2, q1dot, q2dot]^T for two coupled unit masses.
    Matrix A = {
        { 0.0,  0.0, 1.0, 0.0},
        { 0.0,  0.0, 0.0, 1.0},
        {-2.0,  1.0,-0.08,0.0},
        { 1.0, -2.0, 0.0,-0.08}
    };

    // B = [0; 0; G] because M = I in qdd = ... + G u.
    Matrix B_mass1 = {
        {0.0},
        {0.0},
        {1.0},
        {0.0}
    };

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

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

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

    analyzePlacement("force on mass 1 only", A, B_mass1);
    analyzePlacement("force on mass 2 only", A, B_mass2);
    analyzePlacement("same force on both masses", A, B_common);
    analyzePlacement("independent forces on both masses", A, B_both);

    return 0;
}

The Java version mirrors the C++ algorithm and is suitable for students who want to implement reachability computations without relying on a specialized control package.

Chapter10_Lesson4.java

/*
Chapter10_Lesson4.java

Physical Interpretation: Actuator Placement and Authority
Modern Control - Chapter 10, Lesson 4

This Java example computes the reachability matrix and its numerical rank
for several actuator-placement choices in a two-degree-of-freedom system.

Compile and run:
    javac Chapter10_Lesson4.java
    java Chapter10_Lesson4
*/

public class Chapter10_Lesson4 {
    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 kdim = B.length;
        int c = B[0].length;
        double[][] C = zeros(r, c);

        for (int i = 0; i < r; i++) {
            for (int k = 0; k < kdim; 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 col0 = 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][col0 + j] = block[i][j];
                }
            }
            col0 += cols;
        }
        return H;
    }

    static double[][] controllabilityMatrix(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 numericalRank(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 i = rank + 1; i < rows; i++) {
                if (Math.abs(M[i][col]) > Math.abs(M[pivot][col])) pivot = i;
            }

            if (Math.abs(M[pivot][col]) <= tol) continue;

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

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

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

        return rank;
    }

    static double frobeniusNorm(double[][] M) {
        double sum = 0.0;
        for (double[] row : M) {
            for (double value : row) sum += value * value;
        }
        return Math.sqrt(sum);
    }

    static void analyzePlacement(String name, double[][] A, double[][] B) {
        double[][] C = controllabilityMatrix(A, B);
        int rank = numericalRank(C, 1e-9);
        double score = frobeniusNorm(C);

        System.out.println("\nActuator placement: " + name);
        System.out.println("Rank of [B, AB, ..., A^(n-1)B]: " + rank + " out of " + A.length);
        System.out.printf("Frobenius norm authority score of reachability matrix: %.6f%n", score);
    }

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

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

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

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

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

        analyzePlacement("force on mass 1 only", A, BMass1);
        analyzePlacement("force on mass 2 only", A, BMass2);
        analyzePlacement("same force on both masses", A, BCommon);
        analyzePlacement("independent forces on both masses", A, BBoth);
    }
}

MATLAB provides direct control-system functions such as ctrb, while Mathematica is convenient for symbolic and high-precision matrix operations. In Simulink, the same model can be tested using a State-Space block with different candidate \( \mathbf{B} \) matrices.

Chapter10_Lesson4.m

% Chapter10_Lesson4.m
%
% Physical Interpretation: Actuator Placement and Authority
% Modern Control - Chapter 10, Lesson 4
%
% This script compares actuator placements for a two-degree-of-freedom
% mass-spring-damper system using MATLAB Control System Toolbox functions
% and a finite-horizon controllability Gramian.
%
% Related Simulink idea:
%   Use a State-Space block with A, B, C = eye(4), D = zeros(4,m).
%   Feed candidate inputs u(t) through Step, Signal Builder, or From Workspace.
%   Compare state trajectories for different B matrices.

clear; clc;

K = [2 -1; -1 2];
Damp = [0.08 0; 0 0.08];
M = eye(2);

A = [zeros(2) eye(2);
     -M\K     -M\Damp];

placements = {
    'force on mass 1 only', [1; 0];
    'force on mass 2 only', [0; 1];
    'same force on both masses', [1; 1];
    'independent forces on both masses', eye(2)
};

targetNames = {'displace mass 1', 'displace mass 2', 'relative motion', 'common motion'};
targets = [1 0 0 0;
           0 1 0 0;
           1/sqrt(2) -1/sqrt(2) 0 0;
           1/sqrt(2)  1/sqrt(2) 0 0];

T = 5;

for idx = 1:size(placements,1)
    name = placements{idx,1};
    G = placements{idx,2};
    B = [zeros(2,size(G,2)); M\G];

    Ctrb = ctrb(A,B);
    r = rank(Ctrb);

    W = integral(@(tau) expm(A*tau)*B*B'*expm(A'*tau), ...
                 0, T, 'ArrayValued', true);

    fprintf('\n============================================\n');
    fprintf('Actuator placement: %s\n', name);
    fprintf('Rank of ctrb(A,B): %d out of %d\n', r, size(A,1));
    fprintf('Eigenvalues of Wc(%g):\n', T);
    disp(eig(W).');

    fprintf('Directional authority v''*W*v:\n');
    for k = 1:size(targets,1)
        v = targets(k,:).';
        fprintf('  %-18s : %.8f\n', targetNames{k}, v.'*W*v);
    end

    [V,Lambda] = eig(A);
    modalInput = inv(V)*B;
    modalScores = vecnorm(modalInput,2,2);
    fprintf('Modal authority scores ||w_i''*B||_2:\n');
    for k = 1:length(modalScores)
        fprintf('  lambda = %.4f%+.4fi, score = %.8f\n', ...
            real(Lambda(k,k)), imag(Lambda(k,k)), modalScores(k));
    end
end

Chapter10_Lesson4.nb

(* Chapter10_Lesson4.nb *)

(* Physical Interpretation: Actuator Placement and Authority
   Modern Control - Chapter 10, Lesson 4

   Paste this code into a Mathematica notebook. It computes the reachability
   matrix, finite-horizon controllability Gramian, directional authority,
   and modal authority for several actuator placements.
*)

ClearAll["Global`*"];

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

finiteHorizonGramian[A_, B_, T_] := NIntegrate[
  MatrixExp[A tau].B.Transpose[B].MatrixExp[Transpose[A] tau],
  {tau, 0, T}
];

modalAuthority[A_, B_] := Module[{vals, vecs, vinv, scores},
  {vals, vecs} = Eigensystem[A];
  (* Eigensystem returns eigenvectors as rows; transpose makes columns. *)
  vinv = Inverse[Transpose[vecs]];
  scores = Norm /@ (vinv.B);
  Transpose[{vals, scores}]
];

K = { {2, -1}, {-1, 2} };
Damp = { {0.08, 0}, {0, 0.08} };
A = ArrayFlatten[{
   {ConstantArray[0, {2, 2}], IdentityMatrix[2]},
   {-K, -Damp}
}];

placements = {
   {"force on mass 1 only", { {1}, {0} } },
   {"force on mass 2 only", { {0}, {1} } },
   {"same force on both masses", { {1}, {1} } },
   {"independent forces on both masses", IdentityMatrix[2]}
};

targets = {
   {"displace mass 1", {1, 0, 0, 0} },
   {"displace mass 2", {0, 1, 0, 0} },
   {"relative motion", {1/Sqrt[2], -1/Sqrt[2], 0, 0} },
   {"common motion", {1/Sqrt[2], 1/Sqrt[2], 0, 0} }
};

Tfinal = 5;

Table[
  Module[{name, G, B, Ctrb, W, dirScores, modScores},
    name = placement[[1]];
    G = placement[[2]];
    B = ArrayFlatten[{ {ConstantArray[0, {2, Length[Transpose[G]]}]}, {G} }];

    Ctrb = controllabilityMatrix[A, B];
    W = finiteHorizonGramian[A, B, Tfinal];

    dirScores = Table[
      {target[[1]], N[target[[2]].W.target[[2]]]},
      {target, targets}
    ];

    modScores = modalAuthority[A, B];

    <|
      "Placement" -> name,
      "ReachabilityRank" -> MatrixRank[Ctrb],
      "GramianEigenvalues" -> N[Eigenvalues[W]],
      "DirectionalAuthority" -> dirScores,
      "ModalAuthority" -> N[modScores]
    |>
  ],
  {placement, placements}
] // Dataset

Problem 1 (Input Matrix from a Mechanical Model): Consider \( m\ddot{q}+d\dot{q}+kq=gu \) with state \( \mathbf{x}=\begin{bmatrix}q & \dot{q}\end{bmatrix}^T \). Derive \( \mathbf{A} \) and \( \mathbf{B} \). Explain what actuator authority means in this scalar case.

Solution:

\[ \dot{\mathbf{x} } = \begin{bmatrix}0 & 1 \\ -k/m & -d/m\end{bmatrix}\mathbf{x}+ \begin{bmatrix}0 \\ g/m\end{bmatrix}u. \]

Thus the actuator directly affects acceleration, not position. Position becomes reachable through the kinematic relation \( \dot{q}=v \). In this scalar oscillator, larger \( |g|/m \) means stronger acceleration authority per unit input.

Problem 2 (Dynamic Propagation of an Actuator Direction): Let \( \mathbf{A}=\begin{bmatrix}0&1\\0&0\end{bmatrix} \) and \( \mathbf{B}=\begin{bmatrix}0\\1\end{bmatrix} \). Compute \( \mathbf{B} \) and \( \mathbf{A}\mathbf{B} \). Interpret the result.

Solution:

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

The actuator directly changes velocity, but the dynamics propagate velocity into position. Therefore both velocity and position directions are generated. Physically, a force first changes velocity, and then velocity changes position.

Problem 3 (Weak Modal Authority): Suppose \( \mathbf{A} \) has a left eigenvector \( \mathbf{w}_i^T \) and \( \mathbf{w}_i^T\mathbf{B}=\mathbf{0} \). Show that input cannot affect the corresponding modal coordinate.

Solution: The modal coordinate satisfies

\[ \dot{z}_i=\lambda_i z_i+ \mathbf{w}_i^T\mathbf{B}\mathbf{u}= \lambda_i z_i. \]

Its solution is \( z_i(t)=e^{\lambda_i t}z_i(0) \). Since no term containing \( \mathbf{u} \) remains, the input cannot modify that modal coordinate.

Problem 4 (Directional Authority from a Gramian): Let \( \mathbf{W}_c(T)=\operatorname{diag}(10,1,0.01) \). Compare the authority along the three coordinate directions.

Solution: For coordinate direction \( \mathbf{e}_i \),

\[ a_T(\mathbf{e}_i)= \mathbf{e}_i^T\mathbf{W}_c(T)\mathbf{e}_i. \]

Hence the authority values are \( 10 \), \( 1 \), and \( 0.01 \). The third direction is much more difficult to reach energetically because the corresponding Gramian value is small.

Problem 5 (Placement Interpretation in a Two-Mass System): A two-mass system has state \( \mathbf{x}=\begin{bmatrix}q_1&q_2&\dot{q}_1&\dot{q}_2\end{bmatrix}^T \). Compare an actuator force on mass 1 only with two independent force actuators, one on each mass.

Solution: With unit masses, an actuator on mass 1 only gives

\[ \mathbf{B}_1= \begin{bmatrix}0\\0\\1\\0\end{bmatrix}. \]

Two independent force actuators give

\[ \mathbf{B}_2= \begin{bmatrix}0&0\\0&0\\1&0\\0&1\end{bmatrix}. \]

The single actuator may still move both masses if coupling exists in \( \mathbf{A} \), but the two-actuator placement has more direct acceleration authority and usually lower energy cost for independently shaping common and relative motion.

Actuator placement determines the input matrix \( \mathbf{B} \), while system dynamics determine how those input directions spread through state space. The physical authority of an actuator is therefore measured by direct input directions, propagated directions such as \( \mathbf{A}\mathbf{B} \), modal projections \( \mathbf{w}_i^T\mathbf{B} \), and finite-time energy indicators from \( \mathbf{W}_c(T) \). Good actuator placement gives enough geometric authority, strong modal authority, and acceptable energy under real constraints.