In deterministic modeling, sensors and actuators are often represented as ideal algebraic or dynamic blocks. In practice, however, every measurement and every commanded input is corrupted by stochastic effects. The key modeling question is not merely “noise exists,” but where it enters the dynamics and how it transforms through the system. A physically meaningful noise model must distinguish:

\( w(t) \): process disturbances entering plant dynamics (unmodeled forces, torque ripple, etc.), \( n_a(t) \): actuator-channel noise or uncertainty before the plant, \( v(t) \): sensor white/electronic noise, \( b(t) \): sensor bias drift, and \( q(t) \): quantization noise.

For an LTI plant, the same noise source can produce very different output statistics depending on whether it enters through the state equation or the measurement equation. This is why the structure of the stochastic model matters as much as the marginal variance values.

flowchart TD
  UC["Command u_c"] --> ACT["Actuator dynamics"]
  NA["Actuator noise n_a"] --> SUMA["Add at actuator input"]
  ACT --> SUMA
  SUMA --> PLANT["Plant dynamics"]
  WP["Process disturbance w"] --> PLANT
  PLANT --> XS["State x"]
  XS --> SEN["Sensor map Cx"]
  VB["Bias drift b"] --> SUMY["Measurement assembly"]
  VV["Sensor white noise v"] --> SUMY
  QQ["Quantization q"] --> SUMY
  SEN --> SUMY
  SUMY --> YM["Measured output y_m"]
        

A standard continuous-time stochastic state-space model for a dynamic plant with noisy actuator and noisy sensor is:

\[ \dot{\mathbf{x} }(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}\big(\mathbf{u}_c(t) + \mathbf{n}_a(t)\big) + \mathbf{L}\mathbf{w}(t), \qquad \mathbf{y}_m(t) = \mathbf{C}\mathbf{x}(t) + \mathbf{D}\mathbf{u}_c(t) + \mathbf{b}(t) + \mathbf{v}(t) + \mathbf{q}(t) \]

Here \( \mathbf{w}(t) \) and \( \mathbf{v}(t) \) are often idealized as zero-mean white processes (with spectral intensity matrices), while \( \mathbf{n}_a(t) \) and \( \mathbf{b}(t) \) are frequently colored processes with memory.

After sampling with period \( T_s \), a practical discrete-time model is

\[ \mathbf{x}_{k+1} = \mathbf{A}_d\mathbf{x}_k + \mathbf{B}_d\mathbf{u}_{c,k} + \mathbf{B}_d\mathbf{n}_{a,k} + \mathbf{G}_d\mathbf{w}_k, \qquad \mathbf{y}_{m,k} = \mathbf{C}\mathbf{x}_k + \mathbf{b}_k + \mathbf{v}_k + \mathbf{q}_k \]

The decomposition is intentionally explicit: actuator noise enters through the same channel as the control command (scaled by \( \mathbf{B}_d \)), while sensor noise appears additively at the measurement output.

A common bias-drift model is a random walk:

\[ \mathbf{b}_{k+1} = \mathbf{b}_k + \boldsymbol{\eta}_{b,k}, \qquad E[\boldsymbol{\eta}_{b,k}] = \mathbf{0}, \qquad E[\boldsymbol{\eta}_{b,k}\boldsymbol{\eta}_{b,k}^{\top}] = \mathbf{Q}_b \]

and a common actuator colored-noise model is an AR(1) (discrete first-order shaping):

\[ \mathbf{n}_{a,k+1} = \boldsymbol{\rho}_a \mathbf{n}_{a,k} + \boldsymbol{\varepsilon}_{a,k}, \qquad E[\boldsymbol{\varepsilon}_{a,k}] = \mathbf{0}, \qquad E[\boldsymbol{\varepsilon}_{a,k}\boldsymbol{\varepsilon}_{a,k}^{\top}] = \mathbf{Q}_{\varepsilon a} \]

with stability requiring each scalar mode to satisfy \( |\rho_i|<1 \).

Let \( \mathbf{m}_k = E[\mathbf{x}_k] \). Taking expectation of the discrete-time state equation gives

\[ \mathbf{m}_{k+1} = \mathbf{A}_d\mathbf{m}_k + \mathbf{B}_d\mathbf{u}_{c,k} + \mathbf{B}_d E[\mathbf{n}_{a,k}] + \mathbf{G}_d E[\mathbf{w}_k] \]

Therefore, if \( E[\mathbf{n}_{a,k}] = \mathbf{0} \) and \( E[\mathbf{w}_k] = \mathbf{0} \), the state mean evolves exactly as in the deterministic model. However, nonzero actuator bias immediately shifts the state mean.

For measurements:

\[ E[\mathbf{y}_{m,k}] = \mathbf{C}\mathbf{m}_k + E[\mathbf{b}_k] + E[\mathbf{v}_k] + E[\mathbf{q}_k] \]

Even when sensor white noise is zero-mean, a drifting bias \( \mathbf{b}_k \) produces a slowly varying offset. This distinction is crucial in control and estimation: white noise increases uncertainty, but bias drift causes persistent systematic error.

If the bias is modeled as a random walk with zero-mean increments and initial mean \( E[\mathbf{b}_0] \), then

\[ E[\mathbf{b}_k] = E[\mathbf{b}_0] \]

but its covariance grows with time (derived in later sections), which is why bias noise is often called a nonstationary measurement component.

Define the centered state \( \tilde{\mathbf{x} }_k = \mathbf{x}_k - \mathbf{m}_k \). Under the independence assumptions used in basic stochastic system analysis (independent state/noise terms and zero cross-covariances), the centered dynamics become

\[ \tilde{\mathbf{x} }_{k+1} = \mathbf{A}_d\tilde{\mathbf{x} }_k + \mathbf{B}_d\tilde{\mathbf{n} }_{a,k} + \mathbf{G}_d\tilde{\mathbf{w} }_k \]

Let \( \mathbf{P}_k = E[\tilde{\mathbf{x} }_k\tilde{\mathbf{x} }_k^{\top}] \), \( \mathbf{Q}_a = E[\tilde{\mathbf{n} }_{a,k}\tilde{\mathbf{n} }_{a,k}^{\top}] \), and \( \mathbf{Q}_w = E[\tilde{\mathbf{w} }_k\tilde{\mathbf{w} }_k^{\top}] \). Expanding the covariance and dropping zero cross-terms yields

\[ \mathbf{P}_{k+1} = \mathbf{A}_d\mathbf{P}_k\mathbf{A}_d^{\top} + \mathbf{B}_d\mathbf{Q}_a\mathbf{B}_d^{\top} + \mathbf{G}_d\mathbf{Q}_w\mathbf{G}_d^{\top} \]

Proof sketch: Write \( \mathbf{P}_{k+1}=E[(\cdots)(\cdots)^{\top}] \), expand all nine terms, and use \( E[\tilde{\mathbf{x} }_k\tilde{\mathbf{n} }_{a,k}^{\top}] = \mathbf{0} \), \( E[\tilde{\mathbf{x} }_k\tilde{\mathbf{w} }_{k}^{\top}] = \mathbf{0} \), and \( E[\tilde{\mathbf{n} }_{a,k}\tilde{\mathbf{w} }_{k}^{\top}] = \mathbf{0} \). The remaining three terms are exactly the three matrix terms above.

For measurements, define the centered measurement \( \tilde{\mathbf{y} }_{m,k} = \mathbf{y}_{m,k} - E[\mathbf{y}_{m,k}] \). If \( \mathbf{b}_k, \mathbf{v}_k, \mathbf{q}_k \) are independent of \( \mathbf{x}_k \), then the measurement covariance is

\[ \mathbf{S}_{y,k} = E[\tilde{\mathbf{y} }_{m,k}\tilde{\mathbf{y} }_{m,k}^{\top}] = \mathbf{C}\mathbf{P}_k\mathbf{C}^{\top} + \mathbf{R}_{b,k} + \mathbf{R}_v + \mathbf{R}_q \]

where \( \mathbf{R}_{b,k}=E[(\mathbf{b}_k-E[\mathbf{b}_k])(\mathbf{b}_k-E[\mathbf{b}_k])^{\top}] \) may depend on \( k \) (nonstationary bias drift). For a random walk bias:

\[ \mathbf{R}_{b,k} = \mathbf{R}_{b,0} + k\,\mathbf{Q}_b \]

so the measurement variance can grow over time even when the plant itself is stable.

From Lesson 2 and Lesson 3, an LTI system transforms noise spectra through the squared magnitude of transfer functions. For scalar channels, if \( H_w(s) \) maps process noise to output and \( H_a(s) \) maps actuator noise to output, then

\[ S_{yy}(\omega) = |H_w(j\omega)|^2 S_{ww}(\omega) + |H_a(j\omega)|^2 S_{aa}(\omega) + S_{vv}(\omega) + S_{qq}(\omega) \]

This decomposition is powerful because it separates source statistics from system dynamics. The same sensor noise level can be harmless in a low-bandwidth system and harmful in a resonant system with large high-frequency gain.

Colored noise is often generated by passing white noise through a shaping filter. A first-order shaping model is

\[ \dot{\xi}(t) = -\frac{1}{T_n}\xi(t) + \eta(t), \qquad n(t)=\xi(t) \]

with transfer function \( H_n(s)=\dfrac{1}{s + 1/T_n} \). If \( \eta(t) \) is white, then \( n(t) \) is colored with low-pass PSD shape.

The discrete-time analogue is AR(1):

\[ n_{k+1} = \rho n_k + \varepsilon_k, \qquad |\rho|<1 \]

and the stationary variance is obtained by taking variance on both sides:

\[ \sigma_n^2 = \rho^2 \sigma_n^2 + \sigma_{\varepsilon}^2 \quad \Rightarrow \quad \sigma_{\varepsilon}^2 = (1-\rho^2)\sigma_n^2 \]

This formula is used directly in the implementation codes to generate actuator colored noise with a prescribed steady-state variance.

A realistic sensor model usually combines multiple error sources:

\[ y_{m,k} = (1+s_f)\,y_k + b_k + v_k + q_k + n_{j,k} \]

where \( s_f \) is a scale-factor error, \( b_k \) is bias drift, \( v_k \) is white/electronic noise, \( q_k \) is quantization noise, and \( n_{j,k} \) captures sampling-jitter-induced error (if relevant).

6.1 Quantization Noise Model and Variance

If a uniform quantizer with step \( \Delta \) operates in a regime where the signal spans many quantization bins and is not synchronized with the grid, the quantization error is approximated as \( q_k \sim \mathcal{U}(-\Delta/2,\Delta/2) \). Then

\[ E[q_k]=0, \qquad \operatorname{Var}(q_k)=\int_{-\Delta/2}^{\Delta/2} q^2 \frac{1}{\Delta}\,dq = \frac{\Delta^2}{12} \]

This approximation breaks down for slowly varying signals or deterministic periodic signals that lock to quantization levels; in those cases, the “noise” becomes signal-dependent.

6.2 Bias Drift as a Random Walk

Bias drift is frequently modeled as

\[ b_{k+1}=b_k+\eta_{b,k}, \qquad \eta_{b,k}\sim \mathcal{N}(0,Q_b) \]

which implies

\[ \operatorname{Var}(b_k)=\operatorname{Var}(b_0)+kQ_b \]

The growth with \( k \) explains why long-duration measurements can drift significantly even when short-term noise appears small.

6.3 Sampling Jitter Linearization

Suppose a sensor samples at perturbed times \( t_k + \delta t_k \), where \( E[\delta t_k]=0 \) and \( \operatorname{Var}(\delta t_k)=\sigma_t^2 \). First-order Taylor expansion gives

\[ y(t_k+\delta t_k) \approx y(t_k) + \dot{y}(t_k)\,\delta t_k \]

Therefore the jitter-induced variance is approximately

\[ \operatorname{Var}(n_{j,k}) \approx \dot{y}(t_k)^2\,\sigma_t^2 \]

so fast transients are more sensitive to timing jitter than slowly varying regimes.

Actuator imperfections should not be modeled only as “output measurement noise.” They act on the plant through the input channel and can alter the state trajectory directly. A useful structure is

\[ \dot{\mathbf{z} }_a = \mathbf{A}_a \mathbf{z}_a + \mathbf{B}_a \mathbf{u}_c + \mathbf{L}_a \boldsymbol{\eta}_a, \qquad \mathbf{u}_{act} = \mathbf{C}_a \mathbf{z}_a \]

where \( \mathbf{z}_a \) represents actuator internal dynamics (drive amplifier, valve spool, motor current loop, etc.), and \( \boldsymbol{\eta}_a \) is a white excitation that generates colored actuator noise at the output.

The plant then sees

\[ \dot{\mathbf{x} } = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u}_{act} + \mathbf{L}\mathbf{w} \]

which is equivalent to augmenting the state with actuator dynamics. This is often more accurate than inserting ad hoc white noise directly on the plant state.

In addition, actuator saturation and dead-zone (introduced in nonlinear-system chapters) can coexist with stochastic disturbances. Around an operating point, a local linearization gives an effective stochastic gain:

\[ u_{act} \approx \bar{u}_{act} + g_{eff}\,(\delta u_c + \delta n_a) \]

so the input-noise covariance seen by the plant is scaled by \( g_{eff}^2 \).

When a continuous-time model is known, the discrete covariance matrices used in sampled-data simulation should be derived consistently. For

\[ \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{L}\mathbf{w}(t), \qquad E[\mathbf{w}(t)\mathbf{w}(\lambda)^{\top}] = \mathbf{Q}_c\,\delta(t-\lambda) \]

the exact discrete-time covariance over one sample interval is

\[ \mathbf{Q}_d = \int_{0}^{T_s} e^{\mathbf{A}\sigma}\,\mathbf{L}\mathbf{Q}_c\mathbf{L}^{\top}\, e^{\mathbf{A}^{\top}\sigma}\,d\sigma \]

This formula follows from the variation-of-constants solution and the white-noise covariance property. In practice, \( \mathbf{Q}_d \) may be computed numerically, by matrix-fraction methods, or approximated for small \( T_s \) using

\[ \mathbf{Q}_d \approx \mathbf{L}\mathbf{Q}_c\mathbf{L}^{\top}T_s \]

(first-order approximation). The same principle applies to actuator-noise shaping states if they are modeled continuously.

A common mistake is to set \( \mathbf{Q}_d = \mathbf{Q}_c \) directly. This is dimensionally incorrect: continuous-time intensities and discrete-time covariances are not the same object.

In system dynamics and control engineering, noise modeling should be integrated with the physical modeling workflow rather than appended at the end. The recommended sequence is:

1. Build the deterministic plant and actuator/sensor dynamic blocks.
2. Identify noise injection locations (process, actuator, sensor, timing).
3. Choose source models (white, random walk, AR(1), shaping-filter state).
4. Map source statistics to discrete covariances and verify units.
5. Validate the model using sample variance, autocorrelation, and PSD checks.

flowchart TD
  A["Deterministic plant model"] --> B["Add actuator and sensor dynamics"]
  B --> C["Choose noise injection points"]
  C --> D["Select source models: white / bias RW / AR1"]
  D --> E["Map to discrete covariances Q and R"]
  E --> F["Simulate and estimate stats"]
  F --> G["Compare mean, variance, PSD with targets"]
  G --> H["Refine parameters and assumptions"]
        

This workflow is deliberately aligned with the material from Chapters 15 and 16 (simulation and sampling) and Lessons 1–3 of Chapter 17 (stationarity, correlation, PSD, and covariance propagation).

The following implementations simulate the same pedagogical example: a mass–spring–damper plant with colored actuator noise, process noise, sensor white noise, sensor bias random walk, and quantization. Each implementation also computes a covariance recursion prediction

\[ \mathbf{P}_{k+1} = \mathbf{A}_d\mathbf{P}_k\mathbf{A}_d^{\top} + \mathbf{G}_d\mathbf{Q}_w\mathbf{G}_d^{\top} + \mathbf{B}_d\mathbf{Q}_a\mathbf{B}_d^{\top} \]

and compares the predicted measurement variance against empirical variance from simulation.

Chapter17_Lesson4.py

# Chapter17_Lesson4.py
# Noise Modeling in Sensors and Actuators for Dynamic Systems
# Python simulation: mass-spring-damper with actuator colored noise,
# process noise, sensor bias drift, measurement noise, and quantization.

import numpy as np

np.random.seed(17)

# Physical parameters
m = 1.2
c = 0.35
k = 4.0
Ts = 0.01
N = 6000

# Continuous-time matrices (x = [position, velocity]^T)
A_c = np.array([[0.0, 1.0],
                [-k/m, -c/m]])
B_c = np.array([[0.0],
                [1.0/m]])
G_c = np.array([[0.0],
                [1.0/m]])  # process disturbance input channel
C = np.array([[1.0, 0.0]]) # position sensor

# Simple discretization (forward Euler for teaching/demo)
A = np.eye(2) + Ts * A_c
B = Ts * B_c
G = Ts * G_c

# Noise parameters
sigma_w = 0.20             # process force noise std (white)
sigma_v = 0.01             # sensor white noise std
sigma_bias_rw = 0.002      # bias random-walk intensity (per sqrt(s))
delta_q = 0.001            # quantization step
rho_a = 0.97               # actuator colored-noise AR(1) coefficient
sigma_a_ss = 0.08          # actuator noise steady-state std

# Covariances
Qw = np.array([[sigma_w**2]])
Qa = np.array([[sigma_a_ss**2]])
Rv = np.array([[sigma_v**2]])
Rq = np.array([[delta_q**2 / 12.0]])  # uniform quantization model

# Input command
t = np.arange(N) * Ts
u_cmd = 0.8 * np.sin(2.0 * np.pi * 0.7 * t)

# Storage
x = np.zeros((2, N))
y = np.zeros(N)
y_true = np.zeros(N)
bias = 0.0
a_noise = 0.0

# Covariance recursion (mean-subtracted, independent noises)
P = np.zeros((2, 2))
P_hist = np.zeros((N, 2, 2))
Syy_hist = np.zeros(N)

for k_idx in range(N - 1):
    # Colored actuator noise (AR(1))
    a_noise = rho_a * a_noise + np.sqrt(1.0 - rho_a**2) * sigma_a_ss * np.random.randn()

    # Process white force noise
    w_k = sigma_w * np.random.randn()

    # Plant update
    u_actual = u_cmd[k_idx] + a_noise
    x[:, k_idx + 1] = (A @ x[:, k_idx]
                       + (B[:, 0] * u_actual)
                       + (G[:, 0] * w_k))

    # Sensor bias drift (random walk)
    bias = bias + sigma_bias_rw * np.sqrt(Ts) * np.random.randn()

    # Sensor white noise + quantization
    v_k = sigma_v * np.random.randn()
    y_analog = (C @ x[:, k_idx + 1])[0] + bias + v_k
    q_k = delta_q * np.round(y_analog / delta_q) - y_analog
    y[k_idx + 1] = y_analog + q_k
    y_true[k_idx + 1] = (C @ x[:, k_idx + 1])[0]

    # Theoretical covariance propagation for the plant state
    # P_{k+1} = A P_k A^T + G Qw G^T + B Qa B^T
    P = A @ P @ A.T + G @ Qw @ G.T + B @ Qa @ B.T
    P_hist[k_idx + 1] = P

    # Predicted measurement variance excluding bias random walk accumulation
    Syy = C @ P @ C.T + Rv + Rq
    Syy_hist[k_idx + 1] = Syy[0, 0]

# Empirical statistics after transient
burn = 1000
y_centered = y[burn:] - np.mean(y[burn:])
y_true_centered = y_true[burn:] - np.mean(y_true[burn:])

emp_var_y = np.var(y_centered, ddof=1)
emp_var_ytrue = np.var(y_true_centered, ddof=1)
pred_var_y = np.mean(Syy_hist[burn:])  # excludes long-term bias random walk growth

# Estimate quantization variance empirically
q_est = y[burn:] - (np.round(y[burn:] / delta_q) * delta_q)  # residual in [-delta_q/2, delta_q/2]
q_est = q_est - np.mean(q_est)
emp_var_q = np.var(q_est, ddof=1)

print("=== Noise Modeling Demo (Python) ===")
print(f"Sampling time Ts = {Ts:.4f} s, N = {N}")
print(f"Empirical var(y_true)        : {emp_var_ytrue:.6e}")
print(f"Empirical var(y_measured)    : {emp_var_y:.6e}")
print(f"Predicted var(y) (no bias RW): {pred_var_y:.6e}")
print(f"Quantization var theory      : {delta_q**2/12.0:.6e}")
print(f"Quantization var empirical   : {emp_var_q:.6e}")

# Optional simple PSD-like diagnostic using FFT (periodogram)
Y = np.fft.rfft(y_centered)
freqs = np.fft.rfftfreq(len(y_centered), d=Ts)
psd = (np.abs(Y) ** 2) / (len(y_centered) / Ts)

# Print a few dominant frequencies
idx = np.argsort(psd)[-5:][::-1]
print("\nTop spectral peaks (Hz, PSD):")
for i in idx:
    print(f"{freqs[i]:8.3f} Hz, {psd[i]:.6e}")
      

Chapter17_Lesson4.cpp

// Chapter17_Lesson4.cpp
// Noise Modeling in Sensors and Actuators for Dynamic Systems
// C++17 simulation of a mass-spring-damper with actuator and sensor noise.

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <random>
#include <vector>

int main() {
    const double m = 1.2;
    const double c = 0.35;
    const double k = 4.0;
    const double Ts = 0.01;
    const int N = 6000;

    // Noise parameters
    const double sigma_w = 0.20;       // process force noise std
    const double sigma_v = 0.01;       // sensor white noise std
    const double sigma_bias_rw = 0.002;// bias random-walk intensity
    const double delta_q = 0.001;      // quantizer step
    const double rho_a = 0.97;         // actuator colored-noise coefficient
    const double sigma_a_ss = 0.08;    // actuator steady-state std

    const double PI = 3.14159265358979323846;

    std::mt19937 rng(17);
    std::normal_distribution<double> n01(0.0, 1.0);

    // States: x1 = position, x2 = velocity
    double x1 = 0.0, x2 = 0.0;
    double bias = 0.0;
    double aNoise = 0.0;

    std::vector<double> y(N, 0.0), yTrue(N, 0.0);

    // Theoretical covariance recursion for x = [x1, x2]^T
    // A = I + Ts * [[0,1],[-k/m,-c/m]], B = Ts * [[0],[1/m]], G = B
    const double A11 = 1.0;
    const double A12 = Ts;
    const double A21 = -Ts * k / m;
    const double A22 = 1.0 - Ts * c / m;
    const double B1 = 0.0;
    const double B2 = Ts / m;
    const double G1 = 0.0;
    const double G2 = Ts / m;

    // P = [[p11,p12],[p12,p22]]
    double p11 = 0.0, p12 = 0.0, p22 = 0.0;
    const double Qw = sigma_w * sigma_w;
    const double Qa = sigma_a_ss * sigma_a_ss;
    const double Rv = sigma_v * sigma_v;
    const double Rq = (delta_q * delta_q) / 12.0;
    std::vector<double> Syy(N, 0.0);

    for (int kIdx = 0; kIdx < N - 1; ++kIdx) {
        double t = kIdx * Ts;
        double uCmd = 0.8 * std::sin(2.0 * PI * 0.7 * t);

        // Actuator colored noise (AR(1))
        aNoise = rho_a * aNoise + std::sqrt(1.0 - rho_a * rho_a) * sigma_a_ss * n01(rng);

        // Process and sensor noise
        double w = sigma_w * n01(rng);
        bias += sigma_bias_rw * std::sqrt(Ts) * n01(rng);
        double v = sigma_v * n01(rng);

        // Plant update
        double uActual = uCmd + aNoise;
        double x1n = A11 * x1 + A12 * x2 + B1 * uActual + G1 * w;
        double x2n = A21 * x1 + A22 * x2 + B2 * uActual + G2 * w;
        x1 = x1n;
        x2 = x2n;

        // Sensor model
        double yAnalog = x1 + bias + v;
        double q = delta_q * std::round(yAnalog / delta_q) - yAnalog;
        double yMeas = yAnalog + q;

        y[kIdx + 1] = yMeas;
        yTrue[kIdx + 1] = x1;

        // Covariance update P <- A P A^T + GQwG^T + BQaB^T
        double ap11 = A11 * p11 + A12 * p12;
        double ap12 = A11 * p12 + A12 * p22;
        double ap21 = A21 * p11 + A22 * p12;
        double ap22 = A21 * p12 + A22 * p22;

        double p11n = ap11 * A11 + ap12 * A12 + G1 * Qw * G1 + B1 * Qa * B1;
        double p12n = ap11 * A21 + ap12 * A22 + G1 * Qw * G2 + B1 * Qa * B2;
        double p22n = ap21 * A21 + ap22 * A22 + G2 * Qw * G2 + B2 * Qa * B2;

        p11 = p11n;
        p12 = p12n;
        p22 = p22n;

        Syy[kIdx + 1] = p11 + Rv + Rq; // C = [1, 0]
    }

    // Empirical variance after burn-in
    const int burn = 1000;
    int M = N - burn;
    double meanY = 0.0, meanYTrue = 0.0, meanSyy = 0.0;
    for (int i = burn; i < N; ++i) {
        meanY += y[i];
        meanYTrue += yTrue[i];
        meanSyy += Syy[i];
    }
    meanY /= M;
    meanYTrue /= M;
    meanSyy /= M;

    double varY = 0.0, varYTrue = 0.0;
    for (int i = burn; i < N; ++i) {
        varY += (y[i] - meanY) * (y[i] - meanY);
        varYTrue += (yTrue[i] - meanYTrue) * (yTrue[i] - meanYTrue);
    }
    varY /= (M - 1);
    varYTrue /= (M - 1);

    std::cout << "=== Noise Modeling Demo (C++) ===\n";
    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Empirical var(y_true): " << varYTrue << "\n";
    std::cout << "Empirical var(y)    : " << varY << "\n";
    std::cout << "Predicted mean Syy  : " << meanSyy << " (bias RW excluded)\n";
    std::cout << "Quantization theory : " << (delta_q * delta_q / 12.0) << "\n";

    return 0;
}
      

Chapter17_Lesson4.java

// Chapter17_Lesson4.java
// Noise Modeling in Sensors and Actuators for Dynamic Systems
// Java simulation of a noisy mass-spring-damper sensor/actuator chain.

import java.util.Random;

public class Chapter17_Lesson4 {

    static class Gaussian {
        private final Random rng;
        Gaussian(long seed) { rng = new Random(seed); }
        double next() {
            // Box-Muller
            double u1 = Math.max(1e-12, rng.nextDouble());
            double u2 = rng.nextDouble();
            return Math.sqrt(-2.0 * Math.log(u1)) * Math.cos(2.0 * Math.PI * u2);
        }
    }

    public static void main(String[] args) {
        final double m = 1.2, c = 0.35, k = 4.0;
        final double Ts = 0.01;
        final int N = 6000;

        final double sigmaW = 0.20;
        final double sigmaV = 0.01;
        final double sigmaBiasRW = 0.002;
        final double deltaQ = 0.001;
        final double rhoA = 0.97;
        final double sigmaASS = 0.08;

        Gaussian g = new Gaussian(17L);

        double x1 = 0.0, x2 = 0.0;
        double bias = 0.0, aNoise = 0.0;

        double[] y = new double[N];
        double[] yTrue = new double[N];
        double[] syy = new double[N];

        // Discrete model (Euler)
        double A11 = 1.0;
        double A12 = Ts;
        double A21 = -Ts * k / m;
        double A22 = 1.0 - Ts * c / m;
        double B1 = 0.0;
        double B2 = Ts / m;
        double G1 = 0.0;
        double G2 = Ts / m;

        // Covariance recursion
        double p11 = 0.0, p12 = 0.0, p22 = 0.0;
        double Qw = sigmaW * sigmaW;
        double Qa = sigmaASS * sigmaASS;
        double Rv = sigmaV * sigmaV;
        double Rq = deltaQ * deltaQ / 12.0;

        for (int kIdx = 0; kIdx < N - 1; kIdx++) {
            double t = kIdx * Ts;
            double uCmd = 0.8 * Math.sin(2.0 * Math.PI * 0.7 * t);

            aNoise = rhoA * aNoise + Math.sqrt(1.0 - rhoA * rhoA) * sigmaASS * g.next();
            double w = sigmaW * g.next();
            bias += sigmaBiasRW * Math.sqrt(Ts) * g.next();
            double v = sigmaV * g.next();

            double uActual = uCmd + aNoise;
            double x1n = A11 * x1 + A12 * x2 + B1 * uActual + G1 * w;
            double x2n = A21 * x1 + A22 * x2 + B2 * uActual + G2 * w;
            x1 = x1n;
            x2 = x2n;

            double yAnalog = x1 + bias + v;
            double q = deltaQ * Math.rint(yAnalog / deltaQ) - yAnalog;
            y[kIdx + 1] = yAnalog + q;
            yTrue[kIdx + 1] = x1;

            // P <- A P A^T + GQG^T + BQaB^T
            double ap11 = A11 * p11 + A12 * p12;
            double ap12 = A11 * p12 + A12 * p22;
            double ap21 = A21 * p11 + A22 * p12;
            double ap22 = A21 * p12 + A22 * p22;

            double p11n = ap11 * A11 + ap12 * A12 + G1 * Qw * G1 + B1 * Qa * B1;
            double p12n = ap11 * A21 + ap12 * A22 + G1 * Qw * G2 + B1 * Qa * B2;
            double p22n = ap21 * A21 + ap22 * A22 + G2 * Qw * G2 + B2 * Qa * B2;

            p11 = p11n;
            p12 = p12n;
            p22 = p22n;

            syy[kIdx + 1] = p11 + Rv + Rq;
        }

        int burn = 1000;
        int M = N - burn;
        double meanY = 0.0, meanYTrue = 0.0, meanSyy = 0.0;
        for (int i = burn; i < N; i++) {
            meanY += y[i];
            meanYTrue += yTrue[i];
            meanSyy += syy[i];
        }
        meanY /= M;
        meanYTrue /= M;
        meanSyy /= M;

        double varY = 0.0, varYTrue = 0.0;
        for (int i = burn; i < N; i++) {
            varY += (y[i] - meanY) * (y[i] - meanY);
            varYTrue += (yTrue[i] - meanYTrue) * (yTrue[i] - meanYTrue);
        }
        varY /= (M - 1);
        varYTrue /= (M - 1);

        System.out.println("=== Noise Modeling Demo (Java) ===");
        System.out.printf("Empirical var(y_true): %.6e%n", varYTrue);
        System.out.printf("Empirical var(y)    : %.6e%n", varY);
        System.out.printf("Predicted mean Syy  : %.6e (bias RW excluded)%n", meanSyy);
        System.out.printf("Quantization theory : %.6e%n", (deltaQ * deltaQ / 12.0));
    }
}
      

Chapter17_Lesson4.m

The MATLAB script includes a direct simulation and, at the end, a concise Simulink implementation checklist (State-Space block, noise blocks, quantizer, and scope wiring).

% Chapter17_Lesson4.m
% Noise Modeling in Sensors and Actuators for Dynamic Systems
% MATLAB script: state-space simulation with actuator colored noise,
% process noise, sensor bias drift, measurement noise, and quantization.

clear; clc; rng(17);

% Physical parameters
m = 1.2; c = 0.35; k = 4.0;
Ts = 0.01; N = 6000;
t = (0:N-1)' * Ts;

% Continuous-time matrices
A_c = [0 1; -k/m -c/m];
B_c = [0; 1/m];
G_c = [0; 1/m];
C = [1 0];

% Discretization (Euler for teaching/demo)
A = eye(2) + Ts*A_c;
B = Ts*B_c;
G = Ts*G_c;

% Noise parameters
sigma_w = 0.20;
sigma_v = 0.01;
sigma_bias_rw = 0.002;
delta_q = 0.001;
rho_a = 0.97;
sigma_a_ss = 0.08;

Qw = sigma_w^2;
Qa = sigma_a_ss^2;
Rv = sigma_v^2;
Rq = delta_q^2/12;

u_cmd = 0.8*sin(2*pi*0.7*t);

x = zeros(2,N);
y = zeros(N,1);
y_true = zeros(N,1);
Syy = zeros(N,1);

bias = 0;
aNoise = 0;
P = zeros(2,2);

for kIdx = 1:N-1
    % Colored actuator noise (AR(1))
    aNoise = rho_a*aNoise + sqrt(1-rho_a^2)*sigma_a_ss*randn;

    % Process noise
    w = sigma_w * randn;

    % Plant update
    u_actual = u_cmd(kIdx) + aNoise;
    x(:,kIdx+1) = A*x(:,kIdx) + B*u_actual + G*w;

    % Sensor bias random walk
    bias = bias + sigma_bias_rw*sqrt(Ts)*randn;

    % White measurement noise + quantization
    v = sigma_v * randn;
    y_analog = C*x(:,kIdx+1) + bias + v;
    q = delta_q*round(y_analog/delta_q) - y_analog;
    y(kIdx+1) = y_analog + q;
    y_true(kIdx+1) = C*x(:,kIdx+1);

    % Covariance recursion (independent noises)
    P = A*P*A' + G*Qw*G' + B*Qa*B';
    Syy(kIdx+1) = C*P*C' + Rv + Rq;
end

burn = 1000;
idx = burn:N;
empVarY = var(y(idx),1);
empVarYTrue = var(y_true(idx),1);
predVarY = mean(Syy(idx));

fprintf('=== Noise Modeling Demo (MATLAB) ===\n');
fprintf('Empirical var(y_true): %.6e\n', empVarYTrue);
fprintf('Empirical var(y)    : %.6e\n', empVarY);
fprintf('Predicted mean Syy  : %.6e (bias RW excluded)\n', predVarY);
fprintf('Quantization theory : %.6e\n', delta_q^2/12);

% Simple PSD estimate (periodogram)
yc = y(idx) - mean(y(idx));
Y = fft(yc);
L = length(yc);
f = (0:floor(L/2))'/(L*Ts);
Pyy = (abs(Y(1:floor(L/2)+1)).^2)/(L/Ts);

[~, order] = sort(Pyy, 'descend');
fprintf('\nTop spectral peaks (Hz, PSD):\n');
for i = 1:5
    ii = order(i);
    fprintf('%8.3f Hz, %.6e\n', f(ii), Pyy(ii));
end

% Simulink implementation notes (manual setup):
% 1) Use a State-Space block with A_c, B_c, C, D=0.
% 2) Add Band-Limited White Noise for process force and Sum into actuator/plant input.
% 3) Add a Discrete Filter block for actuator AR(1) noise (z-domain: 1/(1-rho_a z^-1)).
% 4) Add a Unit Delay + Random Number to create bias random walk.
% 5) Add Quantizer block with interval delta_q.
% 6) Scope measured y and true position for comparison.
      

Chapter17_Lesson4.nb

(* Chapter17_Lesson4.nb *)
(* Noise Modeling in Sensors and Actuators for Dynamic Systems *)
(* Wolfram Language code: covariance propagation + stochastic simulation *)

SeedRandom[17];

(* Physical and sampling parameters *)
m = 1.2; c = 0.35; k = 4.0;
Ts = 0.01; nSteps = 6000;

(* Continuous-time matrices *)
Ac = { {0., 1.}, {-k/m, -c/m} };
Bc = { {0.}, {1./m} };
Gc = { {0.}, {1./m} };
Cmat = { {1., 0.} };

(* Euler discretization *)
A = IdentityMatrix[2] + Ts Ac;
B = Ts Bc;
G = Ts Gc;

(* Noise parameters *)
sigmaW = 0.20;
sigmaV = 0.01;
sigmaBiasRW = 0.002;
deltaQ = 0.001;
rhoA = 0.97;
sigmaASS = 0.08;

Qw = { {sigmaW^2} };
Qa = { {sigmaASS^2} };
Rv = { {sigmaV^2} };
Rq = { {deltaQ^2/12.0} };

(* Symbolic covariance recursion (discrete Lyapunov step) *)
P = ConstantArray[0., {2, 2}];

x = {0., 0.};
bias = 0.;
aNoise = 0.;

y = ConstantArray[0., nSteps];
yTrue = ConstantArray[0., nSteps];
syy = ConstantArray[0., nSteps];

Do[
  t = (kIdx - 1) Ts;
  uCmd = 0.8 Sin[2 Pi 0.7 t];

  aNoise = rhoA aNoise + Sqrt[1 - rhoA^2] sigmaASS RandomVariate[NormalDistribution[0, 1]];
  w = sigmaW RandomVariate[NormalDistribution[0, 1]];

  x = Flatten[A.x + B (uCmd + aNoise) + G w];

  bias = bias + sigmaBiasRW Sqrt[Ts] RandomVariate[NormalDistribution[0, 1]];
  v = sigmaV RandomVariate[NormalDistribution[0, 1]];

  yAnalog = First[Cmat.x] + bias + v;
  q = deltaQ Round[yAnalog/deltaQ] - yAnalog;
  y[[kIdx + 1]] = yAnalog + q;
  yTrue[[kIdx + 1]] = First[Cmat.x];

  P = A.P.Transpose[A] + G.Qw.Transpose[G] + B.Qa.Transpose[B];
  syy[[kIdx + 1]] = First[First[Cmat.P.Transpose[Cmat] + Rv + Rq]];
  ,
  {kIdx, 1, nSteps - 1}
];

burn = 1000;
ySub = y[[burn ;;]];
yTrueSub = yTrue[[burn ;;]];
syySub = syy[[burn ;;]];

empVarY = Variance[ySub];
empVarYTrue = Variance[yTrueSub];
predVarY = Mean[syySub];

Print["=== Noise Modeling Demo (Mathematica) ==="];
Print["Empirical var(y_true): ", NumberForm[empVarYTrue, {10, 6}]];
Print["Empirical var(y)    : ", NumberForm[empVarY, {10, 6}]];
Print["Predicted mean Syy  : ", NumberForm[predVarY, {10, 6}], " (bias RW excluded)"];
Print["Quantization theory : ", NumberForm[deltaQ^2/12.0, {10, 6}]];

(* Optional exact steady-state covariance for the linear plant (without bias RW) *)
p11 =.; p12 =.; p22 =.;
Psym = { {p11, p12}, {p12, p22} };
Qeff = G.Qw.Transpose[G] + B.Qa.Transpose[B];
sol = Solve[Flatten[Psym == A.Psym.Transpose[A] + Qeff], {p11, p12, p22}];
Print["Steady-state plant covariance (linearized recursion): "];
Print[Psym /. sol[[1]]];
      

Problem 1 (State Covariance with Actuator and Process Noise): Consider the discrete model \( \mathbf{x}_{k+1}=\mathbf{A}_d\mathbf{x}_k+\mathbf{B}_d\mathbf{u}_{c,k}+\mathbf{B}_d\mathbf{n}_{a,k}+\mathbf{G}_d\mathbf{w}_k \) with zero-mean independent \( \mathbf{n}_{a,k} \) and \( \mathbf{w}_k \). Derive the covariance recursion for \( \mathbf{P}_k \).

Solution: Let \( \mathbf{m}_k=E[\mathbf{x}_k] \) and define \( \tilde{\mathbf{x} }_k=\mathbf{x}_k-\mathbf{m}_k \). Then

\[ \tilde{\mathbf{x} }_{k+1}=\mathbf{A}_d\tilde{\mathbf{x} }_k+\mathbf{B}_d\tilde{\mathbf{n} }_{a,k}+\mathbf{G}_d\tilde{\mathbf{w} }_k \]

so

\[ \mathbf{P}_{k+1} = E[\tilde{\mathbf{x} }_{k+1}\tilde{\mathbf{x} }_{k+1}^{\top}] \]

Expanding and using independence (zero cross-covariance terms) gives

\[ \mathbf{P}_{k+1} = \mathbf{A}_d\mathbf{P}_k\mathbf{A}_d^{\top} + \mathbf{B}_d\mathbf{Q}_a\mathbf{B}_d^{\top} + \mathbf{G}_d\mathbf{Q}_w\mathbf{G}_d^{\top} \]

This is the discrete covariance propagation formula used in stochastic simulation and Kalman filtering.

Problem 2 (Quantization Variance): Assume quantization error \( q \) is uniformly distributed over \( [-\Delta/2,\Delta/2] \). Prove that \( \operatorname{Var}(q)=\Delta^2/12 \).

Solution: Since the density is symmetric, \( E[q]=0 \). Therefore

\[ \operatorname{Var}(q)=E[q^2] = \int_{-\Delta/2}^{\Delta/2} q^2 \frac{1}{\Delta}\,dq = \frac{1}{\Delta}\left[\frac{q^3}{3}\right]_{-\Delta/2}^{\Delta/2} = \frac{\Delta^2}{12} \]

The result is exact under the uniform model.

Problem 3 (AR(1) Actuator Noise Parameterization): Let actuator noise obey \( n_{k+1}=\rho n_k+\varepsilon_k \) with \( |\rho|<1 \) and \( \varepsilon_k\sim \mathcal{N}(0,\sigma_{\varepsilon}^2) \). If the desired stationary variance is \( \sigma_n^2 \), derive \( \sigma_{\varepsilon}^2 \).

Solution: In stationarity, \( \operatorname{Var}(n_{k+1})=\operatorname{Var}(n_k)=\sigma_n^2 \). Since \( n_k \) and \( \varepsilon_k \) are independent,

\[ \sigma_n^2 = \rho^2 \sigma_n^2 + \sigma_{\varepsilon}^2 \]

hence

\[ \sigma_{\varepsilon}^2=(1-\rho^2)\sigma_n^2 \]

This is the variance relation used in the code examples.

Problem 4 (Bias Drift Covariance Growth): A scalar sensor bias follows \( b_{k+1}=b_k+\eta_k \) with \( E[\eta_k]=0 \), \( \operatorname{Var}(\eta_k)=Q_b \), and independent increments. Show that \( \operatorname{Var}(b_k)=\operatorname{Var}(b_0)+kQ_b \).

Solution: By recursion,

\[ b_k = b_0 + \sum_{i=0}^{k-1}\eta_i \]

Using independence of the increments,

\[ \operatorname{Var}(b_k) = \operatorname{Var}(b_0) + \sum_{i=0}^{k-1}\operatorname{Var}(\eta_i) = \operatorname{Var}(b_0)+kQ_b \]

Therefore random-walk bias is nonstationary and causes long-term variance growth in the measurement channel.

Problem 5 (Measurement Covariance Assembly): For \( y_{m,k}=Cx_k+b_k+v_k+q_k \), assume all terms are mutually independent and zero-mean except possibly \( x_k \). Show the measurement variance formula for the scalar case.

Solution: Let \( P_k=\operatorname{Var}(x_k) \). Then

\[ \operatorname{Var}(y_{m,k}) = C^2 P_k + \operatorname{Var}(b_k) + \operatorname{Var}(v_k) + \operatorname{Var}(q_k) \]

because the covariance of a sum is the sum of covariances when cross-covariance terms vanish. In vector form, this becomes \( \mathbf{S}_{y,k}=\mathbf{C}\mathbf{P}_k\mathbf{C}^{\top}+\mathbf{R}_{b,k}+\mathbf{R}_v+\mathbf{R}_q \).

We modeled sensor and actuator noise as structured stochastic processes connected to the plant through physically meaningful channels. The key lessons are: (1) actuator noise enters through the input channel and propagates through state dynamics, (2) sensor noise is not a single term but a sum of white noise, bias drift, quantization, and timing effects, (3) covariance propagation depends on where noise enters the model, and (4) PSD-based modeling and shaping filters provide a bridge between measured spectra and state-space simulation models.

These models are the foundation for the next lesson on Monte Carlo simulation of stochastic dynamic systems, where repeated random trials will be used to estimate distributions of time responses and performance metrics.