Consider the continuous-time LTI system already used in the previous lessons: \( \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t) \), \( \mathbf{y}(t)=\mathbf{C}\mathbf{x}(t) \). With unknown initial state \( \mathbf{x}_0 \), the measured output over a finite time interval is

\[ \mathbf{y}(t)=\mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0, \qquad t\in[0,T]. \]

The state-reconstruction problem asks whether \( \mathbf{x}_0 \) can be recovered stably from \( \mathbf{y}(t) \). Observability gives uniqueness; conditioning determines numerical reliability. These are different: a system can be theoretically observable but practically difficult to reconstruct when one direction produces very little output energy.

flowchart TD
  X0["Unknown initial state x0"] --> DYN["State evolution: x(t)=exp(A t) x0"]
  DYN --> OUT["Output: y(t)=C x(t)"]
  OUT --> DATA["Measured signal on time window 0 to T"]
  DATA --> RHS["Compute b = integral exp(A' t) C' y(t) dt"]
  RHS --> SOLVE["Solve W_o xhat = b"]
  SOLVE --> CHECK["Check eigenvalues and condition number"]
  CHECK --> DECIDE["Reliable estimate or noise-amplified estimate?"]
        

The finite-horizon observability Gramian is

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

For any candidate initial state \( \mathbf{z} \), define the output prediction error functional

\[ J(\mathbf{z})=\int_0^T \left\|\mathbf{y}(t)-\mathbf{C}e^{\mathbf{A}t}\mathbf{z} \right\|_2^2\,dt. \]

Expanding the quadratic form gives

\[ J(\mathbf{z})=\int_0^T\|\mathbf{y}(t)\|_2^2dt -2\mathbf{z}^T\mathbf{b}(T) +\mathbf{z}^T\mathbf{W}_o(T)\mathbf{z}, \]

\[ \mathbf{b}(T)=\int_0^T e^{\mathbf{A}^Tt} \mathbf{C}^T\mathbf{y}(t)\,dt. \]

Setting the gradient with respect to \( \mathbf{z} \) equal to zero yields the normal equations for reconstructing the initial state:

\[ \mathbf{W}_o(T)\hat{\mathbf{x} }_0=\mathbf{b}(T). \]

If \( \mathbf{W}_o(T) \) is positive definite, the least-squares minimizer is unique:

\[ \hat{\mathbf{x} }_0=\mathbf{W}_o(T)^{-1}\mathbf{b}(T). \]

In numerical computation, however, this formula should be interpreted as a linear solve, not as an instruction to explicitly compute \( \mathbf{W}_o(T)^{-1} \).

Since \( \mathbf{W}_o(T) \) is symmetric positive definite for an observable system over the chosen interval, its eigenvalue decomposition is

\[ \mathbf{W}_o(T)=\mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^T, \qquad \boldsymbol{\Lambda}=\operatorname{diag}(\lambda_1,\dots,\lambda_n), \qquad 0 < \lambda_1\le\cdots\le\lambda_n. \]

The output energy generated by an initial state is

\[ E_y(\mathbf{x}_0)=\int_0^T\|\mathbf{y}(t)\|_2^2dt =\mathbf{x}_0^T\mathbf{W}_o(T)\mathbf{x}_0. \]

If \( \mathbf{x}_0=\mathbf{q}_i \), where \( \mathbf{q}_i \) is an eigenvector of the Gramian, then

\[ E_y(\mathbf{q}_i)=\lambda_i. \]

Therefore, small eigenvalues represent initial-state directions that produce weak output energy. The spectral condition number is

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

Large \( \kappa_2(\mathbf{W}_o) \) means output energy is distributed very unevenly across state directions. This does not necessarily destroy observability, but it makes the inverse reconstruction problem sensitive.

flowchart TD
  E1["Large eigenvalue direction"] --> O1["Strong output energy"]
  O1 --> R1["Small reconstruction variance"]
  E2["Small eigenvalue direction"] --> O2["Weak output energy"]
  O2 --> R2["Large reconstruction variance"]
  R1 --> C["Condition number = largest eigenvalue / smallest eigenvalue"]
  R2 --> C
        

Suppose the measured output contains additive noise: \( \tilde{\mathbf{y} }(t)=\mathbf{y}(t)+\boldsymbol{\eta}(t) \). Then the right-hand side becomes

\[ \tilde{\mathbf{b} }(T)=\mathbf{b}(T)+\delta\mathbf{b}(T), \qquad \delta\mathbf{b}(T)=\int_0^T e^{\mathbf{A}^Tt} \mathbf{C}^T\boldsymbol{\eta}(t)\,dt. \]

The reconstructed state perturbation satisfies

\[ \delta\mathbf{x}_0 =\tilde{\mathbf{x} }_0-\hat{\mathbf{x} }_0 =\mathbf{W}_o(T)^{-1}\delta\mathbf{b}(T). \]

Hence

\[ \|\delta\mathbf{x}_0\|_2 \le \|\mathbf{W}_o(T)^{-1}\|_2\, \|\delta\mathbf{b}(T)\|_2 =\frac{\|\delta\mathbf{b}(T)\|_2} {\lambda_{\min}(\mathbf{W}_o(T))}. \]

A more classical relative perturbation bound for the linear system \( \mathbf{W}_o\mathbf{x}_0=\mathbf{b} \) is

\[ \frac{\|\delta\mathbf{x}_0\|_2}{\|\mathbf{x}_0\|_2} \le \kappa_2(\mathbf{W}_o) \frac{\|\delta\mathbf{b}\|_2}{\|\mathbf{b}\|_2} +O\!\left(\|\delta\mathbf{b}\|_2^2\right). \]

Thus, a small relative error in the output-derived right-hand side can produce a large relative error in the reconstructed initial state when \( \kappa_2(\mathbf{W}_o) \) is large.

Proof sketch: Since \( \mathbf{W}_o\delta\mathbf{x}_0=\delta\mathbf{b} \), multiplication by \( \mathbf{W}_o^{-1} \) gives \( \delta\mathbf{x}_0=\mathbf{W}_o^{-1}\delta\mathbf{b} \). Taking norms and using submultiplicativity gives the first inequality. The relative bound follows by dividing by \( \|\mathbf{x}_0\|_2 \) and using

\[ \|\mathbf{b}\|_2 =\|\mathbf{W}_o\mathbf{x}_0\|_2 \le \|\mathbf{W}_o\|_2\|\mathbf{x}_0\|_2. \]

If the output noise is idealized as zero-mean white noise with intensity \( \sigma^2 \), the least-squares reconstruction has approximate covariance

\[ \operatorname{Cov}(\hat{\mathbf{x} }_0) \approx \sigma^2\mathbf{W}_o(T)^{-1}. \]

Along the Gramian eigenvector \( \mathbf{q}_i \), the variance is approximately

\[ \operatorname{Var}(\mathbf{q}_i^T\hat{\mathbf{x} }_0) \approx \frac{\sigma^2}{\lambda_i}. \]

This formula explains why small Gramian eigenvalues are dangerous: they correspond to directions in which a small output disturbance produces a large uncertainty in the estimated initial state.

In computation, output data are usually sampled at \( t_k=k\Delta t \). A simple quadrature-based approximation is

\[ \mathbf{W}_N =\sum_{k=0}^{N} e^{\mathbf{A}^Tt_k} \mathbf{C}^T\mathbf{C}e^{\mathbf{A}t_k}\Delta t, \qquad \mathbf{b}_N =\sum_{k=0}^{N} e^{\mathbf{A}^Tt_k} \mathbf{C}^T\mathbf{y}(t_k)\Delta t. \]

The reconstructed estimate is obtained by solving

\[ \mathbf{W}_N\hat{\mathbf{x} }_0=\mathbf{b}_N. \]

For ill-conditioned problems, regularization is often used:

\[ \hat{\mathbf{x} }_{0,\alpha} =\left(\mathbf{W}_N+\alpha\mathbf{I}\right)^{-1} \mathbf{b}_N, \qquad \alpha>0. \]

The parameter \( \alpha \) reduces noise amplification by increasing small eigenvalues, but it also introduces bias. This is the standard bias-variance trade-off in inverse problems.

Let \( \mathbf{x}=\mathbf{T}\mathbf{z} \). Then the transformed output matrix is \( \mathbf{C}_z=\mathbf{C}\mathbf{T} \), and the transformed observability Gramian is

\[ \mathbf{W}_{o,z}(T) =\mathbf{T}^T\mathbf{W}_{o,x}(T)\mathbf{T}. \]

Observability rank is invariant under nonsingular coordinate changes, but numerical condition numbers are not invariant. Therefore, conditioning must be interpreted relative to a physically meaningful scaling of states. A badly scaled state vector may falsely suggest poor observability even when the physical reconstruction problem is acceptable.

The main operations required in this lesson are matrix exponentials, numerical quadrature, eigenvalue computation, linear solves, and condition-number diagnostics.

• Python: NumPy, SciPy, python-control, and optionally Slycot.
• C++: Eigen, Armadillo, Blaze, and SLICOT/LAPACK-backed routines.
• Java: EJML, Apache Commons Math, JBLAS, and ND4J.
• MATLAB/Simulink: Control System Toolbox (ss, obsv, gram, lyap) and the Simulink State-Space block for output simulation.
• Wolfram Mathematica: MatrixExp, Eigenvalues, LinearSolve, and symbolic integration for low-dimensional examples.

Chapter15_Lesson4.py

# Chapter15_Lesson4.py
# Condition Numbers and Sensitivity in State Reconstruction
#
# Required packages:
#   pip install numpy scipy
#
# Optional related libraries for larger Modern Control workflows:
#   python-control: state-space models, observability matrix, gramians
#   slycot: faster Lyapunov/Gramian routines when available

import numpy as np
from scipy.linalg import expm, eigvalsh
from numpy.linalg import solve, cond, norm

def finite_observability_gramian(A, C, T=5.0, N=4000):
    """Compute W_o(T) = integral_0^T exp(A^T t) C^T C exp(A t) dt."""
    n = A.shape[0]
    W = np.zeros((n, n), dtype=float)
    dt = T / N

    for k in range(N):
        # Midpoint quadrature gives better accuracy than left-endpoint quadrature.
        t = (k + 0.5) * dt
        Phi = expm(A * t)
        W += Phi.T @ C.T @ C @ Phi * dt

    return 0.5 * (W + W.T)

def reconstruct_initial_state(A, C, t_grid, y_samples, ridge=0.0):
    """Solve (W + ridge I) x0_hat = b from sampled output data."""
    n = A.shape[0]
    W = np.zeros((n, n), dtype=float)
    b = np.zeros(n, dtype=float)

    for k in range(len(t_grid) - 1):
        t_mid = 0.5 * (t_grid[k] + t_grid[k + 1])
        dt = t_grid[k + 1] - t_grid[k]
        y_mid = 0.5 * (y_samples[k] + y_samples[k + 1])

        Phi = expm(A * t_mid)
        H = C @ Phi
        W += H.T @ H * dt
        b += (H.T @ np.atleast_1d(y_mid)).ravel() * dt

    W = 0.5 * (W + W.T)
    return solve(W + ridge * np.eye(n), b), W, b

def main():
    # A weakly observed second state: epsilon makes the second state visible,
    # but only faintly. The system remains observable for epsilon != 0 when
    # the two modes are distinct, yet the reconstruction can be badly conditioned.
    eps = 0.02
    A = np.array([[-1.0, 0.0],
                  [ 0.0,-2.0]])
    C = np.array([[1.0, eps]])

    T = 5.0
    N = 2500
    t = np.linspace(0.0, T, N + 1)

    x0_true = np.array([1.0, -1.5])

    # Output y(t) = C exp(A t) x0 + measurement noise.
    clean = np.array([(C @ expm(A * ti) @ x0_true).item() for ti in t])
    rng = np.random.default_rng(4)
    noise_std = 1.0e-3
    y_noisy = clean + noise_std * rng.standard_normal(clean.shape)

    x0_hat, W, b = reconstruct_initial_state(A, C, t, y_noisy, ridge=0.0)
    x0_ridge, W_ridge, _ = reconstruct_initial_state(A, C, t, y_noisy, ridge=1.0e-5)

    lam = eigvalsh(W)
    kappa = lam[-1] / lam[0]

    print("Finite-horizon observability Gramian W_o(T):")
    print(W)
    print("\nEigenvalues:", lam)
    print("2-norm condition number:", kappa)
    print("numpy cond(W):", cond(W))

    print("\nTrue x0:                  ", x0_true)
    print("Least-squares x0 estimate:", x0_hat)
    print("Ridge x0 estimate:        ", x0_ridge)

    print("\nRelative LS error:", norm(x0_hat - x0_true) / norm(x0_true))
    print("Relative ridge error:", norm(x0_ridge - x0_true) / norm(x0_true))

    # First-order sensitivity bound:
    # ||delta x|| <= ||W^{-1}|| ||delta b||
    clean_hat, W_clean, b_clean = reconstruct_initial_state(A, C, t, clean, ridge=0.0)
    delta_b = b - b_clean
    delta_x = x0_hat - clean_hat

    bound = norm(solve(W, delta_b))
    print("\nActual perturbation norm ||delta x||:", norm(delta_x))
    print("Bound quantity ||W^{-1} delta b||:   ", bound)

if __name__ == "__main__":
    main()

Chapter15_Lesson4.cpp

// Chapter15_Lesson4.cpp
// Condition Numbers and Sensitivity in State Reconstruction
//
// Build:
//   g++ -std=c++17 -O2 Chapter15_Lesson4.cpp -o Chapter15_Lesson4
//
// This example intentionally uses a small from-scratch 2x2 implementation.
// For larger Modern Control workflows, use Eigen, Armadillo, Blaze, SLICOT,
// or a C++ control library with Lyapunov and state-space routines.

#include <cmath>
#include <iomanip>
#include <iostream>
#include <array>

struct Vec2 {
    double x1;
    double x2;
};

struct Mat2 {
    double a11, a12, a21, a22;
};

double det2(const Mat2& M) {
    return M.a11 * M.a22 - M.a12 * M.a21;
}

Vec2 solve2(const Mat2& M, const Vec2& b) {
    double d = det2(M);
    return {
        ( M.a22 * b.x1 - M.a12 * b.x2) / d,
        (-M.a21 * b.x1 + M.a11 * b.x2) / d
    };
}

std::array<double, 2> eigenvaluesSymmetric2(const Mat2& M) {
    double tr = M.a11 + M.a22;
    double diff = M.a11 - M.a22;
    double rad = std::sqrt(diff * diff + 4.0 * M.a12 * M.a12);
    return {0.5 * (tr - rad), 0.5 * (tr + rad)};
}

int main() {
    const double eps = 0.02;
    const double T = 5.0;
    const int N = 200000;
    const double dt = T / static_cast<double>(N);

    // A = diag(-1, -2), C = [1, eps]
    // True initial condition
    const Vec2 x0_true{1.0, -1.5};

    Mat2 W{0.0, 0.0, 0.0, 0.0};
    Vec2 b{0.0, 0.0};

    for (int k = 0; k < N; ++k) {
        double t = (k + 0.5) * dt;

        double phi1 = std::exp(-t);
        double phi2 = std::exp(-2.0 * t);

        // H(t) = C exp(A t) = [phi1, eps phi2]
        double h1 = phi1;
        double h2 = eps * phi2;

        // Deterministic small measurement disturbance
        double y_clean = h1 * x0_true.x1 + h2 * x0_true.x2;
        double noise = 1.0e-3 * std::sin(37.0 * t);
        double y = y_clean + noise;

        W.a11 += h1 * h1 * dt;
        W.a12 += h1 * h2 * dt;
        W.a21 += h2 * h1 * dt;
        W.a22 += h2 * h2 * dt;

        b.x1 += h1 * y * dt;
        b.x2 += h2 * y * dt;
    }

    Vec2 xhat = solve2(W, b);
    auto evals = eigenvaluesSymmetric2(W);
    double kappa = evals[1] / evals[0];

    std::cout << std::setprecision(12);
    std::cout << "W = [[" << W.a11 << ", " << W.a12 << "], ["
              << W.a21 << ", " << W.a22 << "]]\n";
    std::cout << "lambda_min = " << evals[0] << "\n";
    std::cout << "lambda_max = " << evals[1] << "\n";
    std::cout << "condition number = " << kappa << "\n\n";

    std::cout << "true x0 = [" << x0_true.x1 << ", " << x0_true.x2 << "]\n";
    std::cout << "estimated x0 = [" << xhat.x1 << ", " << xhat.x2 << "]\n";

    // Ridge-regularized reconstruction: (W + alpha I) x = b
    double alpha = 1.0e-5;
    Mat2 Wr{W.a11 + alpha, W.a12, W.a21, W.a22 + alpha};
    Vec2 xr = solve2(Wr, b);
    std::cout << "ridge estimated x0 = [" << xr.x1 << ", " << xr.x2 << "]\n";

    return 0;
}

Chapter15_Lesson4.java

// Chapter15_Lesson4.java
// Condition Numbers and Sensitivity in State Reconstruction
//
// Build:
//   javac Chapter15_Lesson4.java
// Run:
//   java Chapter15_Lesson4
//
// This example uses a minimal from-scratch 2x2 implementation.
// For larger Modern Control workflows in Java, consider EJML, JBLAS,
// Apache Commons Math, or ND4J for matrix algebra.

public class Chapter15_Lesson4 {
    static class Vec2 {
        double x1, x2;
        Vec2(double x1, double x2) {
            this.x1 = x1;
            this.x2 = x2;
        }
    }

    static class Mat2 {
        double a11, a12, a21, a22;
        Mat2(double a11, double a12, double a21, double a22) {
            this.a11 = a11;
            this.a12 = a12;
            this.a21 = a21;
            this.a22 = a22;
        }
    }

    static double det2(Mat2 M) {
        return M.a11 * M.a22 - M.a12 * M.a21;
    }

    static Vec2 solve2(Mat2 M, Vec2 b) {
        double d = det2(M);
        return new Vec2(
            ( M.a22 * b.x1 - M.a12 * b.x2) / d,
            (-M.a21 * b.x1 + M.a11 * b.x2) / d
        );
    }

    static double[] eigenvaluesSymmetric2(Mat2 M) {
        double tr = M.a11 + M.a22;
        double diff = M.a11 - M.a22;
        double rad = Math.sqrt(diff * diff + 4.0 * M.a12 * M.a12);
        return new double[] {0.5 * (tr - rad), 0.5 * (tr + rad)};
    }

    public static void main(String[] args) {
        double eps = 0.02;
        double T = 5.0;
        int N = 200000;
        double dt = T / N;

        Vec2 x0True = new Vec2(1.0, -1.5);

        Mat2 W = new Mat2(0.0, 0.0, 0.0, 0.0);
        Vec2 b = new Vec2(0.0, 0.0);

        for (int k = 0; k < N; k++) {
            double t = (k + 0.5) * dt;

            double phi1 = Math.exp(-t);
            double phi2 = Math.exp(-2.0 * t);

            // H(t) = C exp(A t) = [phi1, eps phi2]
            double h1 = phi1;
            double h2 = eps * phi2;

            double yClean = h1 * x0True.x1 + h2 * x0True.x2;
            double noise = 1.0e-3 * Math.sin(37.0 * t);
            double y = yClean + noise;

            W.a11 += h1 * h1 * dt;
            W.a12 += h1 * h2 * dt;
            W.a21 += h2 * h1 * dt;
            W.a22 += h2 * h2 * dt;

            b.x1 += h1 * y * dt;
            b.x2 += h2 * y * dt;
        }

        Vec2 xhat = solve2(W, b);
        double[] evals = eigenvaluesSymmetric2(W);
        double kappa = evals[1] / evals[0];

        System.out.printf("W = [[%.12f, %.12f], [%.12f, %.12f]]%n",
            W.a11, W.a12, W.a21, W.a22);
        System.out.printf("lambda_min = %.12e%n", evals[0]);
        System.out.printf("lambda_max = %.12e%n", evals[1]);
        System.out.printf("condition number = %.12e%n%n", kappa);

        System.out.printf("true x0 = [%.12f, %.12f]%n", x0True.x1, x0True.x2);
        System.out.printf("estimated x0 = [%.12f, %.12f]%n", xhat.x1, xhat.x2);

        // Ridge-regularized reconstruction: (W + alpha I) x = b
        double alpha = 1.0e-5;
        Mat2 Wr = new Mat2(W.a11 + alpha, W.a12, W.a21, W.a22 + alpha);
        Vec2 xr = solve2(Wr, b);
        System.out.printf("ridge estimated x0 = [%.12f, %.12f]%n", xr.x1, xr.x2);
    }
}

Chapter15_Lesson4.m

% Chapter15_Lesson4.m
% Condition Numbers and Sensitivity in State Reconstruction
%
% Required core functions: expm, eig, cond, mldivide
% Related MATLAB/Simulink Modern Control tools:
%   Control System Toolbox: ss, obsv, gram, lyap
%   Simulink: State-Space block for output simulation

clear; clc;

epsSensor = 0.02;
A = [-1 0; 0 -2];
C = [1 epsSensor];

T = 5;
N = 4000;
t = linspace(0, T, N + 1);
dt = t(2) - t(1);

x0_true = [1; -1.5];

W = zeros(2);
b = zeros(2,1);

for k = 1:N
    tm = 0.5 * (t(k) + t(k + 1));
    Phi = expm(A * tm);
    H = C * Phi;

    y_clean = H * x0_true;
    y_noisy = y_clean + 1.0e-3 * sin(37 * tm);

    W = W + (H' * H) * dt;
    b = b + H' * y_noisy * dt;
end

W = 0.5 * (W + W');
x0_hat = W \ b;

alpha = 1.0e-5;
x0_ridge = (W + alpha * eye(2)) \ b;

lambda = eig(W);
kappa = cond(W);

disp('Finite-horizon observability Gramian W_o(T):');
disp(W);
disp('Eigenvalues:');
disp(lambda.');
fprintf('2-norm condition number: %.12e\n\n', kappa);

disp('True x0:');
disp(x0_true.');
disp('Least-squares estimate:');
disp(x0_hat.');
disp('Ridge estimate:');
disp(x0_ridge.');

% Optional Control System Toolbox checks
if exist('obsv', 'file') == 2
    O = obsv(A, C);
    fprintf('\nRank of Kalman observability matrix: %d\n', rank(O));
end

if exist('ss', 'file') == 2
    sys = ss(A, [], C, []);
    disp('State-space model created with ss(A,[],C,[]).');
end

% Optional Simulink note:
% A corresponding Simulink model can be built with a State-Space block using
% A, B=[], C, D=[] for autonomous response, with initial condition x0_true.
% The reconstruction algorithm above is then applied to the logged output y(t).

Chapter15_Lesson4.nb

(* Chapter15_Lesson4.nb *)
(* Condition Numbers and Sensitivity in State Reconstruction *)
(* This notebook-style Wolfram Language code computes a finite-horizon
   observability Gramian and reconstructs the initial state. *)

ClearAll["Global`*"];

epsSensor = 0.02;
A = { {-1., 0.}, {0., -2.} };
Cmat = { {1., epsSensor} };
T = 5.;
nSteps = 4000;
dt = T/nSteps;

x0True = {1., -1.5};

phi[t_] := MatrixExp[A t];
h[t_] := Cmat . phi[t];

yClean[t_] := First[h[t] . x0True];
yNoisy[t_] := yClean[t] + 10^-3 Sin[37 t];

W = Sum[
   Transpose[h[(k + 1/2) dt]] . h[(k + 1/2) dt] dt,
   {k, 0, nSteps - 1}
];

b = Sum[
   Flatten[Transpose[h[(k + 1/2) dt]] yNoisy[(k + 1/2) dt] dt],
   {k, 0, nSteps - 1}
];

W = (W + Transpose[W])/2;

x0Hat = LinearSolve[W, b];
alpha = 10^-5;
x0Ridge = LinearSolve[W + alpha IdentityMatrix[2], b];

eigs = Eigenvalues[W];
kappa = Max[eigs]/Min[eigs];

Print["Finite-horizon observability Gramian W_o(T):"];
Print[MatrixForm[W]];
Print["Eigenvalues = ", eigs];
Print["Condition number = ", kappa];
Print["True x0 = ", x0True];
Print["Least-squares estimate = ", x0Hat];
Print["Ridge estimate = ", x0Ridge];

(* Optional symbolic diagonal formula for A=diag(a1,a2), C={c1,c2}:
   W_ij(T) = c_i c_j (Exp[(a_i+a_j) T] - 1)/(a_i+a_j)
   when a_i + a_j != 0. *)

Problem 1 (Scalar Reconstruction): Let \( \dot{x}=ax \) and \( y=cx \). Derive the finite-horizon observability Gramian and determine when the scalar state is reconstructible.

Solution:

\[ W_o(T)=\int_0^T c^2e^{2at}dt =\begin{cases} c^2\dfrac{e^{2aT}-1}{2a}, & a\neq 0,\\ c^2T, & a=0. \end{cases} \]

The scalar initial state is reconstructible if and only if \( W_o(T)>0 \), which requires \( c\neq 0 \) and \( T>0 \).

Problem 2 (Weak Sensor Direction): Consider \( \mathbf{A}=\operatorname{diag}(-1,-2) \) and \( \mathbf{C}=\begin{bmatrix}1 & \varepsilon\end{bmatrix} \). Compute the Gramian entries over \( [0,T] \).

Solution:

\[ \mathbf{W}_o(T)= \begin{bmatrix} \dfrac{1-e^{-2T} }{2} & \varepsilon\dfrac{1-e^{-3T} }{3} \\ \varepsilon\dfrac{1-e^{-3T} }{3} & \varepsilon^2\dfrac{1-e^{-4T} }{4} \end{bmatrix}. \]

The second diagonal term is proportional to \( \varepsilon^2 \). Therefore, as \( \varepsilon \) becomes small, the second state remains theoretically visible for nonzero \( \varepsilon \), but its reconstruction becomes increasingly ill-conditioned.

Problem 3 (Sensitivity Bound): Suppose \( \lambda_{\min}(\mathbf{W}_o)=10^{-5} \) and \( \|\delta\mathbf{b}\|_2=10^{-6} \). Bound \( \|\delta\mathbf{x}_0\|_2 \).

Solution:

\[ \|\delta\mathbf{x}_0\|_2 \le \frac{\|\delta\mathbf{b}\|_2} {\lambda_{\min}(\mathbf{W}_o)} =\frac{10^{-6} }{10^{-5} }=10^{-1}. \]

A small right-hand-side perturbation can therefore produce an initial state error of order \( 0.1 \).

Problem 4 (Coordinate Scaling): Let \( \mathbf{x}=\mathbf{T}\mathbf{z} \). Show that \( \mathbf{W}_{o,z}=\mathbf{T}^T \mathbf{W}_{o,x}\mathbf{T} \).

Solution:

Since \( \mathbf{A}_z=\mathbf{T}^{-1}\mathbf{A}_x\mathbf{T} \) and \( \mathbf{C}_z=\mathbf{C}_x\mathbf{T} \), the transformed output satisfies \( \mathbf{y}=\mathbf{C}_x\mathbf{x} =\mathbf{C}_x\mathbf{T}\mathbf{z} \). Substitution into the Gramian integral gives

\[ \mathbf{W}_{o,z} =\int_0^T e^{\mathbf{A}_z^Tt} \mathbf{C}_z^T\mathbf{C}_ze^{\mathbf{A}_zt}dt =\mathbf{T}^T\mathbf{W}_{o,x}\mathbf{T}. \]

Problem 5 (Regularized Reconstruction): If \( \mathbf{W}_o=\mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^T \), express the Tikhonov estimate \( (\mathbf{W}_o+lpha\mathbf{I})^{-1}\mathbf{b} \) in eigen-coordinates.

Solution:

\[ \hat{\mathbf{x} }_{0,\alpha} =\mathbf{Q}(\boldsymbol{\Lambda}+lpha\mathbf{I})^{-1} \mathbf{Q}^T\mathbf{b} =\sum_{i=1}^n \frac{\mathbf{q}_i^T\mathbf{b} }{\lambda_i+lpha} \mathbf{q}_i. \]

Regularization limits amplification in directions with small \( \lambda_i \) by replacing \( 1/\lambda_i \) with \( 1/(\lambda_i+lpha) \).

Observability answers whether an initial state is uniquely determined by output measurements, but the condition number of the observability Gramian answers whether the reconstruction is numerically stable. Small Gramian eigenvalues identify weakly observable directions, cause large error covariance, and motivate careful scaling, stable linear solves, and regularization when reconstructing states from noisy output data.