Suppose a plant model is already available in state-space form and the engineer must decide which state components, linear combinations, or physical locations should be measured. The central question is: which sensors make the initial state easiest to infer from output data? The answer is not determined only by the algebraic rank of the observability matrix. Rank tells us whether reconstruction is possible in exact arithmetic; the observability Gramian tells us how much output energy each initial-state direction produces.

\[ \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t),\qquad \mathbf{y}_S(t)=\mathbf{C}_S\mathbf{x}(t),\qquad 0\le t\le T. \]

The subscript \( S \) denotes the selected set of sensors. A candidate sensor is represented by a row vector \( \mathbf{c}_i \), and the selected output matrix is formed by stacking the chosen rows:

\[ \mathbf{C}_S = \begin{bmatrix} \mathbf{c}_{i_1} \\ \mathbf{c}_{i_2} \\ \vdots \\ \mathbf{c}_{i_p} \end{bmatrix},\qquad S=\{i_1,i_2,\dots,i_p\}. \]

flowchart TD
  A["State matrix A"] --> B["Candidate sensor rows c1, c2, ..., cm"]
  B --> C["Choose sensor set S"]
  C --> D["Build output matrix C_S"]
  D --> E["Compute finite-horizon Gramian Wo(S,T)"]
  E --> F["Inspect trace, logdet, min eigenvalue, condition"]
  F --> G["Qualitative placement decision"]
        

For a selected sensor set \( S \), the finite-horizon observability Gramian is

\[ \mathbf{W}_o(S;T)=\int_0^T e^{\mathbf{A}^Tt}\mathbf{C}_S^T\mathbf{C}_S e^{\mathbf{A}t}\,dt. \]

The associated output energy generated by an initial state \( \mathbf{x}_0 \) is

\[ E_S(\mathbf{x}_0;T)=\int_0^T \mathbf{y}_S^T(t)\mathbf{y}_S(t)\,dt =\mathbf{x}_0^T\mathbf{W}_o(S;T)\mathbf{x}_0. \]

Therefore, sensor placement can be interpreted as shaping the quadratic form that maps initial states to measurable output energy. If \( E_S(\mathbf{x}_0;T) \) is small for some nonzero direction, then that direction is weakly visible over the measurement interval.

Proof of the energy identity.

\[ \begin{aligned} \mathbf{y}_S(t) &=\mathbf{C}_S e^{\mathbf{A}t}\mathbf{x}_0,\\ \mathbf{y}_S^T(t)\mathbf{y}_S(t) &=\mathbf{x}_0^T e^{\mathbf{A}^Tt} \mathbf{C}_S^T\mathbf{C}_S e^{\mathbf{A}t} \mathbf{x}_0. \end{aligned} \]

Integrating the second line from \( 0 \) to \( T \) gives the Gramian expression above.

A useful property for placement intuition is additivity. Since \( \mathbf{C}_S \) is obtained by stacking rows,

\[ \mathbf{C}_S^T\mathbf{C}_S = \sum_{i\in S} \mathbf{c}_i^T\mathbf{c}_i. \]

Hence the selected Gramian decomposes as a sum of individual sensor contributions:

\[ \mathbf{W}_o(S;T)=\sum_{i\in S}\mathbf{W}_i(T),\qquad \mathbf{W}_i(T)=\int_0^T e^{\mathbf{A}^Tt} \mathbf{c}_i^T\mathbf{c}_i e^{\mathbf{A}t}\,dt. \]

This does not mean that sensors are independent in performance. The matrices add, but the scalar quality measure may not. For example, adding a second sensor that observes the same already-visible direction may increase \( \operatorname{tr}(\mathbf{W}_o) \) while doing little for the weakest direction.

Monotonicity in the Loewner order. If \( S_1\subseteq S_2 \), then

\[ \mathbf{W}_o(S_2;T)-\mathbf{W}_o(S_1;T) =\sum_{i\in S_2\setminus S_1}\mathbf{W}_i(T)\succeq 0. \]

Thus adding sensors cannot reduce output energy in any initial-state direction. The practical issue is not whether output energy increases; it is whether it increases in the directions that were previously weak.

Let \( \lambda_1\le\lambda_2\le\cdots\le\lambda_n \) be the eigenvalues of \( \mathbf{W}_o(S;T) \). Several scalar summaries are useful:

\[ \operatorname{tr}(\mathbf{W}_o)=\sum_{i=1}^n\lambda_i. \]

The trace measures total output energy over orthonormal state directions. It is simple and additive, but it can prefer sensors that make already-visible directions even more visible.

\[ \det(\mathbf{W}_o)=\prod_{i=1}^n\lambda_i,\qquad \log\det(\mathbf{W}_o)=\sum_{i=1}^n\log(\lambda_i). \]

The determinant or log-determinant measures an ellipsoid-volume effect. It strongly penalizes near-zero eigenvalues, so it tends to favor more balanced visibility across state directions.

\[ \lambda_{\min}(\mathbf{W}_o)=\lambda_1,\qquad \kappa_2(\mathbf{W}_o)= \frac{\lambda_{\max}(\mathbf{W}_o)}{\lambda_{\min}(\mathbf{W}_o)}. \]

The smallest eigenvalue describes the weakest observable direction. The condition number describes anisotropy: a large value means that reconstruction will be much more sensitive in some directions than in others.

Because the Gramian is symmetric positive semidefinite, it admits an orthonormal eigen-decomposition

\[ \mathbf{W}_o=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^T, \qquad \boldsymbol{\Lambda}=\operatorname{diag}(\lambda_1,\dots,\lambda_n). \]

If \( \mathbf{x}_0=\mathbf{V}\boldsymbol{\alpha} \), then

\[ E_S(\mathbf{x}_0;T) =\boldsymbol{\alpha}^T\boldsymbol{\Lambda}\boldsymbol{\alpha} =\sum_{i=1}^n\lambda_i\alpha_i^2. \]

Thus the eigenvectors of \( \mathbf{W}_o \) are the principal directions of output-energy sensitivity. Small eigenvalues correspond to initial-state combinations that create little output energy and are therefore difficult to distinguish from measurement noise.

\[ \mathcal{E}_\rho(S)= \left\{\mathbf{x}_0: \mathbf{x}_0^T\mathbf{W}_o(S;T)\mathbf{x}_0\le\rho^2 \right\}. \]

The set \( \mathcal{E}_\rho(S) \) contains initial states whose measured output energy is no larger than \( \rho^2 \). Long axes of this ellipsoid indicate weakly visible state directions.

The Gramian suggests several placement heuristics. These rules are not a substitute for a full estimation design, but they are valuable early engineering diagnostics:

• Prefer sensors that reduce the longest axes of the output-energy ellipsoid, not only sensors that increase total energy.
• Compare \( \lambda_{\min} \) and \( \kappa_2 \) when reconstruction sensitivity is important.
• Use \( \operatorname{tr}(\mathbf{W}_o) \) for a quick total-energy score, but check whether it ignores weak modes.
• Use regularized \( \log\det \) when the Gramian is nearly singular or when numerical robustness matters.
• Inspect the eigenvector associated with \( \lambda_{\min} \); it identifies the state combination that remains hardest to observe.

flowchart TD
  A["Compute Wo for each candidate sensor set"] --> B["Is rank full?"]
  B -->|no| C["Reject as exactly unobservable over horizon"]
  B -->|yes| D["Inspect smallest eigenvalue"]
  D --> E["Inspect condition number"]
  E --> F["Compare trace and logdet"]
  F --> G["Choose sensors that balance energy \nacross directions"]
  C --> H["Add sensors that see missing directions"]
  H --> A
        

Proposition 1: positive semidefiniteness. For every selected sensor set \( S \),

\[ \mathbf{z}^T\mathbf{W}_o(S;T)\mathbf{z} =\int_0^T \left\|\mathbf{C}_S e^{\mathbf{A}t}\mathbf{z}\right\|_2^2\,dt \ge 0. \]

Hence \( \mathbf{W}_o(S;T)\succeq 0 \). It is positive definite exactly when no nonzero initial-state direction is invisible over the interval.

Proposition 2: exact observability over the horizon. The pair \( (\mathbf{A},\mathbf{C}_S) \) is observable over \( [0,T] \) if and only if

\[ \operatorname{rank}\mathbf{W}_o(S;T)=n. \]

If the rank is deficient, there exists a nonzero \( \mathbf{z} \) such that \( \mathbf{C}_S e^{\mathbf{A}t}\mathbf{z}=\mathbf{0} \) throughout the interval, so the corresponding initial-state component produces no output.

Proposition 3: sensor redundancy can be invisible to rank. Two sensor sets may both yield full rank but have very different conditioning:

\[ \operatorname{rank}\mathbf{W}_o(S_a;T) =\operatorname{rank}\mathbf{W}_o(S_b;T)=n, \qquad \kappa_2(\mathbf{W}_o(S_a;T))\gg \kappa_2(\mathbf{W}_o(S_b;T)). \]

In this case both sensor sets are theoretically observable, but \( S_b \) is qualitatively better for numerical state reconstruction.

For practical computation, the finite-horizon Gramian may be computed by numerical quadrature or by integrating the matrix differential equation

\[ \dot{\mathbf{W} }(t)= \mathbf{A}^T\mathbf{W}(t)+\mathbf{W}(t)\mathbf{A} +\mathbf{C}^T\mathbf{C},\qquad \mathbf{W}(0)=\mathbf{0}. \]

If \( \mathbf{A} \) is asymptotically stable, the infinite-horizon Gramian is the unique positive semidefinite solution of

\[ \mathbf{A}^T\mathbf{W}_o+\mathbf{W}_o\mathbf{A} +\mathbf{C}^T\mathbf{C}=\mathbf{0}. \]

Useful libraries include NumPy/SciPy and python-control in Python, Eigen or Armadillo in C++, EJML or Apache Commons Math in Java, Control System Toolbox and Simulink in MATLAB, and built-in linear algebra functions in Wolfram Mathematica. The code below deliberately includes from-scratch finite-horizon computation for teaching clarity.

Chapter15_Lesson3.py

# Chapter15_Lesson3.py
# Finite-horizon observability Gramian and qualitative sensor placement.
#
# Libraries:
#   numpy: numerical arrays and eigenvalue routines.
# Optional extensions:
#   scipy.linalg.solve_continuous_lyapunov for infinite-horizon stable systems;
#   python-control for state-space model handling.

from __future__ import annotations

import itertools
import math
from typing import Iterable, List, Sequence, Tuple

import numpy as np


def gramian_rhs(A: np.ndarray, W: np.ndarray, Q: np.ndarray) -> np.ndarray:
    """Right-hand side dW/dt = A.T W + W A + C.T C."""
    return A.T @ W + W @ A + Q


def finite_horizon_observability_gramian(
    A: np.ndarray,
    C: np.ndarray,
    T: float = 6.0,
    steps: int = 4000,
) -> np.ndarray:
    """
    Compute W_o(0,T) = integral_0^T exp(A.T t) C.T C exp(A t) dt
    through the matrix differential equation.

    This implementation uses RK4 and does not require SciPy.
    """
    n = A.shape[0]
    W = np.zeros((n, n), dtype=float)
    Q = C.T @ C
    h = T / steps

    for _ in range(steps):
        k1 = gramian_rhs(A, W, Q)
        k2 = gramian_rhs(A, W + 0.5 * h * k1, Q)
        k3 = gramian_rhs(A, W + 0.5 * h * k2, Q)
        k4 = gramian_rhs(A, W + h * k3, Q)
        W = W + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)

    return 0.5 * (W + W.T)


def row_sensor(index: int, n: int) -> np.ndarray:
    """Return a 1-by-n sensor row measuring x_index."""
    c = np.zeros((1, n), dtype=float)
    c[0, index] = 1.0
    return c


def stack_sensors(sensor_indices: Sequence[int], n: int) -> np.ndarray:
    """Build a C matrix from selected coordinate sensors."""
    return np.vstack([row_sensor(i, n) for i in sensor_indices])


def gramian_metrics(W: np.ndarray, eps: float = 1e-10) -> dict:
    """
    Qualitative scalar summaries:
      trace  -> total output energy averaged over coordinate directions
      logdet -> volume of the output-energy ellipsoid
      lambda_min -> weakest observable direction
      cond   -> anisotropy / numerical sensitivity
    """
    eigvals = np.linalg.eigvalsh(W)
    eigvals_clipped = np.maximum(eigvals, eps)

    rank = int(np.sum(eigvals > 1e-8))
    cond = math.inf if eigvals[0] <= eps else float(eigvals[-1] / eigvals[0])
    logdet = float(np.sum(np.log(eigvals_clipped)))

    return {
        "trace": float(np.trace(W)),
        "logdet_eps": logdet,
        "lambda_min": float(eigvals[0]),
        "lambda_max": float(eigvals[-1]),
        "rank": rank,
        "condition": cond,
        "eigvals": eigvals,
    }


def score_sensor_sets(
    A: np.ndarray,
    candidates: Sequence[int],
    k: int,
    T: float = 6.0,
) -> List[Tuple[Tuple[int, ...], dict]]:
    """Exhaustively score all k-sensor coordinate choices."""
    n = A.shape[0]
    scored = []
    for sensor_set in itertools.combinations(candidates, k):
        C = stack_sensors(sensor_set, n)
        W = finite_horizon_observability_gramian(A, C, T=T)
        scored.append((sensor_set, gramian_metrics(W)))
    scored.sort(key=lambda item: item[1]["logdet_eps"], reverse=True)
    return scored


def greedy_logdet_selection(
    A: np.ndarray,
    candidates: Sequence[int],
    k: int,
    T: float = 6.0,
    regularization: float = 1e-8,
) -> List[int]:
    """
    Greedy qualitative placement rule:
    at each step, add the sensor that gives the largest regularized log-det.
    """
    n = A.shape[0]
    selected: List[int] = []
    remaining = list(candidates)

    for step in range(k):
        best_sensor = None
        best_score = -math.inf

        for sensor in remaining:
            trial = selected + [sensor]
            C_trial = stack_sensors(trial, n)
            W_trial = finite_horizon_observability_gramian(A, C_trial, T=T)
            eigvals = np.linalg.eigvalsh(W_trial + regularization * np.eye(n))
            score = float(np.sum(np.log(eigvals)))

            if score > best_score:
                best_score = score
                best_sensor = sensor

        selected.append(best_sensor)
        remaining.remove(best_sensor)
        print(f"Step {step + 1}: add sensor x{best_sensor + 1}; logdet={best_score:.6f}")

    return selected


def main() -> None:
    # A three-state coupled stable system.
    A = np.array(
        [
            [0.0, 1.0, 0.0],
            [-2.0, -0.45, 0.8],
            [0.0, -0.7, -1.25],
        ],
        dtype=float,
    )

    candidates = [0, 1, 2]
    T = 6.0

    print("Single-sensor qualitative scores:")
    for sensor in candidates:
        C = stack_sensors([sensor], A.shape[0])
        W = finite_horizon_observability_gramian(A, C, T=T)
        metrics = gramian_metrics(W)
        print(
            f"x{sensor + 1}: trace={metrics['trace']:.4f}, "
            f"logdet_eps={metrics['logdet_eps']:.4f}, "
            f"lambda_min={metrics['lambda_min']:.4e}, "
            f"rank={metrics['rank']}, cond={metrics['condition']:.4e}"
        )

    print("\nBest two-sensor sets by regularized log-det:")
    for sensor_set, metrics in score_sensor_sets(A, candidates, k=2, T=T):
        names = ", ".join(f"x{i + 1}" for i in sensor_set)
        print(
            f"{ { {names} } }: trace={metrics['trace']:.4f}, "
            f"logdet_eps={metrics['logdet_eps']:.4f}, "
            f"lambda_min={metrics['lambda_min']:.4e}, "
            f"rank={metrics['rank']}, cond={metrics['condition']:.4e}"
        )

    print("\nGreedy two-sensor choice:")
    greedy = greedy_logdet_selection(A, candidates, k=2, T=T)
    print("Selected:", [f"x{i + 1}" for i in greedy])


if __name__ == "__main__":
    main()

      

Chapter15_Lesson3.cpp

// Chapter15_Lesson3.cpp
// Finite-horizon observability Gramian and qualitative coordinate-sensor scoring.
// From-scratch implementation for small teaching examples.
//
// Compile:
//   g++ -std=c++17 Chapter15_Lesson3.cpp -O2 -o Chapter15_Lesson3

#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#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 transpose(const Matrix& A) {
    int r = static_cast<int>(A.size());
    int c = static_cast<int>(A[0].size());
    Matrix T = zeros(c, r);
    for (int i = 0; i < r; ++i)
        for (int j = 0; j < c; ++j)
            T[j][i] = A[i][j];
    return T;
}

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

Matrix scale(const Matrix& A, double s) {
    Matrix B = A;
    for (auto& row : B)
        for (double& v : row)
            v *= s;
    return B;
}

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

Matrix sensorMatrix(const std::vector<int>& sensors, int n) {
    Matrix C = zeros(static_cast<int>(sensors.size()), n);
    for (int r = 0; r < static_cast<int>(sensors.size()); ++r)
        C[r][sensors[r]] = 1.0;
    return C;
}

Matrix gramianRhs(const Matrix& A, const Matrix& W, const Matrix& Q) {
    Matrix AT = transpose(A);
    return add(add(multiply(AT, W), multiply(W, A)), Q);
}

Matrix finiteHorizonGramian(const Matrix& A, const Matrix& C, double T, int steps) {
    int n = static_cast<int>(A.size());
    Matrix W = zeros(n, n);
    Matrix Q = multiply(transpose(C), C);
    double h = T / static_cast<double>(steps);

    for (int s = 0; s < steps; ++s) {
        Matrix k1 = gramianRhs(A, W, Q);
        Matrix k2 = gramianRhs(A, add(W, scale(k1, 0.5 * h)), Q);
        Matrix k3 = gramianRhs(A, add(W, scale(k2, 0.5 * h)), Q);
        Matrix k4 = gramianRhs(A, add(W, scale(k3, h)), Q);

        Matrix increment = add(add(k1, scale(k2, 2.0)), add(scale(k3, 2.0), k4));
        W = add(W, scale(increment, h / 6.0));
    }

    for (int i = 0; i < n; ++i)
        for (int j = i + 1; j < n; ++j) {
            double avg = 0.5 * (W[i][j] + W[j][i]);
            W[i][j] = avg;
            W[j][i] = avg;
        }

    return W;
}

double trace(const Matrix& A) {
    double s = 0.0;
    for (int i = 0; i < static_cast<int>(A.size()); ++i) s += A[i][i];
    return s;
}

double det3(const Matrix& A) {
    return A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1])
         - A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0])
         + A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]);
}

double regularizedLogDet3(const Matrix& W, double eps) {
    Matrix R = W;
    for (int i = 0; i < 3; ++i) R[i][i] += eps;
    double d = det3(R);
    if (d <= 0.0) return -std::numeric_limits<double>::infinity();
    return std::log(d);
}

void printScore(const std::string& name, const Matrix& W) {
    std::cout << std::setw(8) << name
              << "  trace=" << std::setw(12) << trace(W)
              << "  logdet_eps=" << std::setw(12) << regularizedLogDet3(W, 1e-8)
              << "  det_eps=" << std::setw(12) << std::exp(regularizedLogDet3(W, 1e-8))
              << "\n";
}

int main() {
    Matrix A = {
        {0.0, 1.0, 0.0},
        {-2.0, -0.45, 0.8},
        {0.0, -0.7, -1.25}
    };

    const int n = 3;
    const double T = 6.0;
    const int steps = 4000;

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Single-sensor scores\n";
    for (int i = 0; i < n; ++i) {
        Matrix C = sensorMatrix({i}, n);
        Matrix W = finiteHorizonGramian(A, C, T, steps);
        printScore("x" + std::to_string(i + 1), W);
    }

    std::cout << "\nTwo-sensor scores\n";
    double bestScore = -std::numeric_limits<double>::infinity();
    std::string bestName;
    for (int i = 0; i < n; ++i) {
        for (int j = i + 1; j < n; ++j) {
            Matrix C = sensorMatrix({i, j}, n);
            Matrix W = finiteHorizonGramian(A, C, T, steps);
            std::string name = "x" + std::to_string(i + 1) + ",x" + std::to_string(j + 1);
            double score = regularizedLogDet3(W, 1e-8);
            printScore(name, W);
            if (score > bestScore) {
                bestScore = score;
                bestName = name;
            }
        }
    }

    std::cout << "\nBest qualitative two-sensor set by log-det: " << bestName << "\n";
    return 0;
}

      

Chapter15_Lesson3.java

// Chapter15_Lesson3.java
// Finite-horizon observability Gramian and qualitative coordinate-sensor scoring.
// From-scratch implementation for small teaching examples.
//
// Compile and run:
//   javac Chapter15_Lesson3.java
//   java Chapter15_Lesson3

import java.util.Arrays;

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

    static double[][] transpose(double[][] A) {
        int r = A.length;
        int c = A[0].length;
        double[][] T = zeros(c, r);
        for (int i = 0; i < r; i++)
            for (int j = 0; j < c; j++)
                T[j][i] = A[i][j];
        return T;
    }

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

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

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

    static double[][] sensorMatrix(int[] sensors, int n) {
        double[][] C = zeros(sensors.length, n);
        for (int r = 0; r < sensors.length; r++) {
            C[r][sensors[r]] = 1.0;
        }
        return C;
    }

    static double[][] gramianRhs(double[][] A, double[][] W, double[][] Q) {
        double[][] AT = transpose(A);
        return add(add(multiply(AT, W), multiply(W, A)), Q);
    }

    static double[][] finiteHorizonGramian(double[][] A, double[][] C, double T, int steps) {
        int n = A.length;
        double[][] W = zeros(n, n);
        double[][] Q = multiply(transpose(C), C);
        double h = T / steps;

        for (int s = 0; s < steps; s++) {
            double[][] k1 = gramianRhs(A, W, Q);
            double[][] k2 = gramianRhs(A, add(W, scale(k1, 0.5 * h)), Q);
            double[][] k3 = gramianRhs(A, add(W, scale(k2, 0.5 * h)), Q);
            double[][] k4 = gramianRhs(A, add(W, scale(k3, h)), Q);

            double[][] increment = add(add(k1, scale(k2, 2.0)), add(scale(k3, 2.0), k4));
            W = add(W, scale(increment, h / 6.0));
        }

        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j < n; j++) {
                double avg = 0.5 * (W[i][j] + W[j][i]);
                W[i][j] = avg;
                W[j][i] = avg;
            }
        }

        return W;
    }

    static double trace(double[][] A) {
        double s = 0.0;
        for (int i = 0; i < A.length; i++) s += A[i][i];
        return s;
    }

    static double det3(double[][] A) {
        return A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1])
             - A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0])
             + A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]);
    }

    static double regularizedLogDet3(double[][] W, double eps) {
        double[][] R = new double[3][3];
        for (int i = 0; i < 3; i++) R[i] = Arrays.copyOf(W[i], 3);
        for (int i = 0; i < 3; i++) R[i][i] += eps;
        double d = det3(R);
        if (d <= 0.0) return Double.NEGATIVE_INFINITY;
        return Math.log(d);
    }

    static void printScore(String name, double[][] W) {
        System.out.printf(
            "%8s  trace=%12.6f  logdet_eps=%12.6f  det_eps=%12.6e%n",
            name, trace(W), regularizedLogDet3(W, 1e-8),
            Math.exp(regularizedLogDet3(W, 1e-8))
        );
    }

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

        int n = 3;
        double T = 6.0;
        int steps = 4000;

        System.out.println("Single-sensor scores");
        for (int i = 0; i < n; i++) {
            double[][] C = sensorMatrix(new int[] {i}, n);
            double[][] W = finiteHorizonGramian(A, C, T, steps);
            printScore("x" + (i + 1), W);
        }

        System.out.println("\nTwo-sensor scores");
        double bestScore = Double.NEGATIVE_INFINITY;
        String bestName = "";
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j < n; j++) {
                double[][] C = sensorMatrix(new int[] {i, j}, n);
                double[][] W = finiteHorizonGramian(A, C, T, steps);
                String name = "x" + (i + 1) + ",x" + (j + 1);
                double score = regularizedLogDet3(W, 1e-8);
                printScore(name, W);
                if (score > bestScore) {
                    bestScore = score;
                    bestName = name;
                }
            }
        }

        System.out.println("\nBest qualitative two-sensor set by log-det: " + bestName);
    }
}

      

Chapter15_Lesson3.m

% Chapter15_Lesson3.m
% Finite-horizon observability Gramian and qualitative sensor placement.
%
% Related MATLAB/Simulink tools:
%   Control System Toolbox: ss, obsv, gram, lyap
%   Simulink: State-Space block, To Workspace block, Linear Analysis Tool
%
% This script uses finite-horizon quadrature so that the teaching example
% does not depend on asymptotic stability.

clear; clc;

A = [ 0.0   1.0   0.0;
     -2.0  -0.45  0.8;
      0.0  -0.7  -1.25 ];

n = size(A,1);
T = 6.0;
N = 2000;
t = linspace(0,T,N);
dt = t(2)-t(1);

candidateNames = ["x1","x2","x3"];

fprintf("Single-sensor scores\n");
for i = 1:n
    C = zeros(1,n);
    C(i) = 1.0;
    W = finiteHorizonWo(A,C,t,dt);
    printMetrics(candidateNames(i), W);
end

fprintf("\nTwo-sensor scores\n");
bestScore = -Inf;
bestName = "";
for i = 1:n
    for j = i+1:n
        C = zeros(2,n);
        C(1,i) = 1.0;
        C(2,j) = 1.0;
        W = finiteHorizonWo(A,C,t,dt);
        score = regularizedLogDet(W,1e-8);
        name = candidateNames(i) + "," + candidateNames(j);
        printMetrics(name, W);
        if score > bestScore
            bestScore = score;
            bestName = name;
        end
    end
end
fprintf("\nBest qualitative two-sensor set by log-det: %s\n", bestName);

% Infinite-horizon alternative for asymptotically stable A:
% C = [1 0 0; 0 1 0];
% Wo_inf = lyap(A', C'*C);  % solves A'*Wo + Wo*A + C'*C = 0

% Simulink note:
% Build a model with a State-Space block using matrices A, B, C, D.
% For candidate sensors, change C or route selected states through a Selector
% block, log y(t) with To Workspace, and numerically integrate y(t)'*y(t)
% for several initial conditions to compare output energy.

function W = finiteHorizonWo(A,C,t,dt)
    n = size(A,1);
    W = zeros(n,n);
    Q = C'*C;
    for k = 1:length(t)
        E = expm(A*t(k));
        W = W + E'*Q*E*dt;
    end
    W = 0.5*(W + W');
end

function value = regularizedLogDet(W,epsVal)
    e = eig(W + epsVal*eye(size(W)));
    value = sum(log(max(real(e),epsVal)));
end

function printMetrics(name,W)
    e = eig(0.5*(W+W'));
    tr = trace(W);
    ld = regularizedLogDet(W,1e-8);
    lamMin = min(real(e));
    lamMax = max(real(e));
    if lamMin <= 1e-10
        condVal = Inf;
    else
        condVal = lamMax/lamMin;
    end
    fprintf("%8s  trace=%12.6f  logdet_eps=%12.6f  lambda_min=%12.4e  cond=%12.4e\n", ...
        name, tr, ld, lamMin, condVal);
end

      

Chapter15_Lesson3.nb

(* Chapter15_Lesson3.nb *)
(* Finite-horizon observability Gramian and qualitative sensor placement. *)

ClearAll["Global`*"];

A = { {0.0, 1.0, 0.0},
     {-2.0, -0.45, 0.8},
     {0.0, -0.7, -1.25} };

Tfinal = 6.0;
n = Length[A];

sensorMatrix[sensors_List, n_Integer] := Module[{C},
  C = ConstantArray[0.0, {Length[sensors], n}];
  Do[C[[r, sensors[[r]]]] = 1.0, {r, Length[sensors]}];
  C
];

finiteHorizonWo[A_, C_, T_] := Module[{Q, t, W},
  Q = Transpose[C].C;
  W = NIntegrate[
    Transpose[MatrixExp[A t]].Q.MatrixExp[A t],
    {t, 0, T},
    Method -> "GlobalAdaptive"
  ];
  (W + Transpose[W])/2
];

metrics[W_] := Module[{e, eps = 10^-8},
  e = Eigenvalues[(W + Transpose[W])/2];
  <|
    "Trace" -> Tr[W],
    "LogDetEps" -> Total[Log[Clip[e, {eps, Infinity}]]],
    "LambdaMin" -> Min[e],
    "LambdaMax" -> Max[e],
    "Condition" -> If[Min[e] <= eps, Infinity, Max[e]/Min[e]]
  |>
];

Print["Single-sensor scores"];
Do[
  C = sensorMatrix[{i}, n];
  W = finiteHorizonWo[A, C, Tfinal];
  Print["x", i, " -> ", metrics[W]],
  {i, 1, n}
];

Print["Two-sensor scores"];
sets = Subsets[Range[n], {2}];
scores = Table[
   C = sensorMatrix[s, n];
   W = finiteHorizonWo[A, C, Tfinal];
   {s, metrics[W]},
   {s, sets}
];

SortBy[scores, -#2["LogDetEps"] &] // TableForm

      

Problem 1 (Output Energy from a Candidate Sensor): Let \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \) and \( y_i=\mathbf{c}_i\mathbf{x} \). Show that the output energy from this single sensor is \( \mathbf{x}_0^T\mathbf{W}_i(T)\mathbf{x}_0 \).

Solution: Since \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0 \), the scalar output is \( y_i(t)=\mathbf{c}_i e^{\mathbf{A}t}\mathbf{x}_0 \). Therefore

\[ \int_0^T y_i^2(t)\,dt =\mathbf{x}_0^T \left(\int_0^T e^{\mathbf{A}^Tt} \mathbf{c}_i^T\mathbf{c}_i e^{\mathbf{A}t}\,dt\right) \mathbf{x}_0 =\mathbf{x}_0^T\mathbf{W}_i(T)\mathbf{x}_0. \]

Problem 2 (Additivity): If \( S=\{i,j\} \), prove that \( \mathbf{W}_o(S;T)=\mathbf{W}_i(T)+\mathbf{W}_j(T) \).

Solution: For \( \mathbf{C}_S= \begin{bmatrix}\mathbf{c}_i\\ \mathbf{c}_j\end{bmatrix} \),

\[ \mathbf{C}_S^T\mathbf{C}_S =\mathbf{c}_i^T\mathbf{c}_i+ \mathbf{c}_j^T\mathbf{c}_j. \]

Substituting this into the integral definition of the Gramian and using linearity of integration gives the result.

Problem 3 (Trace Can Be Misleading): Consider two positive definite Gramian matrices with eigenvalues \( \{100,1,1\} \) and \( \{40,40,20\} \). Which has larger trace? Which is better if the goal is balanced visibility?

Solution: The first has trace \( 102 \), and the second has trace \( 100 \). Trace slightly prefers the first. However, the first has condition number \( 100 \), while the second has condition number \( 2 \). The second is better for balanced visibility because its weakest direction is much stronger.

Problem 4 (Smallest Eigenvalue Interpretation): Show that among all unit-norm initial states, the minimum output energy equals \( \lambda_{\min}(\mathbf{W}_o) \).

Solution: By the Rayleigh quotient theorem, for a symmetric matrix \( \mathbf{W}_o \),

\[ \min_{\|\mathbf{x}\|_2=1} \mathbf{x}^T\mathbf{W}_o\mathbf{x} =\lambda_{\min}(\mathbf{W}_o). \]

Thus the eigenvector associated with \( \lambda_{\min} \) is the unit initial state that produces the least output energy.

Problem 5 (Qualitative Sensor Choice): Suppose sensor set A gives eigenvalues \( \{10,0.1,0.1\} \), and sensor set B gives \( \{4,3,2\} \). Which set is preferred under trace, log-determinant, and worst-direction criteria?

Solution: Set A has trace \( 10.2 \), while set B has trace \( 9 \); trace prefers A. The determinant of A is \( 10\cdot0.1\cdot0.1=0.1 \), while the determinant of B is \( 4\cdot3\cdot2=24 \); log-determinant strongly prefers B. The smallest eigenvalue of A is \( 0.1 \), while the smallest eigenvalue of B is \( 2 \); the worst-direction criterion also prefers B. This illustrates why high total energy alone is not sufficient.

Sensor placement can be studied through the observability Gramian because the Gramian maps initial states to output energy. Each candidate sensor contributes a positive semidefinite matrix, and selected sensor sets add these contributions. However, placement quality depends on how the added energy is distributed across state directions. Trace measures total visibility, log-determinant rewards volume and balance, the smallest eigenvalue identifies the weakest direction, and the condition number warns about reconstruction sensitivity. These ideas are qualitative but mathematically grounded and will be used again when studying state reconstruction and observer design.