Propagating Non-Gaussian Correlated Type B Uncertainties through Dynamic Non-Linear State Space Matrices

Flaw

A calibration certificate reporting a temperature transmitter accuracy of ±0.05 °C under a rectangular distribution assigns an identical probability density to zero error and to the extreme boundary. That rectangular assumption, combined with an unstated coverage factor on a second-tier bench multimeter, enters a state estimation filter as an equivalent Gaussian covariance matrix. The filter proceeds to minimize mean squared error under linear assumptions.

The state trajectory calculated across a non-linear process promptly diverges from the physical plant. The error does not originate in the state observer mathematics. It enters directly through the mismatch between the non-Gaussian bounding mechanics of Type B uncertainties and the linear covariance algebra forced upon them.

Static metrology isolates uncertainty evaluation to single scalar outputs governed by the Guide to the Expression of Uncertainty in Measurement (JCGM 100:2008). Dynamic non-linear state space representations require continuous tracking of parameter perturbations through non-linear matrix operators. When input uncertainties originate from calibration certificates, manufacturer limits, digitizer quantization, and physical hysteresis, their underlying distributions are inherently non-Gaussian.

They exhibit hard boundaries, asymmetric limits, or multi-modal forms. Forcing these distributions into equivalent Gaussian standard deviations discards higher-order moments. In non-linear dynamic transformations, skewness and kurtosis dictate the location of mean trajectories and the migration of confidence intervals over time.

Engineers face a structural contradiction when handling dynamic measurement equations. Standard sensor procurement relies on nominal tolerance bands. Yet when sensors feed dynamic state matrices, unmodeled spatial and temporal cross-correlations turn nominal tolerance bounds into runaway state covariances.

JCGM 100:2008 clause 4.3 dictates that systematic Type B components with hard physical boundaries require exact probability density functions rather than premature conversion to symmetric variance.

Correlations among Type B uncertainties exacerbate the distortion. Two transducers calibrated against the same deadweight tester share a systematic calibration bias. Thermal gradients across a common aluminum mounting block couple the resistance changes of adjacent sensor channels.

When dynamic state matrices contain products of these coupled state variables, cross-correlation terms amplify system non-linearities. Ignoring cross-correlations understates output variance. Calculating cross-correlations with simple linear Pearson coefficients when the underlying variables are bounded, asymmetric, and non-Gaussian produces mathematically invalid joint probability densities.

The operational cost lands during plant acceptance. An autonomous process or active flight-control system operating near its operational limits relies on dynamic uncertainty bounds to maintain safety margins. If the propagation routine understates state bounds due to forced Gaussian simplifications, safety interlocks trip unnecessarily or the system breaches allowable boundaries.

Dynamic state space models require a continuous, mathematically rigorous chain linking Type B certificates to time-varying probability density functions.

The consequence of treating non-Gaussian input limits as Gaussian covariances in non-linear dynamic models is the systematic underestimation of tail probabilities, invalidating safety margins before the hardware logs an hour of operational runtime.

Mapping

The dynamic non-linear state space formulation treats both the state evolution and the measurement sequence as non-linear vector functions subjected to parameter perturbations and noise vectors:

x(k+1) = f(x(k), u(k), θ, w(k))

y(k) = h(x(k), θ, v(k))

The state vector x(k) carries the dynamic physical variables. The vector u(k) contains known deterministic controls. The parameter vector θ encapsulates physical dimensions, piezoresistive gauge factors, calibration gains, and thermal coefficients.

The process noise w(k) and measurement noise v(k) contain high-frequency stochastic variations, historically categorized as Type A. The parameter vector θ houses the correlated Type B uncertainties. These parameters do not vary randomly from millisecond to millisecond; they remain fixed or drift slowly over time, constrained by finite intervals defined on calibration reports.

Mathematical Structure of Dynamic State Perturbations

Let the vector θ possess a joint probability density function p(θ). This distribution reflects prior knowledge: rectangular bounds from digital quantization, U-shaped distributions from thermal cycling and hysteresis, and trapezoidal distributions derived from certificate tolerances compounded with environmental drift. The state transition function f(·) is non-linear with respect to both x(k) and θ.

Expanding f(·) in a Taylor series around the estimated state and nominal parameters reveals how the higher moments enter:

f(x, θ) = f(x_nom, θ_nom) + J_x Δx + J_θ Δθ + 0.5 + R_3

Here, J_x and J_θ denote the Jacobian matrices, while H represents the respective Hessian tensors. In linear Gaussian propagation, the second-order terms are discarded, and the covariance propagates via J P J^T. In non-linear systems driven by bounded Type B variables, the expectation of the second-order term 0.5 E is non-zero even when the perturbation Δθ is zero-mean. Non-zero third and fourth moments in p(θ) ~ characteristic of rectangular, triangular, and asymmetric beta distributions ~ directly bias the mean trajectory of the state x(k+1) and distort its variance.

Copula Formulations for Correlated Epistemic Parameters

Correlations among Type B uncertainties cannot be modeled by covariance matrices when marginal distributions depart from the Gaussian profile. A standard correlation matrix describes linear association; it fails to constrain joint probabilities to physical boundaries. If parameter θ_1 is rectangularly distributed across and parameter θ_2 is rectangularly distributed across , assigning a linear correlation coefficient of 0.8 through a multivariate normal distribution generates non-zero probability densities in the corners (-1, 1) and (1, -1), violating physical boundary constraints.

Sklar’s theorem resolves this operational constraint. The joint cumulative distribution function F(θ_1, θ_n) decomposes into marginal cumulative distribution functions F_i(θ_i) and a copula function C ~

F(θ_1, θ_2, θ_n) = C(F_1(θ_1), F_2(θ_2), F_n(θ_n))

To couple non-Gaussian Type B uncertainties without distorting individual calibration bounds, metrology systems implement specific copula structures. The Gaussian copula maps arbitrary marginal distributions to a standard normal space via the probability integral transform:

C_Gauss(u_1, u_2; R) = Φ_R(Φ^-1(u_1), Φ^-1(u_2))

Here, Φ^-1 is the inverse standard normal cumulative distribution, and Φ_R is the joint cumulative distribution of a zero-mean normal vector with correlation matrix R. For sensor populations subject to common thermal shock or common calibration reference failure, tail dependence must be accounted for. Clayton and Gumbel copulas (members of the Archimedean family) model asymmetric clustering at the extremes of the operational window. When two silicon pressure sensors share a common supply rail, voltage sag produces correlated boundary shifts that linear correlation coefficients systematically miss.

The standard Pearson correlation matrix fails whenever marginal distributions depart from the Gaussian bell curve.

Transformation

Propagating the joint density p(θ) through the non-linear state update requires transformation techniques that preserve moments without imposing prohibitive computation costs. Three primary methodologies dominate applied dynamic metrology: Polynomial Chaos Expansion (PCE), unscented deterministic transformations, and Markov Chain Monte Carlo (MCMC) propagation governed by JCGM 101:2008.

Generalized Polynomial Chaos Expansion

Generalized Polynomial Chaos Expansion (gPCE) models the stochastic state vector x(k) as an explicit series of orthogonal polynomial bases evaluated over standard random variables. The dynamic state variable expands as:

x(k, ξ) = Σ_{i=0}^{P} c_i(k) Ψ_i(ξ)

The vector ξ represents independent standard random variables, Ψ_i(ξ) are multidimensional orthogonal polynomials chosen specifically for the distribution families of ξ, and c_i(k) are deterministic coefficient vectors propagated through time. Under the Askey scheme of orthogonal polynomials, uniform distributions map to Legendre polynomials, Gaussian distributions map to Hermite polynomials, and bounded beta distributions map to Jacobi polynomials.

Correlated non-Gaussian Type B variables require a pre-processing transformation before entering the polynomial basis. The Nataf transformation or the Rosenblatt transformation maps the physical correlated vector θ to the independent standard vector ξ:

• Rosenblatt Transformation maps conditional cumulative distribution functions sequentially into uniform random spaces, decoupling non-linear dependence structures across asymmetric calibration limits.
• Nataf Transformation converts marginal distributions into standard Gaussian equivalents while iteratively solving the non-linear integral equation relating the physical correlation coefficient ρ_ij to the transformed correlation coefficient ρ_0,ij.
• Gram-Schmidt Orthogonalization generates numerical polynomial bases directly when empirical calibration histograms refuse to conform to standard analytical families.

Once the orthogonal basis is established, the dynamic state space equations transform into a coupled deterministic system for the coefficients c_i(k). Applying Galerkin projection projects the residual error of the non-linear dynamic model onto the polynomial basis, forcing the error to zero:

⟨ , Ψ_j(ξ) ⟩ = 0, for all j = 0, P

The resulting deterministic state equations propagate the entire probability distribution through time. The mean state vector corresponds exactly to the zeroth coefficient c_0(k). The state covariance matrix derives directly from the sum of squares of the higher-order coefficients multiplied by the inner product of the basis polynomials:

Cov = Σ_{i=1}^{P} c_i(k) c_i^T(k) ⟨Ψ_i^2(ξ)⟩

PCE eliminates the need for repeated time-series numerical integrations. Once the coefficient matrices are solved, confidence intervals, skewness, and probability of boundary exceedance compute algebraically from the vectors c_i(k).

ISO/IEC Guide 98-3 Supplement 1 requires that dynamic numerical transformations validate state probability intervals against an absolute coverage criterion of 95% without assuming distribution symmetry.

Deterministic Sampling and the Transformed Unscented Split

Unscented transformations approximate non-linear state propagation using deterministic sigma points. Standard unscented Kalman filtering selects sigma points assuming a symmetric, Gaussian prior. The unscented weights yield accurate mean and covariance propagation up to the second order, but fail when Type B inputs introduce significant kurtosis.

For non-Gaussian variables with known skewness and kurtosis, the higher-order unscented transformation incorporates third- and fourth-moment constraints into the sigma point algebraic equations:

Σ w_i = 1

Σ w_i (χ_i – μ) (χ_i – μ)^T = P_θ

Σ w_i (χ_i – μ)^3 = Skewness

Σ w_i (χ_i – μ)^4 = Kurtosis

When the number of state variables and uncertain parameters exceeds twelve, solving for deterministic sigma points satisfying higher-moment constraints becomes numerically ill-conditioned. Under these conditions, the engineer switches to the transformed unscented split: the dynamic state space matrices are split into purely Gaussian stochastic states propagated via standard unscented filters, and bounded Type B epistemic parameters propagated via a low-order PCE surrogate model.

How does dimensional scaling impact polynomial chaos execution?

As the count of uncertain parameters n increases, the total number of basis terms P in a full tensor product expansion escalates according to the combinatorial expression (n + d)! / (n! d!), where d is the polynomial degree. At n = 8 and d = 4, the basis requires 495 evaluation terms.

At n = 15 and d = 4, the basis balloons to 3,876 terms. To bypass this combinatorial explosion, dynamic non-linear state matrices require hyperbolic truncation schemes or sparse Smolyak quadrature grids. These sparse grids retain cross-coupling terms up to a specified interaction order while truncating higher-order multi-variable combinations, keeping the evaluation count within real-time compute limits.

Bench

Consider a practical verification scenario: a non-linear flow reactor monitored by two high-temperature piezoresistive pressure sensors and an RTD assembly. The fluid density ρ varies non-linearly with temperature T and pressure P according to a virial equation of state. The dynamic state space model describes the mass flow m_dot(k) and fluid density ρ(k) across dynamic throttling conditions:

m_dot(k+1) = C_d(θ) · A(k) · √ + w(k)

ρ(k+1) = ρ_nom ·

The parameter vector θ = ^T contains correlated Type B uncertainty components derived entirely from supplier documentation and bench checks:

1. Discharge Coefficient C_d is bounded within due to machining tolerances, following an asymmetric triangular distribution with mode at 0.620.
2. Thermal Expansion Parameter α exhibits a rectangular distribution across K^-1, established by the calibration uncertainty of the thermal chamber.
3. Compressibility Parameter β carries a U-shaped distribution across Pa^-2, generated by internal mechanical hysteresis inside the sensor diaphragms.

Sensors P_1 and P_2 share a common regulated 10.0 V excitation bus. A 0.05% shift in excitation voltage shifts both pressure calibrations in parallel, introducing a linear correlation coefficient of r = 0.72 across their Type B calibration errors. The temperature sensor resistance is read by a 16-bit analog-to-digital converter introducing uniform quantization error.

The calculation reveals structural shifts in the state estimation. The mean mass flow calculated through generalized Polynomial Chaos shifts upward by 0.061 kg/s compared to the linearized model. This shift occurs because the expectation of the non-linear square root operator operating on asymmetric input spaces is strictly non-zero.

The standard deviation expands by 24.6%. The true 95% coverage interval extends from 14.071 kg/s to 15.688 kg/s. The linearized model predicts an upper limit of 15.433 kg/s.

A safety relief valve sized on the basis of the linearized model operates with a truncated uncertainty margin, exposing downstream piping to uncalculated overpressure excursions.

The skewness of the calculated state reaches 0.684. Standard Kalman filter tracking assumptions demand zero skewness. The non-linear interaction between the U-shaped hysteresis parameter β and the triangular discharge coefficient C_d generates extended tails in the output probability density function that linear algebra cannot represent.

If you fail to preserve the kurtosis of bounded sensor limits, you miscalculate the probability of exceeding physical plant thresholds by up to an order of magnitude.

Drift

Dynamic state matrices rarely operate under pristine factory conditions. The parameter vector θ changes its uncertainty profile as hardware components age in the field. When sensors endure mechanical vibration, thermal cycling, and chemical exposure, their Type B distributions undergo three physical transformations: boundary expansion, profile flattening, and correlation restructuring.

Silicon strain gauge elements demonstrate zero-point wander caused by the relaxation of residual stresses within the glass-frit die-attach layers. Over three years of operation, this mechanical relaxation shifts the zero-offset uncertainty from a narrow rectangular calibration band into an expanded trapezoidal distribution. The outer limits expand by 0.15% of span per year, while the core flat region widens as hysteresis accumulates.

When this wandering parameter is fed into a dynamic state observer, the state covariance matrix cannot remain stationary; it requires a dynamic Type B time update step.

Temporal Drift Modeling in State Matrices

To capture physical drift within the state formulation, the static parameter vector θ transforms into a dynamic Gauss-Markov or random walk epistemic state:

θ(k+1) = θ(k) + δ_drift(k)

The perturbation vector δ_drift(k) does not represent high-frequency white noise. It represents a non-Gaussian increment. For a platinum resistance thermometer (Pt100) compliant with IEC 60751, exposure to temperatures above 300 °C causes grain boundary migration and trace chemical contamination of the platinum winding.

The drift rate is strictly directional: resistance increases over time, creating a positive skewness in the Type B uncertainty distribution.

Field verification protocols must log drift as an evolving non-Gaussian parameter boundary rather than absorbing calibration offsets into broad, unverified process noise terms.

Correlation restructuring presents an even steeper hurdle. When two identical pressure sensors are deployed in a differential configuration, their initial factory calibration errors are largely independent if calibrated on different dates against different deadweight standards. After two years of continuous operation in a harsh marine environment, common-mode salt crusting, diaphragm oxidation, and ambient temperature exposure couple their drift rates.

Their joint probability density function F(θ_1, θ_2) develops a pronounced diagonal alignment. The copula linking them shifts from an independent product copula to a Clayton copula with strong lower-tail dependence.

The supplier will state that annual drift is fully covered by the standard tolerance band printed on the catalog sheet.

Audit

When procurement teams issue Requests for Quotations (RFQs) for high-reliability sensor channels, they frequently specify accuracy as an isolated line-item tolerance: Accuracy ±0.1% of Full Scale. This phrasing invites lower-tier suppliers to submit certificates derived under highly controlled laboratory conditions with coverage factor k=2, omitting the specific probability density distribution and boundary conditions. A rigorous quality audit enforces a traceable documentation hierarchy before components reach the hardware-in-the-loop dynamic integration bench.

Accredited laboratories operating under ISO/IEC 17025 must report whether stated uncertainty figures include hysteresis, dynamic response lag, and environmental coefficients. If an incoming calibration report states a ±0.2% measurement uncertainty without defining whether the figure represents a 95% Gaussian interval, a rectangular bounding limit, or a three-sigma ceiling, the data is metrologically unusable inside non-linear dynamic state matrices.

Incoming Verification and Compliance Requirements

To ensure incoming components meet the mathematical constraints of dynamic uncertainty propagation models, procurement documentation mandates the following deliverables:

• Certified Distribution Geometry requires the calibration laboratory to explicitly define the underlying probability density function (uniform, triangular, or normal) for each stated Type B uncertainty component.
• Cross-Correlation Calibration Matrix defines the shared physical references used during multi-channel sensor testing, providing empirical correlation coefficients across calibration points.
• Thermal Hysteresis Signatures provide ascending and descending temperature cycle data, enabling the construction of U-shaped distribution widths for dynamic state matrix correction.
• Long-Term Stability Bounds quantify drift per 1,000 hours of continuous operation, bounding the non-Gaussian parameter update model over the required service interval.

Quality engineers review the Scope of Accreditation for any testing facility providing calibration certificates. The accreditation scope must explicitly include the parameter ranges and the specific measurement uncertainty declared. A certificate issued by a laboratory lacking explicit accreditation for dynamic calibration or non-linear sensor characterization provides no legal or technical defense if plant telemetry drifts out of compliance.

Every tolerance claim without an attached distribution geometry forces the integration engineer to guess the tails, invalidating the dynamic state bounds before the system ships.

Expenditure

Tightening sensor specifications to minimize uncertainty propagation through non-linear matrices steepens the procurement cost curve sharply. Sourcing engineers face trade-offs between component-level hardware precision, calibration certificate rigor, and computational software complexity. Eliminating downstream state uncertainty does not always require purchasing higher-grade physical sensors; it frequently requires purchasing complete documentation of the non-Gaussian Type B uncertainty structure.

A standard industrial pressure transmitter carries a baseline procurement cost of approximately 450 USD. This device delivers an accuracy claim of ±0.25% of span, typically backed by a standard factory inspection sheet that assumes a symmetric Gaussian distribution with an unspecified coverage factor. Upgrading to a metrology-grade unit with full characterization ~ individual linearization curves, thermal hysteresis mapping, and an ISO/IEC 17025 accredited calibration certificate specifying rectangular and U-shaped uncertainty bounds ~ escalates the unit procurement price to 2,200 USD.

For an installation requiring sixteen monitored points across a non-linear process plant, this shift increases direct instrument expenditure from 7,200 USD to 35,200 USD.

The economic leverage sits in the mathematical formulation. Implementing generalized Polynomial Chaos Expansion or advanced copula transformations in state-estimation software costs engineering hours upfront during filter design. However, it allows the operational architecture to run on 850 USD characterized industrial sensors rather than 5,800 USD primary reference units.

The advanced mathematics absorbs the bounded non-Gaussian distributions of lower-cost hardware, propagating exact confidence bounds through the plant’s non-linearities without violating safety factors.

The ongoing re-calibration bill dictates long-term life-cycle expenditure. If an organization purchases precision sensors but calibrates them using low-tier field calibrators that inject uncharacterized rectangular errors, the hardware investment is nullified. The dynamic filter immediately defaults to the uncertainty floor of the field calibrator.

Sourcing teams must cost the complete calibration chain ~ including field calibrator certification, specialized laboratory benches, and the technician hours required to verify cross-sensor correlations ~ over the five-year operational lifecycle.

Procuring a tighter tolerance line without budgeting for the traceable reference standards required to re-verify it every twelve months merely defers an unbudgeted instrumentation cost into the operating ledger.