Practical Steady-State Discrete-Time Kalman Filter Design for Uncertain LTI Systems

1. Introduction

The Kalman filter, introduced by Rudolf E. Kálmán in 1960 [1], represents a fundamental tool in modern control theory and signal processing for optimal state estimation of linear dynamic systems in the presence of noise. The filter operates recursively on streams of noisy input data to produce statistically optimal estimates of the underlying system state, making it indispensable in applications ranging from navigation and guidance systems to power converter control [2].

1.1. State Estimation and the Filtering Problem

Consider a discrete-time linear time-invariant (LTI) system described by the state-space representation:

$$\left\{\begin{array}{cc}{x}_{k+1}\hfill & =A{x}_k+B{u}_k+w_k\hfill \\ y_k\hfill & =H{x}_k+v_k\hfill \end{array}\right.$$

(1)

where $x_k\in {\mathbb{R}}^n$ is the state vector, $u_k\in {\mathbb{R}}^q$ is the control input, $y_k\in {\mathbb{R}}^m$ is the measurement (output) vector, and $w_k\in {\mathbb{R}}^n$ and $v_k\in {\mathbb{R}}^m$ represent process and measurement noise, respectively.

The fundamental problem of state estimation is to reconstruct the internal state $x_k$ from the available noisy measurements $y_k$ and known inputs $u_k$. This problem is particularly challenging when:

☐

the system state is not directly measurable,

☐

measurements are corrupted by noise with unknown characteristics,

☐

the system model itself contains uncertainties,

☐

real-time estimation is required for feedback control.

The Kalman filter addresses this problem by computing the minimum variance unbiased estimate of the state, i.e. the estimate that minimizes the mean squared error $\mathbb{E}\left[\left|\right|x_k-{\widehat{x}}_k{\left|\right|}^2\right]$ (${\widehat{x}}_k$ being the state estimate) under the assumption of Gaussian noise statistics [3].

1.2. The Dual Role: Estimation and Filtering

The Kalman filter serves two complementary purposes in control and signal processing applications.

State Estimation

The filter provides optimal estimates of the system’s internal state variables that are either unmeasurable or measured with excessive noise. By fusing information from the dynamic model (predictions) and measurements (corrections), the Kalman filter produces estimates that are more accurate than either source alone. The recursive structure, consisting of prediction and update steps, enables real-time implementation with minimal computational overhead:

$$\left\{\begin{array}{c}\begin{array}{cccc}\hfill {\widehat{x}}_{k+1}^{-}& =A{\widehat{x}}_k^{+}+B{u}_k\hfill & \hfill & (\mathrm{a}\mathrm{priori}\mathrm{estimate}),\hfill \\ \hfill {\widehat{x}}_k^{+}& ={\widehat{x}}_k^{-}+K_k(y_k-H{\widehat{x}}_k^{-})\hfill & \hfill & (\mathrm{a}\mathrm{posteriori}\mathrm{estimate}).\hfill \end{array}\hfill \end{array}\right.$$

(2)

The Kalman gain $K_k$ is computed to optimally weight the relative contributions of the prediction and measurement based on their respective noise characteristics.

Variance Minimization (Filtering)

Beyond point estimation, the Kalman filter explicitly minimizes the covariance of the estimation error. By propagating the error covariance matrix $P_k$ through the prediction and update equations, the filter quantifies the uncertainty in its estimates and automatically adjusts the Kalman gain to minimize this uncertainty. This dual estimation of both the state mean and covariance distinguishes the Kalman filter from simpler estimation techniques and enables:

☐

optimal sensor fusion when multiple measurements with different noise levels are available,

☐

detection of degraded system performance through monitoring of innovation statistics,

☐

robust control design using the filtered state estimates with quantified uncertainty,

☐

adaptive adjustment of controller parameters based on estimation confidence.

The variance reduction property is particularly valuable in power converter applications, where measurement noise from analog-to-digital converters, electromagnetic interference, and quantization effects can significantly degrade control performance if not properly filtered.

1.3. Steady-State Kalman Filter

For time-invariant systems, the Kalman gain converges to a constant matrix, yielding the steady-state Kalman filter. This asymptotic solution, obtained by solving the discrete algebraic Riccati equation (DARE), offers several practical advantages:

☐

reduced computational burden (no recursive covariance propagation),

☐

simplified implementation and parameter tuning,

☐

guaranteed stability under standard detectability and stabilizability assumptions (subject to numerical regularization) [4],

☐

design-time computation of filter parameters,

The steady-state formulation is particularly attractive for embedded control systems where computational resources are limited and deterministic real-time performance is critical.

1.4. Model Uncertainty and Robustness

Classical Kalman filter theory assumes perfect knowledge of the system matrices A, B, and H, as well as Gaussian white noise processes w and v with known covariances. In practice, these assumptions are rarely satisfied:

☐

component tolerances introduce uncertainty in physical parameters (resistance, inductance, capacitance),

☐

operating conditions vary (temperature, ageing, non-linearities),

☐

modelling approximations neglect high-frequency dynamics,

☐

external disturbances exhibit non-Gaussian characteristics.

Non-linearities are expressly addressed by Extended Kalman Filters (EKFs) and not investigated any further in this work; it is however worth mentioning that the standard EKF employs a similar linearization strategy — first-order Taylor expansion combined with second-order moment propagation — which is precisely the approach adopted here. When model uncertainties are significant, the nominal Kalman filter may provide far from optimal or even unstable estimates. This motivates the development of robust Kalman filtering techniques that explicitly account for parametric uncertainty. By incorporating the covariance of model parameters into the filter design, one can synthesize filters that maintain near-optimal performance (such as being the best linear estimator) despite modelling errors.

The literature on robust Kalman filtering has evolved along several distinct directions. The seminal work of De Koning [5] established the theoretical foundation by deriving equivalent noise covariances for systems with stochastic parameters variations modelled as sequences of i.i.d. (independent and identically distributed) zero-mean random variables, demonstrating that model uncertainties can be systematically incorporated into the classical Kalman framework and deriving augmented process and measurement noise statistics.

Building on these foundations, several researchers have developed alternative robust filtering methodologies. Xie et al. [6] proposed a min-max approach based on game theory, designing filters that minimize the worst-case estimation error over all admissible parameter uncertainties. Theodor and Shaked [7] extended this framework to derive robust minimum-variance filters with guaranteed performance bounds. Zhu et al. [8] provided a comprehensive design methodology combining linear matrix inequality (LMI) techniques with parameter-dependent Lyapunov functions. Luo and Bosch [9] analysed the performance robustness of Kalman filters, establishing relationships between parameter uncertainty levels and estimation error degradation.

More recent advances have addressed structured uncertainties as polytopic models: Yu et al. [10] developed robust Kalman filters for systems with convex polytopic uncertainties, where the uncertain system matrices lie within a convex hull of known vertices. Xie et al. [11] proposed improved $H_2$ and $H_{\infty }$ filtering techniques that, although not strictly Kalman filters, offer relevant insights into handling model uncertainties.

While all the approaches mentioned above provide valuable theoretical insights and worst-case performance guarantees, they do not explicitly consider the control input $u_k$ in the uncertainty propagation.

Rocha and Terra [12] presented the most comprehensive treatment to date (to the best of the author’s knowledge), deriving robust filters that handle multiplicative noise, time-correlated uncertainties, and cross-coupled parameter variations through a unified recursive framework; it furthermore includes control input $u_k$. Such an approach is, however, considered difficult to apply operationally for practitioners that are not necessarily experts in robust control as it is based on regularization techniques that need non-trivial tuning.

The present work attempts to extend De Koning’s approach [5] to explicitly account for control inputs; it however departs from it by modelling the parametric uncertainties more rigorously although only approximately. This extension is particularly important for feedback control applications where the system operates under varying inputs, and where the interaction between parametric uncertainty and control effort significantly affects the equivalent noise characteristics. By deriving closed-form expressions for the input-dependent covariance terms, it enables systematic tuning of the steady-state filter robustness through appropriate selection of the nominal operating point. The focus of the present work is mostly on the practicality of the proposed robust Kalman filter.

1.5. Scope and Contribution

This paper presents a steady-state discrete-time Kalman filter design methodology for LTI systems with parametric uncertainties similar to the MVFOSM approach [13,14]. The key contributions are:

☐

derivation of equivalent process and measurement noise covariances that incorporate both stochastic noise and deterministic model uncertainty, with explicit treatment of control input effects,

☐

closed-form expressions for the (approximated) state-dependant covariance matrices Q, R, and M as functions of parameter uncertainty and control input,

☐

extension to time-varying implementations for scenarios where reduced conservatism justifies increased computational cost.

The methodology enables systematic design of robust Kalman filters that explicitly trade off nominal performance and robustness to parameters’ variations. By selecting an appropriate nominal input value (e.g., maximum expected input or RMS value for periodic signals), designers can tune the conservatism of the steady-state solution to match application requirements. This is particularly valuable for power converter control, where component tolerances can significantly affect dynamic behaviour, and where computational constraints favour steady-state over time-varying implementations.

The remainder of this paper is organized as follows: Section 2 develops the system model with parametric uncertainty and establishes the discretization procedure. Section 3 derives the steady-state Kalman filter equations together with detailed, although approximated, calculations of the state-dependent covariance matrices. Section 4 briefly outlines the extension to time-varying implementations. Conclusions and directions for future work are presented in Section 5.

2. Modeling and Discretization

The following LTI continuous-time model is considered:

$$\left\{\begin{array}{c}\dot{x}=\left(A_c+\delta A_c\right)x+\left(B_c+\delta B_c\right)\left(u+\delta u\right)+w^{*},\hfill \\ y=\left(H+\delta H\right)x+v^{*},\hfill \end{array}\right.$$

(3)

where:

$$\left\{\begin{array}{cc}{A}_c\hfill & =A_c\left(\overline{\theta }\right),\hfill \\ B_c\hfill & =B_c\left(\overline{\theta }\right),\hfill \\ H\hfill & =H\left(\overline{\theta }\right),\hfill \end{array}\right.$$

are deterministic matrices and:

$$\left\{\begin{array}{cc}\delta A_c\hfill & =\delta A_c\left(\delta \theta \right),\hfill \\ \delta B_c\hfill & =\delta B_c\left(\delta \theta \right),\hfill \\ \delta H\hfill & =\delta H\left(\delta \theta \right),\hfill \end{array}\right.$$

are random matrices due to the uncertainty of the parameters grouped in the vector $\theta$, which can be expressed as:

$$\theta =\overline{\theta }+\delta \theta .$$

(4)

In eq. (4), $\overline{\theta }$ groups the nominal values of the parameters, and $\delta \theta$ represents their uncertainty such that $\mathbb{E}\left[\delta \theta \right]=0$. Furthermore, it is assumed that the input u is affected by some noise $\delta u$ (which could be quantization noise or coming from un-modelled dynamics), and the state is affected by noise $w^{*}$, while the output is corrupted by measurement noise $v^{*}$. It is also assumed that $\delta \theta$ is statistically independent from $\delta u$, $w^{*}$, and $v^{*}$ and $x_0$, which is a natural assumption. In terms of dimensions: the state $x\in {\mathbb{R}}^n$, the output $y\in {\mathbb{R}}^m$, u and $\delta u\in {\mathbb{R}}^q$, $A_c$ and $\delta A_c\in {\mathbb{R}}^{n\times n}$, $B_c$ and $\delta B_c\in {\mathbb{R}}^{n\times q}$, H and $\delta H\in {\mathbb{R}}^{m\times n}$, $w^{*}\in {\mathbb{R}}^n$ and $v^{*}\in {\mathbb{R}}^m$.

For practical use, this model needs to be discretized (with sampling period T) as follows:

$$\left\{\begin{array}{cc}{x}_{k+1}\hfill & =\left(A+\delta A\right)x_k+\left(B+\delta B\right)\left(u_k+\delta u_k\right)+w_k^{*},\hfill \\ y_k\hfill & =\left(H+\delta H\right)x_k+v_k^{*},\hfill \end{array}\right.$$

(5)

where:

$$\left\{\begin{array}{c}A=e^{A_{c}T},\hfill \\ B={\int }_0^T{e}^{A_c\tau }d\tau B_c=A_c^{-1}\left(e^{A_{c}T}-I_n\right)B_c=A_c^{-1}\left(A-I_n\right)B_c,\hfill \end{array}\right.$$

(6)

and where $I_n$ is the $n\times n$ identity matrix and $A_c$ is assumed invertible. It is convenient to express the last relationship in eq. (6) as follows:

$$B=A_c^{-1}\left(A-I_n\right)B_c=G{B}_c$$

(7)

It can be assumed (standard Kalman filters) that $w_k^{*}$ and $v_k^{*}$ are zero-mean uncorrelated white Gaussian noises. The LTI state equations can be therefore rewritten as:

$$\left\{\begin{array}{cc}{x}_{k+1}\hfill & =A{x}_k+B{u}_k+\left(\delta A{x}_k+\delta B{u}_k+B\delta u_k+\delta B\delta u_k\right)+w_k^{*},\hfill \\ y_k\hfill & =H{x}_k+\left(\delta H{x}_k+v_k^{*}\right),\hfill \end{array}\right.$$

(8)

or equivalently as:

$$\left\{\begin{array}{cc}{x}_{k+1}\hfill & =A{x}_k+B{u}_k+w_k,\hfill \\ y_k\hfill & =H{x}_k+v_k,\hfill \end{array}\right.$$

(9)

where:

$$\left\{\begin{array}{cc}{w}_k\hfill & =\delta A{x}_k+\delta B{u}_k+B\delta u_k+\delta B\delta u_k+w_k^{*},\hfill \\ v_k\hfill & =\delta H{x}_k+v_k^{*}.\hfill \end{array}\right.$$

(10)

Continuous-time process noise $w^{*}\left(t\right)$ and measurement noise $v^{*}\left(t\right)$ are assumed zero-mean, white and Gaussian. In particular $\mathbb{Cov}\left(w^{*}\left(t\right)\right)=\mathbb{E}\left[w^{*}\left(t\right)w^{*T}(t-\tau )\right]=Q_{w^{*}}^c\delta \left(\tau \right)$; the discretized process noise $w_k^{*}$ has a covariance [2]:

$$Q_{w^{*}}=\mathbb{Cov}\left(w_k^{*}\right)={\int }_0^T{e}^{A_c\left(T-\tau \right)}{Q}_{w^{*}}^c{e}^{A_c^T\left(T-\tau \right)}d\tau .$$

(11)

Equation (11) can be used if a knowledge of the covariance of the continuous time process noise is available to numerically calculate the covariance of the discrete time process noise. In particular one might assume $Q_{w^{*}}^c$ to be diagonal; as shown by eq. (11) this is not the case in discrete time unless T is extremely small, in such a case it can be seen that $Q_{w^{*}}\approx Q_{w^{*}}^{c}T$ so diagonality is preserved.

Furthermore [2]:

$$R_{v^{*}}={\displaystyle \frac{R_{v^{*}}^c}{T}},$$

(12)

in this case, diagonality of the measurement noise covariance is preserved irrespective of the value of T.

Equation (10) defines equivalent noises similar to those presented in [5] (which, however, does not discuss the control input $u_k$) where the equivalent noises are state dependent.

However, for the case analysed here, the equivalent noises are not zero-mean:

$$\left\{\begin{array}{cc}\mathbb{E}\left[w_k\right]\hfill & =z_k,\hfill \\ \mathbb{E}\left[v_k\right]\hfill & =t_k.\hfill \end{array}\right.$$

(13)

For the second moments:

$$\left\{\begin{array}{cc}\mathbb{E}\left[w_k{w}_k^T\right]\hfill & =P_{w{w}_k},\hfill \\ \mathbb{E}\left[v_k{v}_k^T\right]\hfill & =P_{v{v}_k},\hfill \\ \mathbb{E}\left[w_k{v}_k^T\right]\hfill & =P_{w{v}_k},\hfill \end{array}\right.$$

(14)

and covariances:

$$\left\{\begin{array}{cc}{Q}_k\hfill & =P_{w{w}_k}-z_k{z}_k^T,\hfill \\ R_k\hfill & =P_{v{v}_k}-t_k{t}_k^T,\hfill \\ M_k\hfill & =P_{w{v}_k}-z_k{t}_k^T.\hfill \end{array}\right.$$

(15)

The expected values of the equivalent noises ($z_k$ and $t_k$ respectively) will be considered as additional states. In terms of dimensions: $Q_k\in {\mathbb{R}}^{n\times n}$ , $R_k\in {\mathbb{R}}^{m\times m}$ , $M_k\in {\mathbb{R}}^{n\times m}$.

2.1. First-Order Modelling of the Uncertainty

At the first order, the continuous-time state-space uncertainties can be written as follows:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \delta A_c& \approx J_{A_c}\delta \theta ,\hfill \\ \hfill \delta B_c& \approx J_{B_c}\delta \theta ,\hfill \\ \hfill \delta H& \approx J_H\delta \theta ,\hfill \end{array}\hfill \end{array}\right.$$

(16)

where the $J_{(.)}$ are the Jacobians (used here as a shorthand for Jacobian matrices) of the state-space matrices with respect to the parameters $\theta$. Actually, the Jacobian is defined for (column) vectors and needs to be generalized for matrices by stacking on top of each other the Jacobians of each column; this is also known as vectorization. In other words, the following should have been written as:

$$\left\{\begin{array}{c}\begin{array}{cccc}\hfill \delta A_c& \approx \mathrm{unvec}\left[{\displaystyle \frac{\partial }{\partial \theta }}\mathrm{vec}\left(A_c\right)\delta \theta \right]\hfill & \hfill & =\mathrm{unvec}\left(J_{A_c}\delta \theta \right),\hfill \\ \hfill \delta H& \approx \mathrm{unvec}\left[{\displaystyle \frac{\partial }{\partial \theta }}\mathrm{vec}\left(H\right)\delta \theta \right]\hfill & \hfill & =\mathrm{unvec}\left(J_H\delta \theta \right),\hfill \end{array}\hfill \end{array}\right.$$

(17)

however the simplified, but abused, notation in eq. (16) is preferred throughout this work. The (partial) derivatives $\partial \theta$ should be considered as evaluated at $\theta =\overline{\theta }$, but here a simplified notation has also been preferred. In terms of dimensions: $\theta \in {\mathbb{R}}^p$, $J_{A_c}\in {\mathbb{R}}^{n^2\times p}$, $J_{B_c}\in {\mathbb{R}}^{nq\times p}$, $J_H\in {\mathbb{R}}^{nm\times p}$.

These uncertainties need to be converted into discrete time, except for $\delta H$, which is identical to the continuous time counterpart.

This is done in eq. (18) where it is considered that $B=G\left(A_c\right)B_c$ is a function of both $B_c$ and $A_c$, i.e. $B=g(B_c,A_c)$:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \delta A& \approx {\displaystyle \frac{\partial \mathrm{vec}\left(A\right)}{\partial \mathrm{vec}{\left(A_c\right)}^T}}\delta A_c=J_A\delta A_c\approx J_A{J}_{A_c}\delta \theta ,\hfill \\ \hfill \delta B& \approx {\displaystyle \frac{\partial \mathrm{vec}\left(g\right)}{\partial \mathrm{vec}{\left(B_c\right)}^T}}\delta B_c+{\displaystyle \frac{\partial \mathrm{vec}\left(g\right)}{\partial \mathrm{vec}{\left(A_c\right)}^T}}\delta A_c=\left(I_q\otimes G\right)\delta B_c+J_B\delta A_c\approx \left(I_q\otimes G\right)J_{B_c}\delta \theta +J_B{J}_{A_c}\delta \theta .\hfill \end{array}\hfill \end{array}\right.$$

(18)

In eq. (18), the symbol ⊗ refers to the Kronecker (or tensor) product used in the vectorization that allowed the different contributions to be expressed as functions of the covariances of the state-space matrices. Furthermore in eqs. (17) and (18) the notation for the expression of the Jacobians strictly follows [15] although it is abused for the expression of the differentials which would require an $\mathrm{unvec}$ operation. This will be clarified only when a disambiguation is strictly needed.

In the following, the definition in eq. (19) is used:

$$V=\left(I_q\otimes G\right)J_{B_c}+J_B{J}_{A_c},$$

(19)

hence it can finally be written, with abused notation (i.e. implicitly assuming vectorization and de-vectorization):

$$\delta B\approx V\delta \theta .$$

(20)

The other newly introduced Jacobians are the following:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill J_A& ={\displaystyle \frac{\partial A}{\partial A_c}}={\displaystyle \frac{\partial \mathrm{vec}\left(A\right)}{\partial \mathrm{vec}{\left(A_c\right)}^T}}={\int }_0^T{e}^{A_c^T\tau }\otimes e^{A_c\left(T-\tau \right)}d\tau ,\hfill \\ \hfill J_B& ={\displaystyle \frac{\partial B}{\partial A_c}}={\displaystyle \frac{\partial \mathrm{vec}\left(B\right)}{\partial \mathrm{vec}{\left(A_c\right)}^T}}=\left(B_c^T\otimes A_c^{-1}\right)J_A-\left(B^T\otimes A_c^{-1}\right).\hfill \end{array}\hfill \end{array}\right.$$

(21)

The derivation of the expression of $J_B$ is detailed in the appendix. In terms of dimensions: $V\in {\mathbb{R}}^{nq\times p}$, $J_A\in {\mathbb{R}}^{n^2\times n^2}$, $J_B\in {\mathbb{R}}^{nq\times n^2}$.

2.1.1. Type B Uncertainty Characterization of Component Tolerances

The robust filter design methodology presented in this work requires knowledge of the parameter uncertainty covariance matrix ${\Sigma }_{\delta \theta }$. In practice, the filter designers have rarely access to statistical population data for individual components. Instead, parameter uncertainties are typically available only as manufacturer-specified tolerances. This section discusses how to construct ${\Sigma }_{\delta \theta }$ from tolerance specifications following the Guide to the Expression of Uncertainty in Measurement (GUM) [16] framework.

The GUM distinguishes between two fundamental approaches to uncertainty quantification:

☐

Type A evaluation: Based on statistical analysis of series of observations. Requires measuring a representative sample of components and computing sample statistics (mean, variance, correlations).

☐

Type B evaluation: Based on means other than statistical analysis, including manufacturer specifications, calibration certificates, handbooks, experience, or scientific judgment. This is the approach typically available to filter designers and is therefore adopted in this work.

Electrical component tolerances are usually specified by manufacturers as symmetric bounds around nominal values. Common examples include:

☐

Resistors: $R=10\mathrm{k}\Omega \pm 1\%$ (meaning $R\in [9.9,10.1]\mathrm{k}\Omega$)

☐

Capacitors: $C=100\mu \mathrm{F}\pm 10\%$ (meaning $C\in [90,110]\mu \mathrm{F}$)

☐

Inductors: $L=1\mathrm{mH}\pm 5\%$ (meaning $L\in [0.95,1.05]\mathrm{mH}$)

These specifications define bounded intervals $[{\overline{\theta }}_i-\Delta {\theta }_i,{\overline{\theta }}_i+\Delta {\theta }_i]$ but provide no information about the probability distribution within these bounds.

When only tolerance bounds are available, the GUM methodology recommends using a uniform distribution:

$$\delta {\theta }_i\sim \mathcal{U}(-\Delta {\theta }_i,\Delta {\theta }_i).$$

(22)

The corresponding variance is:

$${\sigma }_{\delta {\theta }_i}^2={\displaystyle \frac{\Delta {\theta }_i^2}{3}}.$$

(23)

For a system with p uncertain parameters $\theta ={[{\theta }_1,{\theta }_2,\dots ,{\theta }_p]}^T$, the covariance matrix is constructed as:

$${\Sigma }_{\delta \theta }=\left[\begin{array}{cccc}{\sigma }_{\delta {\theta }_1}^2& {\sigma }_{\delta {\theta }_1,\delta {\theta }_2}& \dots & {\sigma }_{\delta {\theta }_1,\delta {\theta }_p}\\ {\sigma }_{\delta {\theta }_2,\delta {\theta }_1}& {\sigma }_{\delta {\theta }_2}^2& \dots & {\sigma }_{\delta {\theta }_2,\delta {\theta }_p}\\ ⋮& ⋮& \ddots & ⋮\\ {\sigma }_{\delta {\theta }_p,\delta {\theta }_1}& {\sigma }_{\delta {\theta }_p,\delta {\theta }_2}& \dots & {\sigma }_{\delta {\theta }_p}^2\end{array}\right].$$

(24)

Diagonal elements are computed from tolerance specifications using, as an example eq. (23); off-diagonal elements (covariances between different parameters) are typically set to zero in the absence of specific information; assuming that component parameters are statistically independent is a sound hypothesis most of the time (there could be exceptions for special arrangement such as current or voltage dividers where two or more components are known to track each other etc...). When assuming: (i) uniform distribution, (ii) diagonal covariance matrix and (iii) symmetry around the mean value, it is straightforward to calculate all higher order moments :

$$\mathbb{E}\left[\delta {\theta }_i^n\right]=\left\{\begin{array}{cc}{\displaystyle \frac{\Delta {\theta }_i^n}{n+1}}\hfill & n\mathrm{even}\hfill \\ 0\hfill & n\mathrm{odd}\hfill \end{array}\right.$$

(25)

What is presented in the rest of the paper does not specifically assume the uniform distribution or the diagonality of the covariance matrix, however the derivation of the proposed robust Kalman filter is driven exactly by the practical case where only tolerances of the components are known to the filter designer.

2.2. First-Order Statistical Modelling

The state can be expressed explicitly as:

$$\left\{\begin{array}{cc}{x}_k\hfill & ={\left(A+\delta A\right)}^k{x}_0+{\displaystyle \sum _{j=0}^{k-1}}{\left(A+\delta A\right)}^{k-1-j}\left[\left(B+\delta B\right)\left(u_j+\delta u_j\right)+w_j^{*}\right],\hfill \\ y_k\hfill & =\left(H+\delta H\right)x_k+v_k^{*}.\hfill \end{array}\right.$$

(26)

Even assuming statistical independence of $x_0$ from $\delta A$ (which is a completely reasonable assumption), the problem as expressed in eq. (26) is very complex in general terms. The evolution of the state is governed by ${(A+\delta A)}^k$ which can be expressed as follows:

$$\begin{array}{cc}\hfill {\left(A+\delta A\right)}^k& =A^k+{\displaystyle \sum _{i=0}^{k-1}}{A}^{k-1-i}\delta A{A}^i+{\alpha }_k\left(\delta A\right),\hfill \end{array}$$

(27)

where ${\alpha }_k\left(\delta A\right)$ contains all terms of order $O\left(\delta A^2\right)$ (i.e. terms like

${\sum }_{0\le i<j,i+j\le k-2}{A}^{k-2-i-j}\delta A{A}^i\delta A{A}^j$) and higher-order ones in the expansion of ${(A+\delta A)}^k$ (note that ${\alpha }_k\left(\delta A\right)=0$ for $k<0$).

From eqs. (26) and (27) the state can be expressed as:

$$\begin{array}{cc}\hfill x_k& =\left(A^k+{\displaystyle \sum _{i=0}^{k-1}}{A}^{k-1-i}\delta A{A}^i+{\alpha }_k\left(\delta A\right)\right)x_0\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{k-1}}\left[A^{k-1-j}(B+\delta B)+\left({\displaystyle \sum _{i=0}^{k-2-j}}{A}^{k-2-j-i}\delta A{A}^i\right)(B+\delta B)+{\alpha }_{k-1-j}\left(\delta A\right)(B+\delta B)\right](u_j+\delta u_j)\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{k-1}}\left(A^{k-1-j}+{\displaystyle \sum _{i=0}^{k-2-j}}{A}^{k-2-j-i}\delta A{A}^i+{\alpha }_{k-1-j}\left(\delta A\right)\right)w_j^{*},\hfill \end{array}$$

(28)

From eq. (28) one can decompose $x_k$ as follows:

$$\begin{array}{cc}\hfill x_k& =x_k^{nom}+\delta x_k+\delta x_k^{ho}=A^k{x}_0+{\displaystyle \sum _{j=0}^{k-1}}{A}^{k-1-j}B{u}_j+{\widetilde{\delta A}}_{k-1}{x}_0+{\displaystyle \sum _{j=0}^{k-1}}{A}^{k-1-j}\delta B{u}_j\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{k-1}}{\widetilde{\delta A}}_{k-2-j}B{u}_j+\delta x_k^{ho},\hfill \end{array}$$

(29)

where:

$$\left\{\begin{array}{cc}{x}_k^{nom}\hfill & =A^k{x}_0+{\displaystyle \sum _{j=0}^{k-1}}{A}^{k-1-j}B{u}_j,\hfill \\ \delta x_k\hfill & ={\widetilde{\delta A}}_{k-1}{x}_0+{\displaystyle \sum _{j=0}^{k-1}}{A}^{k-1-j}\delta B{u}_j+{\displaystyle \sum _{j=0}^{k-1}}{\widetilde{\delta A}}_{k-2-j}B{u}_j,\hfill \end{array}\right.$$

(30)

represent, respectively, the nominal state evolution (i.e. the state evolution assuming perfect knowledge of the state-space model) and the first order perturbation in $\delta A$ and $\delta B$, or equivalently in $\delta \theta$. It is worth noting that $\mathbb{E}\left[\delta x_k\right]=0$ because of the independence of $x_0$ from $\delta \theta$ and $\mathbb{E}\left[\delta A\right]=0$ and $\mathbb{E}\left[\delta B\right]=0$. The fundamental idea, in this work, is to design a robust Kalman filter neglecting the higher order contributions in $\delta x_k^{ho}$. This is similar to the MVFOSM approach (where only the zeroth order is considered for the mean and only contributions up to the second order are considered for the second moments such as the covariances) but slightly more accurate as it also allows to take into account the state and output biases induced by the model uncertainty, by augmenting the state as detailed in the next section.

2.3. Augmented Model

To build the augmented model it can be observed (thanks to the statistical independence assumptions) that:

$$\left\{\begin{array}{c}{z}_k=\mathbb{E}\left[w_k\right]=\mathbb{E}\left[\delta A{x}_k\right],\hfill \\ t_k=\mathbb{E}\left[v_k\right]=\mathbb{E}\left[\delta H{x}_k\right].\hfill \end{array}\right.$$

(31)

For the additional states, the dynamics is as follows:

$$\left\{\begin{array}{cc}{z}_{k+1}\hfill & =\mathbb{E}\left[\delta A{x}_{k+1}\right]=\mathbb{E}\left[\delta AA{x}_k\right]+\mathbb{E}\left[\delta A\delta A{x}_k\right]+\mathbb{E}\left[\delta A\delta B{u}_k\right]\hfill \\ \hfill & \approx \mathbb{E}\left[\delta AA{x}_k\right]+\mathbb{E}\left[\delta A\delta A\right]{\overline{x}}_k^{nom}+\mathbb{E}\left[\delta A\delta B\right]u_k,\hfill \\ t_{k+1}\hfill & =\mathbb{E}\left[\delta H{x}_{k+1}\right]=\mathbb{E}\left[\delta HA{x}_k\right]+\mathbb{E}\left[\delta H\delta A{x}_k\right]+\mathbb{E}\left[\delta H\delta B{u}_k\right]\hfill \\ \hfill & \approx \mathbb{E}\left[\delta HA{x}_k\right]+\mathbb{E}\left[\delta H\delta A\right]{\overline{x}}_k^{nom}+\mathbb{E}\left[\delta H\delta B\right]u_k,\hfill \end{array}\right.$$

(32)

where in (32) all contributions that are higher than second order have been neglected.

In eq. (32) the difficulty comes from the fact that A does not commute with $\delta A$ and $\delta H$ so the terms $\mathbb{E}\left[\delta AA{x}_k\right]$ and $\mathbb{E}\left[\delta HA{x}_k\right]$ are not immediately expressible as functions of $z_k$ and $t_k$ respectively. However, it can be observed that:

$$\left\{\begin{array}{cc}\mathbb{E}\left[\delta AA{x}_k\right]\hfill & =\mathrm{vec}\left(\mathbb{E}\left[\delta AA{x}_k\right]\right)=\mathbb{E}\left[\left(x_k^T\otimes \delta A\right)\mathrm{vec}\left(A\right)\right]=\mathbb{E}[x_k^T\otimes \delta A]\mathrm{vec}\left(A\right),\hfill \\ \mathbb{E}\left[\delta HA{x}_k\right]\hfill & =\mathrm{vec}\left(\mathbb{E}\left[\delta HA{x}_k\right]\right)=\mathbb{E}\left[\left(x_k^T\otimes \delta H\right)\mathrm{vec}\left(A\right)\right]=\mathbb{E}[x_k^T\otimes \delta H]\mathrm{vec}\left(A\right),\hfill \end{array}\right.$$

(33)

and that:

$$\left\{\begin{array}{cc}{z}_{k+1}=\hfill & \mathbb{E}\left[\delta A{x}_{k+1}\right]=\mathrm{vec}\left(\mathbb{E}\left[\delta A{I}_n{x}_{k+1}\right]\right)=\mathbb{E}[x_{k+1}^T\otimes \delta A]\mathrm{vec}\left(I_n\right),\hfill \\ t_{k+1}=\hfill & \mathbb{E}\left[\delta H{x}_{k+1}\right]=\mathrm{vec}\left(\mathbb{E}\left[\delta H{I}_n{x}_{k+1}\right]\right)=\mathbb{E}[x_{k+1}^T\otimes \delta H]\mathrm{vec}\left(I_n\right).\hfill \end{array}\right.$$

(34)

Therefore, the recursion can be written as:

$$\left\{\begin{array}{cc}{z}_{k+1}=\hfill & \mathbb{E}[x_{k+1}^T\otimes \delta A]\mathrm{vec}\left(I_n\right)\approx \mathbb{E}[x_k^T\otimes \delta A]\mathrm{vec}\left(A\right)+\mathbb{E}\left[\delta A\delta A\right]{\overline{x}}_k^{nom}+\mathbb{E}\left[\delta A\delta B\right]u_k,\hfill \\ t_{k+1}=\hfill & \mathbb{E}[x_{k+1}^T\otimes \delta H]\mathrm{vec}\left(I_n\right)\approx \mathbb{E}[x_k^T\otimes \delta H]\mathrm{vec}\left(A\right)+\mathbb{E}\left[\delta H\delta A\right]{\overline{x}}_k^{nom}+\mathbb{E}\left[\delta H\delta B\right]u_k.\hfill \end{array}\right.$$

(35)

The recursion, as written in eq. (35), is still implicit; to be able to actually run it, $z_k$ and $t_k$ need to be expressed explicitly. This can be done as follows; first note, by means of the definition in eq. (34), that:

$$\left\{\begin{array}{cc}{z}_k\hfill & =\mathbb{E}[x_k^T\otimes \delta A]\mathrm{vec}\left(I_n\right),\hfill \\ t_k\hfill & =\mathbb{E}[x_k^T\otimes \delta H]\mathrm{vec}\left(I_n\right),\hfill \end{array}\right.$$

(36)

then define:

$$\left\{\begin{array}{cc}{Z}_k\hfill & =\mathbb{E}[x_k^T\otimes \delta A]\in {\mathbb{R}}^{n\times n^2},\hfill \\ T_k\hfill & =\mathbb{E}[x_k^T\otimes \delta H]\in {\mathbb{R}}^{m\times n^2},\hfill \end{array}\right.$$

(37)

the recursion for these newly introduced matrices is therefore the following:

$$\left\{\begin{array}{cc}{Z}_{k+1}\hfill & \approx Z_k\left(A^T\otimes I_n\right)+\left({\overline{x}}_k^{nom}\otimes I_n\right)\mathbb{E}[\delta A^T\otimes \delta A]+\left(u_k^T\otimes I_n\right)\mathbb{E}[\delta B^T\otimes \delta A]\hfill \\ T_{k+1}\hfill & \approx T_k\left(A^T\otimes I_n\right)+\left({\overline{x}}_k^{nom}\otimes I_m\right)\mathbb{E}[\delta A^T\otimes \delta H]+\left(u_k^T\otimes I_m\right)\mathbb{E}[\delta B^T\otimes \delta H].\hfill \end{array}\right.$$

(38)

The derivation of eq. (38) is detailed in the Appendix.

The following holds true for the expected values of state and output:

$$\left\{\begin{array}{cc}{\overline{x}}_{k+1}\hfill & \approx A{\overline{x}}_k+B{u}_k+z_k\hfill \\ Z_{k+1}\hfill & \approx Z_k\left(A^T\otimes I_n\right)+\left({\overline{x}}_k^{nom}\otimes I_n\right)\mathbb{E}[\delta A^T\otimes \delta A]+\left(u_k^T\otimes I_n\right)\mathbb{E}[\delta B^T\otimes \delta A]\hfill \\ z_k\hfill & =Z_k\mathrm{vec}\left(I_n\right)\hfill \\ T_{k+1}\hfill & \approx T_k\left(A^T\otimes I_n\right)+\left({\overline{x}}_k^{nom}\otimes I_m\right)\mathbb{E}[\delta A^T\otimes \delta H]+\left(u_k^T\otimes I_m\right)\mathbb{E}[\delta B^T\otimes \delta H]\hfill \\ t_k\hfill & =T_k\mathrm{vec}\left(I_n\right)\hfill \\ {\overline{y}}_k\hfill & \approx H{\overline{x}}_k+t_k\hfill \end{array}\right.$$

(39)

Equation (39) can be further simplified as follows by recognizing that ${\overline{x}}_k^{nom}\approx {\overline{x}}_k$ up to $O\left(\delta {\theta }^2\right)$ :

$$\left\{\begin{array}{cc}{\overline{x}}_{k+1}\hfill & \approx A{\overline{x}}_k+B{u}_k+z_k\hfill \\ Z_{k+1}\hfill & \approx Z_k\left(A^T\otimes I_n\right)+\left({\overline{x}}_k\otimes I_n\right)\mathbb{E}[\delta A^T\otimes \delta A]+\left(u_k^T\otimes I_n\right)\mathbb{E}[\delta B^T\otimes \delta A]\hfill \\ z_k\hfill & =Z_k\mathrm{vec}\left(I_n\right)\hfill \\ T_{k+1}\hfill & \approx T_k\left(A^T\otimes I_n\right)+\left({\overline{x}}_k\otimes I_m\right)\mathbb{E}[\delta A^T\otimes \delta H]+\left(u_k^T\otimes I_m\right)\mathbb{E}[\delta B^T\otimes \delta H]\hfill \\ t_k\hfill & =T_k\mathrm{vec}\left(I_n\right)\hfill \\ {\overline{y}}_k\hfill & \approx H{\overline{x}}_k+t_k\hfill \end{array}\right.$$

(40)

this simplification is very important as it allows reducing the computational burden by avoiding the introduction of yet another state $d_k=\mathbb{E}\left[\delta x_k\right]$ to keep track of the difference between ${\overline{x}}_k$ and ${\overline{x}}_k^{nom}$ as detailed in the appendix.

3. Steady-State Kalman Filter

Equations (15) and (40) are sufficient to calculate the steady-state Kalman filter when assuming $Q,R$ and M to be time invariant (which will be discussed in Section 3.3) and that the system in (40) is stable (this is so because the homogeneous dynamics of $Z_k$ and $T_k$ are governed by the matrix $(A^T\otimes I_n)\in {\mathbb{R}}^{n^2\times n^2}$ whose eigenvalues are those of A each with multiplicity n; and the ${\overline{x}}_k$ equation is driven by $z_k$ which decays to a constant for constant input. Since all eigenvalues of A are strictly within the unit circle, all modes are stable.). Let the covariance matrix of the state estimation error, in steady-state, be denoted as P. It can be computed by solving the following generalized DARE (discrete algebraic Riccati equation) as reported in [2]:

$$P=AP{A}^T-A\left(P{H}^T+M\right){\left(HP{H}^T+HM+M^T{H}^T+R\right)}^{-1}\left(HP+M^T\right)A^T+Q$$

(41)

Once P is calculated, the steady-state Kalman gain K can be calculated as follows:

$$K=\left(P{H}^T+M\right){\left(HP{H}^T+HM+M^T{H}^T+R\right)}^{-1}$$

(42)

It is useful to introduce the following Choleski decomposition: $L{L}^T=Q-M{R}^{-1}{M}^T$. A unique positive semi-definite solution P of eq. (41) exists under the following conditions [2]:

$$\left\{\begin{array}{c}Q⪰0,\hfill \\ R\succ 0,\hfill \\ (A,H)\mathrm{detectable},\hfill \\ (A-M{R}^{-1}H,L)\mathrm{stabilizable}.\hfill \end{array}\right.$$

(43)

In such a case the solution leads to a stable steady-state Kalman filter i.e. all eigenvalues of $(I_n-KH)A$ are strictly inside the unit circle. The detectability condition is satisfied when all undetectable modes of $(A,H)$ are stable. For a stable open-loop system (all eigenvalues of A strictly inside the unit circle), detectability is automatically satisfied. The stabilizability of $(A-M{R}^{-1}H,L)$ can be verified once Q, R, M are computed, typically by checking the controllability Gramian or by attempting to solve the DARE numerically. Their calculation is discussed in the remaining of this section.

3.1. Steady-State of the Augmented Model

It is not straightforward to calculate the asymptotic response (i.e. assuming $u_k$ has been kept constant for long enough) from eq. (40) so a different way is presented here. Considering the first order perturbation $\delta x_k$ introduced in eq. (28) the additional states$z_k$ and $t_k$ can be written explicitly as:

$$\left\{\begin{array}{cc}{z}_k\hfill & =\mathbb{E}\left[\delta A{x}_k\right]\approx E\left[\delta A\delta x_k\right]=\mathbb{E}\left[\delta A{\widetilde{\delta A}}_{k-1}\right]\mathbb{E}\left[x_0\right]+\mathbb{E}\left[\delta A{\displaystyle \sum _{j=0}^{k-1}}{A}^{k-1-j}\delta B\right]u_j\hfill \\ \hfill & +\mathbb{E}\left[\delta A{\displaystyle \sum _{j=0}^{k-1}}{\widetilde{\delta A}}_{k-2-j}\right]B{u}_j,\hfill \\ t_k\hfill & =\mathbb{E}\left[\delta H{x}_k\right]\approx E\left[\delta H\delta x_k\right]=\mathbb{E}\left[\delta H{\widetilde{\delta A}}_{k-1}\right]\mathbb{E}\left[x_0\right]+\mathbb{E}\left[\delta H{\displaystyle \sum _{j=0}^{k-1}}{A}^{k-1-j}\delta B\right]u_j\hfill \\ \hfill & +\mathbb{E}\left[\delta H{\displaystyle \sum _{j=0}^{k-1}}{\widetilde{\delta A}}_{k-2-j}\right]B{u}_j.\hfill \end{array}\right.$$

(44)

It is worth introducing the following:

$$W={(I_n-A)}^{-1}$$

(45)

The first term of both expressions in eq. (44), proportional to $\mathbb{E}\left[x_0\right]$, vanishes as $k\to \infty$ (the proof is reported in the appendix).

The steady-state values of the additional states read (with $u_j=u_k=u$):

$$\left\{\begin{array}{cc}{z}_{\infty }\hfill & \approx \mathbb{E}\left[\delta AW\delta B\right]u_k+\mathbb{E}\left[\delta AW\delta A\right]WB{u}_k,\hfill \\ t_{\infty }\hfill & \approx \mathbb{E}\left[\delta HW\delta B\right]u_k+\mathbb{E}\left[\delta HW\delta A\right]WB{u}_k.\hfill \end{array}\right.$$

(46)

The proof of eq. (46) is also given in the appendix.

So the steady-state behaviour of the expected value is described by the following:

$$\left\{\begin{array}{cc}{\overline{x}}_{\infty }\hfill & \approx A{\overline{x}}_{\infty }+B{u}_k+z_{\infty },\hfill \\ z_{\infty }\hfill & \approx \mathbb{E}\left[\delta AW\delta B\right]u_k+\mathbb{E}\left[\delta AW\delta A\right]WB{u}_k,\hfill \\ t_{\infty }\hfill & \approx \mathbb{E}\left[\delta HW\delta B\right]u_k+\mathbb{E}\left[\delta HW\delta A\right]WB{u}_k.\hfill \end{array}\right.$$

(47)

The solution of the algebraic system in eq. (47) can be written as follows:

$$\left\{\begin{array}{cc}{\overline{x}}_{\infty }\hfill & \approx F_x{u}_k\hfill \\ z_{\infty }\hfill & \approx F_z{u}_k\hfill \\ t_{\infty }\hfill & \approx F_t{u}_k\hfill \end{array}\right.$$

(48)

letting $F_{nom}={(I_n-A)}^{-1}B=WB$ be the nominal DC gain, the augmented model gains are:

$$\left\{\begin{array}{cc}{F}_x\hfill & ={(I_n-A)}^{-1}(B+F_z)=F_{nom}+W{F}_z,\hfill \\ F_z\hfill & =\mathbb{E}\left[\delta AW\delta B\right]+\mathbb{E}\left[\delta AW\delta A\right]F_{nom},\hfill \\ F_t\hfill & =\mathbb{E}\left[\delta HW\delta B\right]+\mathbb{E}\left[\delta HW\delta A\right]F_{nom}.\hfill \end{array}\right.$$

(49)

For the actual computation of the above gains, the following relationships are useful:

$$\left\{\begin{array}{cc}\mathbb{E}\left[\delta AW\delta B\right]\hfill & =\mathrm{unvec}\left(\mathbb{E}[\delta B^T\otimes \delta A]\mathrm{vec}\left(W\right)\right)\approx \mathrm{unvec}\left[{\Sigma }_{\left(BA\right)}^{*}\mathrm{vec}\left(W\right)\right],\hfill \\ \mathbb{E}\left[\delta AW\delta A\right]\hfill & =\mathrm{unvec}\left(\mathbb{E}[\delta A^T\otimes \delta A]\mathrm{vec}\left(W\right)\right)\approx \mathrm{unvec}\left[{\Sigma }_{\left(AA\right)}^{*}\mathrm{vec}\left(W\right)\right],\hfill \\ \mathbb{E}\left[\delta HW\delta B\right]\hfill & =\mathrm{unvec}\left(\mathbb{E}[\delta B^T\otimes \delta H]\mathrm{vec}\left(W\right)\right)\approx \mathrm{unvec}\left[{\Sigma }_{\left(BH\right)}^{*}\mathrm{vec}\left(W\right)\right],\hfill \\ \mathbb{E}\left[\delta HW\delta A\right]\hfill & =\mathrm{unvec}\left(\mathbb{E}[\delta A^T\otimes \delta H]\mathrm{vec}\left(W\right)\right)\approx \mathrm{unvec}\left[{\Sigma }_{\left(AH\right)}^{*}\mathrm{vec}\left(W\right)\right].\hfill \end{array}\right.$$

(50)

Finally the DC gains of the augmented model are expressed by the following:

$$\left\{\begin{array}{cc}{F}_x\hfill & =W(B+F_z)=F_{nom}+W{F}_z,\hfill \\ F_z\hfill & \approx \mathrm{unvec}\left[{\Sigma }_{\left(BA\right)}^{*}\mathrm{vec}\left(W\right)\right]+\mathrm{unvec}\left[{\Sigma }_{\left(AA\right)}^{*}\mathrm{vec}\left(W\right)\right]F_{nom},\hfill \\ F_t\hfill & \approx \mathrm{unvec}\left[{\Sigma }_{\left(BH\right)}^{*}\mathrm{vec}\left(W\right)\right]+\mathrm{unvec}\left[{\Sigma }_{\left(AH\right)}^{*}\mathrm{vec}\left(W\right)\right]F_{nom}.\hfill \end{array}\right.$$

(51)

3.2. Calculation of Q, R and M

In accordance with the MVFOSM approach in the calculation of the Q, R and M matrices only contributions up to the second order in $\delta A$, $\delta B$ and $\delta H$, or, in general, $O\left(\delta {\theta }^2\right)$ will be considered. However, the terms $z_k{z}_k^T$, $t_k{t}_k^T$ and $z_k{t}_k^T$ although being $O\left(\delta {\theta }^4\right)$ will be retained to ensure that $Q_k$, $R_k$ and $M_k$ exactly represent the covariances of the equivalent noises: $\mathbb{E}\left[(w_k-{\overline{w}}_k){(w_k-{\overline{w}}_k)}^T\right]$, $\mathbb{E}\left[(v_k-{\overline{v}}_k){(v_k-{\overline{v}}_k)}^T\right]$ and $\mathbb{E}\left[(w_k-{\overline{w}}_k){(v_k-{\overline{v}}_k)}^T\right]$ as required by standard Kalman filter.

The following matrices are going to be useful in the upcoming calculations:

$$J_{A_c}^A=J_A{J}_{A_c}\in {\mathbb{R}}^{n^2\times p},$$

(52)

so

$${\Sigma }_{\delta A}=\mathbb{E}\left[\delta A\delta A^T\right]=\mathbb{E}\left[J_A{J}_{A_c}\delta \theta \delta {\theta }^T{J}_{A_c}^T{J}_A^T\right]\approx J_{A_c}^A\mathbb{E}\left[\delta \theta \delta {\theta }^T\right]J_{A_c}^{A^T}=J_{A_c}^A{\Sigma }_{\delta \theta }{J}_{A_c}^{A^T}.$$

(53)

Furthermore:

$${\Sigma }_{\delta u}=\mathbb{E}\left[\delta u_k\delta u_k^T\right],$$

(54)

where ${\Sigma }_{\delta u}\in {\mathbb{R}}^{q\times q}$.

It is also useful to define the following matrices:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill {\Sigma }_{\left(AA\right)}& ={\Sigma }_{\delta A}=\mathbb{E}\left[\delta A\delta A^T\right]\approx J_{A_c}^A{\Sigma }_{\delta \theta }{J}_{A_c}^{A^T}\in {\mathbb{R}}^{n^2\times n^2},\hfill \\ \hfill {\Sigma }_{\left(BB\right)}& ={\Sigma }_{\delta B}=\mathbb{E}\left[\delta B\delta B^T\right]\approx V{\Sigma }_{\delta \theta }{V}^T\in {\mathbb{R}}^{nq\times nq},\hfill \\ \hfill {\Sigma }_{\left(HH\right)}& ={\Sigma }_{\delta H}=\mathbb{E}\left[\delta H\delta H^T\right]\approx J_H{\Sigma }_{\delta \theta }{J}_H^T\in {\mathbb{R}}^{nm\times nm},\hfill \\ \hfill {\Sigma }_{\left(AB\right)}& =\mathbb{E}\left[\delta A\delta B^T\right]\approx J_{A_c}^A{\Sigma }_{\delta \theta }{V}^T\in {\mathbb{R}}^{n^2\times nq},\hfill \\ \hfill {\Sigma }_{\left(AH\right)}& =\mathbb{E}\left[\delta A\delta H^T\right]\approx J_{A_c}^A{\Sigma }_{\delta \theta }{J}_H^T\in {\mathbb{R}}^{n^2\times nm},\hfill \\ \hfill {\Sigma }_{\left(BH\right)}& =\mathbb{E}\left[\delta B\delta H^T\right]\approx V{\Sigma }_{\delta \theta }{J}_H^T\in {\mathbb{R}}^{nq\times nm},\hfill \\ \hfill {\Sigma }_{\left(AA\right)}^{*}& =\mathbb{E}[\delta A^T\otimes \delta A]\approx r\left({\Sigma }_{\left(AA\right)}\right),\hfill \\ \hfill {\Sigma }_{\left(BA\right)}^{*}& =\mathbb{E}[\delta B^T\otimes \delta A]\approx r\left({\Sigma }_{\left(AB\right)}\right),\hfill \\ \hfill {\Sigma }_{\left(AH\right)}^{*}& =\mathbb{E}[\delta A^T\otimes \delta H]\approx r\left({\Sigma }_{\left(AH\right)}\right),\hfill \\ \hfill {\Sigma }_{\left(BH\right)}^{*}& =\mathbb{E}[\delta B^T\otimes \delta H]\approx r\left({\Sigma }_{\left(BH\right)}\right),\hfill \end{array}\hfill \end{array}\right.$$

(55)

where $r(.)$ is a reshaping function that translates the entries of the un-starred matrices in eq. (55) into the entries of the starred matrices that appear in the augmented state. Detailed expression of $r(.)$ is given in the appendix.

It is important to highlight that the equalities in eq. (55) have to be interpreted taking into account the abused notation that implicitly assumes vectorization and de-vectorization e.g. for ${\Sigma }_{\left(AB\right)}=\mathbb{E}\left[\delta A\delta B^T\right]$ the product $\delta A\delta B^T$ does not ordinarily exist as $\delta A\in {\mathbb{R}}^{n\times n}$ and $\delta B^T\in {\mathbb{R}}^{q\times n}$ so they cannot be multiplied; for the calculation of ${\Sigma }_{\left(AB\right)}$ one would more correctly need to write ${\Sigma }_{\left(AB\right)}=\mathbb{E}\left[\mathrm{vec}\left(\delta A\right)\mathrm{vec}\left(\delta B^T\right)\right]=\mathbb{E}\left[\mathrm{vec}\left[\mathrm{unvec}\left(J_{A_c}^A\delta \theta \right)\right]\mathrm{vec}\left[\mathrm{unvec}\left(\delta {\theta }^T{V}^T\right)\right]\right]=\mathbb{E}\left[J_{A_c}^A\delta \theta \delta {\theta }^T{V}^T\right]=J_{A_c}^A{\Sigma }_{\delta \theta }{V}^T$ and so on.

3.2.1. Calculation of Q

Neglecting the terms higher than quadratic it yields:

$$\begin{array}{cc}\hfill P_{w{w}_k}& \approx \mathbb{E}\left[w_k{w}_k^T\right]\approx \mathbb{E}\left[\delta A{x}_k{x}_k^T\delta A^T\right]+\mathbb{E}\left[\delta A{x}_k{u}_k^T\delta B^T\right]+\mathbb{E}\left[\delta B{u}_k{x}_k^T\delta A^T\right]+\hfill \\ \hfill & +\mathbb{E}\left[\delta B{u}_k{u}_k^T\delta B^T\right]+\mathbb{E}\left[B\delta u_k\delta u_k^T{B}^T\right]+\mathbb{E}\left[w_k^{*}{w}_k^{*T}\right],\hfill \end{array}$$

(56)

where the other terms are exactly zero because of statistical independence or approximately zero due to MVFOSM truncation.

The terms in eq. (56) can be expressed as follows:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \mathbb{E}\left[\delta A{x}_k{x}_k^T\delta A^T\right]& \approx \mathbb{E}\left[\delta A{\overline{x}}_k^{nom}{\overline{x}}_k^{no{m}^T}\delta A^T\right]\hfill \\ \hfill & \approx (u_k^T\otimes I_n)\left(F_{nom}^T\otimes I_n\right){\Sigma }_{\left(AA\right)}(F_{nom}\otimes I_n)(u_k\otimes I_n),\hfill \\ \hfill \mathbb{E}\left[\delta A{x}_k{u}_k^T\delta B^T\right]& \approx \mathbb{E}\left[\delta A{\overline{x}}_k^{nom}{u}_k^T\delta B^T\right]\approx \mathbb{E}\left[\delta A{F}_{nom}{u}_k{u}_k^T\delta B^T\right]\hfill \\ \hfill & \approx (u_k^T\otimes I_n)\left(F_{nom}^T\otimes I_n\right){\Sigma }_{\left(AB\right)}(u_k\otimes I_n),\hfill \\ \hfill \mathbb{E}\left[\delta B{u}_k{x}_k^T\delta A^T\right]& \approx \mathbb{E}\left[\delta B{u}_k{\overline{x}}_k^{no{m}^T}\delta A^T\right]\approx \mathbb{E}\left[\delta B{u}_k{u}_k^T{F}_{nom}^T\delta A^T\right]\hfill \\ \hfill & \approx (u_k^T\otimes I_n){\Sigma }_{\left(AB\right)}^T\left(F_{nom}\otimes I_n\right)(u_k\otimes I_n),\hfill \\ \hfill \mathbb{E}\left[\delta B{u}_k{u}_k^T\delta B^T\right]& =(u_k^T\otimes I_n)\mathbb{E}\left[\delta B\delta B^T\right](u_k\otimes I_n)\approx (u_k^T\otimes I_n){\Sigma }_{\left(BB\right)}(u_k\otimes I_n),\hfill \\ \hfill \mathbb{E}\left[B\delta u_k\delta u_k^T{B}^T\right]& =B\mathbb{E}\left[\delta u_k\delta u_k^T\right]B^T=B{\Sigma }_{\delta u}{B}^T,\hfill \\ \hfill \mathbb{E}\left[w_k^{*}{w}_k^{*T}\right]& ={\Sigma }_{w_k^{*}}=Q_{w^{*}}\hfill \end{array}\hfill \end{array}\right.$$

(57)

where it has been considered that $x_k=x_k^{nom}+\delta x_k^{full}$ and $\delta x_k^{full}$ contains all the terms depending on $\delta A$ and $\delta B$; when calculating the second moments only the term $x_k^{nom}$ contributes up to the second order as $\delta x_k^{full}$ would only contribute to third order or higher. As an example when computing $E[\delta A{x}_k{x}_k^T\delta A^T$], expanding the products yields:

☐

$\mathbb{E}\left[\delta A{x}_k^{nom}{x}_k^{no{m}^T}\delta A^T\right]$: $O\left(\delta {\theta }^2\right)$ - retained,

☐

$\mathbb{E}\left[\delta A{x}_k^{nom}\delta x_k^{ful{l}^T}\delta A^T\right]$ and transpose: $O\left(\delta {\theta }^3\right)$ - neglected

☐

$\mathbb{E}\left[\delta A\delta x_k^{full}\delta x_k^{ful{l}^T}\delta A^T\right]$: $O\left(\delta {\theta }^4\right)$ -neglected.

It is important to note that $x_k^{nom}$ depends only on $x_0$ (and the history of $u_k$ that is however deterministic) which is independent of $\delta \theta$ and therefore of $\delta A$, $\delta B$ and $\delta H$ so $\mathbb{E}\left[\delta A{x}_k^{nom}{x}_k^{no{m}^T}\delta A^T\right]=\mathbb{E}\left[\delta A{\overline{x}}_k^{nom}{\overline{x}}_k^{no{m}^T}\delta A^T\right]$. Thus, retaining only second-order terms, ${\overline{x}}_k^{nom}=\mathbb{E}\left[x_k^{nom}\right]$ is used in place of $x_k$, which is consistent with MVFOSM-like truncation. Furthermore it has also been considered that the terms containing $A^k\mathbb{E}\left[x_0\right]$ vanishes in steady-state (as A has eigenvalues strictly inside the unit circle).

As already done for the augmented state, ${\overline{x}}_k^{nom}$ can be approximated by ${\overline{x}}_k$ leading to:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \mathbb{E}\left[\delta A{x}_k{x}_k^T\delta A^T\right]& \approx \mathbb{E}\left[\delta A{\overline{x}}_k{\overline{x}}_k^T\delta A^T\right]\approx (u_k^T\otimes I_n)\left(F_x^T\otimes I_n\right){\Sigma }_{\left(AA\right)}(F_x\otimes I_n)(u_k\otimes I_n),\hfill \\ \hfill \mathbb{E}\left[\delta A{x}_k{u}_k^T\delta B^T\right]& \approx \mathbb{E}\left[\delta A{\overline{x}}_k{u}_k^T\delta B^T\right]\approx \mathbb{E}\left[\delta A{F}_x{u}_k{u}_k^T\delta B^T\right]\hfill \\ \hfill & \approx (u_k^T\otimes I_n)\left(F_x^T\otimes I_n\right){\Sigma }_{\left(AB\right)}(u_k\otimes I_n),\hfill \\ \hfill \mathbb{E}\left[\delta B{u}_k{x}_k^T\delta A^T\right]& \approx \mathbb{E}\left[\delta B{u}_k{\overline{x}}_k^T\delta A^T\right]\approx \mathbb{E}\left[\delta B{u}_k{u}_k^T{F}_x^T\delta A^T\right]\hfill \\ & \approx (u_k^T\otimes I_n){\Sigma }_{\left(AB\right)}^T\left(F_x\otimes I_n\right)(u_k\otimes I_n),\hfill \\ \hfill \mathbb{E}\left[\delta B{u}_k{u}_k^T\delta B^T\right]& =(u_k^T\otimes I_n)\mathbb{E}\left[\delta B\delta B^T\right](u_k\otimes I_n)\approx (u_k^T\otimes I_n){\Sigma }_{\left(BB\right)}(u_k\otimes I_n),\hfill \\ \hfill \mathbb{E}\left[B\delta u_k\delta u_k^T{B}^T\right]& =B\mathbb{E}\left[\delta u_k\delta u_k^T\right]B^T=B{\Sigma }_{\delta u}{B}^T,\hfill \\ \hfill \mathbb{E}\left[w_k^{*}{w}_k^{*T}\right]& ={\Sigma }_{w_k^{*}}=Q_{w^{*}}.\hfill \end{array}\hfill \end{array}\right.$$

(58)

In the remainder of the paper this approximation is going to be adopted for the covariances computation by using the DC gain $F_x$ instead of $F_{nom}$.

The following time invariant covariance matrix is useful for the upcoming computations:

$$Q_0=Q_{w^{*}}+B{\Sigma }_{\delta u}{B}^T.$$

(59)

Thus, the full expression of $P_{w{w}_k}$ can be written as:

$$\begin{array}{cc}\hfill P_{w{w}_k}& \approx Q_0+(u_k^T\otimes I_n)P_u(u_k\otimes I_n)\hfill \\ \hfill P_u& =\left(F_x^T\otimes I_n\right){\Sigma }_{\left(AA\right)}(F_x\otimes I_n)+(F_x^T\otimes I_n){\Sigma }_{\left(AB\right)}+{\Sigma }_{\left(AB\right)}^T(F_x\otimes I_n)+{\Sigma }_{\left(BB\right)}.\hfill \end{array}$$

(60)

Finally:

$$Q_k^{th}\approx P_{w{w}_k}-z_k{z}_k^T=P_{w{w}_k}-F_z{u}_k{u}_k^T{F}_z^T.$$

(61)

3.2.2. Calculation of R

$$P_{v{v}_k}=\mathbb{E}\left[v_k{v}_k^T\right]=\mathbb{E}\left[\delta H{x}_k{x}_k^T\delta H^T\right]+\mathbb{E}\left[v_k^{*}{v}_k^{*T}\right]\approx \mathbb{E}\left[\delta H{\overline{x}}_k^{nom}{\overline{x}}_k^{no{m}^T}\delta H^T\right]+\mathbb{E}\left[v_k^{*}{v}_k^{*T}\right],$$

(62)

as the expectations $\mathbb{E}\left[\delta H{x}_k{v}_k^{*T}\right]$ and $\mathbb{E}\left[v_k^{*}{x}_k^T\delta H^T\right]$ are both (strictly) zero.

The individual contributions are (again using ${\overline{x}}_k$ instead of ${\overline{x}}_k^{nom}$):

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \mathbb{E}\left[\delta H{\overline{x}}_k{\overline{x}}_k^T\delta H^T\right]& \approx \mathbb{E}\left[\delta H{F}_x{u}_k{u}_k^T{F}_x^T\delta H^T\right]\hfill \\ \hfill & \approx (u_k^T\otimes I_m)\left(F_x^T\otimes I_m\right){\Sigma }_{\left(HH\right)}(F_x\otimes I_m)(u_k\otimes I_m),\hfill \\ \hfill \mathbb{E}\left[v_k^{*}{v}_k^{*T}\right]& ={\Sigma }_{v^{*}}=R_{v^{*}}.\hfill \end{array}\hfill \end{array}\right.$$

(63)

Finally:

$$P_{v{v}_k}\approx R_{v^{*}}+(u_k^T\otimes I_m)\left(F_x^T\otimes I_m\right){\Sigma }_{\left(HH\right)}(F_x\otimes I_m)(u_k\otimes I_m),$$

(64)

and:

$$R_k^{th}\approx P_{v{v}_k}-t_k{t}_k^T=P_{v{v}_k}-F_t{u}_k{u}_k^T{F}_t^T.$$

(65)

3.2.3. Calculation of M

$$\begin{array}{cc}\hfill P_{w{v}_k}& =\mathbb{E}\left[w_k{v}_k^T\right]=\mathbb{E}\left[\delta A{x}_k{x}_k^T\delta H^T\right]+\mathbb{E}\left[\delta B{u}_k{x}_k^T\delta H^T\right]\hfill \\ \hfill & \approx \mathbb{E}\left[\delta A{\overline{x}}_k^{nom}{\overline{x}}_k^{no{m}^T}\delta H^T\right]+\mathbb{E}\left[\delta B{u}_k{\overline{x}}_k^T\delta H^T\right],\hfill \end{array}$$

(66)

as the expectations $\mathbb{E}\left[\delta A{x}_k{v}_k^{*T}\right]$ , $\mathbb{E}\left[\delta B{u}_k{v}_k^{*T}\delta H^T\right]$, $\mathbb{E}\left[B\delta u_k{x}_k^T\delta H^T\right]$ and $\mathbb{E}\left[B\delta u_k{v}_k^{*T}\right]$ are all (strictly) zero.

The non zero contributions are:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \mathbb{E}\left[\delta A{\overline{x}}_k{\overline{x}}_k^T\delta H^T\right]& \approx \mathbb{E}\left[\delta A{F}_x{u}_k{u}_k^T{F}_x^T\delta H^T\right]\approx (u_k^T\otimes I_n)(F_x^T\otimes I_n){\Sigma }_{\left(AH\right)}(F_x\otimes I_m)(u_k\otimes I_m),\hfill \\ \hfill \mathbb{E}\left[\delta B{u}_k{\overline{x}}_k^T\delta H^T\right]& \approx \mathbb{E}\left[\delta B{u}_k{u}_k^T{F}_x^T\delta H^T\right]\approx (u_k^T\otimes I_n){\Sigma }_{\left(BH\right)}(F_x\otimes I_m)(u_k\otimes I_m).\hfill \end{array}\hfill \end{array}\right.$$

(67)

Finally:

$$P_{w{v}_k}\approx (u_k^T\otimes I_n)\left[(F_x^T\otimes I_n){\Sigma }_{\left(AH\right)}(F_x\otimes I_m)+{\Sigma }_{\left(BH\right)}(F_x\otimes I_m)\right](u_k\otimes I_m),$$

(68)

and:

$$M_k\approx P_{w{v}_k}-z_k{t}_k^T=P_{w{v}_k}-F_z{u}_k{u}_k^T{F}_t^T.$$

(69)

3.3. Time Invariant Covariances

As derived in the previous sections Q, R and M have components proportional to $u_k$ (or conversely proportional to $U_{k-1}^{*}$) as such they are time-varying matrices. A time-varying Kalman filter could be straightforwardly implemented if desired and computationally compatible with real time implementation as briefly discussed in Section 4. However, the purpose of this work is to design a steady-state Kalman filter; in order to do that, several choices are possible with different level of robustness.

The degree of robustness can, as an example, be tuned by properly choosing the value of the nominal input u.

For steady-state filter design, time-invariant covariance matrices Q, R and M are sought. Since these depend on $u_k$ (eqs. (60), (64), (68)), a representative constant value must be selected. Two common approaches:

☐

Worst-case robustness: Set $u=u_{max}$ (componentwise) to bound uncertainties for all operating conditions,

☐

Average-case design: For periodic inputs, use $u=u_{rms}$ to match typical operating conditions.

The choice represents a trade-off: larger u increases Q, R, M, making the filter more conservative (higher estimation uncertainty) but more robust to parameter variations. The designer should select u based on the application’s robustness requirements and typical operating regime.

It can happen that, due to the different approximations introduced, $Q_k^{th}$ in eq. (61) is not positive semi-definite and/or $R_k^{th}$ in eq. (65) is not positive definite. If that is the case, one needs to regularize them to numerically solve the DARE equation as follows:

$$\begin{array}{cc}\hfill Q_k& =Q_k^{th}+{ϵ}_Q{I}_n,\mathrm{where}{ϵ}_Q=max\left(0,|{\lambda }_{-}^{max}\left(Q_k^{th}\right)|\right),\hfill \\ \hfill R_k& =R_k^{th}+{ϵ}_R{I}_m,\mathrm{where}{ϵ}_R=max\left(ϵ,|{\lambda }_{-}^{max}\left(R_k^{th}\right)|\right).\hfill \end{array}$$

(70)

and $ϵ$ is a sufficiently small positive number.

3.4. Filter Implementation

It is useful to recall the following property of the a priori and a posteriori estimations:

$$\left\{\begin{array}{c}{\widehat{x}}_k^{+}=\mathbb{E}\left[x_k\right|y_1\dots y_k],\hfill \\ {\widehat{x}}_k^{-}=\mathbb{E}\left[x_k\right|y_1\dots y_{k-1}].\hfill \end{array}\right.$$

(71)

It is important to note that for the "Predict Phase" $x_k$ can be approximated by ${\widehat{x}}_k^{+}$ since in this phase the measurement $y_k$ is available; whereas for the "Update/Correction Phase" $x_k$ can be approximated by ${\widehat{x}}_k^{-}$ as only measurements up to $y_{k-1}$ are available. Furthermore, ${\widehat{z}}_k$ and ${\widehat{t}}_k$ are the estimations of $z_k$ and $t_k$ respectively.

3.4.1. Predict Phase

$$\left\{\begin{array}{cc}{\widehat{x}}_{k+1}^{-}\hfill & =A{\widehat{x}}_k^{+}+B{u}_k+{\widehat{z}}_k,\hfill \\ {\widehat{Z}}_{k+1}\hfill & ={\widehat{Z}}_k\left(A^T\otimes I_n\right)+\left({\widehat{x}}_k^{{+}^T}\otimes I_n\right){\Sigma }_{\left(AA\right)}^{*}+\left(u_k^T\otimes I_n\right){\Sigma }_{\left(BA\right)}^{*},\hfill \\ {\widehat{T}}_{k+1}\hfill & ={\widehat{T}}_k\left(A^T\otimes I_n\right)+\left({\widehat{x}}_k^{{+}^T}\otimes I_m\right){\Sigma }_{\left(AH\right)}^{*}+\left(u_k^T\otimes I_m\right){\Sigma }_{\left(BH\right)}^{*}.\hfill \end{array}\right.$$

(72)

The initial values ${\widehat{Z}}_0$ and ${\widehat{T}}_0$ are both zero.

3.4.2. Update/Correction Phase

$$\left\{\begin{array}{cc}{\widehat{z}}_k\hfill & =Z_k\mathrm{vec}\left(I_n\right),\hfill \\ {\widehat{t}}_k\hfill & =T_k\mathrm{vec}\left(I_n\right),\hfill \\ {\widehat{x}}_k^{+}\hfill & ={\widehat{x}}_k^{-}+K\left[y_k-\left(H{\widehat{x}}_k^{-}+{\widehat{t}}_k\right)\right].\hfill \end{array}\right.$$

(73)

It can be noticed that the correction does not depend only on the (steady-state) Kalman gain K but also on the estimation of the output bias${\widehat{t}}_k$.

3.5. Computational Considerations and Basic MVFOSM Implementation

As it can be seen in eqs. (72) and (73), even if there is no need to update the state estimation covariance and therefore the Kalman gain, to fully profit of the proposed robust filter one needs to update in real time the estimations ${\widehat{z}}_k$ and ${\widehat{t}}_k$ and therefore, during the "Predict Phase", one needs to compute an additional $n\times n^2$ matrix to keep track of ${\widehat{Z}}_k$ and an additional $m\times n^2$ matrix to keep track of ${\widehat{T}}_k$ with respect to the canonical steady-state Kalman filter (the additional computational burden during the "Update/Correction Phase" is negligible 1 ).

A possible simplification can be realized by strictly following the MVFOSM framework i.e. performing:

☐

Zeroth-order approximation for means: $\mathbb{E}\left[x_k\right]={\overline{x}}_k\approx \mathbb{E}\left[x_k^{nom}\right]$ , therefore neglecting $O\left(\delta {\theta }^2\right)$ corrections

☐

Second-order approximation for covariances: retaining all terms up to $O\left(\delta {\theta }^2\right)$ in Q, R and M

This would mean accepting a biased state estimation by considering ${\widehat{z}}_k=0$ and ${\widehat{t}}_k=0$, in this case the resulting robust steady-state Kalman filter would reduce to the canonical Kalman filter described in eq. (2) with $K_k=K$, where K is still the solution of DARE in eq. (42) where one used $Q=P_{w{w}_k}$, $R=P_{v{v}_k}$ and $M=P_{w{v}_k}$; indeed the terms ${\widehat{z}}_k{\widehat{z}}_k^T$, ${\widehat{t}}_k{\widehat{t}}_k^T$ and ${\widehat{z}}_k{\widehat{t}}_k^T$ would be neglected too as they are $O\left(\delta {\theta }^4\right)$.

$$\left\{\begin{array}{c}\begin{array}{cccc}\hfill {\widehat{x}}_{k+1}^{-}& =A{\widehat{x}}_k^{+}+B{u}_k\hfill & \hfill & \mathrm{predict}\mathrm{phase}\hfill \\ \hfill {\widehat{x}}_k^{+}& ={\widehat{x}}_k^{-}+K(y_k-H{\widehat{x}}_k^{-})\hfill & \hfill & \mathrm{correction}/\mathrm{update}\mathrm{phase}\hfill \end{array}\hfill \end{array}\right.$$

(74)

4. Time-Varying Kalman Filter

If computationally feasible a robust time-varying filter can be implemented quite straightforwardly, as briefly sketched in the following. Furthermore, if the original uncertain system is also time-varying but having uncertainties with zero first order moment and constant second order moment, a time-varying filter is the only available option (where one could substitute $A_k$ for A, $B_k$ for B and $H_k$ for H).

4.1. Predict Phase

1.

Calculate $Q_k$ as follows:

$$\begin{array}{cc}\hfill Q_k& \approx Q_0+({\widehat{x}}_k^{{+}^T}\otimes I_n){\Sigma }_{\left(AB\right)}(u_k\otimes I_n)+(u_k^T\otimes I_n){\Sigma }_{\left(AB\right)}^T({\widehat{x}}_k^{+}\otimes I_n)\hfill \\ \hfill & +({\widehat{x}}_k^{{+}^T}\otimes I_n){\Sigma }_{\left(AA\right)}({\widehat{x}}_k^{+}\otimes I_n)+(u_k^T\otimes I_n){\Sigma }_{\left(BB\right)}(u_k\otimes I_n)-{\widehat{z}}_k{\widehat{z}}_k^T.\hfill \end{array}$$

(75)

2.

Perform the prediction step:

$$\left\{\begin{array}{cc}{\widehat{x}}_{k+1}^{-}\hfill & =A{\widehat{x}}_k^{+}+B{u}_k+{\widehat{z}}_k,\hfill \\ {\widehat{Z}}_{k+1}\hfill & ={\widehat{Z}}_k\left(A^T\otimes I_n\right)+\left({\widehat{x}}_k^{{+}^T}\otimes I_n\right){\Sigma }_{\left(AA\right)}^{*}+\left(u_k^T\otimes I_n\right){\Sigma }_{\left(BA\right)}^{*},\hfill \\ {\widehat{T}}_{k+1}\hfill & ={\widehat{T}}_k\left(A^T\otimes I_n\right)+\left({\widehat{x}}_k^{{+}^T}\otimes I_m\right){\Sigma }_{\left(AH\right)}^{*}+\left(u_k^T\otimes I_m\right){\Sigma }_{\left(BH\right)}^{*},\hfill \\ P_{k+1}^{-}\hfill & =A{P}_k^{+}{A}^T+Q_k.\hfill \end{array}\right.$$

(76)

4.2. Update/Correction Phase

1.

Calculate $R_k$ and $M_k$ as follows:

$$R_k\approx R_{v^{*}}+({\widehat{x}}_k^{{-}^T}\otimes I_m){\Sigma }_{\left(HH\right)}({\widehat{x}}_k^{-}\otimes I_m)-{\widehat{t}}_k{\widehat{t}}_k^T,$$

(77)

and

$$M_k\approx ({\widehat{x}}_k^{{-}^T}\otimes I_n){\Sigma }_{\left(AH\right)}({\widehat{x}}_k^{-}\otimes I_m)+(u_k^T\otimes I_n){\Sigma }_{\left(BH\right)}({\widehat{x}}_k^{-}\otimes I_m)-{\widehat{z}}_k{\widehat{t}}_k^T.$$

(78)

2.

Perform the update/correction step (based on the most general formulation of Kalman filter in [2]):

$$\left\{\begin{array}{cc}{S}_k\hfill & =H{P}_k^{-}{H}^T+H{M}_k+M_k^T{H}^T+R_k,\hfill \\ K_k\hfill & =\left(P_k^{-}{H}^T+M_k\right)S_k^{-1},\hfill \\ {\widehat{z}}_k\hfill & =Z_k\mathrm{vec}\left(I_n\right),\hfill \\ {\widehat{t}}_k\hfill & =T_k\mathrm{vec}\left(I_n\right),\hfill \\ {\widehat{x}}_k^{+}\hfill & ={\widehat{x}}_k^{-}+K_k\left[y_k-\left(H{\widehat{x}}_k^{-}+{\widehat{t}}_k\right)\right],\hfill \\ P_k^{+}\hfill & =\left(I_n-K_{k}H\right)P_k^{-}{\left(I_n-K_{k}H\right)}^T+K_k\left(H{M}_k+M_k^T{H}^T+R_k\right)K_k^T-K_k{M}_k^T-M_k{K}_k^T.\hfill \end{array}\right.$$

(79)

5. Conclusions

A practical and systematic methodology for robust steady-state Kalman filter design for uncertain LTI systems has been presented. Building on De Koning’s classical framework, the key contribution is the explicit derivation, within an MVFOSM-like approximation, of input-dependent (and state-dependent) equivalent noise covariance matrices $Q,R$ and M that systematically incorporate second-order statistics of parametric uncertainty alongside the usual process and measurement noise statistics. The resulting steady-state filter is obtained by solving a standard generalized DARE, requiring no iterative robust optimization.

A notable feature of the methodology is that filter robustness is governed by a single physically interpretable design parameter: the nominal input magnitude $\parallel u\parallel$ used to evaluate the input-dependent covariance contributions. Selecting $u=u_{max}$ yields a worst-case robust design; selecting $u=u_{\mathrm{rms}}$ targets average operating conditions. This transparent parametrization distinguishes the proposed approach from LMI- or regularization-based methods that require less intuitive tuning.

Two levels of approximation have been detailed. The simpler formulation (standard MVFOSM truncation, $z_k=0$ and $t_k=0$) yields a biased estimator but retains the canonical Kalman filter prediction-correction structure, making it directly deployable with minimal implementation overhead. The refined formulation tracks second-order bias corrections via an state, improving mean estimate accuracy at the cost of increased memory and computation. An extension to the time-varying case — which replaces the fixed nominal input with the actual time-varying $u_k$ and propagates the error covariance recursively — is also outlined, offering reduced conservatism where computational resources permit.

The methodology is especially relevant to embedded control in power electronics, where component tolerances are often significant, steady-state implementations are strongly preferred, and the physical interpretation of uncertainty parameters (component tolerance bounds) maps naturally onto the required second-order statistics via the GUM framework. Future work will focus on experimental validation in a power converter testbench.