For an LTI system, the matrices \( A \) and \( C \) are constant, so the unmeasured internal modes are fixed directions in state space. In a linear time-varying (LTV) system,

\[ \dot{\mathbf{x} }(t)=\mathbf{A}(t)\mathbf{x}(t)+\mathbf{B}(t)\mathbf{u}(t), \qquad \mathbf{y}(t)=\mathbf{C}(t)\mathbf{x}(t)+\mathbf{D}(t)\mathbf{u}(t) \]

the state dynamics and sensor map may rotate, scale, or lose/gain sensitivity over time. Therefore, a state direction that is invisible at one instant may become visible later. Observability is consequently an interval property: we ask whether \( \mathbf{x}(t_0) \) is uniquely determined from the known input \( \mathbf{u}(t) \) and the measured output \( \mathbf{y}(t) \) for \( t\in[t_0,t_f] \).

The correct replacement for powers \( \mathbf{C}\mathbf{A}^k \) is the time-dependent map \( \mathbf{C}(t)\boldsymbol{\Phi}(t,t_0) \), where \( \boldsymbol{\Phi}(t,t_0) \) is the state-transition matrix of \( \dot{\mathbf{x} }=\mathbf{A}(t)\mathbf{x} \).

flowchart TD
  A["Unknown initial state x(t0)"] --> B["State transition Phi(t,t0)"]
  B --> C["Sensor map C(t)"]
  C --> D["Output history y(t) on interval"]
  D --> E["Remove known input contribution"]
  E --> F["Build information integral W_o"]
  F --> G{"rank(W_o) = n?"}
  G -->|"yes"| H["Initial state is reconstructable"]
  G -->|"no"| I["Some state direction remains invisible"]
        

The LTV state-transition matrix is defined by

\[ \frac{\partial}{\partial t}\boldsymbol{\Phi}(t,t_0) =\mathbf{A}(t)\boldsymbol{\Phi}(t,t_0),\qquad \boldsymbol{\Phi}(t_0,t_0)=\mathbf{I}. \]

For a known input, the solution is

\[ \mathbf{x}(t)=\boldsymbol{\Phi}(t,t_0)\mathbf{x}(t_0) +\int_{t_0}^{t}\boldsymbol{\Phi}(t,s)\mathbf{B}(s)\mathbf{u}(s)\,ds. \]

Substitution into the output equation gives

\[ \mathbf{y}(t)= \mathbf{C}(t)\boldsymbol{\Phi}(t,t_0)\mathbf{x}(t_0) +\mathbf{C}(t)\int_{t_0}^{t}\boldsymbol{\Phi}(t,s)\mathbf{B}(s) \mathbf{u}(s)\,ds+\mathbf{D}(t)\mathbf{u}(t). \]

Since \( \mathbf{u}(t) \), \( \mathbf{B}(t) \), and \( \mathbf{D}(t) \) are assumed known, the forced contribution can be subtracted. Define the homogeneous output component:

\[ \tilde{\mathbf{y} }(t)= \mathbf{y}(t)-\mathbf{C}(t)\int_{t_0}^{t}\boldsymbol{\Phi}(t,s) \mathbf{B}(s)\mathbf{u}(s)\,ds-\mathbf{D}(t)\mathbf{u}(t) =\mathbf{C}(t)\boldsymbol{\Phi}(t,t_0)\mathbf{x}(t_0). \]

Therefore, time-varying observability is fundamentally about the injectivity of the linear operator \( \mathbf{x}(t_0)\mapsto \tilde{\mathbf{y} }(\cdot) \).

Two initial states \( \mathbf{x}_a(t_0) \) and \( \mathbf{x}_b(t_0) \) are indistinguishable on \( [t_0,t_f] \) if they generate the same output for the same known input. Their difference \( \boldsymbol{\eta}=\mathbf{x}_a(t_0)-\mathbf{x}_b(t_0) \) satisfies

\[ \mathbf{C}(t)\boldsymbol{\Phi}(t,t_0)\boldsymbol{\eta} =\mathbf{0},\qquad t\in[t_0,t_f]. \]

Thus the unobservable subspace over the interval is

\[ \mathcal{N}_o[t_0,t_f]= \left\{\boldsymbol{\eta}\in\mathbb{R}^{n}: \mathbf{C}(t)\boldsymbol{\Phi}(t,t_0)\boldsymbol{\eta} =\mathbf{0}\ \text{for all}\ t\in[t_0,t_f]\right\}. \]

The system is observable on \( [t_0,t_f] \) precisely when no nonzero initial-state direction belongs to this set:

\[ \boxed{\text{observable on }[t_0,t_f]\iff \mathcal{N}_o[t_0,t_f]=\{\mathbf{0}\}.} \]

This definition is deliberately finite-interval. A sensor may fail to distinguish states on a short interval but become informative on a longer interval; conversely, a moving sensor may be informative only during a restricted time window.

The finite-interval LTV observability Gramian is

\[ \mathbf{W}_o(t_0,t_f)= \int_{t_0}^{t_f} \boldsymbol{\Phi}^{T}(t,t_0)\mathbf{C}^{T}(t)\mathbf{C}(t) \boldsymbol{\Phi}(t,t_0)\,dt. \]

It is symmetric and positive semidefinite because, for every vector \( \boldsymbol{\eta} \),

\[ \boldsymbol{\eta}^{T}\mathbf{W}_o(t_0,t_f)\boldsymbol{\eta} =\int_{t_0}^{t_f} \left\|\mathbf{C}(t)\boldsymbol{\Phi}(t,t_0)\boldsymbol{\eta} \right\|_2^2\,dt\ge 0. \]

The rank condition is the LTV analog of the Kalman observability test:

\[ \boxed{\text{observable on }[t_0,t_f]\iff \operatorname{rank}\mathbf{W}_o(t_0,t_f)=n \iff \mathbf{W}_o(t_0,t_f)\text{ is positive definite}.} \]

In computational work, the smallest eigenvalue of \( \mathbf{W}_o \) indicates the least visible initial-state direction, while a large condition number indicates sensitivity in state reconstruction.

We prove the equivalence under the usual assumptions that \( \mathbf{A}(t) \) and \( \mathbf{C}(t) \) are piecewise continuous on the finite interval.

Step 1: non-observability implies singular Gramian. If the system is not observable, there exists a nonzero vector \( \boldsymbol{\eta}\ne\mathbf{0} \) such that \( \mathbf{C}(t)\boldsymbol{\Phi}(t,t_0)\boldsymbol{\eta} =\mathbf{0} \) for all \( t\in[t_0,t_f] \). Therefore,

\[ \boldsymbol{\eta}^{T}\mathbf{W}_o\boldsymbol{\eta} =\int_{t_0}^{t_f} \left\|\mathbf{C}(t)\boldsymbol{\Phi}(t,t_0)\boldsymbol{\eta} \right\|_2^2\,dt=0. \]

A symmetric positive semidefinite matrix that maps a nonzero vector to zero quadratic energy cannot be positive definite; hence \( \mathbf{W}_o \) is singular.

Step 2: singular Gramian implies non-observability. If \( \mathbf{W}_o \) is singular, there exists \( \boldsymbol{\eta}\ne\mathbf{0} \) such that \( \boldsymbol{\eta}^{T}\mathbf{W}_o\boldsymbol{\eta}=0 \). Since the integrand is nonnegative,

\[ \int_{t_0}^{t_f} \left\|\mathbf{C}(t)\boldsymbol{\Phi}(t,t_0)\boldsymbol{\eta} \right\|_2^2\,dt=0. \]

Hence the homogeneous output generated by \( \boldsymbol{\eta} \) is zero almost everywhere. By continuity on subintervals where the coefficients are continuous, it is zero throughout the interval. Thus \( \boldsymbol{\eta} \) is an unobservable initial-state direction, and the system is not observable. The two implications prove the theorem.

When \( \mathbf{W}_o(t_0,t_f) \) is nonsingular, the initial state can be reconstructed from the homogeneous output \( \tilde{\mathbf{y} }(t) \). Multiply \( \tilde{\mathbf{y} }(t)=\mathbf{C}(t)\boldsymbol{\Phi}(t,t_0) \mathbf{x}(t_0) \) by \( \boldsymbol{\Phi}^{T}(t,t_0)\mathbf{C}^{T}(t) \) and integrate:

\[ \mathbf{z}= \int_{t_0}^{t_f} \boldsymbol{\Phi}^{T}(t,t_0)\mathbf{C}^{T}(t)\tilde{\mathbf{y} }(t)\,dt = \mathbf{W}_o(t_0,t_f)\mathbf{x}(t_0). \]

Therefore,

\[ \boxed{\mathbf{x}(t_0)= \mathbf{W}_o^{-1}(t_0,t_f) \int_{t_0}^{t_f} \boldsymbol{\Phi}^{T}(t,t_0)\mathbf{C}^{T}(t)\tilde{\mathbf{y} }(t)\,dt.} \]

This formula is mathematically exact for noiseless data, but it is not necessarily numerically safe. If \( \lambda_{\min}(\mathbf{W}_o) \) is very small, the state direction associated with that eigenvector contributes little output energy. Measurement noise in that direction is amplified by \( \mathbf{W}_o^{-1} \). This is why observability should be interpreted both as a rank property and as a conditioning property.

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

The Gramian test is an interval test. For smooth or analytic LTV systems, a local test can be formed from output derivatives. For the homogeneous system \( \dot{\mathbf{x} }=\mathbf{A}(t)\mathbf{x} \), \( \mathbf{y}=\mathbf{C}(t)\mathbf{x} \), define

\[ \mathbf{L}_0(t)=\mathbf{C}(t),\qquad \mathbf{L}_{k+1}(t)=\dot{\mathbf{L} }_k(t)+\mathbf{L}_k(t)\mathbf{A}(t). \]

Then

\[ \frac{d^k\mathbf{y} }{dt^k}(t)=\mathbf{L}_k(t)\mathbf{x}(t). \]

If, at a time \( t^\star \), the stacked matrix

\[ \mathcal{O}_{\mathrm{loc} }(t^\star)= \begin{bmatrix} \mathbf{L}_0(t^\star)\\ \mathbf{L}_1(t^\star)\\ \vdots\\ \mathbf{L}_{n-1}(t^\star) \end{bmatrix} \]

has rank \( n \), then the state is locally distinguishable near \( t^\star \). This resembles the LTI observability matrix, but the derivative term \( \dot{\mathbf{L} }_k(t) \) is essential. Ignoring it is a common error when extending the LTI test to time-varying systems.

The LTI duality \( (\mathbf{A},\mathbf{C}) \) observable \( \iff \) \( (\mathbf{A}^{T},\mathbf{C}^{T}) \) controllable has an LTV counterpart. The natural dual dynamics are

\[ \dot{\mathbf{p} }(t)=-\mathbf{A}^{T}(t)\mathbf{p}(t)+ \mathbf{C}^{T}(t)\mathbf{v}(t). \]

The transition matrix of \( \dot{\mathbf{p} }=-\mathbf{A}^{T}(t)\mathbf{p} \) is related to the original transition matrix by

\[ \boldsymbol{\Phi}_d(t,s)=\boldsymbol{\Phi}^{T}(s,t). \]

The controllability Gramian of the dual system over \( [t_0,t_f] \) satisfies

\[ \mathbf{W}_{c,d}(t_0,t_f)= \boldsymbol{\Phi}^{T}(t_0,t_f)\mathbf{W}_o(t_0,t_f) \boldsymbol{\Phi}(t_0,t_f). \]

Since state-transition matrices are nonsingular, the two Gramians have the same rank. Thus finite-interval observability of \( (\mathbf{A}(t),\mathbf{C}(t)) \) is equivalent to finite-interval controllability of the dual pair \( (-\mathbf{A}^{T}(t),\mathbf{C}^{T}(t)) \) in the appropriate time direction.

flowchart TD
  A["Original LTV pair: A(t), C(t)"] --> B["Output map: C(t) Phi(t,t0)"]
  B --> C["Observability Gramian W_o"]
  C --> D["Rank test for initial-state recovery"]
  A --> E["Dual pair: -A(t)^T, C(t)^T"]
  E --> F["Controllability Gramian of dual"]
  F --> G["Same rank after endpoint transform"]
  G --> D
        

For a discrete-time LTV system

\[ \mathbf{x}_{k+1}=\mathbf{A}_k\mathbf{x}_k+\mathbf{B}_k\mathbf{u}_k, \qquad \mathbf{y}_k=\mathbf{C}_k\mathbf{x}_k+\mathbf{D}_k\mathbf{u}_k, \]

define the transition product

\[ \boldsymbol{\Phi}_{k,k_0}= \mathbf{A}_{k-1}\mathbf{A}_{k-2}\cdots\mathbf{A}_{k_0}, \qquad \boldsymbol{\Phi}_{k_0,k_0}=\mathbf{I}. \]

After subtracting the known forced response, the homogeneous output is \( \tilde{\mathbf{y} }_k=\mathbf{C}_k\boldsymbol{\Phi}_{k,k_0} \mathbf{x}_{k_0} \). The discrete observability Gramian is

\[ \mathbf{W}_o[k_0,k_f]= \sum_{k=k_0}^{k_f} \boldsymbol{\Phi}_{k,k_0}^{T}\mathbf{C}_k^{T}\mathbf{C}_k \boldsymbol{\Phi}_{k,k_0}. \]

The finite sequence is observable if and only if \( \operatorname{rank}\mathbf{W}_o[k_0,k_f]=n \). This form is often easier to compute in sampled-data observers and moving-sensor scheduling problems.

In continuous time, compute the Gramian by integrating the coupled differential equations

\[ \dot{\boldsymbol{\Phi} }(t)=\mathbf{A}(t)\boldsymbol{\Phi}(t), \qquad \dot{\mathbf{W} }(t)= \boldsymbol{\Phi}^{T}(t,t_0)\mathbf{C}^{T}(t)\mathbf{C}(t) \boldsymbol{\Phi}(t,t_0), \]

\[ \boldsymbol{\Phi}(t_0,t_0)=\mathbf{I},\qquad \mathbf{W}(t_0)=\mathbf{0},\qquad \mathbf{W}_o(t_0,t_f)=\mathbf{W}(t_f). \]

A practical procedure is: (i) integrate \( \boldsymbol{\Phi} \); (ii) accumulate the Gramian by quadrature or by the matrix differential equation above; (iii) inspect the rank, eigenvalues, and condition number; (iv) if reconstruction is required, solve \( \mathbf{W}_o\mathbf{x}(t_0)=\mathbf{z} \) rather than explicitly forming \( \mathbf{W}_o^{-1} \).

The following Python script computes the LTV observability Gramian, checks its eigenvalues, and reconstructs the initial state from noiseless homogeneous output data.

Chapter14_Lesson5.py

# Chapter14_Lesson5.py
"""
Time-varying observability demo for Chapter 14, Lesson 5.

The example computes the continuous-time observability Gramian

    W_o = integral Phi(t,t0)^T C(t)^T C(t) Phi(t,t0) dt

for an LTV system xdot = A(t)x, y = C(t)x, then reconstructs x(t0)
from noiseless output samples using x0 = W_o^{-1} z.
"""

import numpy as np
from scipy.integrate import solve_ivp


def A(t: float) -> np.ndarray:
    """Time-varying state matrix."""
    return np.array(
        [
            [0.0, 1.0],
            [-(2.0 + 0.5 * np.sin(1.3 * t)), -(0.15 + 0.05 * np.cos(t))],
        ],
        dtype=float,
    )


def C(t: float) -> np.ndarray:
    """Time-varying output matrix."""
    return np.array([[1.0, 0.2 * np.sin(0.7 * t)]], dtype=float)


def transition_rhs(t: float, phi_flat: np.ndarray) -> np.ndarray:
    """ODE for Phi(t,t0), flattened column-major."""
    n = 2
    phi = phi_flat.reshape((n, n))
    dphi = A(t) @ phi
    return dphi.reshape(-1)


def compute_phi_grid(t0: float, tf: float, samples: int = 1001):
    n = 2
    t_grid = np.linspace(t0, tf, samples)
    phi0 = np.eye(n).reshape(-1)
    sol = solve_ivp(
        transition_rhs,
        (t0, tf),
        phi0,
        t_eval=t_grid,
        rtol=1e-10,
        atol=1e-12,
    )
    if not sol.success:
        raise RuntimeError(sol.message)
    phis = sol.y.T.reshape((-1, n, n))
    return t_grid, phis


def trapezoid_integral(t_grid: np.ndarray, values: np.ndarray) -> np.ndarray:
    """Integrate a vector or matrix-valued array along t using trapezoids."""
    total = np.zeros_like(values[0])
    for k in range(len(t_grid) - 1):
        dt = t_grid[k + 1] - t_grid[k]
        total += 0.5 * dt * (values[k] + values[k + 1])
    return total


def observability_gramian(t_grid: np.ndarray, phis: np.ndarray) -> np.ndarray:
    terms = []
    for t, phi in zip(t_grid, phis):
        ct = C(float(t))
        terms.append(phi.T @ ct.T @ ct @ phi)
    return trapezoid_integral(t_grid, np.array(terms))


def reconstruct_initial_state(t_grid: np.ndarray, phis: np.ndarray, x0_true: np.ndarray):
    y_samples = []
    z_terms = []
    for t, phi in zip(t_grid, phis):
        ct = C(float(t))
        y = ct @ phi @ x0_true
        y_samples.append(float(y[0]))
        z_terms.append(phi.T @ ct.T @ y)
    z = trapezoid_integral(t_grid, np.array(z_terms)).reshape(-1)
    w = observability_gramian(t_grid, phis)
    x0_hat = np.linalg.solve(w, z)
    return np.array(y_samples), z, w, x0_hat


if __name__ == "__main__":
    t0, tf = 0.0, 6.0
    t_grid, phis = compute_phi_grid(t0, tf)

    w = observability_gramian(t_grid, phis)
    eigvals = np.linalg.eigvalsh(w)

    x0_true = np.array([1.0, -0.7])
    y_samples, z, w, x0_hat = reconstruct_initial_state(t_grid, phis, x0_true)

    print("Observability Gramian W_o:")
    print(w)
    print("\nEigenvalues of W_o:", eigvals)
    print("Condition number:", np.linalg.cond(w))
    print("\nTrue x0:", x0_true)
    print("Reconstructed x0:", x0_hat)
    print("Reconstruction error norm:", np.linalg.norm(x0_hat - x0_true))

      

This C++ implementation uses a small 2-by-2 matrix utility and RK4 to approximate \( \boldsymbol{\Phi} \) and \( \mathbf{W}_o \) without external dependencies.

Chapter14_Lesson5.cpp

// Chapter14_Lesson5.cpp
// Continuous-time LTV observability Gramian using a small RK4 integrator.
// No external linear algebra library is required for this 2-state example.

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

using Mat2 = std::array<std::array<double, 2>, 2>;

Mat2 zeros() {
    return { { {0.0, 0.0}, {0.0, 0.0} } };
}

Mat2 identity() {
    return { { {1.0, 0.0}, {0.0, 1.0} } };
}

Mat2 add(const Mat2& A, const Mat2& B) {
    Mat2 C = zeros();
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            C[i][j] = A[i][j] + B[i][j];
    return C;
}

Mat2 scale(const Mat2& A, double s) {
    Mat2 C = zeros();
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            C[i][j] = s * A[i][j];
    return C;
}

Mat2 multiply(const Mat2& A, const Mat2& B) {
    Mat2 C = zeros();
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            for (int k = 0; k < 2; ++k)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}

Mat2 transpose(const Mat2& A) {
    return { { {A[0][0], A[1][0]}, {A[0][1], A[1][1]} } };
}

Mat2 A_of_t(double t) {
    return { { {0.0, 1.0},
             {-(2.0 + 0.5 * std::sin(1.3 * t)), -(0.15 + 0.05 * std::cos(t))} } };
}

Mat2 CtC_of_t(double t) {
    double c1 = 1.0;
    double c2 = 0.2 * std::sin(0.7 * t);
    return { { {c1 * c1, c1 * c2}, {c2 * c1, c2 * c2} } };
}

Mat2 phi_rhs(double t, const Mat2& Phi) {
    return multiply(A_of_t(t), Phi);
}

Mat2 rk4_phi_step(double t, const Mat2& Phi, double h) {
    Mat2 k1 = phi_rhs(t, Phi);
    Mat2 k2 = phi_rhs(t + 0.5 * h, add(Phi, scale(k1, 0.5 * h)));
    Mat2 k3 = phi_rhs(t + 0.5 * h, add(Phi, scale(k2, 0.5 * h)));
    Mat2 k4 = phi_rhs(t + h, add(Phi, scale(k3, h)));

    Mat2 incr = add(add(k1, scale(k2, 2.0)), add(scale(k3, 2.0), k4));
    return add(Phi, scale(incr, h / 6.0));
}

Mat2 gramian_integrand(double t, const Mat2& Phi) {
    return multiply(multiply(transpose(Phi), CtC_of_t(t)), Phi);
}

void print_matrix(const Mat2& M) {
    std::cout << std::fixed << std::setprecision(8);
    for (int i = 0; i < 2; ++i) {
        std::cout << "[ ";
        for (int j = 0; j < 2; ++j) std::cout << std::setw(12) << M[i][j] << " ";
        std::cout << "]\n";
    }
}

int main() {
    const double t0 = 0.0;
    const double tf = 6.0;
    const int steps = 6000;
    const double h = (tf - t0) / steps;

    Mat2 Phi = identity();
    Mat2 W = zeros();

    double t = t0;
    Mat2 F_prev = gramian_integrand(t, Phi);

    for (int k = 0; k < steps; ++k) {
        Mat2 Phi_next = rk4_phi_step(t, Phi, h);
        double t_next = t + h;
        Mat2 F_next = gramian_integrand(t_next, Phi_next);

        W = add(W, scale(add(F_prev, F_next), 0.5 * h));

        Phi = Phi_next;
        F_prev = F_next;
        t = t_next;
    }

    double detW = W[0][0] * W[1][1] - W[0][1] * W[1][0];

    std::cout << "Approximate observability Gramian W_o:\n";
    print_matrix(W);
    std::cout << "det(W_o) = " << std::setprecision(12) << detW << "\n";
    if (detW > 1e-8) {
        std::cout << "The numerical test indicates observability on [0, 6].\n";
    } else {
        std::cout << "The numerical test indicates loss or near-loss of observability.\n";
    }
    return 0;
}

      

The Java version mirrors the C++ algorithm and is useful for students who want to inspect the numerical method without relying on a matrix library.

Chapter14_Lesson5.java

// Chapter14_Lesson5.java
// Continuous-time LTV observability Gramian using RK4 for a 2-state example.

public class Chapter14_Lesson5 {
    static double[][] zeros() {
        return new double[][] { {0.0, 0.0}, {0.0, 0.0} };
    }

    static double[][] identity() {
        return new double[][] { {1.0, 0.0}, {0.0, 1.0} };
    }

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

    static double[][] scale(double[][] A, double s) {
        double[][] C = zeros();
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                C[i][j] = s * A[i][j];
        return C;
    }

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

    static double[][] transpose(double[][] A) {
        return new double[][] { {A[0][0], A[1][0]}, {A[0][1], A[1][1]} };
    }

    static double[][] AofT(double t) {
        return new double[][] {
            {0.0, 1.0},
            {-(2.0 + 0.5 * Math.sin(1.3 * t)), -(0.15 + 0.05 * Math.cos(t))}
        };
    }

    static double[][] CtCofT(double t) {
        double c1 = 1.0;
        double c2 = 0.2 * Math.sin(0.7 * t);
        return new double[][] { {c1 * c1, c1 * c2}, {c2 * c1, c2 * c2} };
    }

    static double[][] phiRhs(double t, double[][] Phi) {
        return multiply(AofT(t), Phi);
    }

    static double[][] rk4PhiStep(double t, double[][] Phi, double h) {
        double[][] k1 = phiRhs(t, Phi);
        double[][] k2 = phiRhs(t + 0.5 * h, add(Phi, scale(k1, 0.5 * h)));
        double[][] k3 = phiRhs(t + 0.5 * h, add(Phi, scale(k2, 0.5 * h)));
        double[][] k4 = phiRhs(t + h, add(Phi, scale(k3, h)));

        double[][] incr = add(add(k1, scale(k2, 2.0)), add(scale(k3, 2.0), k4));
        return add(Phi, scale(incr, h / 6.0));
    }

    static double[][] gramianIntegrand(double t, double[][] Phi) {
        return multiply(multiply(transpose(Phi), CtCofT(t)), Phi);
    }

    static void printMatrix(double[][] M) {
        for (int i = 0; i < 2; i++) {
            System.out.printf("[ %12.8f %12.8f ]%n", M[i][0], M[i][1]);
        }
    }

    public static void main(String[] args) {
        double t0 = 0.0;
        double tf = 6.0;
        int steps = 6000;
        double h = (tf - t0) / steps;

        double[][] Phi = identity();
        double[][] W = zeros();

        double t = t0;
        double[][] Fprev = gramianIntegrand(t, Phi);

        for (int k = 0; k < steps; k++) {
            double[][] PhiNext = rk4PhiStep(t, Phi, h);
            double tNext = t + h;
            double[][] Fnext = gramianIntegrand(tNext, PhiNext);

            W = add(W, scale(add(Fprev, Fnext), 0.5 * h));

            Phi = PhiNext;
            Fprev = Fnext;
            t = tNext;
        }

        double detW = W[0][0] * W[1][1] - W[0][1] * W[1][0];

        System.out.println("Approximate observability Gramian W_o:");
        printMatrix(W);
        System.out.printf("det(W_o) = %.12f%n", detW);
        if (detW > 1e-8) {
            System.out.println("The numerical test indicates observability on [0, 6].");
        } else {
            System.out.println("The numerical test indicates loss or near-loss of observability.");
        }
    }
}

      

The MATLAB script integrates an augmented system containing \( \boldsymbol{\Phi} \), \( \mathbf{W}_o \), and the reconstruction vector \( \mathbf{z} \). In Simulink, the same equations can be implemented with MATLAB Function blocks for \( \mathbf{A}(t) \) and \( \mathbf{C}(t) \), followed by Integrator blocks for \( \dot{\boldsymbol{\Phi} } \), \( \dot{\mathbf{W} } \), and \( \dot{\mathbf{z} } \).

Chapter14_Lesson5.m

% Chapter14_Lesson5.m
% LTV observability Gramian and initial-state reconstruction.
% Simulink note:
%   Implement A(t), C(t), Phi_dot = A(t)Phi, and W_dot = Phi'C'C Phi
%   with MATLAB Function blocks and Integrator blocks. This script gives
%   the numerical reference trajectory for those blocks.

clear; clc;

t0 = 0;
tf = 6;
n = 2;
x0_true = [1; -0.7];

aug0 = [reshape(eye(n), n*n, 1); reshape(zeros(n), n*n, 1); zeros(n,1)];

opts = odeset('RelTol',1e-10,'AbsTol',1e-12);
[t, aug] = ode45(@augmented_rhs, [t0 tf], aug0, opts);

Phi_all = aug(:, 1:n*n);
W_all = aug(:, n*n+1:2*n*n);
z_all = aug(:, 2*n*n+1:end);

W = reshape(W_all(end,:), n, n);
z = z_all(end,:).';
x0_hat = W \ z;

disp('Observability Gramian W_o:');
disp(W);
disp('Eigenvalues of W_o:');
disp(eig(W).');
disp('Condition number:');
disp(cond(W));
disp('True x0:');
disp(x0_true.');
disp('Reconstructed x0:');
disp(x0_hat.');
disp('Reconstruction error norm:');
disp(norm(x0_hat - x0_true));

function d_aug = augmented_rhs(t, aug)
    n = 2;
    Phi = reshape(aug(1:n*n), n, n);

    At = A_of_t(t);
    Ct = C_of_t(t);

    dPhi = At * Phi;
    dW = Phi.' * Ct.' * Ct * Phi;

    y = Ct * Phi * [1; -0.7];  % noiseless homogeneous output
    dz = Phi.' * Ct.' * y;

    d_aug = [reshape(dPhi, n*n, 1); reshape(dW, n*n, 1); dz];
end

function At = A_of_t(t)
    At = [0, 1;
          -(2.0 + 0.5*sin(1.3*t)), -(0.15 + 0.05*cos(t))];
end

function Ct = C_of_t(t)
    Ct = [1.0, 0.2*sin(0.7*t)];
end

      

Mathematica is convenient for symbolic/numeric hybrid workflows. The following notebook-style code computes \( \boldsymbol{\Phi}(t,t_0) \), integrates the Gramian, and solves the reconstruction equation.

Chapter14_Lesson5.nb

(* Chapter14_Lesson5.nb *)
(* Wolfram Mathematica implementation of an LTV observability Gramian. *)

ClearAll["Global`*"];

t0 = 0;
tf = 6;

A[t_] := { {0, 1}, {-(2 + 0.5 Sin[1.3 t]), -(0.15 + 0.05 Cos[t])} };
Cmat[t_] := { {1, 0.2 Sin[0.7 t]} };

phiVars = Array[phi, {2, 2}];
phiMat[t_] := phiVars /. phi[i_, j_] :> phi[i, j][t];

eqnsPhi = Flatten[
   Table[D[phi[i, j][t], t] == (A[t].phiMat[t])[[i, j]], {i, 2}, {j, 2}]
];

icsPhi = Flatten[
   Table[phi[i, j][t0] == IdentityMatrix[2][[i, j]], {i, 2}, {j, 2}]
];

solPhi = NDSolveValue[
   Join[eqnsPhi, icsPhi],
   Flatten[Table[phi[i, j], {i, 2}, {j, 2}]],
   {t, t0, tf},
   WorkingPrecision -> 30
];

Phi[t_?NumericQ] := Partition[solPhi[t], 2];

W = NIntegrate[
   Transpose[Phi[s]].Transpose[Cmat[s]].Cmat[s].Phi[s],
   {s, t0, tf},
   WorkingPrecision -> 30
];

EigenvaluesW = Eigenvalues[W];
ConditionNumberW = Norm[W, 2] Norm[Inverse[W], 2];

x0True = {1, -0.7};

z = NIntegrate[
   Transpose[Phi[s]].Transpose[Cmat[s]].(Cmat[s].Phi[s].x0True),
   {s, t0, tf},
   WorkingPrecision -> 30
];

x0Hat = LinearSolve[W, z];

Print["Observability Gramian W_o = "];
Print[MatrixForm[W]];
Print["Eigenvalues = ", EigenvaluesW];
Print["Condition number = ", ConditionNumberW];
Print["True x0 = ", x0True];
Print["Reconstructed x0 = ", x0Hat];
Print["Error norm = ", Norm[x0Hat - x0True]];

      

Problem 1 (Indistinguishable States): Consider the homogeneous LTV system \( \dot{\mathbf{x} }=\mathbf{A}(t)\mathbf{x} \), \( \mathbf{y}=\mathbf{C}(t)\mathbf{x} \). Prove that two initial states are indistinguishable on \( [t_0,t_f] \) if and only if their difference belongs to \( \mathcal{N}_o[t_0,t_f] \).

Solution: Let \( \boldsymbol{\eta}=\mathbf{x}_a(t_0)-\mathbf{x}_b(t_0) \). The homogeneous outputs differ by

\[ \mathbf{y}_a(t)-\mathbf{y}_b(t)= \mathbf{C}(t)\boldsymbol{\Phi}(t,t_0)\boldsymbol{\eta}. \]

The two outputs are identical for all \( t\in[t_0,t_f] \) exactly when this expression is zero for all times in the interval. This is precisely the defining condition for \( \boldsymbol{\eta}\in\mathcal{N}_o[t_0,t_f] \).

Problem 2 (Gramian for a Rotating Sensor): Let \( \dot{\mathbf{x} }=\mathbf{0} \) and \( y(t)=\begin{bmatrix}\cos t & \sin t\end{bmatrix}\mathbf{x}(0) \) on \( [0,T] \). Compute the observability Gramian and determine when it is nonsingular.

Solution: Since \( \boldsymbol{\Phi}(t,0)=\mathbf{I} \),

\[ \mathbf{W}_o(0,T)= \int_0^T \begin{bmatrix}\cos t\\ \sin t\end{bmatrix} \begin{bmatrix}\cos t & \sin t\end{bmatrix}dt. \]

\[ \mathbf{W}_o(0,T)= \begin{bmatrix} \frac{T}{2}+\frac{\sin 2T}{4} & \frac{1-\cos 2T}{4}\\ \frac{1-\cos 2T}{4} & \frac{T}{2}-\frac{\sin 2T}{4} \end{bmatrix}. \]

Its determinant is

\[ \det\mathbf{W}_o(0,T)= \frac{T^2}{4}-\frac{\sin^2 T}{4}. \]

This determinant is positive for any \( T>0 \) except the limiting case \( T=0 \). Hence any nonzero interval gives two independent sensor directions, although very small \( T \) gives poor conditioning.

Problem 3 (Local Derivative Matrix): For the scalar-output homogeneous system \( \dot{\mathbf{x} }=\mathbf{A}(t)\mathbf{x} \), \( y=\mathbf{C}(t)\mathbf{x} \), derive the first two rows of the local derivative observability matrix.

Solution: The first output equation is \( y=\mathbf{C}(t)\mathbf{x} \), so \( \mathbf{L}_0(t)=\mathbf{C}(t) \). Differentiating,

\[ \dot{y}=\dot{\mathbf{C} }(t)\mathbf{x}(t)+ \mathbf{C}(t)\dot{\mathbf{x} }(t)= \left(\dot{\mathbf{C} }(t)+\mathbf{C}(t)\mathbf{A}(t)\right)\mathbf{x}(t). \]

Therefore \( \mathbf{L}_1(t)=\dot{\mathbf{C} }(t)+\mathbf{C}(t)\mathbf{A}(t) \). The next derivative row is

\[ \mathbf{L}_2(t)=\dot{\mathbf{L} }_1(t)+\mathbf{L}_1(t)\mathbf{A}(t). \]

Problem 4 (Dual Gramian Rank): Show that if \( \mathbf{W}_{c,d}(t_0,t_f)= \boldsymbol{\Phi}^{T}(t_0,t_f)\mathbf{W}_o(t_0,t_f) \boldsymbol{\Phi}(t_0,t_f) \), then \( \operatorname{rank}\mathbf{W}_{c,d}= \operatorname{rank}\mathbf{W}_o \).

Solution: The state-transition matrix is nonsingular, with inverse \( \boldsymbol{\Phi}(t_f,t_0) \). Multiplication by a nonsingular matrix on the left or right does not change matrix rank. Therefore,

\[ \operatorname{rank} \left(\boldsymbol{\Phi}^{T}(t_0,t_f)\mathbf{W}_o \boldsymbol{\Phi}(t_0,t_f)\right) = \operatorname{rank}\mathbf{W}_o. \]

Hence the dual controllability Gramian is full rank exactly when the original observability Gramian is full rank.

Problem 5 (Reconstruction with a Gramian): Suppose \( \mathbf{W}_o(t_0,t_f) \) is nonsingular and the homogeneous output is known. Derive the linear system used to reconstruct \( \mathbf{x}(t_0) \).

Solution: Starting from \( \tilde{\mathbf{y} }(t)=\mathbf{C}(t)\boldsymbol{\Phi}(t,t_0) \mathbf{x}(t_0) \), multiply by \( \boldsymbol{\Phi}^{T}(t,t_0)\mathbf{C}^{T}(t) \) and integrate:

\[ \int_{t_0}^{t_f} \boldsymbol{\Phi}^{T}(t,t_0)\mathbf{C}^{T}(t)\tilde{\mathbf{y} }(t)\,dt = \mathbf{W}_o(t_0,t_f)\mathbf{x}(t_0). \]

Let the left-hand side be \( \mathbf{z} \). Then solve \( \mathbf{W}_o\mathbf{x}(t_0)=\mathbf{z} \). In numerical computation, solving the linear system is preferred over explicitly computing \( \mathbf{W}_o^{-1} \).

Time-varying observability is an interval property: the output history must distinguish all possible initial states. The decisive object is the map \( \mathbf{C}(t)\boldsymbol{\Phi}(t,t_0) \), whose accumulated information is the finite-interval observability Gramian. Full rank of \( \mathbf{W}_o(t_0,t_f) \) is equivalent to unique initial-state reconstruction. The duality with time-varying controllability survives, but the dual dynamics involve \( -\mathbf{A}^{T}(t) \) and endpoint transformations. The next chapter develops observability Gramians and output energy in greater depth.