(1) Time-update phase

Evaluate the cubature points based on the 3-SR rule:

$$\mathfrak{S}\_{k-1|k-1,j} = \mathfrak{S}\_{k-1|k-1} \gamma\_j^{\text{SR3}} + \mathfrak{x}\_{k-1|k-1}, 1 \le j \le 2n \tag{30}$$

where *<sup>S</sup>k*−1|*k*−<sup>1</sup> is the square root of the covariance matrix *<sup>P</sup>k*−1|*k*−<sup>1</sup> so that *<sup>P</sup>k*−1|*k*−<sup>1</sup> <sup>=</sup> *<sup>S</sup>k*−1|*k*−1*S*<sup>T</sup> *k*−1|*k*−1 and could be evaluated using the Cholesky decomposition; and *γ*SR3 *<sup>j</sup>* is formulated as follows:

$$\gamma\_{j}^{\text{SR3}} = \begin{cases} \begin{array}{c} \mathsf{a}\_{j} \\ -\mathsf{a}\_{j-n} \end{array} , \quad 1 \le j \le n \\\ -\mathsf{a}\_{j-n} \quad , \quad n+1 \le j \le 2n \end{cases} \tag{31}$$

$$\begin{array}{lcl} a\_j = \begin{bmatrix} a\_{j,1}, a\_{j,2}, \dots, a\_{j,n} \end{bmatrix}^T\\ a\_{j,i} = \begin{cases} \sqrt{9} & \text{ $i = j$ }\\ 0 & \text{ $i \neq j$ } \end{cases} \end{array} \tag{32}$$

Evaluate the predicted state and the covariance matrix:

$$\begin{aligned} \mathbf{\hat{x}}\_{k|k-1} &= \frac{1}{2\pi} \sum\_{j=1}^{2n} \left( F\_{\mathbf{y}\_j}^x + G\left(\mathbf{\hat{y}}\_j\right) u\_{k-1} + c \right) \\\ \mathbf{P}\_{k|k-1} &= \frac{1}{2\pi} \sum\_{j=1}^{2n} \left( F\_{\mathbf{y}\_j}^x + G\left(\mathbf{\hat{y}}\_j\right) u\_{k-1} + c - \mathbf{\hat{x}}\_{k|k-1} \right) \left( F\_{\mathbf{y}\_j}^x + G\left(\mathbf{\hat{y}}\_j\right) u\_{k-1} + c - \mathbf{\hat{x}}\_{k|k-1} \right)^T + Q \end{aligned} \tag{33}$$

Evaluate the information matrix and the information vector:

$$\begin{aligned} \Upsilon\_{k|k-1} &= \mathcal{P}\_{k|k-1}^{-1} \\ \hat{y}\_{k|k-1} &= \Upsilon\_{k|k-1} \hat{\mathfrak{x}}\_{k|k-1} \end{aligned} \tag{34}$$

(2) Measurement-update phase

Evaluate the cubature points based on the 5-SSR rule

$$\mathbf{f}\_{k|k-1,j}^{\rm m} = \mathbf{S}\_{1,k|k-1} \boldsymbol{\gamma}\_{j}^{\rm SNR} + \boldsymbol{\hat{x}}\_{1,k|k-1} \tag{35}$$

where *γ*SSR5 *<sup>j</sup>* is formulated as follows:

$$\gamma\_{j}^{\text{SSRS}} = \begin{cases} \begin{array}{c} \left[ \begin{array}{c} 0 \ \cdots \ \cdot & 0 \end{array} \right]^{\text{T}} & \text{ } j = 0\\ & \sqrt{n\_{1} + 2c} \ddots \end{array} & \begin{array}{c} j = 0\\ -\sqrt{n\_{1} + 2c} \cdot & \text{ } j = 1, 2, \cdots, n\_{1} + 1\\ -\sqrt{n\_{1} + 2c} \cdot & \text{ } j = n\_{1} + 2, \cdots, 2 \left( n\_{1} + 1 \right) \end{array} & \text{(36)}\\ \begin{array}{c} \sqrt{n\_{1} + 2b} \cdot & \text{ } j = 2n\_{1} + 3, \cdots, \left( n\_{1}^{2} + 5n\_{1} + 4 \right) / 2\\ -\sqrt{n\_{1} + 2b} \cdot & \text{ } j = \left( n\_{1}^{2} + 5n\_{1} + 6 \right) / 2, \cdots, n\_{1}^{2} + 3n\_{1} + 2 \end{array} & \text{(37)} \end{array}$$

where % *cj* & are the vertexes of the *n*1-simplex, and

$$\mathbf{c}\_{j,i} = \begin{cases} \mathbf{c}\_{j,1}, \mathbf{c}\_{j,2}, \dots, \mathbf{c}\_{j,n\_1} \big|^{\mathbf{T}}\\ -\sqrt{\frac{n\_1 + 1}{n\_1(n\_1 - i + 2)(n\_1 - i + 1)}} & \text{, } \quad i < j\\ \sqrt{\frac{(n\_1 + 1)(n\_1 - j + 1)}{n\_1(n\_1 - j + 2)}} & \text{, } \quad i = j \quad 1 \le i \le n\_1, 1 \le j \le n\_1 + 1\\ 0 & \text{, } \quad i > j \end{cases} \tag{37}$$

{*bi*} are the projections of the midpoints of the edges constructed by % *cj* & .

$$\{b\_l\} = \left\{ \sqrt{\frac{n\_1}{2(n\_1 - 1)}} (\mathbf{c}\_l + \mathbf{c}\_m) : l < m, 1 \le l \le n\_1 + 1, 1 \le m \le n\_1 + 1 \right\} \tag{38}$$

The weights for the 5-SSR rule are as follows:

$$w\_j = \begin{cases} \frac{2}{n\_1 + 2} & j = 0\\ \frac{\left(7 - n\_1\right)n\_1^2}{2\left(n\_1 + 1\right)^2 \left(n\_1 + 2\right)^2} & j = 1, \dots, 2\left(n\_1 + 1\right)\\ \frac{2\left(n\_1 - 1\right)^2}{\left(n\_1 + 1\right)^2 \left(n\_1 + 2\right)^2} & j = 2n\_1 + 3, \dots, n\_1^2 + 3n\_1 + 2 \end{cases} \tag{39}$$

Evaluate the predicted measurement ˆ*zk*|*k*−<sup>1</sup> and the cross-covariance matrix *Pxz*:

$$\begin{aligned} \hat{\mathbf{z}}\_{k|k-1} &= \sum\_{j=1}^{n\_1^2 + 3n\_1 + 3} \mathbf{w}\_j h\left(\mathbf{f}\_{k|k-1,j}^{\text{m}}\right) \\ \mathbf{P}\_{\mathbf{X}\mathbf{z}} &= \begin{bmatrix} n\_1^2 + 3n\_1 + 3 \\ \sum\_{j=1}^{n\_1^2 + 3n\_1 + 3} \mathbf{w}\_j \mathbf{f}\_{k|k-1,j}^{\text{m}} h\left(\mathbf{f}\_{k|k-1,j}^{\text{m}}\right)^\text{T} - \hat{\mathbf{x}}\_{1,k|k-1} \hat{\mathbf{z}}\_{k|k-1}^{\text{T}} \\ \mathbf{P}\_{21,k|k-1} \mathbf{P}\_{11,k|k-1}^{-1} \left( \sum\_{j=1}^{n\_1^2 + 3n\_1 + 3} \mathbf{w}\_j \mathbf{f}\_{k|k-1,j}^{\text{m}} h\left(\mathbf{f}\_{k|k-1,j}^{\text{m}}\right)^\text{T} - \hat{\mathbf{x}}\_{1,k|k-1} \hat{\mathbf{z}}\_{k|k-1}^{\text{T}} \right) \end{bmatrix} \end{aligned} \tag{40}$$

As the dimension of the measurement vector is generally large, the inversion of the matrix in Equation (17) is very time-consuming. To omit the inversion, we assumed that the *m* measurements were independent of each other, and the covariance matrix was a diagonal matrix with diagonal elements of *Rj*,*k*(*j* = 1, . . . , *m*). Further assuming that *Rj*,*<sup>k</sup>* was invariant with time and equal to the other diagonal elements, *Rj*,*<sup>k</sup>* was treated as a constant *R*. Thus, the inverse of the *m* × *m* matrix *R* was replaced by the division of the scalar *R*, and the measurement update of the information vector and matrix was obtained as follows.

$$\begin{aligned} \mathfrak{H}\_{k|k} &= \mathfrak{H}\_{k|k-1} + R^{-1} H\_k^\mathrm{T} \left( \mathfrak{z}\_k - h \left( \mathfrak{x}\_{1,k|k-1} \right) + H\_k \mathfrak{x}\_{k|k-1} \right) \\ \mathfrak{Y}\_{k|k} &= \mathfrak{Y}\_{k|k-1} + R^{-1} H\_k^\mathrm{T} H\_k \end{aligned} \tag{41}$$

Finally, the estimated state vector and covariance matrix were recovered based on Equation (18).
