Consider the continuous-time LTI state-space model \( \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t) \), \( \mathbf{y}(t)=\mathbf{C}\mathbf{x}(t)+\mathbf{D}\mathbf{u}(t) \). In observability analysis, we ask whether the initial state \( \mathbf{x}(0) \) can be inferred from the measured output history. Sensor placement determines the rows of \( \mathbf{C} \). Therefore, changing sensors changes the map from internal state coordinates to measured quantities.

If candidate sensors are available as row vectors \( \mathbf{c}_1,\dots,\mathbf{c}_p \in \mathbb{R}^{1\times n} \), a selected sensor set \( \mathcal{S}\subseteq\{1,\dots,p\} \) produces

\[ \mathbf{C}_{\mathcal{S} } = \begin{bmatrix} \mathbf{c}_{i_1}\\ \mathbf{c}_{i_2}\\ \vdots\\ \mathbf{c}_{i_q} \end{bmatrix},\quad \mathcal{S}=\{i_1,i_2,\dots,i_q\}. \]

Direct state sensing is a special case: measuring state \( x_j \) uses the row \( \mathbf{e}_j^T \). However, physical sensors often measure combinations of states: a velocity sensor may measure a velocity coordinate directly; an accelerometer may measure a combination involving dynamics, gravity, and orientation; a voltage sensor may measure only a subset of electrical states. The lesson is: a sensor is a row of \( \mathbf{C} \), but observability depends on both \( \mathbf{C} \) and \( \mathbf{A} \).

flowchart TD
  A["Physical model: xdot = A x + B u"] --> B["Candidate sensors: c1, c2, ..., cp"]
  B --> C["Choose sensor set S"]
  C --> D["Build output matrix C_S"]
  D --> E["Stack dynamic signatures: C_S, C_S A, ..., C_S A^(n-1)"]
  E --> F["Check state directions that appear in output history"]
  F --> G["Observable enough? rank, conditioning, noise sensitivity"]
  G --> H["Accept placement or choose different sensors"]
        

Assume first that the known input effect has been subtracted from the output, so that observability depends on the homogeneous response generated by \( \mathbf{x}(0) \). For zero input,

\[ \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}(0),\quad \mathbf{y}(t)=\mathbf{C}e^{\mathbf{A}t}\mathbf{x}(0). \]

At \( t=0 \), output derivatives reveal repeated dynamic signatures of the initial state:

\[ \mathbf{y}^{(k)}(0)=\mathbf{C}\mathbf{A}^{k}\mathbf{x}(0), \quad k=0,1,\dots,n-1. \]

Thus a candidate sensor placement should not be evaluated only by \( \mathbf{C}\mathbf{x} \). It must be evaluated by the collection of rows \( \mathbf{C},\mathbf{C}\mathbf{A},\dots,\mathbf{C}\mathbf{A}^{n-1} \). A sensor may not measure a state directly, yet that state may become visible through dynamic coupling in \( \mathbf{A} \).

\[ \mathcal{M}(\mathbf{C},\mathbf{A}) = \begin{bmatrix} \mathbf{C}\\ \mathbf{C}\mathbf{A}\\ \mathbf{C}\mathbf{A}^2\\ \vdots\\ \mathbf{C}\mathbf{A}^{n-1} \end{bmatrix}. \]

In the next chapter, this stacked matrix becomes the formal observability matrix. Here, we use it conceptually: it records how a sensor row and its repeated propagation through the dynamics expose hidden state directions.

A nonzero initial state direction \( \mathbf{v}\neq\mathbf{0} \) is invisible to a sensor placement if it generates zero output and zero output derivatives:

\[ \mathbf{C}\mathbf{A}^{k}\mathbf{v}=\mathbf{0}, \quad k=0,1,\dots,n-1. \]

Therefore the unobservable subspace associated with a placement is

\[ \mathcal{N}_o(\mathbf{C},\mathbf{A}) = \bigcap_{k=0}^{n-1}\ker(\mathbf{C}\mathbf{A}^{k}). \]

A good placement shrinks this intersection as much as possible. Full observability corresponds to \( \mathcal{N}_o(\mathbf{C},\mathbf{A})=\{\mathbf{0}\} \). If the intersection contains a nonzero direction, two different initial states separated by that direction produce identical output histories.

Proof idea: Suppose \( \mathbf{v}\in\mathcal{N}_o(\mathbf{C},\mathbf{A}) \). The Taylor expansion of the zero-input output is

\[ \mathbf{C}e^{\mathbf{A}t}\mathbf{v} = \sum_{k=0}^{\infty}\frac{t^k}{k!}\mathbf{C}\mathbf{A}^{k}\mathbf{v}. \]

By the Cayley-Hamilton theorem, all powers \( \mathbf{A}^{k} \) for \( k\ge n \) are linear combinations of \( \mathbf{I},\mathbf{A},\dots,\mathbf{A}^{n-1} \). Hence if the first \( n \) dynamic signatures vanish, all higher signatures vanish too, and \( \mathbf{C}e^{\mathbf{A}t}\mathbf{v}=\mathbf{0} \) for all \( t \). Such a direction cannot be recovered from measurements.

For a selected sensor set \( \mathcal{S} \), define the dynamic measurement matrix \( \mathcal{M}_{\mathcal{S} }=\mathcal{M}(\mathbf{C}_{\mathcal{S} },\mathbf{A}) \). Algebraically, full observability requires enough independent rows to distinguish all \( n \) state directions:

\[ \operatorname{rank}(\mathcal{M}_{\mathcal{S} })=n. \]

A sensor is redundant relative to an existing placement if its rows do not increase the row span of the dynamic signatures. If \( \mathcal{S}_1\subseteq\mathcal{S}_2 \), then \( \mathcal{M}_{\mathcal{S}_1} \) is obtained by deleting rows from \( \mathcal{M}_{\mathcal{S}_2} \). Therefore,

\[ \operatorname{rank}(\mathcal{M}_{\mathcal{S}_1}) \le \operatorname{rank}(\mathcal{M}_{\mathcal{S}_2}). \]

This proves the monotonicity property: adding sensors cannot reduce ideal mathematical observability. It may, however, worsen practical design if extra sensors introduce noise, calibration drift, cost, bandwidth limits, or poor numerical conditioning in the estimator.

A useful numerical distinction is: rank asks whether reconstruction is theoretically possible, while conditioning asks whether reconstruction is reliable under measurement noise. Two placements may both have rank \( n \), but the placement with a larger smallest singular value of \( \mathcal{M}_{\mathcal{S} } \) is usually less sensitive.

\[ \|\Delta\hat{\mathbf{x} }(0)\|_2 \lesssim \frac{\|\Delta\mathbf{z}\|_2} {\sigma_{\min}(\mathcal{M}_{\mathcal{S} })}, \]

where \( \mathbf{z} \) denotes the stacked ideal output derivative data. This bound motivates sensor placement criteria that maximize \( \sigma_{\min} \) or improve the condition number after full rank is achieved.

Suppose sensors can directly measure individual states. Measuring \( x_i \) means choosing \( \mathbf{c}_i=\mathbf{e}_i^T \). If \( \mathbf{A} \) is diagonal, measuring one state gives information only about that same mode:

\[ \mathbf{A}=\operatorname{diag}(\lambda_1,\dots,\lambda_n),\quad \mathbf{e}_i^T\mathbf{A}^k = \lambda_i^k\mathbf{e}_i^T. \]

In that case, one direct sensor on \( x_i \) cannot reveal the other coordinates. However, for coupled dynamics, powers of \( \mathbf{A} \) spread information across states. For example, in a mechanical system, measuring position can reveal velocity because the equation \( \dot{x}_1=x_2 \) makes velocity appear in the time derivative of position.

The following simple second-order model shows the principle:

\[ \dot{\mathbf{x} }= \begin{bmatrix}0&1\\-k/m&-b/m\end{bmatrix}\mathbf{x},\quad y=\begin{bmatrix}1&0\end{bmatrix}\mathbf{x}. \]

Although the sensor measures only position, the first output derivative is

\[ \dot{y}= \begin{bmatrix}1&0\end{bmatrix} \begin{bmatrix}0&1\\-k/m&-b/m\end{bmatrix}\mathbf{x} = \begin{bmatrix}0&1\end{bmatrix}\mathbf{x}. \]

Thus position and velocity become visible through the pair \( y,\dot{y} \). This is the basic mechanism behind state reconstruction from fewer sensors than states.

Practical sensor placement often involves a finite set of candidate sensors and a budget constraint. Let \( s_i\in\{0,1\} \) indicate whether candidate sensor \( i \) is selected. A cardinality-constrained placement problem can be written as

\[ \max_{\mathbf{s}\in\{0,1\}^{p} } \; J(\mathbf{s}) \quad \text{subject to}\quad \sum_{i=1}^{p}s_i\le q. \]

Common choices of \( J \) include:

\[ J_{\mathrm{rank} }(\mathbf{s}) = \operatorname{rank}(\mathcal{M}_{\mathcal{S} }),\quad J_{\sigma}(\mathbf{s}) = \sigma_{\min}(\mathcal{M}_{\mathcal{S} }),\quad J_{\log\det}(\mathbf{s}) = \log\det(\mathcal{M}_{\mathcal{S} }^T\mathcal{M}_{\mathcal{S} }+\epsilon\mathbf{I}). \]

The rank criterion focuses on theoretical observability. The smallest singular value criterion focuses on worst-case sensitivity. The log-determinant criterion measures the volume of the information ellipsoid and is widely used because it rewards balanced information across state directions.

The same idea can be expressed through an output-energy matrix. For the finite derivative-based approximation:

\[ \mathbf{G}_{o,\mathcal{S} } = \mathcal{M}_{\mathcal{S} }^T\mathcal{M}_{\mathcal{S} } = \sum_{k=0}^{n-1}(\mathbf{A}^{k})^T \mathbf{C}_{\mathcal{S} }^T\mathbf{C}_{\mathcal{S} }\mathbf{A}^{k}. \]

If a nonzero state direction \( \mathbf{v} \) has very small quadratic value \( \mathbf{v}^T\mathbf{G}_{o,\mathcal{S} }\mathbf{v} \), then that direction generates weak measured signatures and will be difficult to estimate accurately.

flowchart TD
  S1["Candidate sensor rows"] --> R1["Rank objective: \nreveal all directions"]
  S1 --> R2["Sigma-min objective: \nreduce worst-case noise sensitivity"]
  S1 --> R3["Log-det objective: \nincrease information volume"]
  R1 --> P["Placement decision under budget"]
  R2 --> P
  R3 --> P
  P --> V["Validate with simulation and estimator performance"]
        

Consider a four-state model. Suppose we can place direct sensors on selected state coordinates:

\[ \mathbf{A}= \begin{bmatrix} 0&1&0&0\\ -2&-0.4&0.8&0\\ 0&0&-1&1\\ 0.6&0&-3&-0.5 \end{bmatrix}. \]

If a single sensor measures \( x_1 \), then \( \mathbf{C}=[1\;0\;0\;0] \). The dynamic measurement matrix becomes

\[ \mathcal{M}= \begin{bmatrix} 1&0&0&0\\ 0&1&0&0\\ -2&-0.4&0.8&0\\ 0.8&-1.84&-1.12&0.8 \end{bmatrix}. \]

In this example, the row span already reaches all four state directions, so a single well-placed sensor can be theoretically sufficient. However, if its smallest singular value is small, reconstruction may still be noise-sensitive. A second sensor may not be needed for rank, but it can improve conditioning and robustness.

By contrast, in a block-diagonal system with two uncoupled subsystems, a sensor in one block cannot reveal states in the other block. This is why sensor placement must account for system coupling, not only the number of sensors.

Modern control software typically constructs \( \mathcal{M} \), computes its rank or singular values, and compares alternative sensor sets. Useful libraries include:

• Python: numpy and scipy.linalg for matrix computation; python-control provides obsv(A, C).
• C++: Eigen is commonly used for SVD, QR, matrix powers, and efficient dense/sparse linear algebra.
• Java: EJML or Apache Commons Math can compute rank, SVD, and matrix products.
• MATLAB/Simulink: obsv(A,C), rank, svd, and Control System Toolbox workflows are standard; Simulink can validate estimator behavior with sensor noise blocks.
• Wolfram Mathematica: MatrixPower, MatrixRank, SingularValueList, and symbolic matrix expressions are useful for exact small examples.

The following script enumerates one- and two-state sensor placements, builds \( \mathcal{M} \), computes rank, and reports conditioning indicators.

Chapter13_Lesson4.py

# Chapter13_Lesson4.py
# Relationship of Observability to Sensor Placement
# Enumerates candidate state sensors and ranks their observability matrices.
#
# Recommended scientific stack:
#   pip install numpy scipy control
# This file uses only NumPy so it can run in a minimal environment.

import itertools
import numpy as np


def observability_matrix(A: np.ndarray, C: np.ndarray) -> np.ndarray:
    """Return O = [C; C A; ...; C A^(n-1)]."""
    n = A.shape[0]
    blocks = []
    Ak = np.eye(n)
    for _ in range(n):
        blocks.append(C @ Ak)
        Ak = Ak @ A
    return np.vstack(blocks)


def observability_score(A: np.ndarray, C: np.ndarray, tol: float = 1e-9):
    """Return rank, smallest singular value, and log-det proxy."""
    O = observability_matrix(A, C)
    singular_values = np.linalg.svd(O, compute_uv=False)
    rank = int(np.sum(singular_values > tol))
    gram = O.T @ O
    eigvals = np.linalg.eigvalsh(gram + tol * np.eye(A.shape[0]))
    logdet = float(np.sum(np.log(eigvals)))
    sigma_min = float(singular_values[-1])
    return rank, sigma_min, logdet, O


def build_C_from_state_sensors(sensor_indices, n):
    """Each selected sensor measures one state directly: y_i = x_j."""
    rows = []
    for j in sensor_indices:
        row = np.zeros(n)
        row[j] = 1.0
        rows.append(row)
    return np.vstack(rows)


def enumerate_sensor_sets(A, max_sensors=None):
    """Enumerate direct state-sensor subsets and rank each placement."""
    n = A.shape[0]
    max_sensors = n if max_sensors is None else max_sensors
    results = []
    for r in range(1, max_sensors + 1):
        for subset in itertools.combinations(range(n), r):
            C = build_C_from_state_sensors(subset, n)
            rank, sigma_min, logdet, _ = observability_score(A, C)
            results.append({
                "sensors": tuple(j + 1 for j in subset),
                "rank": rank,
                "sigma_min": sigma_min,
                "logdet": logdet
            })
    results.sort(key=lambda item: (item["rank"], item["sigma_min"], item["logdet"]), reverse=True)
    return results


def main():
    # Four-state cascade-like model.
    # Measuring x1 can reveal downstream dynamic coupling differently from measuring x4.
    A = np.array([
        [0.0, 1.0, 0.0, 0.0],
        [-2.0, -0.4, 0.8, 0.0],
        [0.0, 0.0, -1.0, 1.0],
        [0.6, 0.0, -3.0, -0.5]
    ])

    print("A =")
    print(A)
    print("\nBest direct state-sensor placements:")
    for item in enumerate_sensor_sets(A, max_sensors=2)[:10]:
        print(
            f"sensors={item['sensors']}, "
            f"rank={item['rank']}, "
            f"sigma_min={item['sigma_min']:.3e}, "
            f"logdet={item['logdet']:.3f}"
        )

    # Inspect one placement in detail.
    C = build_C_from_state_sensors([0], A.shape[0])
    rank, sigma_min, logdet, O = observability_score(A, C)
    print("\nPlacement: measure x1")
    print("O =")
    print(O)
    print(f"rank={rank}, sigma_min={sigma_min:.3e}, logdet={logdet:.3f}")


if __name__ == "__main__":
    main()

This C++ version implements the matrix operations from scratch. In a production project, replace the scratch routines with Eigen's matrix classes, QR decomposition, and singular value decomposition.

Chapter13_Lesson4.cpp

// Chapter13_Lesson4.cpp
// Relationship of Observability to Sensor Placement
// Scratch implementation: matrix powers, observability matrix, and rank by Gaussian elimination.
// Compile: g++ -std=c++17 Chapter13_Lesson4.cpp -O2 -o Chapter13_Lesson4

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

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

Matrix eye(int n) {
    Matrix I(n, std::vector<double>(n, 0.0));
    for (int i = 0; i < n; ++i) I[i][i] = 1.0;
    return I;
}

Matrix multiply(const Matrix& A, const Matrix& B) {
    int m = static_cast<int>(A.size());
    int p = static_cast<int>(A[0].size());
    int n = static_cast<int>(B[0].size());
    Matrix C(m, std::vector<double>(n, 0.0));
    for (int i = 0; i < m; ++i) {
        for (int k = 0; k < p; ++k) {
            for (int j = 0; j < n; ++j) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return C;
}

Matrix stackVertical(const std::vector<Matrix>& blocks) {
    int cols = static_cast<int>(blocks[0][0].size());
    Matrix out;
    for (const auto& B : blocks) {
        for (const auto& row : B) {
            if (static_cast<int>(row.size()) != cols) {
                throw std::runtime_error("Inconsistent column count.");
            }
            out.push_back(row);
        }
    }
    return out;
}

Matrix observabilityMatrix(const Matrix& A, const Matrix& C) {
    int n = static_cast<int>(A.size());
    Matrix Ak = eye(n);
    std::vector<Matrix> blocks;
    for (int k = 0; k < n; ++k) {
        blocks.push_back(multiply(C, Ak));
        Ak = multiply(Ak, A);
    }
    return stackVertical(blocks);
}

int rankMatrix(Matrix M, double tol = 1e-9) {
    int rows = static_cast<int>(M.size());
    int cols = static_cast<int>(M[0].size());
    int r = 0;
    for (int c = 0; c < cols && r < rows; ++c) {
        int pivot = r;
        for (int i = r + 1; i < rows; ++i) {
            if (std::fabs(M[i][c]) > std::fabs(M[pivot][c])) pivot = i;
        }
        if (std::fabs(M[pivot][c]) <= tol) continue;
        std::swap(M[r], M[pivot]);
        double piv = M[r][c];
        for (int j = c; j < cols; ++j) M[r][j] /= piv;
        for (int i = 0; i < rows; ++i) {
            if (i == r) continue;
            double factor = M[i][c];
            for (int j = c; j < cols; ++j) M[i][j] -= factor * M[r][j];
        }
        ++r;
    }
    return r;
}

Matrix CFromMask(int mask, int n) {
    Matrix C;
    for (int j = 0; j < n; ++j) {
        if (mask & (1 << j)) {
            std::vector<double> row(n, 0.0);
            row[j] = 1.0;
            C.push_back(row);
        }
    }
    return C;
}

int sensorCount(int mask) {
    int count = 0;
    while (mask) {
        count += mask & 1;
        mask >>= 1;
    }
    return count;
}

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

    int n = static_cast<int>(A.size());
    std::cout << "Sensor placement by rank of O = [C; C A; ...; C A^(n-1)]\n";

    for (int mask = 1; mask < (1 << n); ++mask) {
        if (sensorCount(mask) > 2) continue;
        Matrix C = CFromMask(mask, n);
        Matrix O = observabilityMatrix(A, C);
        int r = rankMatrix(O);
        std::cout << "Sensors {";
        bool first = true;
        for (int j = 0; j < n; ++j) {
            if (mask & (1 << j)) {
                if (!first) std::cout << ",";
                std::cout << (j + 1);
                first = false;
            }
        }
        std::cout << "} rank = " << r << "\n";
    }

    return 0;
}

This Java implementation uses pure arrays. For larger systems, EJML or Apache Commons Math is preferable for stable rank and SVD computations.

Chapter13_Lesson4.java

// Chapter13_Lesson4.java
// Relationship of Observability to Sensor Placement
// Scratch implementation: observability matrix and numerical rank.
// Compile: javac Chapter13_Lesson4.java
// Run:     java Chapter13_Lesson4

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

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

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

    static double[][] stack(List<double[][]> blocks) {
        int cols = blocks.get(0)[0].length;
        int rows = 0;
        for (double[][] block : blocks) rows += block.length;
        double[][] out = new double[rows][cols];
        int r = 0;
        for (double[][] block : blocks) {
            for (double[] row : block) {
                out[r++] = Arrays.copyOf(row, cols);
            }
        }
        return out;
    }

    static double[][] observabilityMatrix(double[][] A, double[][] C) {
        int n = A.length;
        double[][] Ak = eye(n);
        List<double[][]> blocks = new ArrayList<>();
        for (int k = 0; k < n; k++) {
            blocks.add(multiply(C, Ak));
            Ak = multiply(Ak, A);
        }
        return stack(blocks);
    }

    static int rank(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++) M[i] = Arrays.copyOf(input[i], 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(M[i][c]) > Math.abs(M[pivot][c])) pivot = i;
            }
            if (Math.abs(M[pivot][c]) <= tol) continue;

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

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

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

    static double[][] cFromMask(int mask, int n) {
        int count = Integer.bitCount(mask);
        double[][] C = new double[count][n];
        int row = 0;
        for (int j = 0; j < n; j++) {
            if ((mask & (1 << j)) != 0) {
                C[row][j] = 1.0;
                row++;
            }
        }
        return C;
    }

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

        int n = A.length;
        System.out.println("Sensor placement by rank of O = [C; C A; ...; C A^(n-1)]");

        for (int mask = 1; mask < (1 << n); mask++) {
            if (Integer.bitCount(mask) > 2) continue;
            double[][] C = cFromMask(mask, n);
            double[][] O = observabilityMatrix(A, C);
            int r = rank(O, 1e-9);

            System.out.print("Sensors {");
            boolean first = true;
            for (int j = 0; j < n; j++) {
                if ((mask & (1 << j)) != 0) {
                    if (!first) System.out.print(",");
                    System.out.print(j + 1);
                    first = false;
                }
            }
            System.out.println("} rank = " + r);
        }
    }
}

MATLAB users can compute the same quantities with basic matrix operations or with obsv(A,C) from Control System Toolbox. In Simulink, the chosen sensor matrix can be implemented as a matrix gain from the state vector to the measured output, followed by noise and estimator blocks.

Chapter13_Lesson4.m

% Chapter13_Lesson4.m
% Relationship of Observability to Sensor Placement
% Uses MATLAB matrix operations. If Control System Toolbox is available,
% obsv(A,C) gives the same stacked observability matrix.

clear; clc;

A = [ 0.0  1.0  0.0  0.0;
     -2.0 -0.4  0.8  0.0;
      0.0  0.0 -1.0  1.0;
      0.6  0.0 -3.0 -0.5 ];

n = size(A,1);

fprintf("Sensor placement by rank of O = [C; C A; ...; C A^(n-1)]\n");

for r = 1:2
    subsets = nchoosek(1:n, r);
    for idx = 1:size(subsets,1)
        sensors = subsets(idx,:);
        C = zeros(r,n);
        for k = 1:r
            C(k,sensors(k)) = 1;
        end

        O = local_observability_matrix(A,C);
        rk = rank(O);
        G = O' * O;
        sigma_min = min(svd(O));
        logdet_proxy = log(det(G + 1e-9*eye(n)));

        fprintf("Sensors { %s }: rank=%d, sigma_min=%.3e, logdet=%.3f\n", ...
            num2str(sensors), rk, sigma_min, logdet_proxy);
    end
end

% Compare with MATLAB Control System Toolbox, if available:
% C = [1 0 0 0];
% O_toolbox = obsv(A,C);

function O = local_observability_matrix(A,C)
    n = size(A,1);
    O = [];
    Ak = eye(n);
    for k = 1:n
        O = [O; C*Ak]; %#ok<AGROW>
        Ak = Ak*A;
    end
end

Mathematica is useful for symbolic observability analysis when the entries of \( \mathbf{A} \) depend on parameters.

Chapter13_Lesson4.nb

(* Chapter13_Lesson4.nb *)
(* Relationship of Observability to Sensor Placement *)
(* Paste this into a Wolfram Mathematica notebook or save as a .wl script. *)

ClearAll[observabilityMatrix, cFromSensors, scorePlacement];

observabilityMatrix[A_, C_] := Module[
  {n = Length[A]},
  Join @@ Table[C . MatrixPower[A, k], {k, 0, n - 1}]
];

cFromSensors[sensors_, n_] := Module[
  {C = ConstantArray[0, {Length[sensors], n}]},
  Do[C[[i, sensors[[i]]]] = 1, {i, Length[sensors]}];
  C
];

scorePlacement[A_, sensors_] := Module[
  {n = Length[A], C, O, gram, sv},
  C = cFromSensors[sensors, n];
  O = observabilityMatrix[A, C];
  gram = Transpose[O].O;
  sv = SingularValueList[O];
  <|
    "Sensors" -> sensors,
    "Rank" -> MatrixRank[O],
    "SmallestSingularValue" -> Min[sv],
    "LogDetProxy" -> Log[Det[gram + 10^-9 IdentityMatrix[n]]]
  |>
];

A = {
  {0.0, 1.0, 0.0, 0.0},
  {-2.0, -0.4, 0.8, 0.0},
  {0.0, 0.0, -1.0, 1.0},
  {0.6, 0.0, -3.0, -0.5}
};

placements = Join @@ Table[Subsets[Range[Length[A]], {r}], {r, 1, 2}];
scores = scorePlacement[A, #] & /@ placements;

SortBy[scores, {-#Rank, -#SmallestSingularValue, -#LogDetProxy} &] // Dataset

(* Example: inspect the observability matrix for measuring x1 only. *)
C1 = cFromSensors[{1}, Length[A]];
O1 = observabilityMatrix[A, C1];
MatrixForm[O1]
MatrixRank[O1]

Problem 1 (Position Sensor for a Mass-Spring-Damper): Consider \( \dot{\mathbf{x} }=\begin{bmatrix}0&1\\-k/m&-b/m\end{bmatrix}\mathbf{x} \) and \( y=x_1 \). Show that the placement measuring position makes both position and velocity visible.

Solution: The sensor matrix is \( \mathbf{C}=[1\;0] \). The first two dynamic signatures are

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

Stacking them gives the identity matrix, so the two state coordinates are distinguishable from \( y \) and \( \dot{y} \).

Problem 2 (Uncoupled Blocks): Let \( \mathbf{A}=\operatorname{diag}(\mathbf{A}_1,\mathbf{A}_2) \) and \( \mathbf{C}=[\mathbf{C}_1\;\mathbf{0}] \). Prove that states in the second block are invisible.

Solution: Since powers of a block-diagonal matrix remain block diagonal,

\[ \mathbf{A}^{k}= \begin{bmatrix} \mathbf{A}_1^k&\mathbf{0}\\ \mathbf{0}&\mathbf{A}_2^k \end{bmatrix},\quad \mathbf{C}\mathbf{A}^{k} = \begin{bmatrix} \mathbf{C}_1\mathbf{A}_1^k&\mathbf{0} \end{bmatrix}. \]

Therefore any state vector \( \mathbf{v}=[\mathbf{0}^T\;\mathbf{v}_2^T]^T \) produces \( \mathbf{C}\mathbf{A}^{k}\mathbf{v}=\mathbf{0} \) for every \( k \). The second block cannot be reconstructed from sensors placed only on the first block.

Problem 3 (Monotonicity of Sensor Addition): Suppose \( \mathcal{S}_1\subseteq\mathcal{S}_2 \). Prove that adding sensors cannot reduce the rank of the dynamic measurement matrix.

Solution: The matrix \( \mathcal{M}_{\mathcal{S}_1} \) is formed by keeping only some rows of \( \mathcal{M}_{\mathcal{S}_2} \). Row deletion cannot increase rank; equivalently, row addition cannot decrease rank. Hence

\[ \operatorname{rank}(\mathcal{M}_{\mathcal{S}_1}) \le \operatorname{rank}(\mathcal{M}_{\mathcal{S}_2}). \]

Problem 4 (Sensor Redundancy): Let two candidate sensors satisfy \( \mathbf{c}_2=\alpha\mathbf{c}_1 \) for some nonzero scalar \( \alpha \). Show that the second sensor is redundant in the ideal noiseless rank sense.

Solution: For every \( k \), \( \mathbf{c}_2\mathbf{A}^k=\alpha\mathbf{c}_1\mathbf{A}^k \). Therefore all dynamic rows generated by sensor 2 lie in the span of the dynamic rows generated by sensor 1. Adding sensor 2 does not change the row span and cannot increase rank.

Problem 5 (Choosing Between Two Full-Rank Placements): Two placements both satisfy \( \operatorname{rank}(\mathcal{M})=n \). Placement A has \( \sigma_{\min}(\mathcal{M})=10^{-3} \), and placement B has \( \sigma_{\min}(\mathcal{M})=10^{-1} \). Which one is preferable for noisy measurements?

Solution: Placement B is preferable because its smallest singular value is larger. The reconstruction sensitivity bound scales approximately like \( 1/\sigma_{\min}(\mathcal{M}) \). Thus placement A may amplify noise roughly one hundred times more than placement B in the worst direction.

Sensor placement determines the output matrix \( \mathbf{C} \), but observability depends on the interaction between \( \mathbf{C} \) and the state dynamics \( \mathbf{A} \). A sensor can reveal states it does not directly measure when dynamic coupling carries those state directions into output derivatives. The mathematical goal is to choose sensors that eliminate nonzero unobservable directions, while the engineering goal is to do so with good conditioning, acceptable cost, noise robustness, and implementable hardware.