**Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization** Provisional chapter

DOI: 10.5772/intechopen.72465

Efficient Matrix-Free Ensemble Kalman Filter

Elias David Niño Ruiz, Rolando Beltrán Arrieta and Implementations: Accounting for Localization

Alfonso Manuel Mancilla Herrera

Additional information is available at the end of the chapter and Alfonso Manuel Mancilla Herrera

Elias David Niño Ruiz, Rolando Beltrán Arrieta

http://dx.doi.org/10.5772/intechopen.72465 Additional information is available at the end of the chapter

#### Abstract

Signals and Image Processing (IWSSIP), 2016; May 23-25, 2016; Bratislava, Slovakia.

[15] Muller D, Pauli J, Nunn C, Gormer S, Muller-Schneiders S. Time to contact estimation using interest points. In: 12th International Conference on Transportation Systems, 2009. ITSC '09; October 4-7, 2009; St. Louis, MO, USA. IEEE; 2009. pp. 1-6. DOI: 10.1109/ITSC.2009.5309851

[16] Bradski G, Kaehler A. Learning OpenCV. 1st ed. USA: O'Reilly Media, Inc.; 2008. p. 555 [17] Prince S. Computer Vision: Models, Learning and Inference. USA: Cambridge University

[18] Kalman REA. New approach to linear filtering and prediction problems. Transactions of the ASME – Journal of Basic Engineering. 1960;**82**(1):35-45. DOI: 10.1115/1.3662552 [19] Maybeck P. Introduction. In: Stochastic Models, Estimation and Control. 1st ed. London:

IEEE; 2016. pp. 1-6. DOI: 10.1109/IWSSIP.2016.7502702

Press; 2012. p. 580

164 Kalman Filters - Theory for Advanced Applications

Academic Press; 1979. p. 1-16

This chapter discusses efficient and practical matrix-free implementations of the ensemble Kalman filter (EnKF) in order to account for localization during the assimilation of observations. In the EnKF context, an ensemble of model realizations is utilized in order to estimate the moments of its underlying error distribution. Since ensemble members come at high computational costs (owing to current operational model resolutions) ensemble sizes are constrained by the hundreds while, typically, their error distributions range in the order of millions. This induces spurious correlations in estimates of prior error correlations when these are approximated via the ensemble covariance matrix. Localization methods are commonly utilized in order to counteract this effect. EnKF implementations in this context are based on a modified Cholesky decomposition. Different flavours of Choleskybased filters are discussed in this chapter. Furthermore, the computational effort in all formulations is linear with regard to model resolutions. Experimental tests are performed making use of the Lorenz 96 model. The results reveal that, in terms of root-mean-squareerrors, all formulations perform equivalently.

Keywords: ensemble Kalman filter, modified Cholesky decomposition, sampling methods

### 1. Introduction

The ensemble Kalman filter (EnKF) is a sequential Monte Carlo method for parameter and state estimation in highly nonlinear models. The popularity of the EnKF owes to its simple

© The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited.

formulation and relatively easy implementation. In the EnKF, an ensemble of model realizations is utilized in order to, among other things, estimate the moments of the prior error distribution (forecast distribution). During this estimation and the assimilation of observations, many challenges are present in practice. Some of them are addressed in this chapter: the first one is related to the proper estimation of background error correlations, which has a high impact on the quality of the analysis corrections while the last one is related to the proper estimation of the posterior ensemble when the observation operator is nonlinear. In all cases, background error correlations are estimated based on the Bickel and Levina estimator [1].

<sup>B</sup> <sup>≈</sup>P<sup>b</sup> <sup>¼</sup> <sup>1</sup>

<sup>X</sup><sup>a</sup> <sup>¼</sup> <sup>A</sup> � <sup>P</sup><sup>b</sup> � ��<sup>1</sup>

<sup>X</sup><sup>a</sup> <sup>¼</sup> <sup>X</sup><sup>b</sup> <sup>þ</sup> Pb � <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> <sup>þ</sup> <sup>H</sup> � <sup>P</sup><sup>b</sup> � <sup>H</sup><sup>T</sup> � ��<sup>1</sup>

matrix of innovations (on observations) is denoted by <sup>Δ</sup><sup>Y</sup> <sup>¼</sup> <sup>Y</sup><sup>s</sup> � <sup>H</sup> � <sup>X</sup><sup>b</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�N, and the matrix

formed by samples from a normal distribution with moments N ð Þ 0; R . In practice, one of the

The forms of the Cholesky and the modified Cholesky decomposition were discussed by Golub and Van Loan in [5]. If A ∈ R<sup>n</sup>�<sup>n</sup> is symmetric positive definite, then there exist a unique lower triangular <sup>L</sup><sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>n</sup> with positive diagonal entries such that <sup>A</sup> <sup>¼</sup> <sup>L</sup> � <sup>L</sup><sup>T</sup>. If all the leading principal submatrices of A ∈ R<sup>n</sup>�<sup>n</sup> are non-singular, then there exist unique lower triangular matrices L and <sup>M</sup> and a unique diagonal matrix <sup>D</sup> <sup>¼</sup> diagð Þ <sup>d</sup>1; <sup>d</sup>2;…; dn such that <sup>A</sup> <sup>¼</sup> <sup>L</sup> � <sup>D</sup> � <sup>M</sup><sup>T</sup> : If <sup>A</sup> is symmetric, then <sup>L</sup> <sup>¼</sup> <sup>M</sup> and <sup>A</sup> <sup>¼</sup> <sup>L</sup> � <sup>D</sup> � LT . Golub and Van Loan provide proofs of those

where the analysis covariance matrix <sup>A</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>n</sup> is given by <sup>A</sup> <sup>¼</sup> Pb � ��<sup>1</sup>

matrix of member deviations ΔX<sup>b</sup> ∈ R<sup>n</sup>�<sup>N</sup> reads,

tions are done over both error distributions.

of perturbed observations reads <sup>Y</sup><sup>s</sup> <sup>¼</sup> <sup>y</sup> � <sup>1</sup><sup>T</sup>

2.2. Cholesky and modified Cholesky decomposition

most utilized formulations is given by (8).

of the next formulations,

and

<sup>N</sup> � <sup>1</sup> � <sup>Δ</sup>X<sup>b</sup> � <sup>Δ</sup>X<sup>b</sup> � �<sup>T</sup>

Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization

where N is the ensemble size, xb i½ � ∈ R<sup>n</sup>�<sup>1</sup> is the ith ensemble member, for 1 ≤ i ≤ N, x<sup>b</sup> ∈ R<sup>n</sup>�<sup>1</sup> is well known as the background state, B∈ R<sup>n</sup>�<sup>n</sup> stands for the background error covariance matrix, xb is the ensemble mean, and Pb is the ensemble covariance matrix. Likewise, the

A variety of EnKF implementations have been proposed in the current literature, most of them rely on basic assumptions over the moments of two probability distributions: the background (prior) and the analysis (posterior) distributions. In most of EnKF formulations, normal assump-

In general, depending on how the assimilation process is performed, the EnKF implementation will fall in one of two categories: stochastic filter [3] or deterministic filter [4]. In the context of stochastic filters, for instance, the assimilation of observations can be performed by using any

<sup>Δ</sup>X<sup>b</sup> <sup>¼</sup> <sup>X</sup><sup>b</sup> � xb � <sup>1</sup><sup>T</sup>

∈ R<sup>n</sup>�n, (4)

http://dx.doi.org/10.5772/intechopen.72465

167

<sup>N</sup> ∈ R<sup>n</sup>�<sup>N</sup>: (5)

� <sup>Δ</sup><sup>Y</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�N, (8)

<sup>þ</sup> <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> � <sup>H</sup> h i�<sup>1</sup>

, the

� <sup>X</sup><sup>b</sup> <sup>þ</sup> <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> � <sup>Y</sup><sup>s</sup> h i<sup>∈</sup> Rn�N, (6)

<sup>N</sup> <sup>þ</sup> <sup>E</sup><sup>∈</sup> <sup>R</sup><sup>m</sup>�<sup>N</sup>, where the columns of <sup>E</sup><sup>∈</sup> <sup>R</sup><sup>m</sup>�<sup>N</sup> are

<sup>X</sup><sup>a</sup> <sup>¼</sup> <sup>X</sup><sup>b</sup> <sup>þ</sup> <sup>A</sup> � <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> � <sup>Δ</sup><sup>Y</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�N, (7)

### 2. Preliminaries

We want to estimate the state of a dynamical system x<sup>∗</sup> ∈ R<sup>n</sup>�<sup>1</sup> , which (approximately) evolves according to some (imperfect) numerical model operators x<sup>∗</sup> <sup>k</sup> <sup>¼</sup> <sup>M</sup>tk�1!t<sup>k</sup> <sup>x</sup><sup>∗</sup> k�1 � �, for instance, a model which mimics the behavior of the atmosphere and/or the ocean where n is the model dimension (typically associated with the model resolution). The estimation process is based on two sources of information: a prior estimate of <sup>x</sup><sup>∗</sup>, <sup>x</sup><sup>b</sup> <sup>¼</sup> <sup>x</sup><sup>∗</sup> <sup>þ</sup> <sup>ξ</sup> with <sup>ξ</sup> � <sup>N</sup> ð Þ <sup>0</sup>n;<sup>B</sup> , and a noisy observation <sup>y</sup> <sup>¼</sup> <sup>H</sup> <sup>x</sup><sup>∗</sup> ð Þþ <sup>e</sup> with <sup>e</sup> � <sup>N</sup> ð Þ <sup>0</sup>m; <sup>R</sup> , where <sup>m</sup> is the number of observed model components, <sup>R</sup><sup>m</sup>�<sup>m</sup> is the estimated data error covariance matrix, <sup>H</sup> : <sup>R</sup><sup>n</sup>�<sup>1</sup> ! <sup>R</sup><sup>m</sup>�<sup>1</sup> is the (nonlinear) observation operator that maps the information from the model domain to the observation space, 0<sup>n</sup> is the nth dimensional vector whose components are all zeros. In practice, m ≪ n or m < n. The assimilation process can be performed by using sequential data assimilation methods, and one of the most widely used methods is the ensemble Kalman filter, which is detailed in the following paragraphs.

#### 2.1. The ensemble Kalman filter

In the ensemble Kalman filter (EnKF) [2], an ensemble of model realizations

$$\mathbf{x}^{\mathbf{b}} = \left[ \mathbf{x}^{b[\mathbf{1}]}, \mathbf{x}^{b[\mathbf{2}]}, \dots, \mathbf{x}^{b[\mathbf{N}]} \right] \in \mathbb{R}^{n \times N},\tag{1}$$

is utilized in order to estimate the moments of the background error distribution,

$$\mathbf{x}^{b[i]} \sim \mathcal{N}(\mathbf{x}^b, \mathbf{B}),\tag{2}$$

via the empirical moments of the ensemble (1) and therefore,

$$\mathfrak{x}^b \approx \mathfrak{x}^b = \frac{1}{N} \cdot \sum\_{i=1}^{N} \mathfrak{x}^{b[i]} \in \mathbb{R}^{n \times 1},\tag{3}$$

and

Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization http://dx.doi.org/10.5772/intechopen.72465 167

$$\mathbf{B} \approx \mathbf{P}^b = \frac{1}{N-1} \cdot \Delta \mathbf{X}^b \cdot \left[\Delta \mathbf{X}^b\right]^T \in \mathbb{R}^{n \times n},\tag{4}$$

where N is the ensemble size, xb i½ � ∈ R<sup>n</sup>�<sup>1</sup> is the ith ensemble member, for 1 ≤ i ≤ N, x<sup>b</sup> ∈ R<sup>n</sup>�<sup>1</sup> is well known as the background state, B∈ R<sup>n</sup>�<sup>n</sup> stands for the background error covariance matrix, xb is the ensemble mean, and Pb is the ensemble covariance matrix. Likewise, the matrix of member deviations ΔX<sup>b</sup> ∈ R<sup>n</sup>�<sup>N</sup> reads,

$$
\Delta \mathbf{X}^b = \mathbf{X}^b - \overline{\mathbf{x}}^b \cdot \mathbf{1}\_N^T \in \mathbb{R}^{n \times N}. \tag{5}
$$

A variety of EnKF implementations have been proposed in the current literature, most of them rely on basic assumptions over the moments of two probability distributions: the background (prior) and the analysis (posterior) distributions. In most of EnKF formulations, normal assumptions are done over both error distributions.

In general, depending on how the assimilation process is performed, the EnKF implementation will fall in one of two categories: stochastic filter [3] or deterministic filter [4]. In the context of stochastic filters, for instance, the assimilation of observations can be performed by using any of the next formulations,

$$\mathbf{X}^{a} = \mathbf{A} \cdot \left[ \left[ \mathbf{P}^{b} \right]^{-1} \cdot \mathbf{X}^{b} + \mathbf{H}^{T} \cdot \mathbf{R}^{-1} \cdot \mathbf{Y}^{s} \right] \in \mathbf{R}^{n \times N},\tag{6}$$

$$\mathbf{X}^{a} = \mathbf{X}^{b} + \mathbf{A} \cdot \mathbf{H}^{T} \cdot \mathbf{R}^{-1} \cdot \boldsymbol{\Delta} \mathbf{Y} \in \mathbb{R}^{n \times N},\tag{7}$$

and

formulation and relatively easy implementation. In the EnKF, an ensemble of model realizations is utilized in order to, among other things, estimate the moments of the prior error distribution (forecast distribution). During this estimation and the assimilation of observations, many challenges are present in practice. Some of them are addressed in this chapter: the first one is related to the proper estimation of background error correlations, which has a high impact on the quality of the analysis corrections while the last one is related to the proper estimation of the posterior ensemble when the observation operator is nonlinear. In all cases, background error correlations are estimated based on the Bickel and Levina estimator [1].

model which mimics the behavior of the atmosphere and/or the ocean where n is the model dimension (typically associated with the model resolution). The estimation process is based on two sources of information: a prior estimate of <sup>x</sup><sup>∗</sup>, <sup>x</sup><sup>b</sup> <sup>¼</sup> <sup>x</sup><sup>∗</sup> <sup>þ</sup> <sup>ξ</sup> with <sup>ξ</sup> � <sup>N</sup> ð Þ <sup>0</sup>n;<sup>B</sup> , and a noisy observation <sup>y</sup> <sup>¼</sup> <sup>H</sup> <sup>x</sup><sup>∗</sup> ð Þþ <sup>e</sup> with <sup>e</sup> � <sup>N</sup> ð Þ <sup>0</sup>m; <sup>R</sup> , where <sup>m</sup> is the number of observed model components, <sup>R</sup><sup>m</sup>�<sup>m</sup> is the estimated data error covariance matrix, <sup>H</sup> : <sup>R</sup><sup>n</sup>�<sup>1</sup> ! <sup>R</sup><sup>m</sup>�<sup>1</sup> is the (nonlinear) observation operator that maps the information from the model domain to the observation space, 0<sup>n</sup> is the nth dimensional vector whose components are all zeros. In practice, m ≪ n or m < n. The assimilation process can be performed by using sequential data assimilation methods, and one of the most widely used methods is the ensemble Kalman filter,

, which (approximately) evolves

� �, for instance, a

k�1

∈ R<sup>n</sup>�N, (1)

, (3)

<sup>x</sup>b i½ � � <sup>N</sup> <sup>x</sup><sup>b</sup>;<sup>B</sup> � �, (2)

<sup>k</sup> <sup>¼</sup> <sup>M</sup>tk�1!t<sup>k</sup> <sup>x</sup><sup>∗</sup>

2. Preliminaries

166 Kalman Filters - Theory for Advanced Applications

We want to estimate the state of a dynamical system x<sup>∗</sup> ∈ R<sup>n</sup>�<sup>1</sup>

In the ensemble Kalman filter (EnKF) [2], an ensemble of model realizations

<sup>X</sup><sup>b</sup> <sup>¼</sup> <sup>x</sup><sup>b</sup>½ � <sup>1</sup> ; xb½ � <sup>2</sup> ;…; <sup>x</sup>b N½ � h i

is utilized in order to estimate the moments of the background error distribution,

xb <sup>≈</sup> xb <sup>¼</sup> <sup>1</sup>

N � X N

i¼1

xb i½ � ∈ R<sup>n</sup>�<sup>1</sup>

via the empirical moments of the ensemble (1) and therefore,

according to some (imperfect) numerical model operators x<sup>∗</sup>

which is detailed in the following paragraphs.

2.1. The ensemble Kalman filter

and

$$\mathbf{X}^{a} = \mathbf{X}^{b} + \mathbf{P}^{b} \cdot \mathbf{H}^{T} \cdot \left[\mathbf{R}^{-1} + \mathbf{H} \cdot \mathbf{P}^{b} \cdot \mathbf{H}^{T}\right]^{-1} \cdot \Delta \mathbf{Y} \in \mathbb{R}^{n \times N},\tag{8}$$

where the analysis covariance matrix <sup>A</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>n</sup> is given by <sup>A</sup> <sup>¼</sup> Pb � ��<sup>1</sup> <sup>þ</sup> <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> � <sup>H</sup> h i�<sup>1</sup> , the matrix of innovations (on observations) is denoted by <sup>Δ</sup><sup>Y</sup> <sup>¼</sup> <sup>Y</sup><sup>s</sup> � <sup>H</sup> � <sup>X</sup><sup>b</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�N, and the matrix of perturbed observations reads <sup>Y</sup><sup>s</sup> <sup>¼</sup> <sup>y</sup> � <sup>1</sup><sup>T</sup> <sup>N</sup> <sup>þ</sup> <sup>E</sup><sup>∈</sup> <sup>R</sup><sup>m</sup>�<sup>N</sup>, where the columns of <sup>E</sup><sup>∈</sup> <sup>R</sup><sup>m</sup>�<sup>N</sup> are formed by samples from a normal distribution with moments N ð Þ 0; R . In practice, one of the most utilized formulations is given by (8).

#### 2.2. Cholesky and modified Cholesky decomposition

The forms of the Cholesky and the modified Cholesky decomposition were discussed by Golub and Van Loan in [5]. If A ∈ R<sup>n</sup>�<sup>n</sup> is symmetric positive definite, then there exist a unique lower triangular <sup>L</sup><sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>n</sup> with positive diagonal entries such that <sup>A</sup> <sup>¼</sup> <sup>L</sup> � <sup>L</sup><sup>T</sup>. If all the leading principal submatrices of A ∈ R<sup>n</sup>�<sup>n</sup> are non-singular, then there exist unique lower triangular matrices L and <sup>M</sup> and a unique diagonal matrix <sup>D</sup> <sup>¼</sup> diagð Þ <sup>d</sup>1; <sup>d</sup>2;…; dn such that <sup>A</sup> <sup>¼</sup> <sup>L</sup> � <sup>D</sup> � <sup>M</sup><sup>T</sup> : If <sup>A</sup> is symmetric, then <sup>L</sup> <sup>¼</sup> <sup>M</sup> and <sup>A</sup> <sup>¼</sup> <sup>L</sup> � <sup>D</sup> � LT . Golub and Van Loan provide proofs of those decompositions, as well as various ways to compute them. We can use the Cholesky factor of a covariance matrix to quickly solve linear systems that involve a covariance matrix. Note that if we compute the Cholesky factorization and solve the triangular system <sup>L</sup> � <sup>y</sup> <sup>¼</sup> <sup>b</sup> and <sup>L</sup><sup>T</sup> � <sup>x</sup> <sup>¼</sup> <sup>y</sup>, then

$$\mathbf{A} \cdot \mathbf{x} = \left(\mathbf{L} \cdot \mathbf{L}^T\right) \cdot \mathbf{x} = \mathbf{L} \cdot \left(\mathbf{L}^T \cdot \mathbf{x}\right) = \mathbf{L} \cdot \mathbf{y} = \mathbf{b}.\tag{9}$$

or equivalently

estimation of <sup>B</sup><sup>b</sup> �<sup>1</sup>

<sup>B</sup><sup>≈</sup> <sup>B</sup><sup>b</sup> <sup>¼</sup> <sup>L</sup>�<sup>1</sup> � <sup>D</sup> � <sup>L</sup>�<sup>T</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�n, (14)

http://dx.doi.org/10.5772/intechopen.72465

169

<sup>A</sup>�<sup>1</sup> <sup>¼</sup> <sup>B</sup>�<sup>1</sup> <sup>þ</sup> <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> � <sup>H</sup> � <sup>R</sup><sup>n</sup>�n, (15)

<sup>x</sup> � <sup>N</sup> xa ð Þ ; <sup>A</sup> , (16)

<sup>þ</sup> <sup>X</sup> � <sup>X</sup><sup>T</sup> , (17)

, sparse estimators of the precision background can be obtained. Further-

Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization

where L∈ R<sup>n</sup>�<sup>n</sup> is a lower-triangular matrix whose diagonal elements are all ones, and D ∈ R<sup>n</sup>�<sup>n</sup> is a diagonal matrix. Even more, when only local effects are considered during the

more, the matrix L can be sparse with only a few non-zero elements per row. Typically, the number of non-zero elements depends on the radius of influence during the estimation of background error correlations. For instance, in the one-dimensional case, the radius of influence denotes the maximum number of non-zero elements, per row, in L. The EnKF-MC is then obtained by plugging in the estimator (3) in (6). Given the structure of the Cholesky factors, the

moreover, since <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> � <sup>H</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>n</sup> can be written as a sum of <sup>m</sup> rank-one matrices, the factors (3) can be updated in order to obtain an estimate of the inverse analysis covariance matrix. The

In the stochastic posterior ensemble Kalman filter (PEnKF-S) [7, 8], we want to estimate the

based on the background ensemble X<sup>b</sup> � �, where xa is the analysis state and A ∈ R<sup>n</sup>�<sup>n</sup> is the analysis covariance matrix. Consider the estimate of the inverse background error covariance matrix Eq. (3), the precision analysis covariance matrix Eq. (4) can be approximated as follows:

<sup>¼</sup> <sup>B</sup><sup>b</sup> �<sup>1</sup>

<sup>2</sup> ∈ R<sup>n</sup>�<sup>m</sup>. The matrix (17) can be written as follows:

<sup>¼</sup> <sup>L</sup><sup>T</sup> � <sup>D</sup>�<sup>1</sup> � <sup>L</sup> <sup>þ</sup>X<sup>m</sup>

where xi denotes the ith column of matrix X, for 1 ≤ i ≤ m. Consider the sequence of factor

i¼1

xi � xT i ,

<sup>A</sup>�<sup>1</sup> <sup>≈</sup> <sup>A</sup><sup>b</sup> �<sup>1</sup>

Ab �1

EnKF-MC can be seen as a matrix-free implementation of the EnKF.

Recall that the precision analysis covariance matrix reads,

next section discussed such approximation in detail.

moments of the analysis distribution,

where <sup>X</sup> <sup>¼</sup> <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup>

updates,

3.1.1. Posterior ensemble Kalman filter stochastic (PEnKF-S)

For a symmetric matrix A ∈ R<sup>n</sup>�<sup>n</sup>, the modified Cholesky decomposition involves finding a non-negative diagonal matrix E such that A þ E is positive definite. In particular, E ¼ 0 whenever A is positive definite, otherwise k kE should be sufficiently small. This allows applying the usual Cholesky factorization to <sup>A</sup> <sup>þ</sup> <sup>E</sup>, i.e. finding matrices <sup>L</sup>, <sup>D</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>n</sup> such that <sup>A</sup> <sup>þ</sup> <sup>E</sup> <sup>¼</sup> <sup>L</sup> � <sup>D</sup> � LT .

Recall <sup>Δ</sup>X<sup>b</sup> <sup>¼</sup> <sup>X</sup><sup>b</sup> � xb <sup>⊗</sup> <sup>1</sup><sup>T</sup> <sup>N</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>N</sup>. Thus, <sup>Δ</sup>xb i½ � � <sup>N</sup> ð Þ <sup>0</sup>n;<sup>B</sup> for <sup>1</sup> <sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>n</sup>. Let xb i½ � <sup>∈</sup> <sup>R</sup><sup>N</sup>�<sup>1</sup> , the vector holding the ith row across all columns of ΔX<sup>b</sup> for 1 ≤ i ≤ n. Bickel and Levina in [1] discussed a modified Cholesky decomposition for the estimation of precision covariance matrix, and this allows in our context to estimate of B�<sup>1</sup> by fitting regressions of the form:

$$\mathbf{x}^{b[i]} = \sum\_{j=1}^{i-1} \mathbf{x}^{b[j]} \cdot \boldsymbol{\beta}\_{i,j} + \boldsymbol{\xi}^{[i]} \in \mathbb{R}^{N \times 1},\tag{10}$$

then, f gL ij ¼ �βi,j , for 1 ≤ i < j ≤ n.

#### 3. Filters based on modified Cholesky decomposition

#### 3.1. Ensemble Kalman filter based on modified Cholesky

In the ensemble Kalman filter based on a modified Cholesky decomposition (EnKF-MC) [6], the analysis ensemble is estimated by using the following equations,

$$\mathbf{X}^{a} = \mathbf{X}^{b} + \left[\mathbf{B}^{-1} + \mathbf{H}^{T}\mathbf{R}^{-1} \cdot \mathbf{H}\right]^{-1} \cdot \mathbf{H}^{T} \cdot \mathbf{R}^{-1} \cdot \boldsymbol{\Delta} \mathbf{Y},\tag{11}$$

and

$$
\Delta \mathbf{Y} = \left[ \mathbf{Y}^{\mathcal{S}} - \mathbf{H} \mathbf{X}^{\mathcal{b}} \right] \in \mathbb{R}^{n \times N}. \tag{12}
$$

In this context, error correlations are estimated via a modified Cholesky decomposition. This provides an estimate of the inverse background error covariance matrix of the form,

$$\mathbf{B}^{-1} \approx \widehat{\mathbf{B}}^{-1} = \mathbf{L}^{T} \cdot \mathbf{D}^{-1} \cdot \mathbf{L} \in \mathbb{R}^{n \times n},\tag{13}$$

or equivalently

decompositions, as well as various ways to compute them. We can use the Cholesky factor of a covariance matrix to quickly solve linear systems that involve a covariance matrix. Note that if we compute the Cholesky factorization and solve the triangular system <sup>L</sup> � <sup>y</sup> <sup>¼</sup> <sup>b</sup> and <sup>L</sup><sup>T</sup> � <sup>x</sup> <sup>¼</sup> <sup>y</sup>,

For a symmetric matrix A ∈ R<sup>n</sup>�<sup>n</sup>, the modified Cholesky decomposition involves finding a non-negative diagonal matrix E such that A þ E is positive definite. In particular, E ¼ 0 whenever A is positive definite, otherwise k kE should be sufficiently small. This allows applying the usual Cholesky factorization to <sup>A</sup> <sup>þ</sup> <sup>E</sup>, i.e. finding matrices <sup>L</sup>, <sup>D</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>n</sup> such that <sup>A</sup> <sup>þ</sup> <sup>E</sup>

vector holding the ith row across all columns of ΔX<sup>b</sup> for 1 ≤ i ≤ n. Bickel and Levina in [1] discussed a modified Cholesky decomposition for the estimation of precision covariance matrix, and this allows in our context to estimate of B�<sup>1</sup> by fitting regressions of the form:

In the ensemble Kalman filter based on a modified Cholesky decomposition (EnKF-MC) [6],

In this context, error correlations are estimated via a modified Cholesky decomposition. This

provides an estimate of the inverse background error covariance matrix of the form,

<sup>B</sup>�<sup>1</sup> <sup>≈</sup> <sup>B</sup><sup>b</sup> �<sup>1</sup>

xb j½ � � <sup>β</sup>i,j <sup>þ</sup> <sup>ξ</sup>½ �<sup>i</sup> <sup>∈</sup> <sup>R</sup><sup>N</sup>�<sup>1</sup>

xb i½ � <sup>¼</sup> <sup>X</sup> i�1

3. Filters based on modified Cholesky decomposition

the analysis ensemble is estimated by using the following equations,

<sup>X</sup><sup>a</sup> <sup>¼</sup> <sup>X</sup><sup>b</sup> <sup>þ</sup> <sup>B</sup>�<sup>1</sup> <sup>þ</sup> <sup>H</sup>TR�<sup>1</sup> � <sup>H</sup> � ��<sup>1</sup>

3.1. Ensemble Kalman filter based on modified Cholesky

, for 1 ≤ i < j ≤ n.

j¼1

<sup>A</sup> � <sup>x</sup> <sup>¼</sup> <sup>L</sup> � LT � � � <sup>x</sup> <sup>¼</sup> <sup>L</sup> � <sup>L</sup><sup>T</sup> � <sup>x</sup> � � <sup>¼</sup> <sup>L</sup> � <sup>y</sup> <sup>¼</sup> <sup>b</sup>: (9)

<sup>N</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>N</sup>. Thus, <sup>Δ</sup>xb i½ � � <sup>N</sup> ð Þ <sup>0</sup>n;<sup>B</sup> for <sup>1</sup> <sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>n</sup>. Let xb i½ � <sup>∈</sup> <sup>R</sup><sup>N</sup>�<sup>1</sup>

, the

, (10)

� <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> � <sup>Δ</sup>Y, (11)

<sup>Δ</sup><sup>Y</sup> <sup>¼</sup> <sup>Y</sup><sup>S</sup> � <sup>H</sup>:X<sup>b</sup> � � <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>N</sup>: (12)

<sup>¼</sup> <sup>L</sup><sup>T</sup> � <sup>D</sup>�<sup>1</sup> � <sup>L</sup><sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>n</sup>, (13)

then

<sup>¼</sup> <sup>L</sup> � <sup>D</sup> � LT .

then, f gL ij ¼ �βi,j

and

Recall <sup>Δ</sup>X<sup>b</sup> <sup>¼</sup> <sup>X</sup><sup>b</sup> � xb <sup>⊗</sup> <sup>1</sup><sup>T</sup>

168 Kalman Filters - Theory for Advanced Applications

$$\mathbf{B} \approx \widehat{\mathbf{B}} = \mathbf{L}^{-1} \cdot \mathbf{D} \cdot \mathbf{L}^{-T} \in \mathbb{R}^{n \times n},\tag{14}$$

where L∈ R<sup>n</sup>�<sup>n</sup> is a lower-triangular matrix whose diagonal elements are all ones, and D ∈ R<sup>n</sup>�<sup>n</sup> is a diagonal matrix. Even more, when only local effects are considered during the estimation of <sup>B</sup><sup>b</sup> �<sup>1</sup> , sparse estimators of the precision background can be obtained. Furthermore, the matrix L can be sparse with only a few non-zero elements per row. Typically, the number of non-zero elements depends on the radius of influence during the estimation of background error correlations. For instance, in the one-dimensional case, the radius of influence denotes the maximum number of non-zero elements, per row, in L. The EnKF-MC is then obtained by plugging in the estimator (3) in (6). Given the structure of the Cholesky factors, the EnKF-MC can be seen as a matrix-free implementation of the EnKF.

Recall that the precision analysis covariance matrix reads,

$$A^{-1} = \mathbf{B}^{-1} + \mathbf{H}^T \cdot \mathbf{R}^{-1} \cdot \mathbf{H} \cdot \mathbb{R}^{n \times n},\tag{15}$$

moreover, since <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> � <sup>H</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>n</sup> can be written as a sum of <sup>m</sup> rank-one matrices, the factors (3) can be updated in order to obtain an estimate of the inverse analysis covariance matrix. The next section discussed such approximation in detail.

#### 3.1.1. Posterior ensemble Kalman filter stochastic (PEnKF-S)

In the stochastic posterior ensemble Kalman filter (PEnKF-S) [7, 8], we want to estimate the moments of the analysis distribution,

$$\mathfrak{x} \sim \mathcal{N}(\mathfrak{x}^a, A), \tag{16}$$

based on the background ensemble X<sup>b</sup> � �, where xa is the analysis state and A ∈ R<sup>n</sup>�<sup>n</sup> is the analysis covariance matrix. Consider the estimate of the inverse background error covariance matrix Eq. (3), the precision analysis covariance matrix Eq. (4) can be approximated as follows:

$$\mathbf{A}^{-1} \approx \hat{\mathbf{A}}^{-1} = \hat{\mathbf{B}}^{-1} + \mathbf{X} \cdot \mathbf{X}^{T} \,\tag{17}$$

where <sup>X</sup> <sup>¼</sup> <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> <sup>2</sup> ∈ R<sup>n</sup>�<sup>m</sup>. The matrix (17) can be written as follows:

$$\widehat{A}^{-1} = L^T \cdot \mathcal{D}^{-1} \cdot L + \sum\_{i=1}^{m} \mathbf{x}\_i \cdot \mathbf{x}\_i^T \cdot \mathbf{x}\_i$$

where xi denotes the ith column of matrix X, for 1 ≤ i ≤ m. Consider the sequence of factor updates,

$$\left[\boldsymbol{L}^{i}\right]^{T}\cdot\mathbf{D}^{i}\cdot\boldsymbol{L}^{i} = \left[\boldsymbol{L}^{(i-1)}\right]^{T}\cdot\mathbf{D}^{i-1}\cdot\boldsymbol{L}^{i-1} + \boldsymbol{x}\_{i}\cdot\boldsymbol{x}\_{i}^{T}$$

$$\left[\boldsymbol{L}^{i}\right]^{T}\cdot\mathbf{D}^{i}\cdot\boldsymbol{L}^{i} = \left[\boldsymbol{L}^{(i-1)}\right]^{T}\cdot\left[\boldsymbol{D}^{(i-1)} + p\_{i}\cdot p\_{i}^{T}\right]\cdot\boldsymbol{L}^{(i-1)}$$

$$\left[\boldsymbol{L}^{i}\right]^{T}\cdot\boldsymbol{D}^{i}\cdot\boldsymbol{L}^{i} = \left[\boldsymbol{L}^{(i-1)}\cdot\boldsymbol{L}^{(i-1)}\right]^{T}\cdot\boldsymbol{\hat{D}}^{(i-1)}\cdot\left[\boldsymbol{L}^{(i-1)}\cdot\boldsymbol{L}^{(i-1)}\right],$$

$$\text{where }\boldsymbol{L}^{i-1}\cdot p\_{i} = \boldsymbol{x}\_{i}\in\mathbb{R}^{n\times 1}, \text{ for }1\leq i\leq m, \stackrel{\frown}{\mathbf{B}}^{-1} = \left[\boldsymbol{L}^{(0)}\right]^{T}\cdot\mathbf{D}^{(0)}\cdot\boldsymbol{L}^{(0)} \text{ and }$$

$$\mathbf{D}^{(i)} + \boldsymbol{\mathcal{p}}\_{i} \cdot \boldsymbol{\mathcal{p}}\_{i}^{\mathrm{T}} = \left[\boldsymbol{\tilde{L}}^{(i)}\right]^{\mathrm{T}} \cdot \boldsymbol{\tilde{D}}^{(i)} \cdot \boldsymbol{\tilde{L}}^{(i)} \in \mathbb{R}^{n \times n},\tag{18}$$

for 1 ≤ k < n and 1 ≤ j ≤ k � 1. The set of Eqs. (19) and (20) can be used in order to derive an algorithm for rank-one update of Cholesky factors, and the updating process is shown in

� �

Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization

http://dx.doi.org/10.5772/intechopen.72465

171

for all column vectors in X, and

� <sup>D</sup>ð Þ <sup>m</sup> � <sup>L</sup>ð Þ <sup>m</sup>

� �

� �

Algorithm 1. Rank-one update for the factors Lð Þ <sup>i</sup>�<sup>1</sup> and Dð Þ <sup>i</sup>�<sup>1</sup>

1. function UPD\_CHOLESKY\_FACTORS Lð Þ <sup>i</sup>�<sup>1</sup> ; Dð Þ <sup>i</sup>�<sup>1</sup> ; x<sup>i</sup>

2. Compute <sup>p</sup><sup>i</sup> from <sup>L</sup><sup>i</sup> � �<sup>T</sup> � pi <sup>¼</sup> xi

7. Compute ~lkk according to (11).

10. Set <sup>L</sup>ð Þ<sup>i</sup> <sup>L</sup>ð Þ <sup>i</sup>�<sup>1</sup> � <sup>L</sup>ð Þ <sup>i</sup>�<sup>1</sup> and <sup>D</sup>ð Þ<sup>i</sup> <sup>D</sup><sup>~</sup> ð Þ<sup>i</sup> .

this process is detailed in Algorithm 2.

2.

Algorithm 1 can be used in order to update the factors of <sup>B</sup><sup>b</sup> �<sup>1</sup>

Algorithm 2. Computing the factors <sup>L</sup>ð Þ <sup>m</sup> and <sup>D</sup>ð Þ <sup>m</sup> of <sup>L</sup>ð Þ <sup>m</sup> h i<sup>T</sup>

13. function COMPUTE\_ANALYSIS\_FACTOR Lð Þ<sup>0</sup> ; Dð Þ<sup>0</sup> ; H; R

16. Set <sup>L</sup>ð Þ<sup>i</sup> ; <sup>D</sup>ð Þ<sup>i</sup> h i UPD\_CHOLESKY\_FACTORS <sup>L</sup>ð Þ <sup>i</sup>�<sup>1</sup> ; <sup>D</sup>ð Þ <sup>i</sup>�<sup>1</sup> ; <sup>x</sup><sup>i</sup>

4. Compute ~dk via Eq. (10).

6. for j ¼ 1 ! k � 1 do

3. for k ¼ n � 1 do

5. Set lkk 1.

8. end for 9. end for

11. return Lð Þ<sup>i</sup> , Dð Þ<sup>i</sup> 12. end function

14. Set <sup>X</sup> <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup>

15. for i ¼ 1 ! m do

18. return Lð Þ <sup>m</sup> , Dð Þ <sup>m</sup>

17. end for

19. end function

Algorithm 1.

We can use of Dolittle's method in order to compute the factors D~ ð Þ<sup>i</sup> and L~ð Þ<sup>i</sup> in (9), and it is enough to note that,

1 ~l<sup>21</sup> ~l<sup>31</sup> ⋯ ~ln<sup>1</sup> 0 1 ~l<sup>32</sup> ⋯ ~ln<sup>2</sup> 00 1 ⋯ ~ln<sup>3</sup> ⋮⋮ ⋮ ⋱⋮ 00 0 ⋯ 1 2 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 5 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} <sup>L</sup>~ð Þ<sup>i</sup> ½ �<sup>T</sup> ~d<sup>1</sup> 0000 0 ~d<sup>2</sup> 000 0 0 ~d<sup>3</sup> 0 0 ⋮⋮⋮⋱ 0 0000 ~dn 2 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 5 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} <sup>D</sup><sup>~</sup> ð Þ<sup>i</sup> ½ � 100 … 0 ~l<sup>21</sup> 1 0 … 0 ~l<sup>31</sup> ~l<sup>21</sup> 1 … 0 ⋮⋮⋮⋱ 0 ~ln<sup>1</sup> ~ln<sup>2</sup> ~ln<sup>3</sup> … 1 2 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 5 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} <sup>L</sup>~ð Þ<sup>i</sup> ½ � ¼ <sup>d</sup><sup>1</sup> <sup>þ</sup> <sup>p</sup><sup>2</sup> <sup>1</sup> p<sup>1</sup> � p<sup>2</sup> p<sup>1</sup> � p<sup>3</sup> … p<sup>1</sup> � pn <sup>p</sup><sup>2</sup> � <sup>p</sup><sup>1</sup> <sup>d</sup><sup>2</sup> <sup>þ</sup> <sup>p</sup><sup>2</sup> <sup>2</sup> p<sup>2</sup> � p<sup>3</sup> … p<sup>2</sup> � pn <sup>p</sup><sup>3</sup> � <sup>p</sup><sup>1</sup> <sup>p</sup><sup>3</sup> � <sup>p</sup><sup>2</sup> <sup>d</sup><sup>3</sup> <sup>þ</sup> <sup>p</sup><sup>2</sup> <sup>3</sup> … p<sup>3</sup> � pn ⋮ ⋮ ⋮ ⋱⋮ pn � <sup>p</sup><sup>1</sup> pn � <sup>p</sup><sup>2</sup> pn � <sup>p</sup><sup>3</sup> … dn <sup>þ</sup> <sup>p</sup><sup>2</sup> n 2 6 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 7 5 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} <sup>D</sup>ð Þ<sup>i</sup> <sup>þ</sup>pi�pT <sup>i</sup> ½ �

After some math simplifications, the next equations are obtained,

$$
\tilde{d}\_k = p\_k^2 + d\_k - \sum\_{q=k+1}^n \tilde{d}\_q \cdot \tilde{l}\_{q i'}^2 \tag{19}
$$

and

$$\tilde{I}\_{k\dot{j}} = \frac{1}{\tilde{d}\_k} \cdot \left[ p\_k \cdot p\_j - \sum\_{q=k+1}^n \tilde{d}\_q \cdot \tilde{I}\_{q\dot{i}} \cdot \tilde{I}\_{q\dot{j}} \right],\tag{20}$$

for 1 ≤ k < n and 1 ≤ j ≤ k � 1. The set of Eqs. (19) and (20) can be used in order to derive an algorithm for rank-one update of Cholesky factors, and the updating process is shown in Algorithm 1.

Algorithm 1. Rank-one update for the factors Lð Þ <sup>i</sup>�<sup>1</sup> and Dð Þ <sup>i</sup>�<sup>1</sup>


<sup>L</sup><sup>i</sup> � �<sup>T</sup> � <sup>D</sup><sup>i</sup> � Li <sup>¼</sup> <sup>L</sup>ð Þ <sup>i</sup>�<sup>1</sup> h i<sup>T</sup>

<sup>L</sup><sup>i</sup> � �<sup>T</sup> � <sup>D</sup><sup>i</sup> � Li <sup>¼</sup> <sup>L</sup>ð Þ <sup>i</sup>�<sup>1</sup> h i<sup>T</sup>

where Li�<sup>1</sup> � pi <sup>¼</sup> xi <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>1</sup>

170 Kalman Filters - Theory for Advanced Applications

enough to note that,

and

1 ~l<sup>21</sup> ~l<sup>31</sup> ⋯ ~ln<sup>1</sup>

0 1 ~l<sup>32</sup> ⋯ ~ln<sup>2</sup>

00 1 ⋯ ~ln<sup>3</sup>

⋮⋮ ⋮ ⋱⋮

00 0 ⋯ 1


¼

<sup>L</sup><sup>i</sup> � �<sup>T</sup> � <sup>D</sup><sup>i</sup> � Li <sup>¼</sup> <sup>L</sup>~ð Þ <sup>i</sup>�<sup>1</sup> � <sup>L</sup>ð Þ <sup>i</sup>�<sup>1</sup> h i<sup>T</sup>

, for 1 <sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>m</sup>, <sup>B</sup><sup>b</sup> �<sup>1</sup>

<sup>D</sup>ð Þ<sup>i</sup> <sup>þ</sup> pi � pT

<sup>d</sup><sup>1</sup> <sup>þ</sup> <sup>p</sup><sup>2</sup>

After some math simplifications, the next equations are obtained,

<sup>~</sup>lkj <sup>¼</sup> <sup>1</sup> ~dk

<sup>~</sup>dk <sup>¼</sup> <sup>p</sup><sup>2</sup>

2 4

<sup>p</sup><sup>2</sup> � <sup>p</sup><sup>1</sup> <sup>d</sup><sup>2</sup> <sup>þ</sup> <sup>p</sup><sup>2</sup>

<sup>p</sup><sup>3</sup> � <sup>p</sup><sup>1</sup> <sup>p</sup><sup>3</sup> � <sup>p</sup><sup>2</sup> <sup>d</sup><sup>3</sup> <sup>þ</sup> <sup>p</sup><sup>2</sup>

⋮ ⋮ ⋮ ⋱⋮

pn � <sup>p</sup><sup>1</sup> pn � <sup>p</sup><sup>2</sup> pn � <sup>p</sup><sup>3</sup> … dn <sup>þ</sup> <sup>p</sup><sup>2</sup>


<sup>k</sup> <sup>þ</sup> dk � <sup>X</sup><sup>n</sup>

� pk � pj � <sup>X</sup><sup>n</sup>

q¼kþ1

q¼kþ1

<sup>~</sup>dq � <sup>~</sup><sup>l</sup> 2

<sup>~</sup>dq � <sup>~</sup>lqi � <sup>~</sup>lqj

3

� <sup>D</sup><sup>i</sup>�<sup>1</sup> � Li�<sup>1</sup> <sup>þ</sup> xi � xT

h i

� <sup>D</sup>ð Þ <sup>i</sup>�<sup>1</sup> <sup>þ</sup> pi � <sup>p</sup><sup>T</sup>

<sup>¼</sup> <sup>L</sup>ð Þ<sup>0</sup> h i<sup>T</sup>

We can use of Dolittle's method in order to compute the factors D~ ð Þ<sup>i</sup> and L~ð Þ<sup>i</sup> in (9), and it is

~d<sup>1</sup> 0000

0 ~d<sup>2</sup> 000

0 0 ~d<sup>3</sup> 0 0

⋮⋮⋮⋱ 0 0000 ~dn


<sup>1</sup> p<sup>1</sup> � p<sup>2</sup> p<sup>1</sup> � p<sup>3</sup> … p<sup>1</sup> � pn

<sup>2</sup> p<sup>2</sup> � p<sup>3</sup> … p<sup>2</sup> � pn

<sup>3</sup> … p<sup>3</sup> � pn

n

qi, (19)

5, (20)

<sup>i</sup> <sup>¼</sup> <sup>L</sup>~ð Þ<sup>i</sup> h i<sup>T</sup>

i

� <sup>D</sup><sup>~</sup> ð Þ <sup>i</sup>�<sup>1</sup> � <sup>L</sup>~ð Þ <sup>i</sup>�<sup>1</sup> � <sup>L</sup>ð Þ <sup>i</sup>�<sup>1</sup> h i

� <sup>D</sup>ð Þ<sup>0</sup> � <sup>L</sup>ð Þ<sup>0</sup> and

� <sup>L</sup>ð Þ <sup>i</sup>�<sup>1</sup>

,

� <sup>D</sup><sup>~</sup> ð Þ<sup>i</sup> � <sup>L</sup>~ð Þ<sup>i</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>n</sup>, (18)

100 … 0 ~l<sup>21</sup> 1 0 … 0

~l<sup>31</sup> ~l<sup>21</sup> 1 … 0

⋮⋮⋮⋱ 0 ~ln<sup>1</sup> ~ln<sup>2</sup> ~ln<sup>3</sup> … 1


i


Algorithm 1 can be used in order to update the factors of <sup>B</sup><sup>b</sup> �<sup>1</sup> for all column vectors in X, and this process is detailed in Algorithm 2.

Algorithm 2. Computing the factors <sup>L</sup>ð Þ <sup>m</sup> and <sup>D</sup>ð Þ <sup>m</sup> of <sup>L</sup>ð Þ <sup>m</sup> h i<sup>T</sup> � <sup>D</sup>ð Þ <sup>m</sup> � <sup>L</sup>ð Þ <sup>m</sup>

13. function COMPUTE\_ANALYSIS\_FACTOR Lð Þ<sup>0</sup> ; Dð Þ<sup>0</sup> ; H; R � �


$$\text{16.} \qquad \mathsf{Set}\left[\mathbf{L}^{(i)}, \mathbf{D}^{(i)}\right] \leftarrow \mathsf{UPD\\_CHOLESKY\\_FACTORS}\left(\mathbf{L}^{(i-1)}, \mathbf{D}^{(i-1)}, \mathbf{x}\_i\right).$$


Once the updating process has been performed, the resulting factors form an estimate of the inverse analysis covariance matrix,

$$\widehat{\mathbf{A}}^{-1} = \left[\mathbf{L}^{(m)}\right]^T \cdot \mathbf{D}^{(m)} \cdot \mathbf{L}^{(m)} \in \mathbb{R}^{n \times n},\tag{21}$$

From this covariance matrix, the posterior mode of the distribution can be approximated as follows:

$$
\overline{\mathfrak{X}}^{\mathfrak{a}} = \overline{\mathfrak{x}}^{\mathfrak{b}} + \mathfrak{z} \in \mathbb{R}^{n \times 1}, \tag{22}
$$

Bb �1

<sup>X</sup><sup>a</sup> <sup>¼</sup> xa � <sup>1</sup><sup>T</sup>

A1 <sup>2</sup> ≈ Ab 1 2

> B�<sup>1</sup> <sup>2</sup> <sup>≈</sup>B<sup>b</sup> �<sup>1</sup>

• The prediction is obtained from the posterior distribution: Pð Þ xjy .

• The posterior distribution is calculated: <sup>P</sup>ð Þ <sup>x</sup>j<sup>y</sup> <sup>∝</sup> <sup>P</sup><sup>b</sup>ð Þ� <sup>x</sup> <sup>P</sup>ð Þ <sup>y</sup>j<sup>x</sup> .

3.2. Markov Chain Monte Carlo-based filters

situation can be summarized as follows:

mean of Pð Þ yjx .

equation of the PEnKF-D

where

and

process.

<sup>2</sup> � <sup>X</sup><sup>b</sup> � xb � <sup>1</sup><sup>T</sup>

<sup>N</sup> þ Ab 1 2 � <sup>B</sup><sup>b</sup> �<sup>1</sup>

N

which are consistent with the dynamics of the numerical model, for instance, samples from (27) are driven by the physics and the dynamics of the numerical model. Finally, the analysis

<sup>¼</sup> <sup>L</sup>~<sup>T</sup> � <sup>D</sup><sup>~</sup> <sup>1</sup>

<sup>2</sup> <sup>¼</sup> <sup>L</sup> � <sup>D</sup><sup>1</sup>

Since the moments of the posterior distribution are unchanged, we expect both methods PEnKF-S and PEnKF-D to perform equivalently in terms of errors during the estimation

Definitely, the EnKF represents a breakthrough in the data assimilation context, perhaps its biggest appeal is that we can obtain a closed form expression for the analysis members. Nevertheless, as can be noted, during the derivation of the analysis equations, some constrains are imposed, for instance, the observation operator is assumed linear, and therefore, the likelihood Pð Þ yjx obeys a Gaussian distribution. Of course, in practice, this assumption can be easily broken, and therefore, bias can be introduced on the posterior prediction; this is notoriously significant if one considers the fields in which data assimilation lives are submerged. From a statistical point of view, more specifically, under the Bayesian framework, this

• The information from numerical model is known like the prior of background: <sup>P</sup><sup>b</sup>ð Þ<sup>x</sup> :

• The information from the observations is incorporated through the likelihood: Pð Þ yjx .

• If <sup>P</sup><sup>b</sup>ð Þ<sup>x</sup> and <sup>P</sup>ð Þ <sup>y</sup>j<sup>x</sup> are Gaussian, <sup>P</sup>ð Þ <sup>y</sup>j<sup>x</sup> is Gaussian too, and the analysis is equal to the

2 h i�<sup>1</sup>

� � � <sup>N</sup> ð Þ <sup>0</sup>;<sup>I</sup> , (27)

Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization

<sup>x</sup><sup>a</sup> <sup>¼</sup> xb <sup>þ</sup> <sup>A</sup> � <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> � <sup>y</sup> � <sup>H</sup> � xb � �, (29)

<sup>2</sup> � <sup>Δ</sup>X<sup>b</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>N</sup>, (28)

http://dx.doi.org/10.5772/intechopen.72465

173

, (30)

<sup>2</sup>, (31)

where

$$\mathbf{L}\left[\mathbf{L}^{(m)}\right]^T \cdot \mathbf{D}^{(m)} \cdot \mathbf{L}^{(m)} \cdot \mathbf{z} = \mathbf{q},\tag{23}$$

with <sup>q</sup> <sup>¼</sup> <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> � <sup>y</sup> � <sup>H</sup> � xb � �<sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>1</sup> . Note that the linear system (23) involves lower and upper triangular matrices, and therefore, xa can be estimated without the need of matrix inversion. Once the posterior mode is computed, the analysis ensemble is built about it. Note that <sup>A</sup><sup>b</sup> reads <sup>A</sup><sup>b</sup> <sup>¼</sup> <sup>L</sup>ð Þ <sup>m</sup> h i�<sup>1</sup> � <sup>D</sup>ð Þ <sup>m</sup> h i�<sup>1</sup> � <sup>L</sup>ð Þ <sup>m</sup> h i�<sup>T</sup> , and therefore, a square root of Ab can be approximated as follows:

$$
\widehat{\boldsymbol{A}}^{\frac{1}{2}} = \left[\boldsymbol{L}^{(m)}\right]^{\mathbf{1}} \cdot \left[\boldsymbol{D}^{(m)}\right]^{-\frac{1}{2}} \in \mathbb{R}^{n \times n},\tag{24}
$$

which can be utilized in order to build the analysis ensemble,

$$\mathbf{X}^{a} = \overline{\mathbf{x}}^{a} \cdot \mathbf{1}\_{N}^{T} + \Delta \mathbf{X}^{a},\tag{25}$$

where ΔX<sup>a</sup> ∈ R<sup>n</sup>�<sup>N</sup> is given by the solution of the linear system,

$$\mathbf{L}^{(m)} \cdot \left[\mathbf{D}^{(m)}\right]^{-\frac{1}{2}} \cdot \Delta \mathbf{X}^{\mathfrak{a}} = \mathcal{W} \in \mathbb{R}^{n \times N},\tag{26}$$

In addition, the columns ofW ∈ R<sup>n</sup>�<sup>N</sup> are formed by samples from a multivariate standard normal distribution. Again, since Lð Þ <sup>m</sup> is lower triangular, the solution of (26) can be obtained readily.

#### 3.1.2. Posterior ensemble Kalman filter deterministic (PEnKF-D)

The deterministic posterior ensemble Kalman filter (PEnKF-D) is a square root formulation of the PEnKF-S. The main difference between them lies on the use of synthetic data. In both methods, the matrices A�<sup>1</sup> and B�<sup>1</sup> are computed similarly (i.e., by using Dolittle's method and the algorithms described in the previous section).

In the PEnKF-D, the computation of <sup>B</sup><sup>b</sup> �<sup>1</sup> is performed (13) based on the empirical moments of the ensemble. The analysis update is performed by using the perturbations

Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization http://dx.doi.org/10.5772/intechopen.72465 173

$$
\widehat{\mathbf{B}}^{-\frac{1}{2}} \cdot \left( \mathbf{X}^{b} - \overline{\mathbf{x}}^{b} \cdot \mathbf{1}\_{N}^{T} \right) \sim \mathcal{N}(\mathbf{0}, \mathbf{I}),
\tag{27}
$$

which are consistent with the dynamics of the numerical model, for instance, samples from (27) are driven by the physics and the dynamics of the numerical model. Finally, the analysis equation of the PEnKF-D

$$\mathbf{X}^{a} = \overline{\mathbf{x}}^{a} \cdot \mathbf{1}\_{N}^{T} + \widehat{\mathbf{A}}^{\frac{1}{2}} \cdot \widehat{\mathbf{B}}^{-\frac{1}{2}} \cdot \Delta \mathbf{X}^{b} \in \mathbb{R}^{n \times N},\tag{28}$$

where

Once the updating process has been performed, the resulting factors form an estimate of the

From this covariance matrix, the posterior mode of the distribution can be approximated as

xa <sup>¼</sup> <sup>x</sup><sup>b</sup> <sup>þ</sup> <sup>z</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>1</sup>

upper triangular matrices, and therefore, xa can be estimated without the need of matrix inversion. Once the posterior mode is computed, the analysis ensemble is built about it. Note

> � <sup>D</sup>ð Þ <sup>m</sup> h i�<sup>1</sup> 2

� <sup>L</sup>ð Þ <sup>m</sup> h i�<sup>T</sup>

<sup>X</sup><sup>a</sup> <sup>¼</sup> <sup>x</sup><sup>a</sup> � <sup>1</sup><sup>T</sup>

2

In addition, the columns ofW ∈ R<sup>n</sup>�<sup>N</sup> are formed by samples from a multivariate standard normal distribution. Again, since Lð Þ <sup>m</sup> is lower triangular, the solution of (26) can be obtained readily.

The deterministic posterior ensemble Kalman filter (PEnKF-D) is a square root formulation of the PEnKF-S. The main difference between them lies on the use of synthetic data. In both methods, the matrices A�<sup>1</sup> and B�<sup>1</sup> are computed similarly (i.e., by using Dolittle's method

<sup>¼</sup> <sup>L</sup>ð Þ <sup>m</sup> h i-<sup>1</sup>

� <sup>D</sup>ð Þ <sup>m</sup> � <sup>L</sup>ð Þ <sup>m</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>n</sup>, (21)

� <sup>D</sup>ð Þ <sup>m</sup> � <sup>L</sup>ð Þ <sup>m</sup> � <sup>z</sup> <sup>¼</sup> <sup>q</sup>, (23)

. Note that the linear system (23) involves lower and

, and therefore, a square root of Ab can be

<sup>N</sup> <sup>þ</sup> <sup>Δ</sup>X<sup>a</sup>, (25)

� <sup>Δ</sup>X<sup>a</sup> <sup>¼</sup> <sup>W</sup> <sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>N</sup>, (26)

is performed (13) based on the empirical moments of

∈ R<sup>n</sup>�<sup>n</sup>, (24)

, (22)

inverse analysis covariance matrix,

172 Kalman Filters - Theory for Advanced Applications

with <sup>q</sup> <sup>¼</sup> <sup>H</sup><sup>T</sup> � <sup>R</sup>�<sup>1</sup> � <sup>y</sup> � <sup>H</sup> � xb � �<sup>∈</sup> <sup>R</sup><sup>n</sup>�<sup>1</sup>

that <sup>A</sup><sup>b</sup> reads <sup>A</sup><sup>b</sup> <sup>¼</sup> <sup>L</sup>ð Þ <sup>m</sup> h i�<sup>1</sup>

approximated as follows:

follows:

where

Ab �1

<sup>¼</sup> <sup>L</sup>ð Þ <sup>m</sup> h i<sup>T</sup>

<sup>L</sup>ð Þ <sup>m</sup> h i<sup>T</sup>

� <sup>D</sup>ð Þ <sup>m</sup> h i�<sup>1</sup>

Ab 1 2

which can be utilized in order to build the analysis ensemble,

where ΔX<sup>a</sup> ∈ R<sup>n</sup>�<sup>N</sup> is given by the solution of the linear system,

3.1.2. Posterior ensemble Kalman filter deterministic (PEnKF-D)

and the algorithms described in the previous section).

In the PEnKF-D, the computation of <sup>B</sup><sup>b</sup> �<sup>1</sup>

<sup>L</sup>ð Þ <sup>m</sup> � <sup>D</sup>ð Þ <sup>m</sup> h i�<sup>1</sup>

the ensemble. The analysis update is performed by using the perturbations

$$
\overline{\mathfrak{X}}^a = \overline{\mathfrak{x}}^b + A \cdot H^T \cdot \mathbb{R}^{-1} \cdot [y - H \cdot \overline{\mathfrak{X}}^b],
\tag{29}
$$

$$\boldsymbol{\dot{A}^{\frac{1}{2}}} \approx \boldsymbol{\hat{\dot{A}}^{\frac{1}{2}}} = \left[\boldsymbol{\tilde{L}}^{T} \cdot \boldsymbol{\tilde{D}}^{\frac{1}{2}}\right]^{-1} \text{.} \tag{30}$$

and

$$\mathbf{B}^{-\frac{1}{2}} \approx \widehat{\mathbf{B}}^{-\frac{1}{2}} = \mathbf{L} \cdot \mathbf{D}^{\frac{1}{2}},\tag{31}$$

Since the moments of the posterior distribution are unchanged, we expect both methods PEnKF-S and PEnKF-D to perform equivalently in terms of errors during the estimation process.

#### 3.2. Markov Chain Monte Carlo-based filters

Definitely, the EnKF represents a breakthrough in the data assimilation context, perhaps its biggest appeal is that we can obtain a closed form expression for the analysis members. Nevertheless, as can be noted, during the derivation of the analysis equations, some constrains are imposed, for instance, the observation operator is assumed linear, and therefore, the likelihood Pð Þ yjx obeys a Gaussian distribution. Of course, in practice, this assumption can be easily broken, and therefore, bias can be introduced on the posterior prediction; this is notoriously significant if one considers the fields in which data assimilation lives are submerged. From a statistical point of view, more specifically, under the Bayesian framework, this situation can be summarized as follows:


A significant number of approaches have been proposed in the current literature in order to deal with these constraints; for instance, the particle filters (PFs) and the maximum likelihood ensemble filter (MLEF) are ensemble-based methods, which deal with nonlinear observation operators. Unfortunately, in the context of PF, its use is questionable under realistic weather forecast scenarios since the number of particles (ensemble members) increases exponentially regarding the number of components in the model state. Anyways, an extended analysis of these methods exceeds the scope of this document, but its exploration is highly recommended.

The terms that involve Qð Þ� can be canceled if a symmetric distribution is chosen as the proposal one, for instance, in the Gaussian distribution case. The procedure for the calculation

Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization

http://dx.doi.org/10.5772/intechopen.72465

175

1. Initialize the Markov Chain, <sup>C</sup>, assigning the background value xb to <sup>C</sup>ð Þ<sup>0</sup> .

xa <sup>¼</sup> <sup>1</sup>

ð Þ <sup>k</sup> � <sup>p</sup> �

MCMC methods are straightforward to implement; when an enough number of iterations is performed, the behavior of the target function is captured on the chain [10]. This is true for even, complex density functions such as those with multiple modes. Briefly, let us focus on the fact that, generally, simulation-based methods such as MCMC explore a discretized grid, and as the mesh is refined, a huge number of iterations are needed before a high probability zone of the posterior distribution is reached [11]. Concretely, Hu et al. [12] proposed a family of modified MCMC dimension-independent algorithms under the name of preconditioned Crank-Nicolson (pCN) MCMC. These methods are robust regarding the curse of dimensionality in the statistical context. Initially, the Crank-Nicolson discretization is applied to a stochastic partial differential equation (SPDE) in order to obtain a new expression for the proposal

δKL xð Þ <sup>t</sup>�<sup>1</sup> þ xð Þ<sup>c</sup>

where L is the precision matrix of B, K is a preconditioning matrix, e<sup>0</sup> is white noise, if K ¼ B

δKL xð Þ <sup>t</sup>�<sup>1</sup> þ xð Þ<sup>c</sup>

� � <sup>þ</sup> ffiffiffiffiffiffiffiffiffi

� � <sup>þ</sup> ffiffiffiffiffiffiffiffiffi

. Of course, we use P<sup>b</sup> like a computationally efficient

<sup>α</sup> <sup>¼</sup> min <sup>1</sup>; exp J xð Þ� <sup>t</sup>�<sup>1</sup>; <sup>y</sup> J x<sup>c</sup> ð Þ ð Þ ; <sup>y</sup> � �, (35)

<sup>2</sup>K<sup>δ</sup> <sup>p</sup> <sup>e</sup>0, (33)

<sup>2</sup>K<sup>δ</sup> <sup>p</sup> , (34)

X k

Cð Þ<sup>i</sup> ,

i¼k�p

of the analysis applying MH is described below:

2. Generate a candidate vector state xð Þ<sup>c</sup> , fromQ. 3. ObtainU from a uniform (0, 1) distribution.

4. If U ≤ α Then: Cð Þ <sup>t</sup>þ<sup>1</sup> ¼ xð Þ<sup>c</sup> Else: Cð Þ <sup>t</sup>þ<sup>1</sup> ¼ Cð Þ<sup>t</sup>

The analysis is computed over the sample:

distribution:

5. Repeat Steps 2 through 4, for t ¼ 1 until t ¼ k � 1

6. Remove the first p vectors of the chain, the burned ones.

<sup>x</sup>ð Þ<sup>c</sup> <sup>¼</sup> <sup>x</sup>ð Þ <sup>t</sup>�<sup>1</sup> � <sup>1</sup>

and δ∈ ½ � 0; 2 in (33), we get the pCN proposal described in (34):

estimation of B. The acceptance probability is defined in (35):

where <sup>ω</sup>e<sup>N</sup> ð Þ <sup>0</sup>;<sup>B</sup> and <sup>β</sup> <sup>¼</sup> <sup>8</sup>δ=ð Þ <sup>2</sup> <sup>þ</sup> <sup>δ</sup> <sup>2</sup>

<sup>x</sup>ð Þ<sup>c</sup> <sup>¼</sup> <sup>x</sup>ð Þ <sup>t</sup>�<sup>1</sup> � <sup>1</sup>

2

2

Taking samples directly from the posterior distribution is a strategy that can help to remove the bias induced by wrong assumptions on the posterior distribution (i.e., normal assumptions). We do not put the sights on finding the mode of the posterior distribution; instead of this, we want a set of state vectors that allow us to create a satisfactory representation of the posterior error distribution, then, based on these samples, it is possible to estimate moments of the posterior distribution from which the analysis ensemble can be obtained [9].

Hereafter, we will construct a modification of the sequential scheme of data assimilation; first, we will describe how to compute the analysis members by using variations in Markov Chain Monte Carlo (MCMC)-based methods, and then, we will include the modified Cholesky decomposition in order to estimate a precision background covariance.

An overview of the proposed method is as follows:


In order to be more specific in the explanation, let us define Metropolis-Hastings (MH) as the selected algorithm from the MCMC family. Now let us focus on MCMC in the most intuitive way possible. Let us define J xð Þ ; y equal to the posterior pdf obtained from the expression <sup>P</sup>ð Þ <sup>x</sup>j<sup>y</sup> <sup>∝</sup> <sup>P</sup><sup>b</sup>ð Þ� <sup>x</sup> <sup>P</sup>ð Þ <sup>y</sup>j<sup>x</sup> . MH explores the state space in order to include model states in a socalled Markov Chain. The selection of candidates is based on the condition <sup>x</sup>ð Þ<sup>c</sup> ; <sup>y</sup> <sup>&</sup>gt; J xð Þ <sup>t</sup>�<sup>1</sup> ; <sup>y</sup> , that is, if the value of <sup>J</sup>ð Þ� for a candidate <sup>x</sup>ð Þ<sup>c</sup> is greater than that for a previous vector xð Þ <sup>t</sup>�<sup>1</sup> . Candidates are generated from a pre-defined proposal distribution Qð Þ� ; generally, a Gaussian multivariate distribution with mean xð Þ <sup>t</sup>�<sup>1</sup> is chosen, and a covariance matrix is defined in order to handle the step sizes. Concisely, MCMC methods proceed as follows: at first iteration, the chain is started with xb vector. At end, the sample is obtained by extracting the last t vectors on the chain, and the other vectors are dismissed or "burned". In (32), we can see the expression to calculate the acceptance probability that relates the target pdf, Jð Þ� , with the proposal Qð Þ�

$$\alpha = \frac{\mathbf{J}(\mathbf{x}\_{(c)}|\mathbf{y})\mathbf{Q}(\mathbf{x}\_{(t-1)}|\mathbf{x}\_{(c)})}{\mathbf{J}(\mathbf{x}\_{(t-1)}|\mathbf{y})\mathbf{Q}(\mathbf{x}\_{(c)}|\mathbf{x}\_{(t-1)})},\tag{32}$$

The terms that involve Qð Þ� can be canceled if a symmetric distribution is chosen as the proposal one, for instance, in the Gaussian distribution case. The procedure for the calculation of the analysis applying MH is described below:


A significant number of approaches have been proposed in the current literature in order to deal with these constraints; for instance, the particle filters (PFs) and the maximum likelihood ensemble filter (MLEF) are ensemble-based methods, which deal with nonlinear observation operators. Unfortunately, in the context of PF, its use is questionable under realistic weather forecast scenarios since the number of particles (ensemble members) increases exponentially regarding the number of components in the model state. Anyways, an extended analysis of these methods exceeds the scope of this document, but its exploration is highly recommended. Taking samples directly from the posterior distribution is a strategy that can help to remove the bias induced by wrong assumptions on the posterior distribution (i.e., normal assumptions). We do not put the sights on finding the mode of the posterior distribution; instead of this, we want a set of state vectors that allow us to create a satisfactory representation of the posterior error distribution, then, based on these samples, it is possible to estimate moments of

the posterior distribution from which the analysis ensemble can be obtained [9].

decomposition in order to estimate a precision background covariance.

to the ensemble of states of the previous assimilation cycle.

An overview of the proposed method is as follows:

174 Kalman Filters - Theory for Advanced Applications

pdf, Jð Þ� , with the proposal Qð Þ�

Hereafter, we will construct a modification of the sequential scheme of data assimilation; first, we will describe how to compute the analysis members by using variations in Markov Chain Monte Carlo (MCMC)-based methods, and then, we will include the modified Cholesky

1. The forecast step is unchanged; the forecast ensemble is obtained by applying the model

2. The analysis step is modified, so that the analysis is not obtained anymore by, for instance (28), but we perform k iterations of an algorithm from the Markov Chain Monte Carlo

In order to be more specific in the explanation, let us define Metropolis-Hastings (MH) as the selected algorithm from the MCMC family. Now let us focus on MCMC in the most intuitive way possible. Let us define J xð Þ ; y equal to the posterior pdf obtained from the expression <sup>P</sup>ð Þ <sup>x</sup>j<sup>y</sup> <sup>∝</sup> <sup>P</sup><sup>b</sup>ð Þ� <sup>x</sup> <sup>P</sup>ð Þ <sup>y</sup>j<sup>x</sup> . MH explores the state space in order to include model states in a socalled Markov Chain. The selection of candidates is based on the condition <sup>x</sup>ð Þ<sup>c</sup> ; <sup>y</sup> <sup>&</sup>gt; J xð Þ <sup>t</sup>�<sup>1</sup> ; <sup>y</sup> , that is, if the value of <sup>J</sup>ð Þ� for a candidate <sup>x</sup>ð Þ<sup>c</sup> is greater than that for a previous vector xð Þ <sup>t</sup>�<sup>1</sup> . Candidates are generated from a pre-defined proposal distribution Qð Þ� ; generally, a Gaussian multivariate distribution with mean xð Þ <sup>t</sup>�<sup>1</sup> is chosen, and a covariance matrix is defined in order to handle the step sizes. Concisely, MCMC methods proceed as follows: at first iteration, the chain is started with xb vector. At end, the sample is obtained by extracting the last t vectors on the chain, and the other vectors are dismissed or "burned". In (32), we can see the expression to calculate the acceptance probability that relates the target

<sup>α</sup> <sup>¼</sup> J xð Þ<sup>c</sup> <sup>j</sup><sup>y</sup> Q xð Þ <sup>t</sup>�<sup>1</sup> <sup>j</sup>xð Þ<sup>c</sup>

J xð Þ <sup>t</sup>�<sup>1</sup> <sup>j</sup><sup>y</sup> Q xð Þ<sup>c</sup> <sup>j</sup>xð Þ <sup>t</sup>�<sup>1</sup>

, (32)

(MCMC) family in order to obtain samples from the posterior distribution.


The analysis is computed over the sample:

$$\mathfrak{x}^a = \frac{1}{(k-p)} \cdot \sum\_{i=k-p}^{k} C\_{(i)\prime}$$

MCMC methods are straightforward to implement; when an enough number of iterations is performed, the behavior of the target function is captured on the chain [10]. This is true for even, complex density functions such as those with multiple modes. Briefly, let us focus on the fact that, generally, simulation-based methods such as MCMC explore a discretized grid, and as the mesh is refined, a huge number of iterations are needed before a high probability zone of the posterior distribution is reached [11]. Concretely, Hu et al. [12] proposed a family of modified MCMC dimension-independent algorithms under the name of preconditioned Crank-Nicolson (pCN) MCMC. These methods are robust regarding the curse of dimensionality in the statistical context. Initially, the Crank-Nicolson discretization is applied to a stochastic partial differential equation (SPDE) in order to obtain a new expression for the proposal distribution:

$$\mathbf{x}\_{(c)} = \mathbf{x}\_{(t-1)} - \frac{1}{2}\delta\mathsf{KC}\mathcal{L}\left(\mathbf{x}\_{(t-1)} + \mathbf{x}\_{(c)}\right) + \sqrt{2\mathsf{KC}\delta}\mathbf{e}\_{0\prime} \tag{33}$$

where L is the precision matrix of B, K is a preconditioning matrix, e<sup>0</sup> is white noise, if K ¼ B and δ∈ ½ � 0; 2 in (33), we get the pCN proposal described in (34):

$$\mathfrak{x}\_{(\varepsilon)} = \mathfrak{x}\_{(t-1)} - \frac{1}{2}\delta\mathfrak{K}\mathcal{L}\left(\mathfrak{x}\_{(t-1)} + \mathfrak{x}\_{(\varepsilon)}\right) + \sqrt{2\mathfrak{K}\delta},\tag{34}$$

where <sup>ω</sup>e<sup>N</sup> ð Þ <sup>0</sup>;<sup>B</sup> and <sup>β</sup> <sup>¼</sup> <sup>8</sup>δ=ð Þ <sup>2</sup> <sup>þ</sup> <sup>δ</sup> <sup>2</sup> . Of course, we use P<sup>b</sup> like a computationally efficient estimation of B. The acceptance probability is defined in (35):

$$\alpha = \min\{1, \exp(\mathbf{J}(\mathbf{x}\_{t-1}, \mathbf{y}) - \mathbf{J}(\mathbf{x}\_t, \mathbf{y}))\},\tag{35}$$

The procedure for the calculation of the analysis applying pCN Metropolis-Hastings is described as follows:

2. The reference solution is perturbed by using samples from a normal distribution with parameters N ð Þ 0n; σ<sup>B</sup> � I . Three different values for σ<sup>B</sup> are considered during the numerical experiments σ<sup>B</sup> ∈f g 0:05; 0:10; 0:15 . This perturbed state is propagated in time in order to make it consistent with the physics and dynamics of the numerical model. From here, an

Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization

http://dx.doi.org/10.5772/intechopen.72465

3. A similar procedure is performed in order to build a perturbed ensemble about xb

ensemble members are propagated in time from where an ensemble of model realizations

4. The assimilation windows consist of 15 equidistant observations. The frequency of obser-

5. The dimension of the vector state is n ¼ 40. The external force of the numerical model is

8. As a measure of quality, the L<sup>2</sup> norm of the analysis state and the reference solution are

9. A total of 100 runs are performed for each pair ð Þ N; σ<sup>B</sup> . For each run, a different initial random vector is utilized in order to build the initial perturbed reference solution x<sup>þ</sup>

The average of the error norms of each pair ð Þ N; σ<sup>B</sup> for the LETKF and the PEnKF implementations is shown in Table 1. As can be seen, on average across 100 runs, the

> 40 20,5671 18,2548 60 20,0567 17,8824

> 40 20,9513 18,3542 60 18,5048 17,8240

> 40 21,1314 18,1731 60 20,8487 17,7590

�<sup>2</sup>). This yields to different initial ensembles as well as

�<sup>1</sup>. The

177

�3

�<sup>1</sup> is obtained.

vations is 0.5 time units, which represents 3.5 days in the atmosphere.

7. Three ensemble sizes are tried during the experiments ∈f g 20; 40; 60 .

synthetic data for the different runs of each configuration (pair).

σ<sup>B</sup> N LETKF PEnKF-S 0.05 20 22,6166 21,2591

0.10 20 23,1742 21,0725

0.15 20 24,8201 20,9059

Table 1. Average of L-2 norm of errors for 100 runs of each configuration ð Þ σB; N for the compared filter

6. The number of observed components is 50% of the dimension of the vector state.

initial background state xb

of model is obtained.

computed across assimilation steps.

(before the model is applied x<sup>∗</sup>

x b i½ � 0 n o<sup>N</sup>

implementations.

i¼1

set to F ¼ 80.


Finally, in the data assimilation context, by using the modified Cholesky decomposition described in the earlier sections, the pCN-MH filter reads:


Now we are ready to numerically test the methods discussed in this chapter.

### 4. Experimental results

#### 4.1. Stochastic filter PEnKF-S

We assess the accuracy of the PEnKF-S and compare it against that of the LETFK implementation proposed by Hunt [13]. The numerical model is the Lorenz 96 model [11], which mimics the behavior of the atmosphere:

$$\frac{d\mathbf{x}\_k}{dt} = -\mathbf{x}\_{k-1} \cdot (\mathbf{x}\_{k-2} - \mathbf{x}\_{k+1}) - \mathbf{x}\_k + F\_\prime f \text{ for } 1 \le k \le n,\tag{36}$$

where n is the number of components in the model and F is an external force; with the value of the parameter, F ¼ 8, the model presents a great entropy.

The experimental settings are described below:

1. An initial random solution x<sup>þ</sup> �<sup>3</sup> is propagated for a while in time by using Lorenz 96 model and a fourth-order Runge Kutta method in order to obtain a vector state x<sup>þ</sup> �<sup>2</sup> whose physics are consistent with the dynamics of such numerical model. This vector state serves as our reference solution.


The procedure for the calculation of the analysis applying pCN Metropolis-Hastings is

Finally, in the data assimilation context, by using the modified Cholesky decomposition

1. The forecast step remains unchanged; the forecast ensemble is obtained by applying the model to the ensemble of states of the previous iteration, but the estimation of background

2. The update step is modified, so that the analysis is obtained by run k iterations of pCN-

We assess the accuracy of the PEnKF-S and compare it against that of the LETFK implementation proposed by Hunt [13]. The numerical model is the Lorenz 96 model [11], which mimics

where n is the number of components in the model and F is an external force; with the value of

and a fourth-order Runge Kutta method in order to obtain a vector state x<sup>þ</sup>

physics are consistent with the dynamics of such numerical model. This vector state serves

2

ð Þ <sup>k</sup>�<sup>p</sup> �

P<sup>k</sup>

dt ¼ �xk�<sup>1</sup> � ð Þ� xk�<sup>2</sup> � xkþ<sup>1</sup> xk <sup>þ</sup> F, for <sup>1</sup> <sup>≤</sup> <sup>k</sup> <sup>≤</sup> n, (36)

�<sup>3</sup> is propagated for a while in time by using Lorenz 96 model

�<sup>2</sup> whose

<sup>i</sup>¼k�<sup>p</sup> <sup>C</sup>ð Þ<sup>i</sup> .

xð Þ <sup>t</sup>�<sup>1</sup> þ βω.

1. Initialize the Markov Chain, <sup>C</sup>, assigning the background value xb to <sup>C</sup>ð Þ<sup>0</sup> .

2. Generate a candidate vector state using (3): <sup>x</sup>ð Þ<sup>c</sup> <sup>¼</sup> <sup>1</sup> � <sup>β</sup><sup>2</sup> � �<sup>1</sup>

3. Obtain U from a uniform (0, 1) distribution.

4. If U ≤ α, then Cð Þ <sup>t</sup>þ<sup>1</sup> ¼ xð Þ<sup>c</sup> . Else: Cð Þ <sup>t</sup>þ<sup>1</sup> ¼ Cð Þ<sup>t</sup> :

5. Repeat Steps 2 through 4, for t ¼ 1 until t ¼ k � 1.

7. The analysis is calculated over the sample: xa <sup>¼</sup> <sup>1</sup>

described in the earlier sections, the pCN-MH filter reads:

covariance matrix is calculated by modified Cholesky.

MH, to obtain a sample of the posterior error distribution.

Now we are ready to numerically test the methods discussed in this chapter.

6. Remove the first p vectors of the chain, the burned ones.

described as follows:

176 Kalman Filters - Theory for Advanced Applications

4. Experimental results

4.1. Stochastic filter PEnKF-S

the behavior of the atmosphere:

1. An initial random solution x<sup>þ</sup>

as our reference solution.

dxk

The experimental settings are described below:

the parameter, F ¼ 8, the model presents a great entropy.


The average of the error norms of each pair ð Þ N; σ<sup>B</sup> for the LETKF and the PEnKF implementations is shown in Table 1. As can be seen, on average across 100 runs, the


Table 1. Average of L-2 norm of errors for 100 runs of each configuration ð Þ σB; N for the compared filter implementations.

Figure 1. Local ensemble transform Kalmar filter (LETKF).

performance of the proposed EnKF implementation outperforms that of the LETKF in terms of L<sup>2</sup> norm of the error. Even more, the PEnKF-S seems to be invariant to the initial background error σ<sup>B</sup> since, in all cases, when the ensemble size is increased a better estimation of the reference state x<sup>∗</sup> at different observation times is obtained. This can also obey to the estimation of background error correlations via the modified Cholesky decomposition since it is drastically improved whenever the ensemble size is increased as pointed out by Bickel and Levina in [1]. In such case, the error decreases by Oð Þ logð Þ n =N . This is crucial in the PEnKF-S formulation since estimates of the precision analysis covariance matrix are obtained by rankone updates on the inverse background error covariance matrix. On the other hand, in the LETKF context, increasing the ensemble size can improve the accuracy of the method but that is not better than the one shown by the PEnKF.

Some plots of the L<sup>2</sup> norm of error for the PEnKF and the LETKF across different configurations and runs are shown in Figure 1. Note that the error of the PEnKF decreases aggressively since the earlier iterations. In the LETKF context, the accuracy is similar to that of the PEnKF only at the end of the assimilation window.

model is the Lorenz 96 model. In this section, we are mainly interested on describing the

Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization

http://dx.doi.org/10.5772/intechopen.72465

179

<sup>i</sup> , where xi denotes the ith compo-

behavior of the method, especially when the observational operator is nonlinear.

The experimental settings are as follows:

2. The vector of states has n ¼ 40 components.

Figure 2. Posterior ensemble Kalmar filter stochastic (PEnKF-S).

4. The length of the Markov Chain is 1000.

5. The proposal is pCN with β ¼ 0:09.

nent of the vector state x.

3. The ensemble size is N ¼ 20.

1. The observational operator is quadratic: f g <sup>H</sup>ð Þ<sup>x</sup> <sup>i</sup> <sup>¼</sup> <sup>x</sup><sup>2</sup>

Figure 3. Posterior ensemble Kalmar filter deterministic (PEnKF-D).

#### 4.2. Deterministic filter PEnKF-D

Holding the settings from the previous section, experiments are performed by using the PEnKF-D. The results have similar behavior to that of the LETFK and the PEnKF-S implementation proposed by Niño et al. in [8]. This can be appreciated in Figures 1–3. As can be seen, the error distributions reveal similar behavior across all compared filters.

#### 4.3. MCMC filter based on modified Cholesky decomposition

We will describe experiments of our proposed filter based on a modified Cholesky decomposition and the preconditioned Crank-Nicolson Metropolis-Hastings. Again, the numerical Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization http://dx.doi.org/10.5772/intechopen.72465 179

Figure 2. Posterior ensemble Kalmar filter stochastic (PEnKF-S).

performance of the proposed EnKF implementation outperforms that of the LETKF in terms of L<sup>2</sup> norm of the error. Even more, the PEnKF-S seems to be invariant to the initial background error σ<sup>B</sup> since, in all cases, when the ensemble size is increased a better estimation of the reference state x<sup>∗</sup> at different observation times is obtained. This can also obey to the estimation of background error correlations via the modified Cholesky decomposition since it is drastically improved whenever the ensemble size is increased as pointed out by Bickel and Levina in [1]. In such case, the error decreases by Oð Þ logð Þ n =N . This is crucial in the PEnKF-S formulation since estimates of the precision analysis covariance matrix are obtained by rankone updates on the inverse background error covariance matrix. On the other hand, in the LETKF context, increasing the ensemble size can improve the accuracy of the method but that

Some plots of the L<sup>2</sup> norm of error for the PEnKF and the LETKF across different configurations and runs are shown in Figure 1. Note that the error of the PEnKF decreases aggressively since the earlier iterations. In the LETKF context, the accuracy is similar to that of the PEnKF

Holding the settings from the previous section, experiments are performed by using the PEnKF-D. The results have similar behavior to that of the LETFK and the PEnKF-S implementation proposed by Niño et al. in [8]. This can be appreciated in Figures 1–3. As can be seen, the

We will describe experiments of our proposed filter based on a modified Cholesky decomposition and the preconditioned Crank-Nicolson Metropolis-Hastings. Again, the numerical

error distributions reveal similar behavior across all compared filters.

4.3. MCMC filter based on modified Cholesky decomposition

is not better than the one shown by the PEnKF.

Figure 1. Local ensemble transform Kalmar filter (LETKF).

178 Kalman Filters - Theory for Advanced Applications

only at the end of the assimilation window.

4.2. Deterministic filter PEnKF-D

Figure 3. Posterior ensemble Kalmar filter deterministic (PEnKF-D).

model is the Lorenz 96 model. In this section, we are mainly interested on describing the behavior of the method, especially when the observational operator is nonlinear.

The experimental settings are as follows:


Acknowledgements

Author details

References

2017

Universidad del Norte in Barranquilla, Colombia.

\*Address all correspondence to: enino@uninorte.edu.co

2008;36(6):2577-2604. DOI: 10.1214/08-AOS600

filter. Monthly Weather Review 1998;126(6):1719-1724

tive to ensemble square root filters. Tellus A. 2008;60(2):361-371

a modified Cholesky decomposition. Atmosphere. 2017;8(7):125

tation. Ocean Dynamics. 2003;53(4):343-367

Applied Mathematics, 1988

tion. AIMS Geosciences. 2015;1:41-78

Universidad Del Norte, Barranquilla, Colombia

This work was supported by the Applied Math and Computational Science Lab (AML-CS) at

Efficient Matrix-Free Ensemble Kalman Filter Implementations: Accounting for Localization

http://dx.doi.org/10.5772/intechopen.72465

181

Elias David Niño Ruiz\*, Rolando Beltrán Arrieta and Alfonso Manuel Mancilla Herrera

Applied Math and Computer Science Laboratory, Department of Computer Science,

[1] Bickel PJ, Levina E. Covariance regularization by thresholding. The Annals of Statistics.

[2] Evensen G. The ensemble Kalman filter: Theoretical formulation and practical implemen-

[3] Burgers G, Jan van Leeuwen P, Evensen G. Analysis scheme in the ensemble Kalman

[4] Sakov P, Oke PR. A deterministic formulation of the ensemble Kalman filter: An alterna-

[5] Coleman TF, Van Loan C. Handbook for matrix computations. Society for Industrial and

[6] Nino-Ruiz ED, Sandu A, Deng X. A parallel implementation of the ensemble Kalman filter based on modified Cholesky decomposition. Journal of Computational Science.

[7] Nino-Ruiz ED, Mancilla A, Calabria JC. A posterior ensemble Kalman filter based on a modified Cholesky decomposition. Procedia Computer Science. 2017;108:2049-2058 [8] Nino-Ruiz ED. A matrix-free posterior ensemble Kalman filter implementation based on

[9] Attia A, Sandu A. A hybrid Monte Carlo sampling filter for non-gaussian data assimila-

[10] David H, Shane Reese C, David Moulton J, Vrugt JA, Colin F. Posterior exploration for computationally intensive. In: Steve B, Andrew G, Jones GL, Meng X-L, editors. Hand-

book of Markov Chain Monte Carlo. Chapman & Hall; 2011. p. 401-418

Figure 4. Histogram of six components of the state vector after 1000 iterations of the algorithm, the dashed lines indicate the position of the true value for the respective component.

Figure 4 describes the distribution of six of the 40 components of the state vector after 1000 iterations of the algorithm, attempting to visualize if the actual value is within or near the region of the highest probability density described by the sample contained in the Markov Chain.

Note that, not in all cases, the actual values of the model components were located at the peaks of the histogram, but most of them were within or near to zones of high probability described by the sample. This result is important if we take into account that probably the posterior distribution is not normal as is the case of quadratic observation operators.

### 5. Conclusions

In this chapter, efficient EnKF implementations were discussed. All of them are based on a modified Cholesky decomposition wherein a precision background covariance is obtained in terms of Cholesky factors. In the first filter, the PEnKF-S, synthetic data are utilized in order to compute the posterior members, as done in stochastic formulations of the filter. Even more, a sequence of rank-one updates can be applied over the factors of the prior precision matrix in order to estimate those of the posterior precision. In the second filter, the PEnKF-D, synthetic data are avoided by using perturbations obtained from the physics and the dynamics of the numerical model. Finally, a MCMC-based filter is obtained in order to reduce the impact of bias when Gaussian assumptions are broken during the assimilation of observations, for instance, when the observation operator is nonlinear. Numerical experiments with the Lorenz 96 model reveal that the proposed filters are comparable to filters from the specialized literature.

### Acknowledgements

This work was supported by the Applied Math and Computational Science Lab (AML-CS) at Universidad del Norte in Barranquilla, Colombia.

### Author details

Elias David Niño Ruiz\*, Rolando Beltrán Arrieta and Alfonso Manuel Mancilla Herrera

\*Address all correspondence to: enino@uninorte.edu.co

Applied Math and Computer Science Laboratory, Department of Computer Science, Universidad Del Norte, Barranquilla, Colombia

### References

Figure 4 describes the distribution of six of the 40 components of the state vector after 1000 iterations of the algorithm, attempting to visualize if the actual value is within or near the region of the highest probability density described by the sample contained in the Markov Chain.

Figure 4. Histogram of six components of the state vector after 1000 iterations of the algorithm, the dashed lines indicate

Note that, not in all cases, the actual values of the model components were located at the peaks of the histogram, but most of them were within or near to zones of high probability described by the sample. This result is important if we take into account that probably the posterior

In this chapter, efficient EnKF implementations were discussed. All of them are based on a modified Cholesky decomposition wherein a precision background covariance is obtained in terms of Cholesky factors. In the first filter, the PEnKF-S, synthetic data are utilized in order to compute the posterior members, as done in stochastic formulations of the filter. Even more, a sequence of rank-one updates can be applied over the factors of the prior precision matrix in order to estimate those of the posterior precision. In the second filter, the PEnKF-D, synthetic data are avoided by using perturbations obtained from the physics and the dynamics of the numerical model. Finally, a MCMC-based filter is obtained in order to reduce the impact of bias when Gaussian assumptions are broken during the assimilation of observations, for instance, when the observation operator is nonlinear. Numerical experiments with the Lorenz 96 model

reveal that the proposed filters are comparable to filters from the specialized literature.

distribution is not normal as is the case of quadratic observation operators.

the position of the true value for the respective component.

180 Kalman Filters - Theory for Advanced Applications

5. Conclusions


[11] Cotter SL, Roberts GO, Stuart AM, White D, et al. MCMC methods for functions: Modifying old algorithms to make them faster. Statistical Science. 2013;28(3):424-446

**Chapter 10**

**Provisional chapter**

**Kalman Filters for Reference Current Generation in**

Shunt active power filter (APF) method have been used by many researchers as a solution in reducing the harmonics creating by the non-liner loads. Therefore, this research is targeted to design and implement a three-phase shunt APF employing Kalman filter estimator. Conventionally, low-pass filter (LPF) is used to filter out the unwanted DC component of the non-linear load to produce the sinusoidal waveform called the reference current. However, when applying LPF it contributes with the phase shift and high transient at the supply current. Therefore, to reduce these problems, the digital Kalman filter estimator is used to replace the LPF for generating the reference current. Details on the investigation between conventional and proposed methods under simulation based on Matlab Simulink platform and experimental that are made for two types of load, namely, three-phase rectifier with RC-load and three-phase induction motor, are presented. The performance criteria of the shunt APF are determined by the supply current waveform, total harmonic distortion (THD), harmonic spectrum and power quality measurements, which were also obtained by simulation and experimental. In conclusion, by employing Kalman filter estimator for generating the reference current, it reduces the time delay and high transient current at the power supply and, thus, improved the overall THD from 0.1 to 0.42% compared to the LPF. **Keywords:** three-phase system, harmonic reduction, active power filter (APF), reference

**Kalman Filters for Reference Current Generation in** 

DOI: 10.5772/intechopen.72467

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution,

© 2018 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

and reproduction in any medium, provided the original work is properly cited.

Electrical power is essential to people's modern lifestyle. In recent five decades, due to the development of the industry contributed to the increase of the types and capacity of the

**Shunt Active Power Filter (APF)**

Syed Mohd Fairuz Bin Syed Mohd Dardin

**Shunt Active Power Filter (APF)**

Syed Mohd Fairuz Bin Syed Mohd Dardin and

Additional information is available at the end of the chapter

Additional information is available at the end of the chapter

Ahmad Shukri Bin Abu Hasim,

and Zulkifilie Bin Ibrahim

Ahmad Shukri Bin Abu Hasim,

http://dx.doi.org/10.5772/intechopen.72467

current generation, Kalman filter

Zulkifilie Bin Ibrahim

**Abstract**

**1. Introduction**


**Provisional chapter**

### **Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF) Shunt Active Power Filter (APF)**

**Kalman Filters for Reference Current Generation in** 

DOI: 10.5772/intechopen.72467

[11] Cotter SL, Roberts GO, Stuart AM, White D, et al. MCMC methods for functions: Modifying old algorithms to make them faster. Statistical Science. 2013;28(3):424-446

182 Kalman Filters - Theory for Advanced Applications

[12] Hu Z, Yao Z, Li J. On an adaptive preconditioned Crank–Nicolson MCMC algorithm for infinite dimensional Bayesian inference. Journal of Computational Physics. 2017;332:492-503

[13] Fatkullin I, Vanden-Eijnden E. A computational strategy for multiscale systems with applications to Lorenz 96 model. Journal of Computational Physics. 2004;200(2):605

> Ahmad Shukri Bin Abu Hasim, Syed Mohd Fairuz Bin Syed Mohd Dardin and Zulkifilie Bin Ibrahim Syed Mohd Fairuz Bin Syed Mohd Dardin and Zulkifilie Bin Ibrahim Additional information is available at the end of the chapter

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/intechopen.72467

Ahmad Shukri Bin Abu Hasim,

### **Abstract**

Shunt active power filter (APF) method have been used by many researchers as a solution in reducing the harmonics creating by the non-liner loads. Therefore, this research is targeted to design and implement a three-phase shunt APF employing Kalman filter estimator. Conventionally, low-pass filter (LPF) is used to filter out the unwanted DC component of the non-linear load to produce the sinusoidal waveform called the reference current. However, when applying LPF it contributes with the phase shift and high transient at the supply current. Therefore, to reduce these problems, the digital Kalman filter estimator is used to replace the LPF for generating the reference current. Details on the investigation between conventional and proposed methods under simulation based on Matlab Simulink platform and experimental that are made for two types of load, namely, three-phase rectifier with RC-load and three-phase induction motor, are presented. The performance criteria of the shunt APF are determined by the supply current waveform, total harmonic distortion (THD), harmonic spectrum and power quality measurements, which were also obtained by simulation and experimental. In conclusion, by employing Kalman filter estimator for generating the reference current, it reduces the time delay and high transient current at the power supply and, thus, improved the overall THD from 0.1 to 0.42% compared to the LPF.

**Keywords:** three-phase system, harmonic reduction, active power filter (APF), reference current generation, Kalman filter

### **1. Introduction**

Electrical power is essential to people's modern lifestyle. In recent five decades, due to the development of the industry contributed to the increase of the types and capacity of the

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2018 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

grid connected load drastically. For that reason, all the electrical consumers at all levels of usage have facing an issue of power quality problems. Both industrial/commercial sector and domestic environment commonly use sensitive equipment and non-linear loads (NLL). Inadvertently, these result in a non-sinusoidal current being drawn from the supply, which contains the harmful harmonic component and fed back to the supply system on the same point of common coupling (PCC). Passive filter is one of the common methods that have been used to overcome this problem. The passive filter is connected in parallel between the supply and NLL for improvement of power factor and harmonic suppression and thus exhibits lower impedance at tuned harmonic frequency. However, this approach does not solve the problem effectively due to its inability to compensate random frequency variation in the current, tuning problem and parallel resonant. Among the techniques, the *d-q* algorithm has been widely used to eliminate the harmonics due to its simplicity of control design relative to the rest. Commonly, the *d-q* algorithm is using LPF to generate the reference current. However, time delay introduced when applying the LPF will contribute to the phase shift in harmonics and higher transient current. Therefore, a new proposed technique of the current reference generator embedded with Kalman filter for shunt APF system is proposed where it reduces the time delay, thus producing improvement of the overall total harmonic distortion (THD) in the system.

words, it allows only fundamental component of the current to flow in the system, suppressing other higher-frequency components. It can also be used to regulate the negative sequence voltage at the load. The series active filter is ideal for eliminating and/or maintaining the output voltage while balancing three-phase voltages [20, 21, 23–26]. On the other hand, shunt APF has been widely used to mitigate the harmonics. It cancels the load-current harmonics and provides reactive compensation to the supply, through the act of injecting equal but opposite harmonic compensating current to the supply [27–38]. Shunt APF has the advantage of carrying only the compensation current plus a small amount of active fundamental current to compensate for system losses. It is also possible to connect several filters in parallel to cater for higher currents, making this type of circuit suitable for a wide range of power ratings [26, 39, 40]. The most common configuration of shunt APF is the inverter type where the role of the filter inductor is to suppress the high frequency at tuned current generated at tuned frequency, while the converter provides complementary filtering on others that includes any random variations through switching techniques [28, 34, 41, 42]. The shunt APF controller can be used in direct or indirect connection. Hybrid APF can be characterized by a combination of passive filter and APF in series or parallel. The combination between series APF with parallel passive filter is the most popular arrangement because the solid-state devices used in active series part help in reducing the size and cost, to about 60–80% of load size [43, 44]. Furthermore, the passive parallel LC filter is used to eliminate lower-order harmonics at reasonable cost [26, 37, 44–49]. Another arrangement is the combination of active filter in series with a parallel passive filter, which is used especially for mediumand high-voltage applications [26]. Further arrangements also include a combination of parallel active and passive filters where the APF part is designed to eliminate the lower order of harmonics, while the passive filter works to eliminate the bulk load-current harmonic [26]. The combination of series active and parallel APF will produce unified power quality conditioner (also known as universal AF). The DC-link element of either inductor or capacitor is shared between two current sources or voltage-source bridges operating as active series and active parallel compensator [50, 51]. This universal AF is considered as an ideal AF, which eliminates voltage and current harmonics, thus capable of providing clean power to critical and harmonic-prone loads, such as computer, medical equipment and others. The main drawbacks are large costs and complex control due to dependency on the number of solid-state devices involved [26, 50, 51].

Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF)

http://dx.doi.org/10.5772/intechopen.72467

185

Many control approaches have been developed to extract and estimate the harmonics in the system. Instantaneous reactive power theory (*p-q* theory), modified *p-q* theory [52–54], *p-q-r* theory [55, 56], vectorial theory [57] and *d-q* theory [58–60] are the techniques that fall into the extraction technique. Due to its simplicity of control design relative to the rest, for that reason this *d-q* algorithm has been widely used to eliminate the harmonics [61]. On the other hand, estimation approach is used to estimate harmonics of frequency component present in the signal and measurement or estimation of the amplitude and phases of those frequencies [62]. This approach can be divided into two classes, non-parametric and parametric methods. The non-parametric methods are based on transformation of the given time-series data sequence. During the estimation process, these methods are not capable of incorporating with any available information about the system. Frequency domain approach using Fourier transform is most commonly used for spectrum analysis in this harmonic estimation [62]. In addition, parametric methods use an appropriate model to represent the signal and then estimate the

### **1.1. State of the art**

The active power filter (APF) technology is now mature in providing compensation for harmonics, reactive power and neutral current in AC networks. It has evolved for the past quarter century of development with varying configurations, control strategies and solid-state devices. Commonly, the APFs are used to eliminate the voltage harmonics, regulate terminal voltage, suppress voltage flicker and improve voltage balance in three-phase systems. This wide range of objectives can be achieved either individually or in combination depending upon the requirements, control strategy and configuration, which have to be selected appropriately. This section describes the history of development and the present status of the APF technology.

With the proliferation of power electronics in energy conversions, power quality is fast becoming an issue of an increasingly important aspect of electrical consumers at all levels of usage. A large number of publications have been covering the power quality survey, measurements, analysis, cause and effects of harmonics and reactive power in the electric networks [1–9]. APFs can be categorized into three types, namely, two-wire (single-phase), three-wire and four-wire threephase configurations, to meet the requirements of the three types of NLL on supply systems. Domestic lights and ovens, TVs, computer power supplies, air conditioners, laser printers and Xerox machines behave as NLL and cause power quality problems for single-phase loads. For this type of load, the APFs are investigated in varying configurations and control strategies [10– 19]. Starting in 1971, many configurations of APF have been developed for improving the power system quality. It can be categorized into four basic types, namely, series, parallel (shunt), hybrid APFs and unified power quality conditioner (universal AF). The series APF operates mainly as a voltage regulator and a harmonic isolator between NLL and utility system [20–23]. In other words, it allows only fundamental component of the current to flow in the system, suppressing other higher-frequency components. It can also be used to regulate the negative sequence voltage at the load. The series active filter is ideal for eliminating and/or maintaining the output voltage while balancing three-phase voltages [20, 21, 23–26]. On the other hand, shunt APF has been widely used to mitigate the harmonics. It cancels the load-current harmonics and provides reactive compensation to the supply, through the act of injecting equal but opposite harmonic compensating current to the supply [27–38]. Shunt APF has the advantage of carrying only the compensation current plus a small amount of active fundamental current to compensate for system losses. It is also possible to connect several filters in parallel to cater for higher currents, making this type of circuit suitable for a wide range of power ratings [26, 39, 40]. The most common configuration of shunt APF is the inverter type where the role of the filter inductor is to suppress the high frequency at tuned current generated at tuned frequency, while the converter provides complementary filtering on others that includes any random variations through switching techniques [28, 34, 41, 42]. The shunt APF controller can be used in direct or indirect connection. Hybrid APF can be characterized by a combination of passive filter and APF in series or parallel. The combination between series APF with parallel passive filter is the most popular arrangement because the solid-state devices used in active series part help in reducing the size and cost, to about 60–80% of load size [43, 44]. Furthermore, the passive parallel LC filter is used to eliminate lower-order harmonics at reasonable cost [26, 37, 44–49]. Another arrangement is the combination of active filter in series with a parallel passive filter, which is used especially for mediumand high-voltage applications [26]. Further arrangements also include a combination of parallel active and passive filters where the APF part is designed to eliminate the lower order of harmonics, while the passive filter works to eliminate the bulk load-current harmonic [26]. The combination of series active and parallel APF will produce unified power quality conditioner (also known as universal AF). The DC-link element of either inductor or capacitor is shared between two current sources or voltage-source bridges operating as active series and active parallel compensator [50, 51]. This universal AF is considered as an ideal AF, which eliminates voltage and current harmonics, thus capable of providing clean power to critical and harmonic-prone loads, such as computer, medical equipment and others. The main drawbacks are large costs and complex control due to dependency on the number of solid-state devices involved [26, 50, 51].

grid connected load drastically. For that reason, all the electrical consumers at all levels of usage have facing an issue of power quality problems. Both industrial/commercial sector and domestic environment commonly use sensitive equipment and non-linear loads (NLL). Inadvertently, these result in a non-sinusoidal current being drawn from the supply, which contains the harmful harmonic component and fed back to the supply system on the same point of common coupling (PCC). Passive filter is one of the common methods that have been used to overcome this problem. The passive filter is connected in parallel between the supply and NLL for improvement of power factor and harmonic suppression and thus exhibits lower impedance at tuned harmonic frequency. However, this approach does not solve the problem effectively due to its inability to compensate random frequency variation in the current, tuning problem and parallel resonant. Among the techniques, the *d-q* algorithm has been widely used to eliminate the harmonics due to its simplicity of control design relative to the rest. Commonly, the *d-q* algorithm is using LPF to generate the reference current. However, time delay introduced when applying the LPF will contribute to the phase shift in harmonics and higher transient current. Therefore, a new proposed technique of the current reference generator embedded with Kalman filter for shunt APF system is proposed where it reduces the time delay, thus producing improvement of the overall total harmonic

The active power filter (APF) technology is now mature in providing compensation for harmonics, reactive power and neutral current in AC networks. It has evolved for the past quarter century of development with varying configurations, control strategies and solid-state devices. Commonly, the APFs are used to eliminate the voltage harmonics, regulate terminal voltage, suppress voltage flicker and improve voltage balance in three-phase systems. This wide range of objectives can be achieved either individually or in combination depending upon the requirements, control strategy and configuration, which have to be selected appropriately. This section describes the history of development and the present status of the APF

With the proliferation of power electronics in energy conversions, power quality is fast becoming an issue of an increasingly important aspect of electrical consumers at all levels of usage. A large number of publications have been covering the power quality survey, measurements, analysis, cause and effects of harmonics and reactive power in the electric networks [1–9]. APFs can be categorized into three types, namely, two-wire (single-phase), three-wire and four-wire threephase configurations, to meet the requirements of the three types of NLL on supply systems. Domestic lights and ovens, TVs, computer power supplies, air conditioners, laser printers and Xerox machines behave as NLL and cause power quality problems for single-phase loads. For this type of load, the APFs are investigated in varying configurations and control strategies [10– 19]. Starting in 1971, many configurations of APF have been developed for improving the power system quality. It can be categorized into four basic types, namely, series, parallel (shunt), hybrid APFs and unified power quality conditioner (universal AF). The series APF operates mainly as a voltage regulator and a harmonic isolator between NLL and utility system [20–23]. In other

distortion (THD) in the system.

184 Kalman Filters - Theory for Advanced Applications

**1.1. State of the art**

technology.

Many control approaches have been developed to extract and estimate the harmonics in the system. Instantaneous reactive power theory (*p-q* theory), modified *p-q* theory [52–54], *p-q-r* theory [55, 56], vectorial theory [57] and *d-q* theory [58–60] are the techniques that fall into the extraction technique. Due to its simplicity of control design relative to the rest, for that reason this *d-q* algorithm has been widely used to eliminate the harmonics [61]. On the other hand, estimation approach is used to estimate harmonics of frequency component present in the signal and measurement or estimation of the amplitude and phases of those frequencies [62]. This approach can be divided into two classes, non-parametric and parametric methods. The non-parametric methods are based on transformation of the given time-series data sequence. During the estimation process, these methods are not capable of incorporating with any available information about the system. Frequency domain approach using Fourier transform is most commonly used for spectrum analysis in this harmonic estimation [62]. In addition, parametric methods use an appropriate model to represent the signal and then estimate the parameters of the model from the available data points. Estimated parameters are applied to the selected model to determine harmonic contents in the signal. This parametric methods offer higher resolution and better accuracy than the non-parametric methods [62]. Kalman filter (KF) estimator is one of the methods that fall into the parametric method category which have been widely studied and used for different applications [62–69].

### **1.2. Main contribution**

There are three main contributions regarding with this research:

	- The new design of shunt APF for generating the reference currents using Kalman filter estimator was proposed to reduce the delay time and high transient current when applying the conventional technique. In addition, the developed system was tested for two different types of loads such as three-phase rectifier with RC-load and three-phase induction motor.
	- The investigation criteria are on the harmonic spectrum, THD and power quality for different three types of load because these criteria affect directly the performance of system that used active power filter.

as shown in **Figure 2**. The KF used a form of feedback control in which the filter estimates the process at any time and then obtains feedback in the form of noisy measurements. These noisy measurements can be further exploited to improve the next estimates in which KF is able to perform because it has both time update and measurement update equations. The time update also known as predictor equation is responsible for projecting forward (in time) current state and error covariance estimate to obtain the estimation in the next time step, while the measurement update equation also called corrector equation is responsible for the feedback such as for incorporating a new measurement into the estimator to improve the

Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF)

http://dx.doi.org/10.5772/intechopen.72467

187

**Figure 2.** New technique of three-phase reference current generator employing Kalman filter estimator.

**Figure 1.** Overall system of shunt APF.


### **1.3. Proposal of the research**

Power quality problems have becoming a critical issue when dealing with power electronic converter and NLL due to the effects of the harmonic contamination in power system. Many techniques have been proposed to overcome these problems such as passive filter which contribute to improve the power factor and harmonic suppression and exhibit lower impedance at a tuned harmonic frequency. However this approach provides incomplete solutions particularly when compensating random frequency variations in the current, tuning and parallel resonant problems. Therefore, various active power filter (APF) configurations with their respective control strategies have been proposed and have been recognized as a viable solution to the problem created by harmonics. Among the technique, the *d-q* algorithm has been widely used to eliminate the harmonics due to its simplicity of control design relative to the rest. Commonly, the *d-q* algorithm is using LPF to generate the reference current. However, time delay introduced when applying the LPF will contribute to the phase shift in harmonics and higher transient current. A new proposed technique of the current reference generator embedded with Kalman filter (KF) for shunt APF system is proposed as shown in **Figure 1**. The KF in the system was used as a LPF to produce a reference current in three-phase system Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF) http://dx.doi.org/10.5772/intechopen.72467 187

**Figure 1.** Overall system of shunt APF.

parameters of the model from the available data points. Estimated parameters are applied to the selected model to determine harmonic contents in the signal. This parametric methods offer higher resolution and better accuracy than the non-parametric methods [62]. Kalman filter (KF) estimator is one of the methods that fall into the parametric method category which

• The new design of shunt APF for generating the reference currents using Kalman filter estimator was proposed to reduce the delay time and high transient current when applying the conventional technique. In addition, the developed system was tested for two different types of loads such as three-phase rectifier with RC-load and three-phase induction motor.

**b.** Investigation of the details of performance based on simulation and experimental for con-

• The investigation criteria are on the harmonic spectrum, THD and power quality for different three types of load because these criteria affect directly the performance of system

**c.** Comparative studies between the conventional and the proposed technique upon experi-

**d.** An analysis is carried out in terms of harmonic spectrum, THD and power quality as well to validate the advantages offered by employing the new techniques relative to the com-

Power quality problems have becoming a critical issue when dealing with power electronic converter and NLL due to the effects of the harmonic contamination in power system. Many techniques have been proposed to overcome these problems such as passive filter which contribute to improve the power factor and harmonic suppression and exhibit lower impedance at a tuned harmonic frequency. However this approach provides incomplete solutions particularly when compensating random frequency variations in the current, tuning and parallel resonant problems. Therefore, various active power filter (APF) configurations with their respective control strategies have been proposed and have been recognized as a viable solution to the problem created by harmonics. Among the technique, the *d-q* algorithm has been widely used to eliminate the harmonics due to its simplicity of control design relative to the rest. Commonly, the *d-q* algorithm is using LPF to generate the reference current. However, time delay introduced when applying the LPF will contribute to the phase shift in harmonics and higher transient current. A new proposed technique of the current reference generator embedded with Kalman filter (KF) for shunt APF system is proposed as shown in **Figure 1**. The KF in the system was used as a LPF to produce a reference current in three-phase system

have been widely studied and used for different applications [62–69].

**a.** Developed a new design of shunt APF employing Kalman filter estimator.

There are three main contributions regarding with this research:

ventional and the proposed technique.

mon implementation of an active power filter.

that used active power filter.

mental implementation.

**1.3. Proposal of the research**

**1.2. Main contribution**

186 Kalman Filters - Theory for Advanced Applications

as shown in **Figure 2**. The KF used a form of feedback control in which the filter estimates the process at any time and then obtains feedback in the form of noisy measurements. These noisy measurements can be further exploited to improve the next estimates in which KF is able to perform because it has both time update and measurement update equations. The time update also known as predictor equation is responsible for projecting forward (in time) current state and error covariance estimate to obtain the estimation in the next time step, while the measurement update equation also called corrector equation is responsible for the feedback such as for incorporating a new measurement into the estimator to improve the

**Figure 2.** New technique of three-phase reference current generator employing Kalman filter estimator.

estimation. Therefore the estimation resembles the combination of predictor-corrector algorithm, which is used in the system. By applying KF, it improves the overall performance and also reduces the time delay and transient current which occurs in the conventional technique. Furthermore, this technique uses every measurement that the system has to further improve on the results by giving a better estimate at each time epoch. The significant improvement can be observed at the total harmonic distortion (THD) reduction at 2.38% compared to when the shunt APF is not implemented at all which performs at 168.39%. The TDH of the source current after the compensation is at 2.18% which is way below the IEEE 519 Standard which imposed a limit at less than 5% of the overall harmonics. In fact, for comparison, the use of KF also performed better than the established low-pass filter, which performs at 2.8% of the THD.

### **2. Mathematical formulation**

There are three elements that involved in generating the required current reference that is used to compensate the undesirable load current components as shown in **Figure 2**. These elements are stationary reference frame, Kalman filter (KF) and DC voltage regulator. The mathematical formulation for each element is further explained in the next subtopics.

#### **2.1. Stationary reference frame**

Stationary reference frame also known as *d-q* algorithm was developed based on Park transformation. This method transforms three-phase into *d-q* coordinates (rotating reference frame with fundamental frequency) using Park transformations. In this case, the load currents are measured and transformed into *d-q* coordinates. The equations to transform *a-b-c* coordinate into *α-β-0* coordinate is presented in Eq. (1):

$$
\begin{bmatrix}
\dot{i}\_o \\
\dot{i}\_\alpha \\
\dot{i}\_\beta
\end{bmatrix} = \sqrt{\frac{2}{3}} \begin{bmatrix}
\sqrt{\varepsilon\_2} & \sqrt{\varepsilon\_2} & \sqrt{\varepsilon\_2} \\
1 & ^{-1/2}\_{-1/2} & ^{-1/2}\_{-1/2} \\
0 & ^{\sqrt{3}/2}\_{\sqrt{2}} & ^{-\sqrt{3}/2}\_{-\sqrt{3}}
\end{bmatrix} \begin{bmatrix}
\dot{i}\_{La} \\
\dot{i}\_{Lb} \\
\dot{i}\_{Lc}
\end{bmatrix} \tag{1}
$$

**2.2. DC voltage regulator**

**2.3. Kalman filter**

and

system and *yk*

defined as.

**Figure 3.** Kalman filter.

KF as depicted in **Figure 3**.

*xk* <sup>=</sup> [

PI controller with two gains: *Kp*

The DC voltage regulator is controlled with a traditional PI controller. The DC voltage, *Vdc*, is measured and then compared with a constant reference value *Vdc\**. The error is processed by a

The use of Kalman filter (KF) provides an efficient computational means to estimate the state of a process which is able to minimize the means of the squared error. This is achieved by keeping tracks of the estimated state of the system as well as the variance of the estimates via

*xk* = *A xk*−1 + *B uk* + *wk* (3)

*zk* = *H xk* + *vk* (4)

cess noise vector and observation noise vector, respectively, and it is assumed to be mutually independent and normally distributed. Relative to the system, since the fundamental positive sequence components of the non-linear load current appears as DC quantities of the synchronous reference frame rotating at 50 Hz, it can then be separated from the load currents using

In this case, the state transition matrix is the differential equation that relates the state at

*i d*(*k*) *i*

the measurement or sometimes called observation vector. *wk*

where *A* is the state transition matrix, *B* is the control matrix that is applied to *uk*

control vector of the system, and *H* is defined as observation matrix with *xk*

the previous time step *k* − 1 to the current step *k*. Therefore, the state vector *xk*

. Both gains are calculated and tuned accordingly to the

Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF)

and 91 for *Ki*

http://dx.doi.org/10.5772/intechopen.72467

.

189

, which is the

are the pro-

can be further

the state of the

and *vk*

*<sup>q</sup>*(*k*)] (5)

and *Ki*

dynamic response in which the values of both gain are set to 4 for *Kp*

two distinct phases: predict and update. The basic KF can be defined as.

By employing Park transformation, the *α-β-0* coordinate is transformed into *d-q* coordinate as shown in Eq. (1):

$$
\begin{bmatrix}
\dot{i}\_d\\\dot{i}\_q
\end{bmatrix} = \begin{bmatrix}
\cos\theta & \sin\theta\\-\sin\theta & \cos\theta
\end{bmatrix} \begin{bmatrix}
\dot{i}\_a\\\dot{i}\_\beta
\end{bmatrix} \tag{2}
$$

where *θ* = *tan*−1 ( *i β*\_\_ *i α* ).

The phase angle, *θ*, in *d-q* frame is the same with fundamental frequency which makes the DC fundamental current component (*<sup>i</sup> d*¯ , *i q* ¯) and harmonic AC component (*<sup>i</sup> d* <sup>~</sup>, *i q* ~ ) arise due to harmonics at the load [61]. Low-pass filter (LPF) is normally used to determine the DC component. Nevertheless, for such system, phase shift in harmonics and high transient response is unavoidable before attaining its steady state. This is where KF estimator is used to replace the LPF and further improve the overall performance of the system. In order to stabilize the voltage on the DC side of the VSI, the measurement voltage, *Vdc*, measure must follow the reference voltage, *Vdc ref*. Therefore, DC voltage regulator loop is designed by integrating a suitable PI controller.

#### **2.2. DC voltage regulator**

The DC voltage regulator is controlled with a traditional PI controller. The DC voltage, *Vdc*, is measured and then compared with a constant reference value *Vdc\**. The error is processed by a PI controller with two gains: *Kp* and *Ki* . Both gains are calculated and tuned accordingly to the dynamic response in which the values of both gain are set to 4 for *Kp* and 91 for *Ki* .

#### **2.3. Kalman filter**

The use of Kalman filter (KF) provides an efficient computational means to estimate the state of a process which is able to minimize the means of the squared error. This is achieved by keeping tracks of the estimated state of the system as well as the variance of the estimates via two distinct phases: predict and update. The basic KF can be defined as.

$$\mathbf{x}\_k = A\mathbf{x}\_{k-1} + B\boldsymbol{u}\_k + \boldsymbol{w}\_k \tag{3}$$

and

estimation. Therefore the estimation resembles the combination of predictor-corrector algorithm, which is used in the system. By applying KF, it improves the overall performance and also reduces the time delay and transient current which occurs in the conventional technique. Furthermore, this technique uses every measurement that the system has to further improve on the results by giving a better estimate at each time epoch. The significant improvement can be observed at the total harmonic distortion (THD) reduction at 2.38% compared to when the shunt APF is not implemented at all which performs at 168.39%. The TDH of the source current after the compensation is at 2.18% which is way below the IEEE 519 Standard which imposed a limit at less than 5% of the overall harmonics. In fact, for comparison, the use of KF also performed better than the established low-pass filter, which performs at 2.8% of the THD.

There are three elements that involved in generating the required current reference that is used to compensate the undesirable load current components as shown in **Figure 2**. These elements are stationary reference frame, Kalman filter (KF) and DC voltage regulator. The

Stationary reference frame also known as *d-q* algorithm was developed based on Park transformation. This method transforms three-phase into *d-q* coordinates (rotating reference frame with fundamental frequency) using Park transformations. In this case, the load currents are measured and transformed into *d-q* coordinates. The equations to transform *a-b-c* coordinate

> 1 −1⁄2 −1⁄2 0 <sup>√</sup> \_\_ 3⁄2 −√ \_\_ 3⁄2 ]

By employing Park transformation, the *α-β-0* coordinate is transformed into *d-q* coordinate as

The phase angle, *θ*, in *d-q* frame is the same with fundamental frequency which makes the DC

monics at the load [61]. Low-pass filter (LPF) is normally used to determine the DC component. Nevertheless, for such system, phase shift in harmonics and high transient response is unavoidable before attaining its steady state. This is where KF estimator is used to replace the LPF and further improve the overall performance of the system. In order to stabilize the voltage on the DC side of the VSI, the measurement voltage, *Vdc*, measure must follow the reference voltage, *Vdc ref*. Therefore, DC voltage regulator loop is designed by integrating a suitable PI controller.

*cos sin* −*sin cos*][ ⎡ ⎢ ⎣

*i α i*

¯) and harmonic AC component (*<sup>i</sup>*

*i La i Lb i Lc*

⎤ ⎥ ⎦

*<sup>β</sup>*] (2)

) arise due to har-

*d* <sup>~</sup>, *i q* ~ (1)

mathematical formulation for each element is further explained in the next subtopics.

**2. Mathematical formulation**

188 Kalman Filters - Theory for Advanced Applications

**2.1. Stationary reference frame**

shown in Eq. (1):

where *θ* = *tan*−1

into *α-β-0* coordinate is presented in Eq. (1):

[

fundamental current component (*<sup>i</sup>*

( *i β*\_\_ *i α* ). ⎡ ⎢ ⎣

*i o i α i β*

⎤ ⎥ ⎦ = √ \_\_ \_\_2 3[ 1⁄√ \_\_ <sup>2</sup> 1⁄√ \_\_ <sup>2</sup> 1⁄√ \_\_ 2

*i d i <sup>q</sup>*] <sup>=</sup> [

*d*¯ , *i q*

$$\mathbf{z}\_k = H \,\mathbf{x}\_k + \mathbf{v}\_k \tag{4}$$

where *A* is the state transition matrix, *B* is the control matrix that is applied to *uk* , which is the control vector of the system, and *H* is defined as observation matrix with *xk* the state of the system and *yk* the measurement or sometimes called observation vector. *wk* and *vk* are the process noise vector and observation noise vector, respectively, and it is assumed to be mutually independent and normally distributed. Relative to the system, since the fundamental positive sequence components of the non-linear load current appears as DC quantities of the synchronous reference frame rotating at 50 Hz, it can then be separated from the load currents using KF as depicted in **Figure 3**.

In this case, the state transition matrix is the differential equation that relates the state at the previous time step *k* − 1 to the current step *k*. Therefore, the state vector *xk* can be further defined as.

$$\mathbf{x}\_k = \begin{bmatrix} \dot{i}\_{d\psi} \\ \dot{i}\_{q\dot{\psi}} \end{bmatrix} \tag{5}$$

**Figure 3.** Kalman filter.

Furthermore, the optional control input which defined the control matrix B can be neglected. Since the system only measured two parameters, respectively, therefore the measurement matrix can be simply represented by two-by-two identity matrix. Therefore, the implementation of the KF can be rewritten as.

\*\*1001\*\*: \*\*C\*\*: \*\*1001\*\*: \*\*C\*\*: \*\*1001\*\*: \*\*n $-n$ .

$$\begin{aligned} x\_k &= A\_k x\_{k-1} + w\_k \\ P\_k^- &= A\_k P\_{k-1} A^T + Q \end{aligned} \tag{6}$$

**3. Simulation and experimental result**

**3.1. Non-linear load**

by using the formula defined as.

%*THD* <sup>=</sup> <sup>√</sup>

**Figure 6.** Harmonic spectrum before the compensation.

respectively.

The proposed simulation and experimental results designed for the three-phase reference current generation employing KF estimator for three-phase shunt APF are presented. The

Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF)

http://dx.doi.org/10.5772/intechopen.72467

The results for the APF before and after compensation are simulated using Matlab Simulink, while Fluke Power Quality Analyzer captures the results for the experimental. **Figure 5(a)** and **(b)** shows the supply current waveform before the compensation for simulation and experimental result; thus, the harmonic spectrum of both simulation and experimental is shown in **Figure 6**,

From the harmonic spectrum results, the total harmonic distortion (THD) can be determined

\_\_\_\_\_ ∑ *n*=2 ∞ *I h* 2 \_\_\_\_\_ *I f*

(9)

191

work is simulated and implemented using Matlab Simulink and dSPACE.

**Figure 5.** Simulation and experimental result without shunt APF (a) simulation and (b) Experiment.

for the predictor:

$$\begin{aligned} \mathbf{z}\_k &= H \hat{\mathbf{x}}\_k^- + \mathbf{v}\_k \\ \mathbf{S}\_k &= H P\_k^- H^T + \mathbf{R}\_k \\ \hat{\mathbf{x}}\_k &= \hat{\mathbf{x}}\_k^- + \mathbf{K}\_k \{ \mathbf{z}\_k - H \hat{\mathbf{x}}\_k^- \} \end{aligned} \tag{7}$$

The measurement update equation *x*̂ *k* is the estimate reference current of *i d* and *i q* , *x*̂ *k* − is the predicted state, *zk* is the measurement of actual current, *Pk* is the estimate error covariance, *Rk* is the observation covariance matrix and *Kk* is the Kalman gain. In this representation, matrix *P* is the variance matrix of the error *xk* <sup>−</sup> *<sup>x</sup>*̂ *k* where the goal is to minimize this value. Here the Kalman gain calculation will be based on the conventional calculation defined in Eq. (8):

ganın uchununon wının  $\mathbf{e}$  uchun on üre uchunun  $\mathbf{e}$  uchunun  $\mathbf{e}$  uchunun  $\mathbf{u}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$  uchun  $\mathbf{e}$ 

The process noise covariance matrix *Q* and observation noise covariance matrix *R* are tuned manually in order to achieve the optimal performance of the filter. **Figure 4** shows the cycle of KF.

**Figure 4.** Cycle of discrete Kalman filter.

### **3. Simulation and experimental result**

The proposed simulation and experimental results designed for the three-phase reference current generation employing KF estimator for three-phase shunt APF are presented. The work is simulated and implemented using Matlab Simulink and dSPACE.

### **3.1. Non-linear load**

Furthermore, the optional control input which defined the control matrix B can be neglected. Since the system only measured two parameters, respectively, therefore the measurement matrix can be simply represented by two-by-two identity matrix. Therefore, the implementa-

*Pk*

*zk* = *H x*̂ *k* <sup>−</sup> + *vk Sk* <sup>=</sup> *<sup>H</sup> Pk* <sup>−</sup> *HT* <sup>+</sup> *Rk*

<sup>−</sup> + *Kk*(*zk* − *H x*̂

<sup>−</sup> *HT S*−1 *Pk* <sup>=</sup> (*<sup>I</sup>* <sup>−</sup> *Kk <sup>H</sup>*) *Pk*

The process noise covariance matrix *Q* and observation noise covariance matrix *R* are tuned manually in order to achieve the optimal performance of the filter. **Figure 4** shows the cycle

*k* − )

is the estimate reference current of *i*

*x*̂ *<sup>k</sup>* = *x*̂ *k*

*k*

gain calculation will be based on the conventional calculation defined in Eq. (8):

is the measurement of actual current, *Pk*

*k*

<sup>−</sup> <sup>=</sup> *<sup>A</sup> Pk*−1 *AT* <sup>+</sup> *<sup>Q</sup>* (6)

*d* and *i q* , *x*̂ *k* −

<sup>−</sup> (8)

is the estimate error covariance, *Rk*

is the Kalman gain. In this representation, matrix *P* is the

where the goal is to minimize this value. Here the Kalman

(7)

is the pre-

is the

tion of the KF can be rewritten as.

190 Kalman Filters - Theory for Advanced Applications

The measurement update equation *x*̂

observation covariance matrix and *Kk*

*Kk* <sup>=</sup> *Pk*

variance matrix of the error *xk* <sup>−</sup> *<sup>x</sup>*̂

**Figure 4.** Cycle of discrete Kalman filter.

for the predictor:

dicted state, *zk*

of KF.

*xk* <sup>=</sup> *<sup>A</sup> xk*−1 <sup>+</sup> *wk*

The results for the APF before and after compensation are simulated using Matlab Simulink, while Fluke Power Quality Analyzer captures the results for the experimental. **Figure 5(a)** and **(b)** shows the supply current waveform before the compensation for simulation and experimental result; thus, the harmonic spectrum of both simulation and experimental is shown in **Figure 6**, respectively.

From the harmonic spectrum results, the total harmonic distortion (THD) can be determined by using the formula defined as.

$$\text{\textbullet\textbullet}\_{\text{\textbullet}} \text{\textbullet}\_{\text{\textbullet}} = \frac{\sqrt{\sum\_{h} l\_{h}^{2}}}{l\_{f}} \tag{9}$$

**Figure 5.** Simulation and experimental result without shunt APF (a) simulation and (b) Experiment.

**Figure 6.** Harmonic spectrum before the compensation.

**Figure 7.** Simulation Butterworth low-pass filter.

where *Ih* is the harmonic component, *I f* is the fundamental component, and *n* is the harmonic number: 2, 3, 4, etc.

performance of the active power filter (APF). But, time delay which contributes to the phase shift in harmonics and high transient current is the common effect when applying the LPF. **Figure 7** shows the shunt APF when applying Butterworth LPF. It is clearly shown that from the figure, the time delay is recorded at 0.02 s with 43.36% of the THD. On the other hand, there is no time delay when applying the shunt APF using KF estimator which is shown in **Figure 8**. Therefore, the THD produced by the KF estimator is 98% improvement compared to LPF. On the other hand, the experimental results for low-pass and KF are shown in **Figures 9**

Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF)

http://dx.doi.org/10.5772/intechopen.72467

193

and **10**, respectively.

**Figure 10.** Experiment of Kalman filter.

**Figure 9.** Experimental Butterworth low-pass filter.

Therefore, the THD of the line current obtained by the simulation is 56.14%, while the experimental obtains about 47.26%. There are slightly different between the simulation and experimental results because the simulation is simulated at an ideal condition.

#### **3.2. Kalman filter estimator result versus low-pass filter**

Commonly, a Butterworth low-pass filter (LPF) was applied to filter out the unwanted DC component for *d-q* algorithm to ensure that the correct reference currents are generated in the system. Failure to obtain the correct reference current resolves reduction of the overall

**Figure 8.** Simulation of Kalman filter.

**Figure 9.** Experimental Butterworth low-pass filter.

where *Ih*

number: 2, 3, 4, etc.

**Figure 8.** Simulation of Kalman filter.

is the harmonic component, *I*

**Figure 7.** Simulation Butterworth low-pass filter.

192 Kalman Filters - Theory for Advanced Applications

*f*

mental results because the simulation is simulated at an ideal condition.

**3.2. Kalman filter estimator result versus low-pass filter**

Therefore, the THD of the line current obtained by the simulation is 56.14%, while the experimental obtains about 47.26%. There are slightly different between the simulation and experi-

Commonly, a Butterworth low-pass filter (LPF) was applied to filter out the unwanted DC component for *d-q* algorithm to ensure that the correct reference currents are generated in the system. Failure to obtain the correct reference current resolves reduction of the overall

is the fundamental component, and *n* is the harmonic

performance of the active power filter (APF). But, time delay which contributes to the phase shift in harmonics and high transient current is the common effect when applying the LPF. **Figure 7** shows the shunt APF when applying Butterworth LPF. It is clearly shown that from the figure, the time delay is recorded at 0.02 s with 43.36% of the THD. On the other hand, there is no time delay when applying the shunt APF using KF estimator which is shown in **Figure 8**. Therefore, the THD produced by the KF estimator is 98% improvement compared to LPF. On the other hand, the experimental results for low-pass and KF are shown in **Figures 9** and **10**, respectively.

**Figure 10.** Experiment of Kalman filter.

#### **3.3. Three-phase shunt active power filter**

In this shunt active power filter feeding a non-linear load, the results between Butterworth LPF and KF estimator are compared between simulation and experimental, which are shown in **Figures 11** and **12**, respectively.

From the results obtained, it can be concluded that almost the same waveform was produced for both simulation and experimental approaches. Furthermore, the harmonic spectrum form the experimental is shown in **Figure 13**.

Both methods have demonstrated a harmonic reduction with almost identical fundamental current between simulation and experimental.

The THD results obtained show that the new technique shunt APF abides the regulation of IEEE 519–1992 standard. **Tables 1** and **2** show the THD after simulation and experimental results, respectively.

From the observation, the shunt APF using KF estimator technique produces about 0.1% better THD than LPF either in simulation or experimental.

**Figure 11.** Simulation result for shunt APF (a) low-pass filter and (b) Kalman filter.

**3.4. Operation with three-phase induction motor speed drive**

**Figure 13.** Harmonics spectrum.

**Table 1.** Simulation results.

**Table 2.** Experimental results.

A 1.5 kW, 380 V variable speed induction motor (IM) drive is connected in parallel to the APF and the three-phase supply voltages. The motor is operated as a non-linear load and starts to accelerate from standstill at time, *t* = 0.06 s, until it reached the required reference speed which is set at 1400 rpm. **Figures 14** and **15** show the supply current when the IM starts to accelerate without and with shunt APF, respectively. Fluke Power Quality Analyzer was used to measure

**Types of reference current generation THD before (%) THD after (%)**

Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF)

http://dx.doi.org/10.5772/intechopen.72467

195

**Types of reference current generation THD before (%) THD after (%)**

Low-pass filter 55.88 2.09 Kalman filter estimator 1.99

Low-pass filter 47.27 2.30 Kalman filter estimator 2.18

**Figure 12.** Experimental results for shunt APF: (a) low-pass filter and (b) Kalman filter.

**Figure 13.** Harmonics spectrum.

**3.3. Three-phase shunt active power filter**

in **Figures 11** and **12**, respectively.

194 Kalman Filters - Theory for Advanced Applications

results, respectively.

the experimental is shown in **Figure 13**.

current between simulation and experimental.

ter THD than LPF either in simulation or experimental.

**Figure 12.** Experimental results for shunt APF: (a) low-pass filter and (b) Kalman filter.

**Figure 11.** Simulation result for shunt APF (a) low-pass filter and (b) Kalman filter.

In this shunt active power filter feeding a non-linear load, the results between Butterworth LPF and KF estimator are compared between simulation and experimental, which are shown

From the results obtained, it can be concluded that almost the same waveform was produced for both simulation and experimental approaches. Furthermore, the harmonic spectrum form

Both methods have demonstrated a harmonic reduction with almost identical fundamental

The THD results obtained show that the new technique shunt APF abides the regulation of IEEE 519–1992 standard. **Tables 1** and **2** show the THD after simulation and experimental

From the observation, the shunt APF using KF estimator technique produces about 0.1% bet-


**Table 1.** Simulation results.


**Table 2.** Experimental results.

#### **3.4. Operation with three-phase induction motor speed drive**

A 1.5 kW, 380 V variable speed induction motor (IM) drive is connected in parallel to the APF and the three-phase supply voltages. The motor is operated as a non-linear load and starts to accelerate from standstill at time, *t* = 0.06 s, until it reached the required reference speed which is set at 1400 rpm. **Figures 14** and **15** show the supply current when the IM starts to accelerate without and with shunt APF, respectively. Fluke Power Quality Analyzer was used to measure the THD at steady-state condition (*t =* 0.28 s). THD obtained before applying shunt APF is 168.39%, whereas when applying the shunt APF using both KF and LFP, the THD reduced to 2.38 and 2.80%. Furthermore, the harmonic spectrum with or without shunt APF for both KF and LPF is shown in **Figures 16**–**18**, respectively. It can conclude that from the results, the shunt APF employing KF-based estimator produced lower THD than LPF for an induction motor drive application.

The overall total harmonic distortion with or without shunt APF is shown tabulated in **Table 3**, which shows that the KF estimator produces lower THD than LPF for three-phase

Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF)

http://dx.doi.org/10.5772/intechopen.72467

197

induction motor.

**Figure 16.** Harmonic spectrum without shunt APF.

**Figure 17.** Harmonics spectrum after applying shunt APF.

**Figure 14.** Supply current waveform without shunt APF.

**Figure 15.** Supply current when applying shunt APF with Kalman filter estimator.

**Figure 16.** Harmonic spectrum without shunt APF.

the THD at steady-state condition (*t =* 0.28 s). THD obtained before applying shunt APF is 168.39%, whereas when applying the shunt APF using both KF and LFP, the THD reduced to 2.38 and 2.80%. Furthermore, the harmonic spectrum with or without shunt APF for both KF and LPF is shown in **Figures 16**–**18**, respectively. It can conclude that from the results, the shunt APF employing KF-based estimator produced lower THD than LPF for an induction motor

**Figure 15.** Supply current when applying shunt APF with Kalman filter estimator.

**Figure 14.** Supply current waveform without shunt APF.

drive application.

196 Kalman Filters - Theory for Advanced Applications

The overall total harmonic distortion with or without shunt APF is shown tabulated in **Table 3**, which shows that the KF estimator produces lower THD than LPF for three-phase induction motor.

**1.** The hysteresis band plays a significant influence to the THD reduction. With the lower hysteresis band, more accurate PWM generated, thereby improving the THD. In this thesis, the hysteresis band is set to ±0.08 while in simulation at ±0.001. Therefore, in order to have faster and more accurate result, the combination of dSpace and FPGA can be implemented

Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF)

http://dx.doi.org/10.5772/intechopen.72467

199

**2.** Further combination between shunt APF and passive filter (hybrid APF) can be used to improve the performance of the APF. Where the passive filter used to filter the higher order

**3.** Apply the Kalman filter to the instantaneous real and reactive algorithm. Conventionally, the technique used *p-q* algorithm that combined with high-pass filter. However, Kalman filter can be used to replace it, and the performance of the system can be further

1 Department of Electrical and Electronics Engineering, Faculty of Engineering, National

2 Faculty of Electrical Engineering, Universiti Teknikal Malaysia Melaka, Melaka, Malaysia

[1] Duffey CK, Stratford RP. Update of harmonic standard IEEE-519-IEEE recommended practices and requirements for harmonic control in electric power systems. In: Conference Record of the IEEE Industry Applications Society Annual Meeting. 1989;**2**:1618-1624 [2] Shuter TC, Vollkommer HT, Kirkpatrick TL. Survey of harmonic levels on the American electric power distribution system. IEEE Transactions on Power Delivery.

[3] Subjak JS Jr, McQuilkin JS. Harmonics-causes, effects, measurements, and analysis: An

[4] Beides HM, Heydt GT. Dynamic state estimation of power system harmonics using Kalman filter methodology. IEEE Transactions on Power Delivery. 1991;**6**:1663-1670 [5] Gruzs TM. Uncertainties in compliance with harmonic current distortion limits in electric power systems. IEEE Transactions on Industry Applications. 1991;**27**:680-685

[6] Emanuel AE, Orr JA, Cyganski D, Gulachenski EM. A survey of harmonic voltages and currents at the customer's bus. IEEE Transactions on Power Delivery. 1993;**8**:411-421

update. IEEE Transactions on Industry Applications. 1990;**26**:1034-1042

\*, Syed Mohd Fairuz Bin Syed Mohd Dardin1

and

to reduce the computational time from the dSpace.

\*Address all correspondence to: shukri@upnm.edu.my

Defence University of Malaysia, Kuala Lumpur, Malaysia

investigated.

**Author details**

**References**

1989;**4**:2204-2213

Zulkifilie Bin Ibrahim2

Ahmad Shukri Bin Abu Hasim1

harmonic while shunt APF filter the lower order harmonic.

**Figure 18.** Harmonics spectrum without fundamental.


**Table 3.** THD of supply current before and after applying shunt APF.

### **4. Conclusion**

The new techniques of reference current generator by employing KF estimator for shunt APF technique have been presented. The results of the proposed technique in generating the three-phase reference current towards reducing the THD are established using simulation and experimental. For the three-phase rectifier connected with RC load, the performance of the proposed technique is comparable with those based on the LPF reference current generation. The THD of the source current from the experimental result after the compensation is 2.18% which is less than 5% of the harmonic limit imposed by the IEEE 519 standard. Furthermore, nearly 0.1% THD improvement was gained by the proposed techniques compared to LPF. Thus, the comparison of different reference current grid generations for shunt APF is also presented. The performance of KF estimator reference current generation was also studied for induction motor variable speed drive. In induction motor, almost 0.42% improvement of THD was gathered when applying KF estimator compared to LPF.

### **5. Future works**

Although the proposed technique improved the overall performance of shunt APF, there is a room of improvement and suggestion for further research work such as:


### **Author details**

Ahmad Shukri Bin Abu Hasim1 \*, Syed Mohd Fairuz Bin Syed Mohd Dardin1 and Zulkifilie Bin Ibrahim2

\*Address all correspondence to: shukri@upnm.edu.my

1 Department of Electrical and Electronics Engineering, Faculty of Engineering, National Defence University of Malaysia, Kuala Lumpur, Malaysia

2 Faculty of Electrical Engineering, Universiti Teknikal Malaysia Melaka, Melaka, Malaysia

### **References**

**4. Conclusion**

**5. Future works**

The new techniques of reference current generator by employing KF estimator for shunt APF technique have been presented. The results of the proposed technique in generating the three-phase reference current towards reducing the THD are established using simulation and experimental. For the three-phase rectifier connected with RC load, the performance of the proposed technique is comparable with those based on the LPF reference current generation. The THD of the source current from the experimental result after the compensation is 2.18% which is less than 5% of the harmonic limit imposed by the IEEE 519 standard. Furthermore, nearly 0.1% THD improvement was gained by the proposed techniques compared to LPF. Thus, the comparison of different reference current grid generations for shunt APF is also presented. The performance of KF estimator reference current generation was also studied for induction motor variable speed drive. In induction motor, almost 0.42% improve-

**Reference current generation THD before (%) THD after (%)**

Low-pass filter 168.39 2.80 Kalman filter estimator 2.38

**Table 3.** THD of supply current before and after applying shunt APF.

**Figure 18.** Harmonics spectrum without fundamental.

198 Kalman Filters - Theory for Advanced Applications

Although the proposed technique improved the overall performance of shunt APF, there is a

ment of THD was gathered when applying KF estimator compared to LPF.

room of improvement and suggestion for further research work such as:


[7] Henderson RD, Rose PJ. Harmonics: The effects on power quality and transformers. IEEE Transactions on Industry Applications. 1994;**30**:528-532

[21] Rahim NA, Mekhilef S, Zahrul I. A single-phase active power filter for harmonic compensation. In: IEEE International Conference on Industrial Technology, 2005. ICIT 2005.

Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF)

http://dx.doi.org/10.5772/intechopen.72467

201

[22] Wang X, Liu J, Yuan C, Wang Z. A comparative study on voltage-source control and current-source control of series active power filter. In: Applied Power Electronics Conference and Exposition, 2006. APEC '06. Twenty-First Annual IEEE. 2006. 6 pp

[23] Chang GW, Chen SK, Chin YC, Chen WC. An a-b-c reference frame-based compensation strategy for series active power filter control. In: 2006 1ST IEEE Conference on Industrial

[24] Hurng-Liahng J, Jinn-Chang W, Yao-Jen C, Ya-Tsung F. A novel active power filter for harmonic suppression. IEEE Transactions on Power Delivery. 2005;**20**:1507-1513

[25] Moran LT, Mahomar JJ, Dixon JR. Careful connection. IEEE Industry Applications

[26] El-Habrouk M, Darwish MK, Mehta P. Active power filters: A review. IEE Proceedings—

[27] Golwala H, Chudamani R. Comparative study of switching signal generation techniques for three phase four wire shunt active power filter. In: Electric Machines & Drives

[28] Zeng FP, Tan GH, Wang JZ, Ji YC. Novel single-phase five-level voltage-source inverter

[29] Rahmani S, Mendalek N, Al-Haddad K. Experimental design of a nonlinear control technique for three-phase shunt active power filter. IEEE Transactions on Industrial

[30] Golwala H, Chudamani R. Simulation of three-phase four-wire shunt active power filter using novel switching technique. In: 2010 Joint International Conference on Power Electronics, Drives and Energy Systems (PEDES) & 2010 Power India. 2010. pp. 1-7

[31] Abdalla II, Rao KSR, Perumal N. Harmonics mitigation and power factor correction with a modern three-phase four-leg shunt active power filter. In: 2010 IEEE International

[32] Vodyakho O, Mi CC. Three-level inverter-based shunt active power filter in threephase three-wire and four-wire systems. IEEE Transactions on Power Electronics.

[33] Vodyakho O, Kim T, Kwak S, Edrington CS. Comparison of the space vector current controls for shunt active power filters. IET Power Electronics. 2009;**2**:653-664

[34] Vodyakho O, Kim T. Shunt active filter based on three-level inverter for three-phase

Conference (IEMDC), 2011 IEEE International. 2011. pp. 1409-1414

Conference on Power and Energy (PECon). 2010. pp. 156-161

four-wire systems. IET Power Electronics. 2009;**2**:216-226

for the shunt active power filter. IET Power Electronics. 2010;**3**:480-489

2005. pp. 1075-1079

Magazine. 2004;**10**:43-50

Electronics. 2010;**57**:3364-3375

2009;**24**:1350-1363

Electronics and Applications. 2006. pp. 1-4

Electric Power Applications. 2000;**147**:403-413


[21] Rahim NA, Mekhilef S, Zahrul I. A single-phase active power filter for harmonic compensation. In: IEEE International Conference on Industrial Technology, 2005. ICIT 2005. 2005. pp. 1075-1079

[7] Henderson RD, Rose PJ. Harmonics: The effects on power quality and transformers.

[8] Mansoor A, Grady WM, Staats PT, Thallam RS, Doyle MT, Samotyj MJ. Predicting the net harmonic currents produced by large numbers of distributed single-phase computer

[9] Clark SL, Famouri P, Cooley WL. Elimination of supply harmonics. IEEE Industry

[10] Kazerani M, Ziogas PD, Joos G. A novel active current waveshaping technique for solid-state input power factor conditioners. IEEE Transactions on Industrial Electronics.

[11] Enjeti P, Shireen W, Pitel I. Analysis and design of an active power filter to cancel harmonic currents in low voltage electric power distribution systems. In: Proceedings of the 1992 International Conference on Industrial Electronics, Control, Instrumentation, and

[12] Duke RM, Round SD. The steady-state performance of a controlled current active filter.

[13] Round SD, Mohan N. Comparison of frequency and time domain neural network controllers for an active power filter. In: Proceedings of the IECON '93, International Conference

[14] Jou HL, Wu JC, Chu HY. New single-phase active power filter. IEE Proceedings—Electric

[15] Nastran J, Cajhen R, Seliger M, Jereb P. Active power filter for nonlinear AC loads. IEEE

[16] Jae-Ho C, Ga-Woo P, Dewan SB. Standby power supply with active power filter ability using digital controller. In: Applied Power Electronics Conference and Exposition, 1995.

[17] Torrey DA, Al-Zamel AMAM. Single-phase active power filters for multiple nonlinear

[18] Kim YJ, Kim JS, Kim YS. Single-phase active power filter based on rotating reference frame method. In: Proceedings of the Eighth International Conference on Electrical

[19] Khan MM, Feng JF, Chen C, Zhiming W. Single-phase dynamically decoupled active power filter for system-integrated application. IEEE Proceedings—Electric Power

[20] Rudnick H, Dixon J, Moran L. Delivering clean and pure power. IEEE Power and Energy

APEC '95. Conference Proceedings 1995, Tenth Annual. 1995;**2**:783-789

loads. IEEE Transactions on Power Electronics. 1995;**10**:263-272

Machines and Systems, 2005. ICEMS 2005. 2005. pp. 1428-1431

on Industrial Electronics, Control, and Instrumentation; 1993. 1993;**2**:1099-1104

IEEE Transactions on Industry Applications. 1994;**30**:528-532

loads. IEEE Transactions on Power Delivery. 1995;**10**:2001-2006

IEEE Transactions on Power Electronics. 1993;**8**:140-146

Applications Magazine. 1997;**3**:62-67

200 Kalman Filters - Theory for Advanced Applications

Automation. 1992;**1**:368-373

Power Applications. 1994;**141**:129-134

Applications. 2006;**153**:625-631

Magazine. 2003;**1**:32-40

Transactions on Power Electronics. 1994;**9**:92-96

1991;**38**:72-78


[35] Uyyuru KR, Mishra MK, Ghosh A. An optimization-based algorithm for shunt active filter under distorted supply voltages. IEEE Transactions on Power Electronics. 2009;**24**:1223-1232

[50] Pal Y, Swarup A, Singh B. A review of compensating type custom power devices for power quality improvement. In: Joint International Conference on Power System Technology and IEEE Power India Conference, 2008. POWERCON 2008. 2008. pp. 1-8

Kalman Filters for Reference Current Generation in Shunt Active Power Filter (APF)

http://dx.doi.org/10.5772/intechopen.72467

203

[51] Kazemi A, Sarlak M, Barkhordary M. An adaptive noise canceling method for singlephase unified power quality conditioner. In: 2006 1ST IEEE Conference on Industrial

[52] Herrera RS, Salmeron P, Vazquez JR, Litran SP, Perez A. GENERALIZED instantaneous reactive power theory in poly-phase power systems. In: 13th European Conference on

[53] Kelesidis K, Adamidis G, Tsengenes G. Investigation of a control scheme based on modified p-q theory for single phase single stage grid connected PV system. In: 2011

[54] Mulla MA, Rajagopalan C, Chowdhury A. Hardware implementation of series hybrid active power filter using a novel control strategy based on generalised instantaneous

[55] Esfandiari A, Parniani M, Mokhtari H. A new control strategy of shunt active filters for power quality improvement of highly and randomly varying loads. In: 2004 IEEE

[56] Sawant RR, Chandorkar MC. Methods for multi-functional converter control in three-

[57] Salmeron P, Herrera RS, Vazquez JR. Mapping matrices against vectorial frame in the instantaneous reactive power compensation. IET Electric Power Applications. 2007;**1**:727-736 [58] Bhattacharya A, Chakraborty C, Bhattacharya S. Parallel-connected shunt hybrid active power filters operating at different switching frequencies for improved performance.

[59] Suresh Y, Panda AK, Suresh M. Real-time implementation of adaptive fuzzy hysteresis-band current control technique for shunt active power filter. IET Power Electronics.

[60] Quoc-Nam T, Hong-Hee L. An advanced current control strategy for three-phase shunt active power filters. IEEE Transactions on Industrial Electronics. 2013;**60**:5400-5410 [61] Asiminoaei L, Blaabjerg F, Hansen S. Detection is key—harmonic detection methods for active power filter applications. IEEE Industry Applications Magazine. 2007;**13**:22-33

[62] Jain SK, Singh SN. Harmonics estimation in emerging power system: Key issues and

[63] Moreno VM, Lopez AP, Garcias RID. Reference current estimation under distorted line voltage for control of shunt active power filters. IEEE Transactions on Power Electronics.

International Conference on Clean Electrical Power (ICCEP). 2011. pp. 535-540

Power Electronics and Applications, 2009. EPE '09. 2009. pp. 1-10

International Symposium on Industrial Electronics. 2004;**2**:1297-1302

phase four-wire systems. IET Power Electronics. 2009;**2**:52-66

IEEE Transactions on Industrial Electronics. 2012;**59**:4007-4019

challenges. Electric Power Systems Research. 2011;**81**:1754-1766

2012;**5**:1188-1195

2004;**19**:988-994

power theory. IET Power Electronics. 2013;**6**:592-600

Electronics and Applications. 2006. pp. 1-6


[50] Pal Y, Swarup A, Singh B. A review of compensating type custom power devices for power quality improvement. In: Joint International Conference on Power System Technology and IEEE Power India Conference, 2008. POWERCON 2008. 2008. pp. 1-8

[35] Uyyuru KR, Mishra MK, Ghosh A. An optimization-based algorithm for shunt active filter under distorted supply voltages. IEEE Transactions on Power Electronics.

[36] Singh B, Solanki J. An implementation of an adaptive control algorithm for a three-phase shunt active filter. IEEE Transactions on Industrial Electronics. 2009;**56**:2811-2820

[37] Rahmani S, Hamadi A, Mendalek N, Al-Haddad K. A new control technique for three-phase shunt hybrid power filter. IEEE Transactions on Industrial Electronics.

[38] Peng X, Venayagamoorthy GK, Corzine KA. Seven-level shunt active power filter for high-power drive systems. IEEE Transactions on Power Electronics. 2009;**24**:6-13

[39] Akagi H. The state-of-the-art of active filters for power conditioning. In: 2005 European

[41] Pini SH, Barbi I. A single-phase high-power-factor rectifier, based on a two-quadrant shunt active filter. IEEE Transactions on Power Electronics. 2011;**26**:3131-3143

[42] Karuppanan P, Mahapatra K. PLL with PI, PID and fuzzy logic controllers based shunt active power line conditioners. In: 2010 Joint International Conference on Power Electronics, Drives and Energy Systems (PEDES) & 2010 Power India. 2010. pp. 1-6 [43] Singh B, Verma V, Chandra A, Al-Haddad K. Hybrid filters for power quality improvement. IEE Proceedings—Generation, Transmission and Distribution. 2005;**152**:365-378

[44] Zhao W, Luo A, Shen ZJ, Wu C. injection-type hybrid active power filter in high-power grid with background harmonic voltage. IET Power Electronics. 2011;**4**:63-71

[45] A. Varschavsky, J. Dixon, M. Rotella, Moran L. Cascaded nine-level inverter for hybridseries active power filter, using industrial controller, IEEE Transactions on Industrial

[46] Salmeron P, Litran SP. Improvement of the electric power quality using series active and shunt passive filters. IEEE Transactions on Power Delivery. 2010;**25**:1058-1067

[47] Luo A, Shuai Z, Zhu W, Shen ZJ, Tu C. Design and application of a hybrid active power

[48] Litran SP, Salmeron P, Vazquez JR, Herrera RS, Perez A. Control strategy for hybrid power filter to compensate unbalanced and non-linear, three-phase loads. In: 13th European Conference on Power Electronics and Applications, 2009. EPE '09. 2009.

[49] An L, Ci T, Zhi Kang S, Wei Z, Fei R, Ke Z. A novel three-phase hybrid active power filter with a series resonance circuit tuned at the fundamental frequency. IEEE Transactions

filter with injection circuit. IET Power Electronics. 2010;**3**:54-64

on Industrial Electronics. 2009;**56**:2431-2440

Conference on Power Electronics and Applications. 2005. pp. 15

[40] Akagi H. Active harmonic filters. Proceedings of the IEEE. 2005;**93**:2128-2141

2009;**24**:1223-1232

202 Kalman Filters - Theory for Advanced Applications

2009;**56**:2904-2915

Electronics; 2010;**57**:2761-2767

pp. 1-10


[64] Petit JF, Robles G, Amaris H. Predictive algorithm for harmonic mitigation in non-linear loads based on active filters. In: Power Tech, 2005 IEEE Russia. 2005. pp. 1-6

**Chapter 11**

Provisional chapter

**Applications of Kalman Filters for Coherent Optical**

DOI: 10.5772/intechopen.71617

Applications of Kalman Filters for Coherent Optical

In this chapter, we review various applications of Kalman filtering for coherent optical communication systems. First, we briefly discuss the principles of Kalman filter and its variations including extended Kalman filter (EKF) and adaptive Kalman filter (AKF). Later on, we illustrate the applicability of Kalman filters for joint tracking of several optical transmission impairments, simultaneously, by formulating the state space model (SSM) and detailing the principles. A detailed methodology is presented for the joint tracking of linear and nonlinear phase noise along with amplitude noise using EKF. Also, approaches to enhance the performance obtained by EKF by combining with other existing digital signal processing (DSP) techniques are presented. Frequency and phase offset estimation using a two stage linear Kalman filter (LKF)/EKF is also discussed. A cascaded structure of LKF and EKF by splitting the SSM to jointly mitigate the effects of polarization, phase and amplitude noise is also presented. The numerical analysis concludes that the Kalman filter based approaches outperform the conventional methods with better tracking capability and faster convergence besides offering more feasibility

Keywords: optical communication systems, coherent optical transmission, digital signal processing, nonlinear mitigation, phase noise, amplitude noise, QAM

In order to meet the yearning demands of bandwidth and capacity due to ever increasing data traffic, the contemporary research in the field of optical transmission, is focused on developing 400 Gbps and above, Ethernet transmission [1–5]. The achievable information rates using optical fiber as communication channel have been rapidly increased over the past few decades. Some of the technology breakthroughs behind this rapid increase in the transmission capacity, can be listed as the invention and development of the erbium doped fiber amplifiers (EDFA),

> © The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited.

© 2018 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Communication Systems**

Communication Systems

http://dx.doi.org/10.5772/intechopen.71617

for real-time implementations.

Abstract

1. Introduction

Lalitha Pakala and Bernhard Schmauss

Lalitha Pakala and Bernhard Schmauss

Additional information is available at the end of the chapter

Additional information is available at the end of the chapter


Provisional chapter

### **Applications of Kalman Filters for Coherent Optical Communication Systems** Applications of Kalman Filters for Coherent Optical

DOI: 10.5772/intechopen.71617

Lalitha Pakala and Bernhard Schmauss

Additional information is available at the end of the chapter Lalitha Pakala and Bernhard Schmauss

http://dx.doi.org/10.5772/intechopen.71617 Additional information is available at the end of the chapter

Communication Systems

#### Abstract

[64] Petit JF, Robles G, Amaris H. Predictive algorithm for harmonic mitigation in non-linear loads based on active filters. In: Power Tech, 2005 IEEE Russia. 2005. pp. 1-6

[65] Rosendo JA, Bachiller A, Gomez A. Application of self-tuned Kalman filters to control of active power filters. In: Power Tech, 2007 IEEE Lausanne. 2007. pp. 1262-1265

[66] Cardoso R, de Camargo RF, Pinheiro H, Grundling HA. Kalman filter based synchroni-

[67] Panigrahi R, Panda PC, Subudhi BD. Comparison of performances of hysteresis and dead beat controllers in active power filtering. In: 2012 IEEE Third International Conference

[68] Regulski P, Terzija V. Estimation of frequency and fundamental power components using an unscented Kalman filter. IEEE Transactions on Instrumentation and Measurement.

[69] Kanieski JM, Cardoso R, Pinheiro H, Grundling HA. Kalman filter-based control system for power quality conditioning devices. IEEE Transactions on Industrial Electronics.

sation methods. IET *Generation, Transmission & Distribution*. 2008;**2**:542-555

on Sustainable Energy Technologies (ICSET). 2012. pp. 287-292

2012;**61**:952-962

204 Kalman Filters - Theory for Advanced Applications

2013;**60**:5214-5227

In this chapter, we review various applications of Kalman filtering for coherent optical communication systems. First, we briefly discuss the principles of Kalman filter and its variations including extended Kalman filter (EKF) and adaptive Kalman filter (AKF). Later on, we illustrate the applicability of Kalman filters for joint tracking of several optical transmission impairments, simultaneously, by formulating the state space model (SSM) and detailing the principles. A detailed methodology is presented for the joint tracking of linear and nonlinear phase noise along with amplitude noise using EKF. Also, approaches to enhance the performance obtained by EKF by combining with other existing digital signal processing (DSP) techniques are presented. Frequency and phase offset estimation using a two stage linear Kalman filter (LKF)/EKF is also discussed. A cascaded structure of LKF and EKF by splitting the SSM to jointly mitigate the effects of polarization, phase and amplitude noise is also presented. The numerical analysis concludes that the Kalman filter based approaches outperform the conventional methods with better tracking capability and faster convergence besides offering more feasibility for real-time implementations.

Keywords: optical communication systems, coherent optical transmission, digital signal processing, nonlinear mitigation, phase noise, amplitude noise, QAM

### 1. Introduction

In order to meet the yearning demands of bandwidth and capacity due to ever increasing data traffic, the contemporary research in the field of optical transmission, is focused on developing 400 Gbps and above, Ethernet transmission [1–5]. The achievable information rates using optical fiber as communication channel have been rapidly increased over the past few decades. Some of the technology breakthroughs behind this rapid increase in the transmission capacity, can be listed as the invention and development of the erbium doped fiber amplifiers (EDFA),

© 2018 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited.

wavelength division multiplexing (WDM) systems, coherent detection, digital signal processing (DSP) techniques and forward error correction (FEC) schemes ensuring reliable transmission. The advent of coherent detection along with subsequent DSP made it possible to deploy spectrally efficient higher order modulation formats and multiplexing techniques [6, 7]. Moreover, it has also made feasible to digitally equalize the optical fiber transmission impairments [8], which are the main hurdle to increase the bandwidth-distance product. The transmission capacity can be increased several times by employing complex modulation formats like m-ary quadrature amplitude modulation (with m = 4, 16, 64 and so on), and multiplexing techniques like polarization division multiplexing (PDM) and WDM. However, they are more vulnerable to the optical transmission impairments as well as to the carrier phase and frequency offset (FO). Hence, effective DSP algorithms for combatting with the channel impairments were under active research over the past decade [8–23]. Consequently, coherent optical receivers are well developed and employ digital filters that allow for effective equalization of fiber linear impairments like chromatic dispersion (CD) and polarization mode dispersion (PMD) in the electric domain [9]. Typically, CD can be compensated by either frequency or time domain filters using finite impulse response (FIR) or infinite impulse response (IIR) design. Optical receivers exploiting polarization diversity should also compensate for the random fluctuations of the polarization state caused by the stochastic change of fiber birefringence. PMD compensation is widely performed using constant modulus algorithm (CMA) [15] or multi modulus algorithm (MMA) [16]. Attributed to these well-developed linear equalization techniques, fiber nonlinearity still remains a bottleneck for increasing the capacity and transmission reach [24].

step size is sufficiently small when solving the inverse nonlinear Schrodinger equation (NLSE). Several strategies have been proposed in the literature to reduce the number of required DBP steps and enhancing the performance either by including temporal correlations [27, 28] or optimizing the parameters [29]. Nevertheless, the performance of DBP significantly deteriorates in the presence of stochastic impairments like laser phase noise and NLPN [30–32]. Moreover, its

Applications of Kalman Filters for Coherent Optical Communication Systems

http://dx.doi.org/10.5772/intechopen.71617

207

Apart from fiber linear and nonlinear impairment compensation, digital carrier synchronization has also become an essential component of the coherent receivers, for synchronizing the phase and frequency offsets between the transmitter laser and the local oscillator (LO), eliminating the necessity of a phase locked loop. Several carrier phase estimation (CPE) techniques have been proposed for suppressing the laser phase noise [35–44]. CPE being a low complex technique is also under wide investigation to compensate the nonlinear phase shift owing to Kerr effect, besides laser phase noise [30, 44–47]. Investigations were also carried on the combined performance of DBP and CPE, in order to reduce the number of DBP steps per span and there by its complexity [31, 48]. It was reported in [30, 45], that the considered CPE methods outperform the DBP technique implemented using asymmetric split step Fourier method (SSFM) with one step per span and without any parameter optimization. However, the accumulated ASE noise and the NLPN at high signal powers pose a challenging constraint on the conventional CPE limiting its nonlinear mitigation capability. Moreover, a phase unwrapping function [35] is typically required by CPE, which might increase the probability of cycle slips [49] and error propagation. Furthermore, the commonly employed CPE techniques have low tolerance towards the frequency offset (FO) between the transmitter laser and

In the recent era, Kalman filtering has gained huge attention in the field of optical communication systems, owing to its potential capability to mitigate several optical transmission impairments simultaneously. The Kalman filter was first developed by R. E. Kalman in 1960. In [51], he presented a new approach to the linear filtering and prediction problems by introducing the state space notation, where the random processes/signals to be estimated are represented as the output of a linear dynamic system perturbed by uncorrelated noise. This approach facilitates recursive computation of the optimal solution and highly reduces the computational effort as compared to the conventional Wiener filter besides eliminating the memory problems. The so called Kalman filter computes the optimal solution recursively in the minimum mean square error (MMSE) sense. While the applicability of the Wiener filter is limited to stationary processes, the Kalman filters can be also applied to the non-stationary processes. An added advantage of the Kalman filter is its extended applicability also to the nonlinear systems through an approximate linearization and the so called filter is known as extended Kalman filter (EKF). This has attracted the Kalman filters for numerous real-time applications in the fields of navigation, radar, mobile communications, speech signal processing and so forth. Currently, EKF is under wide investigation in coherent optical communication systems for tracking and mitigating linear and nonlinear phase noise, amplitude noise, phase and frequency offsets as well as polarization de-multiplexing [12, 13, 31, 52–66].

applicability is limited to single channel systems [33, 34].

the LO. Therefore, a separate FO estimation module is necessary [50].

Although, multiple information bits being encoded in a single symbol significantly increase the spectral efficiency, the signal becomes more sensitive to the amplified spontaneous emission (ASE) noise that is added in the optical amplifiers along the transmission link. Therefore, a reliable transmission over long distance demands the signal to be launched into the optical fiber at higher power, to ensure a sufficiently high optical signal to noise ratio (OSNR) at the receiver. However, the maximum transmittable launch power per fiber span is constrained by the Kerr nonlinear effects, including self-phase modulation (SPM) and cross phase modulation (XPM) in case of WDM systems, which results in signal degradation [25]. This degrading impact of Kerr nonlinearity is much more severe in multi-channel systems with increasing number of channels [26]. Moreover, the nonlinear phase noise (NLPN) resulting from the signal and ASE noise interactions at high launch powers, deteriorates the signal quality further. On the other hand, signal transmission at low launch powers is limited by the ASE noise. Therefore, mitigation of fiber nonlinearity is vital to enhance the capacity ensuring reliable transmission. Consequently, several nonlinear mitigation techniques have been proposed in the recent era, of which, digital backward propagation (DBP) [10], maximum likelihood sequence estimation (MLSE) based nonlinearity mitigation [17], spectral inversion [18–20], phase conjugated twin waves [21, 22] and perturbation based approaches [23] gained considerable attention. However, the real time implementation of these algorithms is extremely challenging owing to either the high required computational effort or the higher bandwidth consumption. Although computationally complex, DBP has drawn significant attention owing to its capability of mitigating linear and nonlinear impairments simultaneously, provided the channel characteristics are known and the step size is sufficiently small when solving the inverse nonlinear Schrodinger equation (NLSE). Several strategies have been proposed in the literature to reduce the number of required DBP steps and enhancing the performance either by including temporal correlations [27, 28] or optimizing the parameters [29]. Nevertheless, the performance of DBP significantly deteriorates in the presence of stochastic impairments like laser phase noise and NLPN [30–32]. Moreover, its applicability is limited to single channel systems [33, 34].

wavelength division multiplexing (WDM) systems, coherent detection, digital signal processing (DSP) techniques and forward error correction (FEC) schemes ensuring reliable transmission. The advent of coherent detection along with subsequent DSP made it possible to deploy spectrally efficient higher order modulation formats and multiplexing techniques [6, 7]. Moreover, it has also made feasible to digitally equalize the optical fiber transmission impairments [8], which are the main hurdle to increase the bandwidth-distance product. The transmission capacity can be increased several times by employing complex modulation formats like m-ary quadrature amplitude modulation (with m = 4, 16, 64 and so on), and multiplexing techniques like polarization division multiplexing (PDM) and WDM. However, they are more vulnerable to the optical transmission impairments as well as to the carrier phase and frequency offset (FO). Hence, effective DSP algorithms for combatting with the channel impairments were under active research over the past decade [8–23]. Consequently, coherent optical receivers are well developed and employ digital filters that allow for effective equalization of fiber linear impairments like chromatic dispersion (CD) and polarization mode dispersion (PMD) in the electric domain [9]. Typically, CD can be compensated by either frequency or time domain filters using finite impulse response (FIR) or infinite impulse response (IIR) design. Optical receivers exploiting polarization diversity should also compensate for the random fluctuations of the polarization state caused by the stochastic change of fiber birefringence. PMD compensation is widely performed using constant modulus algorithm (CMA) [15] or multi modulus algorithm (MMA) [16]. Attributed to these well-developed linear equalization techniques, fiber nonlinearity still remains a bottleneck

Although, multiple information bits being encoded in a single symbol significantly increase the spectral efficiency, the signal becomes more sensitive to the amplified spontaneous emission (ASE) noise that is added in the optical amplifiers along the transmission link. Therefore, a reliable transmission over long distance demands the signal to be launched into the optical fiber at higher power, to ensure a sufficiently high optical signal to noise ratio (OSNR) at the receiver. However, the maximum transmittable launch power per fiber span is constrained by the Kerr nonlinear effects, including self-phase modulation (SPM) and cross phase modulation (XPM) in case of WDM systems, which results in signal degradation [25]. This degrading impact of Kerr nonlinearity is much more severe in multi-channel systems with increasing number of channels [26]. Moreover, the nonlinear phase noise (NLPN) resulting from the signal and ASE noise interactions at high launch powers, deteriorates the signal quality further. On the other hand, signal transmission at low launch powers is limited by the ASE noise. Therefore, mitigation of fiber nonlinearity is vital to enhance the capacity ensuring reliable transmission. Consequently, several nonlinear mitigation techniques have been proposed in the recent era, of which, digital backward propagation (DBP) [10], maximum likelihood sequence estimation (MLSE) based nonlinearity mitigation [17], spectral inversion [18–20], phase conjugated twin waves [21, 22] and perturbation based approaches [23] gained considerable attention. However, the real time implementation of these algorithms is extremely challenging owing to either the high required computational effort or the higher bandwidth consumption. Although computationally complex, DBP has drawn significant attention owing to its capability of mitigating linear and nonlinear impairments simultaneously, provided the channel characteristics are known and the

for increasing the capacity and transmission reach [24].

206 Kalman Filters - Theory for Advanced Applications

Apart from fiber linear and nonlinear impairment compensation, digital carrier synchronization has also become an essential component of the coherent receivers, for synchronizing the phase and frequency offsets between the transmitter laser and the local oscillator (LO), eliminating the necessity of a phase locked loop. Several carrier phase estimation (CPE) techniques have been proposed for suppressing the laser phase noise [35–44]. CPE being a low complex technique is also under wide investigation to compensate the nonlinear phase shift owing to Kerr effect, besides laser phase noise [30, 44–47]. Investigations were also carried on the combined performance of DBP and CPE, in order to reduce the number of DBP steps per span and there by its complexity [31, 48]. It was reported in [30, 45], that the considered CPE methods outperform the DBP technique implemented using asymmetric split step Fourier method (SSFM) with one step per span and without any parameter optimization. However, the accumulated ASE noise and the NLPN at high signal powers pose a challenging constraint on the conventional CPE limiting its nonlinear mitigation capability. Moreover, a phase unwrapping function [35] is typically required by CPE, which might increase the probability of cycle slips [49] and error propagation. Furthermore, the commonly employed CPE techniques have low tolerance towards the frequency offset (FO) between the transmitter laser and the LO. Therefore, a separate FO estimation module is necessary [50].

In the recent era, Kalman filtering has gained huge attention in the field of optical communication systems, owing to its potential capability to mitigate several optical transmission impairments simultaneously. The Kalman filter was first developed by R. E. Kalman in 1960. In [51], he presented a new approach to the linear filtering and prediction problems by introducing the state space notation, where the random processes/signals to be estimated are represented as the output of a linear dynamic system perturbed by uncorrelated noise. This approach facilitates recursive computation of the optimal solution and highly reduces the computational effort as compared to the conventional Wiener filter besides eliminating the memory problems. The so called Kalman filter computes the optimal solution recursively in the minimum mean square error (MMSE) sense. While the applicability of the Wiener filter is limited to stationary processes, the Kalman filters can be also applied to the non-stationary processes. An added advantage of the Kalman filter is its extended applicability also to the nonlinear systems through an approximate linearization and the so called filter is known as extended Kalman filter (EKF). This has attracted the Kalman filters for numerous real-time applications in the fields of navigation, radar, mobile communications, speech signal processing and so forth. Currently, EKF is under wide investigation in coherent optical communication systems for tracking and mitigating linear and nonlinear phase noise, amplitude noise, phase and frequency offsets as well as polarization de-multiplexing [12, 13, 31, 52–66]. Moreover, enhancing the performance obtained from EKF by incorporating with other existing techniques like DBP has also been studied in [31, 58, 59]. EKF requires only a few complex multiplications to recover one data symbol. Besides the advantage from low complexity, it further offers other benefits including faster convergence, joint tracking and compensation of fiber impairments. Therefore, it is worth discussing and reviewing the applications of Kalman filters for coherent optical communications in a nutshell.

corresponds to the conditional mean [67] as given in Eq. (3). Here, E[�] denotes the expectation

The Kalman filter computes the optimal state recursively, following a predictor-corrector structure, where a prediction is computed prior to the availability of the observation at current time instant k and updates the prediction when the observation at time instant k is available. Throughout this Chapter, we follow the typical notation convention for the Kalman filter equations: any variable with subscript k|k – 1 denotes prediction or apriori estimate, and any variable with subscript k|k or simply, denotes the updated or aposteriori estimate. During the prediction step, the Kalman filter makes the best guess about the system's state based on its dynamics, prior to the availability of the current observation. The state prediction denoted by x^<sup>k</sup>jk�<sup>1</sup> is given in Eq. (4). The uncertainty associated with the prediction is given by the apriori error covariance matrix Pk|k–1, as in Eq. (5). Under the given assumptions and initial conditions, the conditional probability density function (pdf) p(xk|y1,y2,…,yk–1) is also Gaussian, where the apriori state estimate x^<sup>k</sup>jk�<sup>1</sup> and the apriori error covariance Pk|k–1, reflects the mean and variance of the

distribution, as given in Eq. (6). Here, N denotes normal or Gaussian distribution.

P<sup>k</sup>jk�<sup>1</sup> ¼ E½ðx<sup>k</sup> � x^<sup>k</sup>jk�<sup>1</sup>Þðx<sup>k</sup> � x^<sup>k</sup>jk�<sup>1</sup>Þ

computing the optimal state estimate.

During the update step, when the new observation at time k is available, the optimal estimate is computed as a linear combination of the prediction and the new information available from the current measurement weighted by an optimal weighting matrix known as Kalman gain. The update equations can be summarized in Eqs. (7)–(10). The innovation denoted by vk, can be interpreted as the new information that is available in the observation y<sup>k</sup> relative to all the past observations up to time instant k – 1. It is computed as the difference between the actual and the predicted observation <sup>y</sup>^<sup>k</sup>jk�<sup>1</sup>, and is given in Eq. (7). The Kalman gain, denoted by <sup>K</sup>k, determines the extent up to which the innovation should be taken into account in updating the apriori state estimate and is computed according to Eq. (8). Here, H denotes the Hermitian operator. The updated or aposteriori state estimate x^<sup>k</sup>j<sup>k</sup>, and the aposteriori error covariance Pk|k, are computed as given in Eqs. (9) and (10), respectively. The aposteriori pdf p(xk|y1,y2,…,yk) is also Gaussian distributed with mean and variance given by the aposteriori state estimate x^<sup>k</sup>j<sup>k</sup> and the aposteriori error covariance Pk|k, respectively, as given in Eq. (11). Thus, the Kalman filter propagates the first and second order moments of the state distribution recursively for

x<sup>k</sup> ¼ Fkx<sup>k</sup>�<sup>1</sup> þ w<sup>k</sup> (1)

Applications of Kalman Filters for Coherent Optical Communication Systems

http://dx.doi.org/10.5772/intechopen.71617

209

y<sup>k</sup> ¼ Hkx<sup>k</sup> þ n<sup>k</sup> (2)

x^<sup>k</sup> ¼ E½xkjy1, y2, …, yk� (3)

<sup>x</sup>^<sup>k</sup>jk�<sup>1</sup> <sup>¼</sup> <sup>E</sup>½xkjy1, <sup>y</sup>2, …, <sup>y</sup><sup>k</sup>�<sup>1</sup>� ¼ <sup>F</sup>kx^<sup>k</sup>�1jk�<sup>1</sup> (4)

<sup>H</sup>� ¼ <sup>F</sup>kPkFH

<sup>p</sup>ðxkjy1, <sup>y</sup>2, …, <sup>y</sup><sup>k</sup>�<sup>1</sup>Þ � <sup>N</sup>ðx^<sup>k</sup>jk�<sup>1</sup>, <sup>P</sup><sup>k</sup>jk�<sup>1</sup> <sup>Þ</sup> (6)

<sup>k</sup> þ Q<sup>k</sup> (5)

operator.

This chapter is organized as follows: In Section 2, we discuss the principles of Kalman filter by describing the state space notation and the recursive equations. We further present some variations of the Kalman filter, namely, EKF and adaptive Kalman filter (AKF). Section 3 details the applications of EKF for coherent optical communications. We illustrate how to employ Kalman filtering for the joint tracking of several optical transmission impairments by formulating the state space model (SSM) and detailing the working principles. We also describe our numerical model and present the results to justify the theoretical findings. Finally, the chapter is concluded with a note on the key points, in Section 4.

### 2. The Kalman filter

A Kalman filter is an optimal recursive linear MMSE estimator that estimates the state of a linear dynamic perturbed by noise. Since the true state of the system is not observable, instead we obtain the measurements or observations that are corrupted by noise. Now, the goal of the Kalman filter is to obtain an optimal estimate of the unknown state from the noisy observations recursively. The stochastic process under estimation is modeled by a state space model (SSM) which facilitates the recursive nature of the Kalman filter. In the following, we present the general framework of the Kalman filtering and also discuss briefly the principles of the EKF and AKF.

#### 2.1. Principles of Kalman filter

Consider a discrete-time, linear, time varying system in the state space notation, given by Eqs. (1) and (2). Eq. (1) describes how the true state of the system evolves over time and is known as the state or the process equation. Eq. (2) describes how the measurements are related to the states and is known as measurement or observation equation. Here, k denotes the time instant, x<sup>k</sup> and y<sup>k</sup> denote the state vector and the measurement vector, respectively. F<sup>k</sup> denotes the state transition matrix that relates the states at the time instances k and k – 1, in the absence of process noise wk. H<sup>k</sup> denotes the measurement matrix that relates the states to the measurements in the absence of measurement noise nk. The process and measurement noise vectors w<sup>k</sup> and n<sup>k</sup> are assumed to be zero mean white Gaussian noise processes with co-variance matrices Q<sup>k</sup> and Rk, respectively. It is also assumed that the initial state x<sup>0</sup> at time instant 0, is a Gaussian random vector. Given the SSM and these assumptions, the objective of the Kalman filter is to obtain a linear MMSE estimate of x<sup>k</sup> based on the observations {y1, y2,…yk}. The solution corresponds to the conditional mean [67] as given in Eq. (3). Here, E[�] denotes the expectation operator.

Moreover, enhancing the performance obtained from EKF by incorporating with other existing techniques like DBP has also been studied in [31, 58, 59]. EKF requires only a few complex multiplications to recover one data symbol. Besides the advantage from low complexity, it further offers other benefits including faster convergence, joint tracking and compensation of fiber impairments. Therefore, it is worth discussing and reviewing the applications of Kalman

This chapter is organized as follows: In Section 2, we discuss the principles of Kalman filter by describing the state space notation and the recursive equations. We further present some variations of the Kalman filter, namely, EKF and adaptive Kalman filter (AKF). Section 3 details the applications of EKF for coherent optical communications. We illustrate how to employ Kalman filtering for the joint tracking of several optical transmission impairments by formulating the state space model (SSM) and detailing the working principles. We also describe our numerical model and present the results to justify the theoretical findings. Finally,

A Kalman filter is an optimal recursive linear MMSE estimator that estimates the state of a linear dynamic perturbed by noise. Since the true state of the system is not observable, instead we obtain the measurements or observations that are corrupted by noise. Now, the goal of the Kalman filter is to obtain an optimal estimate of the unknown state from the noisy observations recursively. The stochastic process under estimation is modeled by a state space model (SSM) which facilitates the recursive nature of the Kalman filter. In the following, we present the general framework of the Kalman filtering and also discuss briefly the principles of the

Consider a discrete-time, linear, time varying system in the state space notation, given by Eqs. (1) and (2). Eq. (1) describes how the true state of the system evolves over time and is known as the state or the process equation. Eq. (2) describes how the measurements are related to the states and is known as measurement or observation equation. Here, k denotes the time instant, x<sup>k</sup> and y<sup>k</sup> denote the state vector and the measurement vector, respectively. F<sup>k</sup> denotes the state transition matrix that relates the states at the time instances k and k – 1, in the absence of process noise wk. H<sup>k</sup> denotes the measurement matrix that relates the states to the measurements in the absence of measurement noise nk. The process and measurement noise vectors w<sup>k</sup> and n<sup>k</sup> are assumed to be zero mean white Gaussian noise processes with co-variance matrices Q<sup>k</sup> and Rk, respectively. It is also assumed that the initial state x<sup>0</sup> at time instant 0, is a Gaussian random vector. Given the SSM and these assumptions, the objective of the Kalman filter is to obtain a linear MMSE estimate of x<sup>k</sup> based on the observations {y1, y2,…yk}. The solution

filters for coherent optical communications in a nutshell.

208 Kalman Filters - Theory for Advanced Applications

2. The Kalman filter

EKF and AKF.

2.1. Principles of Kalman filter

the chapter is concluded with a note on the key points, in Section 4.

$$\mathbf{x}\_k = \mathbf{F}\_k \mathbf{x}\_{k-1} + \mathbf{w}\_k \tag{1}$$

$$\mathbf{y}\_k = \mathbf{H}\_k \mathbf{x}\_k + \mathbf{n}\_k \tag{2}$$

$$
\hat{\mathbf{x}}\_k = \mathbb{E}[\mathbf{x}\_k | \mathbf{y}\_{1'}, \mathbf{y}\_{2'}, \dots, \mathbf{y}\_k] \tag{3}
$$

The Kalman filter computes the optimal state recursively, following a predictor-corrector structure, where a prediction is computed prior to the availability of the observation at current time instant k and updates the prediction when the observation at time instant k is available. Throughout this Chapter, we follow the typical notation convention for the Kalman filter equations: any variable with subscript k|k – 1 denotes prediction or apriori estimate, and any variable with subscript k|k or simply, denotes the updated or aposteriori estimate. During the prediction step, the Kalman filter makes the best guess about the system's state based on its dynamics, prior to the availability of the current observation. The state prediction denoted by x^<sup>k</sup>jk�<sup>1</sup> is given in Eq. (4). The uncertainty associated with the prediction is given by the apriori error covariance matrix Pk|k–1, as in Eq. (5). Under the given assumptions and initial conditions, the conditional probability density function (pdf) p(xk|y1,y2,…,yk–1) is also Gaussian, where the apriori state estimate x^<sup>k</sup>jk�<sup>1</sup> and the apriori error covariance Pk|k–1, reflects the mean and variance of the distribution, as given in Eq. (6). Here, N denotes normal or Gaussian distribution.

$$
\hat{\mathbf{x}}\_{k|k-1} = \mathbb{E}[\mathbf{x}\_k | \mathbf{y}\_{1'} \mathbf{y}\_{2'} \dots \mathbf{y}\_{k-1}] = \mathbf{F}\_k \hat{\mathbf{x}}\_{k-1|k-1} \tag{4}
$$

$$\mathbf{P}\_{k|k-1} = \mathbf{E}[(\mathbf{x}\_k - \hat{\mathbf{x}}\_{k|k-1})(\mathbf{x}\_k - \hat{\mathbf{x}}\_{k|k-1})^\mathbf{H}] = \mathbf{F}\_k \mathbf{P}\_k \mathbf{F}\_k^\mathbf{H} + \mathbf{Q}\_k \tag{5}$$

$$p(\mathbf{x}\_k | \mathbf{y}\_{1'}, \mathbf{y}\_{2'}, \dots, \mathbf{y}\_{k-1}) \sim \mathcal{N}(\hat{\mathbf{x}}\_{k|k-1}, \mathbf{P}\_{k|k-1}) \tag{6}$$

During the update step, when the new observation at time k is available, the optimal estimate is computed as a linear combination of the prediction and the new information available from the current measurement weighted by an optimal weighting matrix known as Kalman gain. The update equations can be summarized in Eqs. (7)–(10). The innovation denoted by vk, can be interpreted as the new information that is available in the observation y<sup>k</sup> relative to all the past observations up to time instant k – 1. It is computed as the difference between the actual and the predicted observation <sup>y</sup>^<sup>k</sup>jk�<sup>1</sup>, and is given in Eq. (7). The Kalman gain, denoted by <sup>K</sup>k, determines the extent up to which the innovation should be taken into account in updating the apriori state estimate and is computed according to Eq. (8). Here, H denotes the Hermitian operator. The updated or aposteriori state estimate x^<sup>k</sup>j<sup>k</sup>, and the aposteriori error covariance Pk|k, are computed as given in Eqs. (9) and (10), respectively. The aposteriori pdf p(xk|y1,y2,…,yk) is also Gaussian distributed with mean and variance given by the aposteriori state estimate x^<sup>k</sup>j<sup>k</sup> and the aposteriori error covariance Pk|k, respectively, as given in Eq. (11). Thus, the Kalman filter propagates the first and second order moments of the state distribution recursively for computing the optimal state estimate.

$$\mathbf{v}\_k = \mathbf{y}\_k - \hat{\mathbf{y}}\_{k|k-1} \tag{7}$$

2.3. Adaptive Kalman filtering

statistical smoothing.

FkPkF<sup>H</sup>

The Kalman filter computes the optimal solution, provided the process noise and measurement noise covariances, Q<sup>k</sup> and Rk, respectively, are known apriori. However, in practice, a precise knowledge about the noise statistics might not be available. The Kalman gain K<sup>k</sup> takes into account the noise covariances, Q<sup>k</sup> and Rk, to determine the extent of reliability between the predicted state x^<sup>k</sup>jk�<sup>1</sup> and the innovation vk. Therefore, a poor knowledge of the noise statistics might significantly degrade the filter performance and even leads to divergence. To overcome these difficulties, an adaptive approach can be followed to adaptively estimate the noise covariances from the noise samples, (for example, the innovation sequence) that are generated during the Kalman recursions at each time instant. This leads to the adaptive Kalman filtering. The different approaches for adaptive filtering are classified into four types: Bayesian, maximum likelihood, correlation and covariance matching methods [68]. Here, we discuss the approach based on covariance matching [68, 69] for adaptive estimation of noise statistics. The basic idea behind this approach lies on the fact that for an optimal filter, the theoretical covariance of the innovation vk, denoted by Sk, given in Eq. (22) should be consistent with the empirically estimated covariance given in Eq. (23). Here, m denotes the window size to provide

Applications of Kalman Filters for Coherent Optical Communication Systems

http://dx.doi.org/10.5772/intechopen.71617

211

<sup>S</sup><sup>k</sup> <sup>¼</sup> <sup>H</sup>kP<sup>k</sup>jk�<sup>1</sup>HH

<sup>k</sup> Þ ¼ <sup>1</sup> m Xm�1 i¼0

<sup>k</sup> <sup>þ</sup> <sup>R</sup><sup>k</sup> <sup>¼</sup> <sup>1</sup>

Since the Kalman gain K<sup>k</sup> depends on the ratio of the process and measurement noise covariances Qk/Rk, rather than on their individual values, if either of Q<sup>k</sup> or Rk, is known, the other can be adaptively estimated by satisfying the condition for covariance matching, given in Eq. (24). When Q<sup>k</sup> is known, R<sup>k</sup> can be directly estimated from Eq. (24). Alternatively, when R<sup>k</sup> is given, Q<sup>k</sup> can be estimated by a scaling procedure to improve the robustness of the filter. The basic idea behind this scaling method is that if the estimated covariance of vk, on the right hand side of Eq. (24), is much larger than the theoretical covariance, then Q<sup>k</sup> (please note that P<sup>k</sup>jk�<sup>1</sup> ¼

<sup>k</sup> þ Qk) should be increased to bring the theoretical covariance closer to the estimated one and vice-versa. Therefore, Q<sup>k</sup> can be adaptively updated in order to balance any deviations between the theoretical and estimated innovation covariance by considering a scaling factor αk,

<sup>i</sup>¼<sup>0</sup> <sup>v</sup><sup>k</sup>�<sup>i</sup>v<sup>H</sup>

traceðHkP<sup>k</sup>jk�<sup>1</sup>H<sup>H</sup>

<sup>k</sup>�<sup>i</sup> � <sup>R</sup>k<sup>Þ</sup>

<sup>v</sup><sup>k</sup>�<sup>i</sup>v<sup>H</sup>

<sup>v</sup><sup>k</sup>�<sup>i</sup>v<sup>H</sup>

m Xm�1 i¼0

<sup>E</sup>ðvkv<sup>H</sup>

<sup>H</sup>kP<sup>k</sup>jk�<sup>1</sup>HH

as given in Eq. (25). The estimate of Qk, denoted by Q^ <sup>k</sup>, is given in Eq. (26).

<sup>α</sup><sup>k</sup> <sup>¼</sup> trace<sup>ð</sup> <sup>1</sup>

m X<sup>m</sup>�<sup>1</sup> <sup>k</sup> þ R<sup>k</sup> (22)

<sup>k</sup>�<sup>i</sup> (23)

<sup>k</sup> <sup>Þ</sup> (25)

<sup>Q</sup>^ <sup>k</sup> <sup>¼</sup> <sup>α</sup>kQ^ <sup>k</sup>�<sup>1</sup> (26)

<sup>k</sup>�<sup>i</sup> (24)

$$\mathbf{K}\_{k} = \,^{\mathbf{P}}\_{k|k-1} \, \mathbf{H}\_{k}^{\text{H}} (\mathbf{H}\_{k} \mathbf{P}\_{k|k-1} \mathbf{H}\_{k}^{\text{H}} + \,^{\mathbf{R}} \mathbf{k})^{-1} \tag{8}$$

$$
\hat{\mathbf{x}}\_{k|k} = \hat{\mathbf{x}}\_{k|k-1} + \mathbf{K}\_k \mathbf{v}\_k \tag{9}
$$

$$\mathbf{P}\_{k|k} = \mathbf{P}\_{k|k-1} - \mathbf{P}\_{k|k-1} \mathbf{K}\_k \mathbf{H}\_k \tag{10}$$

$$p(\mathbf{x}\_k | \mathbf{y}\_{1'}, \mathbf{y}\_{2'}, \dots, \mathbf{y}\_k) \sim N(\hat{\mathbf{x}}\_{k|k}, \mathbf{P}\_{k|k}) \tag{11}$$

#### 2.2. Extended Kalman filtering

In Section 2.1, we addressed the problem of estimating the unknown state of a linear dynamic system from noisy observations. Now, we consider the filtering problem for nonlinear system dynamics (either the process or observation model or both being nonlinear). The Kalman filter solution can be adopted for the nonlinear dynamic systems through an approximate linearization procedure and the resulting filter is known as EKF. Consider a nonlinear dynamic system described by the SSM given in Eqs. (12) and (13). Here, fk(�) and hk(�) denote the nonlinear state transition function and the measurement function, respectively.

$$\mathbf{x}\_k = \mathbf{f}\_k(\mathbf{x}\_{k-1}) + \mathbf{w}\_k \tag{12}$$

$$\mathbf{y}\_k = \mathbf{h}\_k(\mathbf{x}\_k) + \mathbf{n}\_k \tag{13}$$

The nonlinear system dynamics can be linearized through a first order Taylor approximation at each time instant, around the most recent state estimate. This forms the basic idea of EKF. Let, A<sup>k</sup> and B<sup>k</sup> be the Jacobian matrices of fk(�) and hk(�), respectively, and are computed according to Eqs. (14) and (15). Under the given assumptions and, the initial conditions as discussed in the earlier section, the EKF recursive equations can be summarized in Eqs. (16)– (20).

$$\mathbf{A}\_{k} = \frac{\partial \mathbf{f}\_{k}(\mathbf{x})}{\partial \mathbf{x}} \text{ at } \mathbf{x} = \,\, \hat{\mathbf{x}}\_{k-1|k-1} \tag{14}$$

$$\mathbf{B}\_{k} = \frac{\partial \mathbf{h}\_{k}(\mathbf{x})}{\partial \mathbf{x}} \text{ at } \mathbf{x} = \,\,\hat{\mathbf{x}}\_{k|k-1} \tag{15}$$

$$
\hat{\mathbf{x}}\_{k|k-1} = \mathbf{f}\_k(\hat{\mathbf{x}}\_{k-1|k-1}) \tag{16}
$$

$$\mathbf{P}\_{k|k-1} = \mathbf{A}\_k \mathbf{P}\_k \mathbf{A}\_k^H + \mathbf{Q}\_k \tag{17}$$

$$\mathbf{v}\_{k} = \mathbf{y}\_{k} - \hat{\mathbf{y}}\_{k|k-1} = \mathbf{y}\_{k} - \mathbf{h}\_{k}(\hat{\mathbf{x}}\_{k|k-1}) \tag{18}$$

$$\mathbf{K}\_{k} = \,^{\mathbf{P}}\_{k|k-1} \,^{\mathbf{B}}\_{k}^{\mathbf{H}} (\mathbf{B}\_{k} \mathbf{P}\_{k|k-1} \mathbf{B}\_{k}^{\mathbf{H}} + \,^{\mathbf{R}}\_{k})^{-1} \tag{19}$$

$$
\hat{\mathbf{x}}\_{k|k} = \hat{\mathbf{x}}\_{k|k-1} + \mathbf{K}\_k \mathbf{v}\_k \tag{20}
$$

$$\mathbf{P}\_{k|k} = \mathbf{P}\_{k|k-1} - \mathbf{P}\_{k|k-1} \mathbf{K}\_k \mathbf{B}\_k \tag{21}$$

#### 2.3. Adaptive Kalman filtering

<sup>v</sup><sup>k</sup> <sup>¼</sup> <sup>y</sup><sup>k</sup> � <sup>y</sup>^<sup>k</sup>jk�<sup>1</sup> (7)

x^<sup>k</sup>j<sup>k</sup> ¼ x^<sup>k</sup>jk�<sup>1</sup> þ Kkv<sup>k</sup> (9)

x<sup>k</sup> ¼ fkðx<sup>k</sup>�<sup>1</sup>Þ þ w<sup>k</sup> (12)

y<sup>k</sup> ¼ hkðxkÞ þ n<sup>k</sup> (13)

at x ¼ x^<sup>k</sup>�1jk�<sup>1</sup> (14)

at x ¼ x^<sup>k</sup>jk�<sup>1</sup> (15)

<sup>k</sup> þ Q<sup>k</sup> (17)

�<sup>1</sup> (19)

x^<sup>k</sup>jk�<sup>1</sup> ¼ fkðx^<sup>k</sup>�1jk�<sup>1</sup> Þ (16)

x^<sup>k</sup>j<sup>k</sup> ¼ x^<sup>k</sup>jk�<sup>1</sup> þ Kkv<sup>k</sup> (20)

P<sup>k</sup>j<sup>k</sup> ¼ P<sup>k</sup>jk�<sup>1</sup> � P<sup>k</sup>jk�<sup>1</sup>KkB<sup>k</sup> (21)

<sup>v</sup><sup>k</sup> <sup>¼</sup> <sup>y</sup><sup>k</sup> � <sup>y</sup>^<sup>k</sup>jk�<sup>1</sup> <sup>¼</sup> <sup>y</sup><sup>k</sup> � <sup>h</sup>kðx^<sup>k</sup>jk�<sup>1</sup> <sup>Þ</sup> (18)

<sup>k</sup> þ RkÞ

P<sup>k</sup>j<sup>k</sup> ¼ P<sup>k</sup>jk�<sup>1</sup> � P<sup>k</sup>jk�<sup>1</sup>KkH<sup>k</sup> (10)

pðxkjy1, y2, …, ykÞ � Nðx^<sup>k</sup>jk, P<sup>k</sup>j<sup>k</sup> Þ (11)

�<sup>1</sup> (8)

<sup>k</sup> þ RkÞ

<sup>K</sup><sup>k</sup> <sup>¼</sup> <sup>P</sup><sup>k</sup>jk�<sup>1</sup> <sup>H</sup><sup>H</sup>

transition function and the measurement function, respectively.

2.2. Extended Kalman filtering

210 Kalman Filters - Theory for Advanced Applications

(20).

<sup>k</sup> <sup>ð</sup>HkP<sup>k</sup>jk�<sup>1</sup>H<sup>H</sup>

In Section 2.1, we addressed the problem of estimating the unknown state of a linear dynamic system from noisy observations. Now, we consider the filtering problem for nonlinear system dynamics (either the process or observation model or both being nonlinear). The Kalman filter solution can be adopted for the nonlinear dynamic systems through an approximate linearization procedure and the resulting filter is known as EKF. Consider a nonlinear dynamic system described by the SSM given in Eqs. (12) and (13). Here, fk(�) and hk(�) denote the nonlinear state

The nonlinear system dynamics can be linearized through a first order Taylor approximation at each time instant, around the most recent state estimate. This forms the basic idea of EKF. Let, A<sup>k</sup> and B<sup>k</sup> be the Jacobian matrices of fk(�) and hk(�), respectively, and are computed according to Eqs. (14) and (15). Under the given assumptions and, the initial conditions as discussed in the earlier section, the EKF recursive equations can be summarized in Eqs. (16)–

> <sup>A</sup><sup>k</sup> <sup>¼</sup> <sup>∂</sup>fkðx<sup>Þ</sup> ∂x

> > <sup>B</sup><sup>k</sup> <sup>¼</sup> <sup>∂</sup>hkðx<sup>Þ</sup> ∂x

<sup>K</sup><sup>k</sup> <sup>¼</sup> <sup>P</sup><sup>k</sup>jk�<sup>1</sup> <sup>B</sup><sup>H</sup>

<sup>P</sup><sup>k</sup>jk�<sup>1</sup> <sup>¼</sup> <sup>A</sup>kPkA<sup>H</sup>

<sup>k</sup> <sup>ð</sup>BkP<sup>k</sup>jk�<sup>1</sup>B<sup>H</sup>

The Kalman filter computes the optimal solution, provided the process noise and measurement noise covariances, Q<sup>k</sup> and Rk, respectively, are known apriori. However, in practice, a precise knowledge about the noise statistics might not be available. The Kalman gain K<sup>k</sup> takes into account the noise covariances, Q<sup>k</sup> and Rk, to determine the extent of reliability between the predicted state x^<sup>k</sup>jk�<sup>1</sup> and the innovation vk. Therefore, a poor knowledge of the noise statistics might significantly degrade the filter performance and even leads to divergence. To overcome these difficulties, an adaptive approach can be followed to adaptively estimate the noise covariances from the noise samples, (for example, the innovation sequence) that are generated during the Kalman recursions at each time instant. This leads to the adaptive Kalman filtering. The different approaches for adaptive filtering are classified into four types: Bayesian, maximum likelihood, correlation and covariance matching methods [68]. Here, we discuss the approach based on covariance matching [68, 69] for adaptive estimation of noise statistics. The basic idea behind this approach lies on the fact that for an optimal filter, the theoretical covariance of the innovation vk, denoted by Sk, given in Eq. (22) should be consistent with the empirically estimated covariance given in Eq. (23). Here, m denotes the window size to provide statistical smoothing.

$$\mathbf{S}\_{k} = \mathbf{H}\_{k}\mathbf{P}\_{k|k-1}\mathbf{H}\_{k}^{\mathbf{H}} + \mathbf{R}\_{k} \tag{22}$$

$$\mathbb{E}(\mathbf{v}\boldsymbol{\omega}\mathbf{v}\_k^H) = \frac{1}{m} \sum\_{i=0}^{m-1} \mathbf{v}\_{k-i} \mathbf{v}\_{k-i}^H \tag{23}$$

$$\mathbf{H}\_k \mathbf{P}\_{k|k-1} \mathbf{H}\_{\mathbf{k}}^{\mathbf{H}} + \mathbf{R}\_k = \frac{1}{m} \sum\_{i=0}^{m-1} \mathbf{v}\_{k-i} \mathbf{v}\_{k-i}^{\mathbf{H}} \tag{24}$$

Since the Kalman gain K<sup>k</sup> depends on the ratio of the process and measurement noise covariances Qk/Rk, rather than on their individual values, if either of Q<sup>k</sup> or Rk, is known, the other can be adaptively estimated by satisfying the condition for covariance matching, given in Eq. (24). When Q<sup>k</sup> is known, R<sup>k</sup> can be directly estimated from Eq. (24). Alternatively, when R<sup>k</sup> is given, Q<sup>k</sup> can be estimated by a scaling procedure to improve the robustness of the filter. The basic idea behind this scaling method is that if the estimated covariance of vk, on the right hand side of Eq. (24), is much larger than the theoretical covariance, then Q<sup>k</sup> (please note that P<sup>k</sup>jk�<sup>1</sup> ¼ FkPkF<sup>H</sup> <sup>k</sup> þ Qk) should be increased to bring the theoretical covariance closer to the estimated one and vice-versa. Therefore, Q<sup>k</sup> can be adaptively updated in order to balance any deviations between the theoretical and estimated innovation covariance by considering a scaling factor αk, as given in Eq. (25). The estimate of Qk, denoted by Q^ <sup>k</sup>, is given in Eq. (26).

$$\alpha\_{k} = \frac{\text{trace}(\frac{1}{m} \sum\_{i=0}^{m-1} \mathbf{v}\_{k-i} \mathbf{v}\_{k-i}^{H} - \mathbf{R}\_{k})}{\text{trace}(\mathbf{H}\_{k} \mathbf{P}\_{k|k-1} \mathbf{H}\_{k}^{H})} \tag{25}$$

$$
\hat{\mathbf{Q}}\_k = \,\,\alpha\_k \hat{\mathbf{Q}}\_{k-1} \tag{26}
$$

In the case of EKF, the same procedure can be followed for adaptive estimation of noise covariances, by replacing the measurement matrix with the Jacobian matrix.

following, we first describe the general principles of CPE. Later, we explain our proposed

Consider an m-ary QAM received signal on single polarization, which is sampled and compensated for linear impairments. Assuming perfect linear equalization, the k-th input signal to the CPE can be written as in Eq. (27). Here, rk denotes the k-th input signal to CPE, ak denotes the transmitted symbol, and nk denotes the collective amplified spontaneous emission (ASE) noise which is assumed to be white Gaussian process. θ<sup>k</sup> denotes the phase noise arising from the laser linewidth effects and fiber nonlinearity, which is typically modeled as a Wiener process and is given in Eq. (28). Figure 1 (a) describes the input signal model to CPE. It can be seen that after ak is rotated by phase noise θk, nk further adds additional phase noise n<sup>0</sup>

amplitude noise n~k. The objective of CPE is to estimate the phase noise θk, and derotate the received signal rk, in order to recover the transmitted symbol ak, as given in Eq. (29) and Figure 1 (b). However, since the CPE targets at estimating an accurate θ^k, the recovered transmitted symbol ^ak, still suffers from the residual phase noise or amplitude noise or both.

The effects of nk, as discussed in Section 3.1.1, can be taken into account by reformulating Eq. (27) as given in Eq. (30) which forms the input signal to CPANE. Here, rk is modeled as the transmitted symbol ak being rotated by a complex quantity ψk, that considers the effects of both phase and amplitude noise in its real and imaginary parts, respectively, as given in Eq. (30).

rk <sup>¼</sup> ake<sup>j</sup>θ<sup>k</sup> <sup>þ</sup> nk (27)

Applications of Kalman Filters for Coherent Optical Communication Systems

http://dx.doi.org/10.5772/intechopen.71617

θ<sup>k</sup> ¼ θ<sup>k</sup>�<sup>1</sup> þ wk (28)

^ak <sup>¼</sup> rke�<sup>j</sup>θ^<sup>k</sup> (29)

I

aˆk(CPE) aˆk(CPANE)

˜nk

rk

θk n

<sup>k</sup> a<sup>k</sup>

(b)

Q

<sup>k</sup> and

213

CPANE algorithm illustrating its principles and implementation details using EKF.

3.1.1. Principles of CPE

For more details, please refer [12].

3.1.2. Principles of CPANE

˜nk

θk n k

(a)

n<sup>k</sup>

rk

I

Figure 1. (a) Input signal model to CPE (b) recovered symbols using CPE and CPANE [54].

ak

akejθ*<sup>k</sup>*

Q

### 3. Kalman filtering for coherent optical communications

### 3.1. Kalman filtering for carrier phase and amplitude noise estimation (CPANE)

Digital carrier phase estimation (CPE) has become an essential component of coherent optical receivers to recover the carrier phase perturbed by laser phase noise arising from the transmitter laser or LO [35–43]. Several CPE techniques have been developed in the literature based on feedback [39, 40] or feed forward loops [35–37]. Depending on how the data phase is wiped off, they can be further classified into decision directed (DD) [35, 40, 42, 46] or non-decision directed (NDD) methods [39, 41, 43]. NDD methods like Viterbi-Viterbi [41] CPE has gained high attention due to its ease of implementation. However, it employs m-th power scheme to remove the data modulation and therefore, are only better suited for QPSK systems. However, for higher QAM systems, DD-CPE methods exhibit better performance compared to NDD CPE methods [35, 42].

Apart from tracking the carrier phase, CPE being a low complex technique, can also be employed for compensating the nonlinear phase shift stemming from the Kerr nonlinear effects [30, 44–47]. However, the nonlinear mitigation performance of CPE is limited in the presence of ASE noise and at high launch powers. Moreover, a phase unwrapping function is typically required for CPE that might increase the probability of cycle slips [35, 49]. Addressing these problems, we have proposed a CPANE algorithm using EKF in [12, 53] for the joint mitigation of linear and nonlinear phase noise as well as ASE induced phase and amplitude distortions. Unlike CPE, EKF-CPANE estimates a complex quantity, and therefore, no argument function is required which eliminates the ambiguity associated with multiples of 2π and consequent cycle slips.

Kalman filter based CPE has been introduced and numerically verified in [52]. From the numerical results, it was reported that the Kalman based phase estimation combined with DD equalizer in a feedback configuration outperforms the conventional CMA based approach [52]. CPE based on EKF was demonstrated and verified experimentally for QPSK and 16-QAM systems in [57]. In [55], EKF has been investigated for characterizing the laser phase and amplitude noise. EKF based carrier synchronization has also been verified experimentally, in combination with expectation maximization (EM). A carrier recovery scheme based on block estimation process with Kalman filter has been demonstrated in [56]. This approach was verified experimentally for 16 and 64-QAM signals. However, these Kalman filter based approaches estimate an argument which involves sine and cosine functions, computation of the Jacobian matrix and also require matrix multiplications and inversions, which increases the computational complexity. The proposed method in [12, 53], estimates a complex quantity accounting also for the phase and amplitude distortions arising from the ASE noise in addition to the carrier phase. The variables in the SSM reduce to scalars and therefore, the vectors and matrices are reduced to scalars which will ease the computational effort. In the following, we first describe the general principles of CPE. Later, we explain our proposed CPANE algorithm illustrating its principles and implementation details using EKF.

#### 3.1.1. Principles of CPE

In the case of EKF, the same procedure can be followed for adaptive estimation of noise

Digital carrier phase estimation (CPE) has become an essential component of coherent optical receivers to recover the carrier phase perturbed by laser phase noise arising from the transmitter laser or LO [35–43]. Several CPE techniques have been developed in the literature based on feedback [39, 40] or feed forward loops [35–37]. Depending on how the data phase is wiped off, they can be further classified into decision directed (DD) [35, 40, 42, 46] or non-decision directed (NDD) methods [39, 41, 43]. NDD methods like Viterbi-Viterbi [41] CPE has gained high attention due to its ease of implementation. However, it employs m-th power scheme to remove the data modulation and therefore, are only better suited for QPSK systems. However, for higher QAM systems, DD-CPE methods exhibit better performance compared to NDD

Apart from tracking the carrier phase, CPE being a low complex technique, can also be employed for compensating the nonlinear phase shift stemming from the Kerr nonlinear effects [30, 44–47]. However, the nonlinear mitigation performance of CPE is limited in the presence of ASE noise and at high launch powers. Moreover, a phase unwrapping function is typically required for CPE that might increase the probability of cycle slips [35, 49]. Addressing these problems, we have proposed a CPANE algorithm using EKF in [12, 53] for the joint mitigation of linear and nonlinear phase noise as well as ASE induced phase and amplitude distortions. Unlike CPE, EKF-CPANE estimates a complex quantity, and therefore, no argument function is required which eliminates the ambiguity associated with multiples of 2π and

Kalman filter based CPE has been introduced and numerically verified in [52]. From the numerical results, it was reported that the Kalman based phase estimation combined with DD equalizer in a feedback configuration outperforms the conventional CMA based approach [52]. CPE based on EKF was demonstrated and verified experimentally for QPSK and 16-QAM systems in [57]. In [55], EKF has been investigated for characterizing the laser phase and amplitude noise. EKF based carrier synchronization has also been verified experimentally, in combination with expectation maximization (EM). A carrier recovery scheme based on block estimation process with Kalman filter has been demonstrated in [56]. This approach was verified experimentally for 16 and 64-QAM signals. However, these Kalman filter based approaches estimate an argument which involves sine and cosine functions, computation of the Jacobian matrix and also require matrix multiplications and inversions, which increases the computational complexity. The proposed method in [12, 53], estimates a complex quantity accounting also for the phase and amplitude distortions arising from the ASE noise in addition to the carrier phase. The variables in the SSM reduce to scalars and therefore, the vectors and matrices are reduced to scalars which will ease the computational effort. In the

covariances, by replacing the measurement matrix with the Jacobian matrix.

3.1. Kalman filtering for carrier phase and amplitude noise estimation (CPANE)

3. Kalman filtering for coherent optical communications

CPE methods [35, 42].

212 Kalman Filters - Theory for Advanced Applications

consequent cycle slips.

Consider an m-ary QAM received signal on single polarization, which is sampled and compensated for linear impairments. Assuming perfect linear equalization, the k-th input signal to the CPE can be written as in Eq. (27). Here, rk denotes the k-th input signal to CPE, ak denotes the transmitted symbol, and nk denotes the collective amplified spontaneous emission (ASE) noise which is assumed to be white Gaussian process. θ<sup>k</sup> denotes the phase noise arising from the laser linewidth effects and fiber nonlinearity, which is typically modeled as a Wiener process and is given in Eq. (28). Figure 1 (a) describes the input signal model to CPE. It can be seen that after ak is rotated by phase noise θk, nk further adds additional phase noise n<sup>0</sup> <sup>k</sup> and amplitude noise n~k. The objective of CPE is to estimate the phase noise θk, and derotate the received signal rk, in order to recover the transmitted symbol ak, as given in Eq. (29) and Figure 1 (b). However, since the CPE targets at estimating an accurate θ^k, the recovered transmitted symbol ^ak, still suffers from the residual phase noise or amplitude noise or both. For more details, please refer [12].

$$r\_k = a\_k \mathbf{e}^{i\theta\_k} + n\_k \tag{27}$$

$$
\theta\_k = \, \theta\_{k-1} + \varpi\_k \tag{28}
$$

$$
\hat{a}\_k = r\_k \mathbf{e}^{-j\hat{\Theta}\_k} \tag{29}
$$

#### 3.1.2. Principles of CPANE

The effects of nk, as discussed in Section 3.1.1, can be taken into account by reformulating Eq. (27) as given in Eq. (30) which forms the input signal to CPANE. Here, rk is modeled as the transmitted symbol ak being rotated by a complex quantity ψk, that considers the effects of both phase and amplitude noise in its real and imaginary parts, respectively, as given in Eq. (30).

Figure 1. (a) Input signal model to CPE (b) recovered symbols using CPE and CPANE [54].

The objective of CPANE is to recover θ<sup>k</sup> more accurately, by estimating the complex quantity ψk. The recovered transmitted symbol ^ak is given in Eq. (32). Since ψ<sup>k</sup> takes into account, both the phase and amplitude distortions, ^ak can be recovered more accurately by employing CPANE compared to CPE, as depicted in Figure 1 (b). Moreover, unlike CPE, CPANE eliminates the necessity of phase unwrapping function.

$$r\_k = a\_k \mathbf{e}^{i\psi\_k} \tag{30}$$

simulations on single channel systems in [12, 53] and multi-channel systems in [54]. Here, we briefly discuss the numerical model and present a few simulation results reproduced from [12], so that the flow of the readers is not interrupted. The numerical model of polarization multiplexed (PM) m-QAM coherent transmission system including a DSP module at the receiver, is depicted in Figure 3. Here, we consider the PM-m-QAM transmitter with m = 16 and 64, operated at 28 and 18.667 GBaud, respectively. These signals are transmitted through a standard single mode fiber (SSMF) link at different launch powers. The SSMF has the following parameters: attenuation coefficient (α) = 0.2 dB/km, dispersion coefficient (D) = 16 ps/nm-km, and nonlinearity coefficient (γ) = 1.2/W-km. The span length of SSMF is 80 km and a number of 12 and 6 spans have been considered for 16 and 64 QAM, respectively, yielding a total transmission distance of 960 and 480 km. The span losses are compensated by an EDFA with a gain of 16 dB and noise figure (NF) of 4 dB. For simplicity, PMD has been neglected in this study. At the receive end, we employ a dual polarization coherent receiver which is followed by a DSP module. The laser linewidth of the LO has been set to 100 kHz. After coherent detection, the signals are re-sampled to twice the symbol rate and are followed by linear compensation. Then the signals are further down sampled to the symbol rate and are further processed by the EKF-CPANE for mitigating linear and nonlinear phase noise besides amplitude noise. The performance of EKF-CPANE is compared to

Applications of Kalman Filters for Coherent Optical Communication Systems

http://dx.doi.org/10.5772/intechopen.71617

215

$$
\psi\_k = (\theta\_k + n'\_k) + \mathrm{j}\tilde{n}\_k \tag{31}
$$

$$
\hat{a}\_k = r\_k \mathbf{e}^{-j\ddot{\psi}\_k} \tag{32}
$$

#### 3.1.3. EKF-CPANE for joint mitigation of phase and amplitude noise

As discussed earlier, CPANE algorithm can be employed for the joint mitigation of phase and amplitude noise. However, it requires a reliable tracking of the complex quantity ψk, which can be accomplished by an EKF. The required SSM for the EKF can be formulated using Eqs. (33) and (34). Eq. (33) represents the state or process equation that describes the time evolution of ψk. Eq. (34) represents the observation equation that describes the relation of the states ψ<sup>k</sup> to the observations rk. Eq. (34) is similar to Eq. (30), however, for consistency of the filter, the measurement noise mk has been taken into account. Here, all the variables in the SSM are scalar quantities. Comparing to the standard SSM for EKF described in Section 2.2, it can be noted that the state transition is identity and the measurement matrix is the transmitted symbol, ak, for simplicity, we call it measurement weight (MW), since it is a scalar. The EKF recursive equations can be derived analogously by relating the SSM to the standard SSM of EKF discussed in Section 2.2. Since the MW ak is required to compute the update equations, which is not known apriori, EKF-CPANE is DD. The required decisions of ak, denoted by dk are obtained by de-rotating rk with an average of the past updated estimates <sup>ψ</sup>^<sup>k</sup> over a window length of N, as given in Eq. (35). For more details on the prediction and update equations of EKF-CPANE, please refer [12]. Figure 2 depicts the schematic of the EKF-CPANE algorithm, illustrating that the prediction <sup>ψ</sup>^<sup>k</sup>jk�<sup>1</sup> is the delayed version of the past updated estimate and the current updated state <sup>ψ</sup>^<sup>k</sup>j<sup>k</sup> is the linear combination of the prediction <sup>ψ</sup>^<sup>k</sup>jk�<sup>1</sup> and the innovation vk weighted by the Kalman gain Kk. The process of making the required decisions for the update step has also been illustrated in Figure 2.

$$
\psi\_k = \,\,\psi\_{k-1} + \,\text{w}\_k \tag{33}
$$

$$r\_k = a\_k \mathbf{e}^{i\psi\_k} + m\_k \tag{34}$$

224 Gb/s

PM-m-QAM Tx 80km

*rk*

*r*ˆ*k|k−*<sup>1</sup>

SSMF α = 0.2 dB/km *D* = 16 ps/nm-km γ = 1.2 /km-W

Figure 2. Block diagram of EKF-CPANE algorithm [54].

EDFA PBS

Figure 3. Numerical model of PM-m-QAM coherent transmission system with DSP module [12].

x N

*− Kk* +

e *<sup>j</sup>*(*.*)

e*−j*(*.*)

*vk*

Decision

*dk*

*tk*

LO

Dual Polarization Coherent

Receiver

Delay

1 *<sup>N</sup>* ∑ *N* ψˆ *<sup>k</sup>−<sup>N</sup>*

<sup>ψ</sup><sup>ˆ</sup> *<sup>k</sup>|k−*<sup>1</sup>

I*<sup>x</sup>* +jQ*<sup>x</sup>*

I*<sup>y</sup>* +jQ*<sup>y</sup>*

Linear compensation *a*ˆ*k*

<sup>ψ</sup><sup>ˆ</sup> *<sup>k</sup>|<sup>k</sup>* <sup>e</sup>*−j*(*.*)

CPANE/CPE

BER Evaluation

$$d\_k = \text{decision}(t\_k) \text{ where } t\_k = r\_k \mathbf{e}^{-\frac{1}{N} \sum\_N \hat{\psi}\_{k-N}} \tag{35}$$

#### 3.1.4. Numerical analysis of EKF-CPANE

The performance of EKF-CPANE algorithm for mitigating the laser phase noise, fiber nonlinearity besides the ASE induced phase and amplitude distortions has been verified through numerical simulations on single channel systems in [12, 53] and multi-channel systems in [54]. Here, we briefly discuss the numerical model and present a few simulation results reproduced from [12], so that the flow of the readers is not interrupted. The numerical model of polarization multiplexed (PM) m-QAM coherent transmission system including a DSP module at the receiver, is depicted in Figure 3. Here, we consider the PM-m-QAM transmitter with m = 16 and 64, operated at 28 and 18.667 GBaud, respectively. These signals are transmitted through a standard single mode fiber (SSMF) link at different launch powers. The SSMF has the following parameters: attenuation coefficient (α) = 0.2 dB/km, dispersion coefficient (D) = 16 ps/nm-km, and nonlinearity coefficient (γ) = 1.2/W-km. The span length of SSMF is 80 km and a number of 12 and 6 spans have been considered for 16 and 64 QAM, respectively, yielding a total transmission distance of 960 and 480 km. The span losses are compensated by an EDFA with a gain of 16 dB and noise figure (NF) of 4 dB. For simplicity, PMD has been neglected in this study. At the receive end, we employ a dual polarization coherent receiver which is followed by a DSP module. The laser linewidth of the LO has been set to 100 kHz. After coherent detection, the signals are re-sampled to twice the symbol rate and are followed by linear compensation. Then the signals are further down sampled to the symbol rate and are further processed by the EKF-CPANE for mitigating linear and nonlinear phase noise besides amplitude noise. The performance of EKF-CPANE is compared to

Figure 2. Block diagram of EKF-CPANE algorithm [54].

The objective of CPANE is to recover θ<sup>k</sup> more accurately, by estimating the complex quantity ψk. The recovered transmitted symbol ^ak is given in Eq. (32). Since ψ<sup>k</sup> takes into account, both the phase and amplitude distortions, ^ak can be recovered more accurately by employing CPANE compared to CPE, as depicted in Figure 1 (b). Moreover, unlike CPE, CPANE elimi-

ψ<sup>k</sup> ¼ ðθ<sup>k</sup> þ n<sup>0</sup>

As discussed earlier, CPANE algorithm can be employed for the joint mitigation of phase and amplitude noise. However, it requires a reliable tracking of the complex quantity ψk, which can be accomplished by an EKF. The required SSM for the EKF can be formulated using Eqs. (33) and (34). Eq. (33) represents the state or process equation that describes the time evolution of ψk. Eq. (34) represents the observation equation that describes the relation of the states ψ<sup>k</sup> to the observations rk. Eq. (34) is similar to Eq. (30), however, for consistency of the filter, the measurement noise mk has been taken into account. Here, all the variables in the SSM are scalar quantities. Comparing to the standard SSM for EKF described in Section 2.2, it can be noted that the state transition is identity and the measurement matrix is the transmitted symbol, ak, for simplicity, we call it measurement weight (MW), since it is a scalar. The EKF recursive equations can be derived analogously by relating the SSM to the standard SSM of EKF discussed in Section 2.2. Since the MW ak is required to compute the update equations, which is not known apriori, EKF-CPANE is DD. The required decisions of ak, denoted by dk are obtained by de-rotating rk with an average of the past updated estimates <sup>ψ</sup>^<sup>k</sup> over a window length of N, as given in Eq. (35). For more details on the prediction and update equations of EKF-CPANE, please refer [12]. Figure 2 depicts the schematic of the EKF-CPANE algorithm, illustrating that the prediction <sup>ψ</sup>^<sup>k</sup>jk�<sup>1</sup> is the delayed version of the past updated estimate and the current updated state <sup>ψ</sup>^<sup>k</sup>j<sup>k</sup> is the linear combination of the prediction <sup>ψ</sup>^<sup>k</sup>jk�<sup>1</sup> and the innovation vk weighted by the Kalman gain Kk. The process of making the required decisions

rk <sup>¼</sup> ake<sup>j</sup>ψ<sup>k</sup> (30)

^ak <sup>¼</sup> rke�<sup>j</sup>ψ^<sup>k</sup> (32)

<sup>ψ</sup><sup>k</sup> <sup>¼</sup> <sup>ψ</sup><sup>k</sup>�<sup>1</sup> <sup>þ</sup> wk (33)

rk <sup>¼</sup> ake<sup>j</sup>ψ<sup>k</sup> <sup>þ</sup> mk (34)

<sup>N</sup>ψ^ <sup>k</sup>�<sup>N</sup> (35)

�<sup>j</sup> <sup>1</sup> N X

<sup>k</sup>Þ þ jn~<sup>k</sup> (31)

nates the necessity of phase unwrapping function.

214 Kalman Filters - Theory for Advanced Applications

3.1.3. EKF-CPANE for joint mitigation of phase and amplitude noise

for the update step has also been illustrated in Figure 2.

3.1.4. Numerical analysis of EKF-CPANE

dk ¼ decisionðtkÞ where tk ¼ rke

The performance of EKF-CPANE algorithm for mitigating the laser phase noise, fiber nonlinearity besides the ASE induced phase and amplitude distortions has been verified through numerical

Figure 3. Numerical model of PM-m-QAM coherent transmission system with DSP module [12].

feedforward DD-CPE [46], feedback DD phase locked loop (DD-PLL) [36] and a NDD universal CPE (U-CPE) [39]. The noise covariances for EKF-CPANE, the tap length or step size for DD-CPE, DD-PLL and U-CPE are set to optimize the performance.

and nonlinear phase noise as well as amplitude noise simultaneously besides less computational effort. Although, EKF-CPANE outperforms several other considered CPE methods, the effectiveness of EKF-CPANE in mitigating fiber nonlinear effects can be further enhanced by reducing the number of errors in the pre-decisions dk. We have proposed a weighted innovation approach (WIA) in [12], where the innovation is computed as a weighted combination of the two nearest likely decisions. Although a gain of ≈ 0.3 dB in the Q-factor can be obtained compared to conventional EKF-CPANE, in the linear regime, no notable improvement can be seen in the nonlinear regime. On the other hand, DBP has emerged to be an effective technique in mitigating linear and nonlinear impairments simultaneously, provided the channel parameters are known a-priori and the step size is sufficiently small. However, DBP can compensate only the deterministic impairments of self-phase modulation and its performance deteriorates significantly in the presence of stochastic impairments like laser phase noise, ASE and NLPN. Moreover, the required huge computational effort keeps it far away from real-time implementation. Nevertheless, by employing a few DBP steps prior to EKF-CPANE would yield an enhanced tolerance towards nonlinearities since DBP is well capable of mitigating deterministic impairments and EKF takes into account the stochastic nature of ASE noise and NLPN. By partially compensating fiber nonlinear effects employing few DBP steps prior to EKF, would result in improved pre-decisions and thereby facilitates the residual compensation of nonlinearities along with amplitude and phase noise effectively. These theoretical findings are verified through numerical simulations on both single [31] and multichannel systems [58].

Applications of Kalman Filters for Coherent Optical Communication Systems

http://dx.doi.org/10.5772/intechopen.71617

217

In [31], it was reported that the EKF-CPANE outperforms the asymmetric split step Fourier method (ASSFM) based one step per span (OSPS) DBP with optimized nonlinear co-efficient γ (ODBP), for single channel systems, for transmission on both SSFM and non-zero dispersion shifted fiber (NZ-DSF). A detailed investigation has also been carried out on the combined performance of DBP and EKF-CPANE with an analysis on the influence of the nonlinear coefficient and the step size of DBP when employed prior to EKF-CPANE. The numerical model employed in this study is similar to the one discussed in Section 3.1.4, with a few changes in the parameters of NF being 5 dB and the linewidth of LO being 500 kHz. The influence of DBP step size on the combined performance of DBP and EKF-CPANE for both SSMF as well as NZ-DSF transmission is illustrated in Figure 5(a) [31]. Here, OCDBP denotes the optimized DBP which has a nonlinear coefficient different from ODBP when employed prior to EKF. A worth noting result is that at a launch power of 3 dBm and a transmission distance of 960 km, a gain of 1 dB in the Q-factor can be obtained by employing 0.3 DBP steps per span prior to EKF-CPANE, for both SSMF and NZ-DSF transmission. At the expense of additional computational effort, the deployment of a few DBP steps prior to EKF-CPANE

For the case of multi-channel systems, also, a detailed analysis has been performed in [58], on the combined performance of DBP and EKF for mitigation of inter and intra channel nonlinearities besides phase and amplitude noise. Here, the DBP is employed by considering the temporal correlations between the neighboring signal samples and is termed as correlated DBP (CDBP) [27, 28]. This approach will improve the accuracy in computing the nonlinear phase shift and there by enhances the nonlinear mitigation performance. Since the optimization of nonlinear coefficient plays a vital role on the performance of DBP, we proposed an amplitude

further enhances its performance trading off to complexity.

The bit error rate (BER) performance of the considered algorithms is evaluated and a Q-factor is computed as 20log10 erfcinv(2\*BER). The Q-factor vs. launch power curves for 16-QAM and 64-QAM are depicted in Figure 4(a) and (b), respectively. It can be seen that EKF-CPANE exhibits better performance compared to DD-CPE, DD-PLL and U-CPE in both linear and nonlinear regimes. This performance enhancement is better visible compared to the DD-CPE method. For PM-64-QAM, it can also be seen that the DD-CPE experiences cycle slips occurring through the error propagation of wrong decisions which can be seen in Figure 4(b) at launch powers ranging from 2 to 1 dBm [12]. Since the performance of DD algorithms strongly depends on the pre-decisions made by the algorithm, we study the impact of ideal error free decisions on their performance by replacing the pre-decisions dk with the true data symbols ak. The algorithms with the ideal case are denoted by IEKF-CPANE, IDD-CPE and IDD-PLL. It can be seen from Figure 4(a) and (b), that the IEKF-CPANE shows significant performance enhancement and better tolerance towards linear and nonlinear phase noise as well as amplitude noise, compared to IDD-CPE and IDD-PLL. Unlike EKF-CAPNE, no notable improvement can be obtained for the DD-CPE and DD-PLL between their practical and ideal cases. Although, the ideal case, where the true symbols ak are already known, is not possible in practice, it should be noted that the performance of EKF-CPANE can be further improved by reducing the number of decision errors, which will be further discussed in the next Section 3.2.

#### 3.2. EKF and DBP for fiber nonlinear mitigation

In Section 3.1, we have described how the EKF can be employed for the joint mitigation of phase and amplitude noise. From the numerical results discussed in Section 3.1.4, it can be concluded that the EKF-CAPNE algorithm shows promising results in mitigating the linear

Figure 4. Q-factor vs. launch power curves for the considered algorithms (a) PM-16-QAM transmission over 960 km of SSMF transmission (b) PM-64-QAM transmission over 480 km of SSMF transmission [12].

and nonlinear phase noise as well as amplitude noise simultaneously besides less computational effort. Although, EKF-CPANE outperforms several other considered CPE methods, the effectiveness of EKF-CPANE in mitigating fiber nonlinear effects can be further enhanced by reducing the number of errors in the pre-decisions dk. We have proposed a weighted innovation approach (WIA) in [12], where the innovation is computed as a weighted combination of the two nearest likely decisions. Although a gain of ≈ 0.3 dB in the Q-factor can be obtained compared to conventional EKF-CPANE, in the linear regime, no notable improvement can be seen in the nonlinear regime. On the other hand, DBP has emerged to be an effective technique in mitigating linear and nonlinear impairments simultaneously, provided the channel parameters are known a-priori and the step size is sufficiently small. However, DBP can compensate only the deterministic impairments of self-phase modulation and its performance deteriorates significantly in the presence of stochastic impairments like laser phase noise, ASE and NLPN. Moreover, the required huge computational effort keeps it far away from real-time implementation. Nevertheless, by employing a few DBP steps prior to EKF-CPANE would yield an enhanced tolerance towards nonlinearities since DBP is well capable of mitigating deterministic impairments and EKF takes into account the stochastic nature of ASE noise and NLPN. By partially compensating fiber nonlinear effects employing few DBP steps prior to EKF, would result in improved pre-decisions and thereby facilitates the residual compensation of nonlinearities along with amplitude and phase noise effectively. These theoretical findings are verified through numerical simulations on both single [31] and multichannel systems [58].

feedforward DD-CPE [46], feedback DD phase locked loop (DD-PLL) [36] and a NDD universal CPE (U-CPE) [39]. The noise covariances for EKF-CPANE, the tap length or step size for DD-CPE,

The bit error rate (BER) performance of the considered algorithms is evaluated and a Q-factor is computed as 20log10 erfcinv(2\*BER). The Q-factor vs. launch power curves for 16-QAM and 64-QAM are depicted in Figure 4(a) and (b), respectively. It can be seen that EKF-CPANE exhibits better performance compared to DD-CPE, DD-PLL and U-CPE in both linear and nonlinear regimes. This performance enhancement is better visible compared to the DD-CPE method. For PM-64-QAM, it can also be seen that the DD-CPE experiences cycle slips occurring through the error propagation of wrong decisions which can be seen in Figure 4(b) at launch powers ranging from 2 to 1 dBm [12]. Since the performance of DD algorithms strongly depends on the pre-decisions made by the algorithm, we study the impact of ideal error free decisions on their performance by replacing the pre-decisions dk with the true data symbols ak. The algorithms with the ideal case are denoted by IEKF-CPANE, IDD-CPE and IDD-PLL. It can be seen from Figure 4(a) and (b), that the IEKF-CPANE shows significant performance enhancement and better tolerance towards linear and nonlinear phase noise as well as amplitude noise, compared to IDD-CPE and IDD-PLL. Unlike EKF-CAPNE, no notable improvement can be obtained for the DD-CPE and DD-PLL between their practical and ideal cases. Although, the ideal case, where the true symbols ak are already known, is not possible in practice, it should be noted that the performance of EKF-CPANE can be further improved by reducing the number of decision errors, which will be further discussed in the next Section 3.2.

In Section 3.1, we have described how the EKF can be employed for the joint mitigation of phase and amplitude noise. From the numerical results discussed in Section 3.1.4, it can be concluded that the EKF-CAPNE algorithm shows promising results in mitigating the linear

> *−*4 *−*2

Q-factor (dB)

Figure 4. Q-factor vs. launch power curves for the considered algorithms (a) PM-16-QAM transmission over 960 km of

SSMF transmission (b) PM-64-QAM transmission over 480 km of SSMF transmission [12].

*−*6 *−*4 *−*20 2 4 6

DD-CPE IDD-CPE DD-PLL IDD-PLL U-CPE EKF-CPANE

IEKF-CPANE

Launch Power (dBm)

(b)

DD-PLL and U-CPE are set to optimize the performance.

216 Kalman Filters - Theory for Advanced Applications

3.2. EKF and DBP for fiber nonlinear mitigation

*−*6 *−*4 *−*20 2 4 6

IEKF-CPANE

DD-CPE IDD-CPE DD-PLL IDD-PLL U-CPE EKF-CPANE

Launch Power (dBm)

(a)

4

6

Q-factor (dB) 8

10

In [31], it was reported that the EKF-CPANE outperforms the asymmetric split step Fourier method (ASSFM) based one step per span (OSPS) DBP with optimized nonlinear co-efficient γ (ODBP), for single channel systems, for transmission on both SSFM and non-zero dispersion shifted fiber (NZ-DSF). A detailed investigation has also been carried out on the combined performance of DBP and EKF-CPANE with an analysis on the influence of the nonlinear coefficient and the step size of DBP when employed prior to EKF-CPANE. The numerical model employed in this study is similar to the one discussed in Section 3.1.4, with a few changes in the parameters of NF being 5 dB and the linewidth of LO being 500 kHz. The influence of DBP step size on the combined performance of DBP and EKF-CPANE for both SSMF as well as NZ-DSF transmission is illustrated in Figure 5(a) [31]. Here, OCDBP denotes the optimized DBP which has a nonlinear coefficient different from ODBP when employed prior to EKF. A worth noting result is that at a launch power of 3 dBm and a transmission distance of 960 km, a gain of 1 dB in the Q-factor can be obtained by employing 0.3 DBP steps per span prior to EKF-CPANE, for both SSMF and NZ-DSF transmission. At the expense of additional computational effort, the deployment of a few DBP steps prior to EKF-CPANE further enhances its performance trading off to complexity.

For the case of multi-channel systems, also, a detailed analysis has been performed in [58], on the combined performance of DBP and EKF for mitigation of inter and intra channel nonlinearities besides phase and amplitude noise. Here, the DBP is employed by considering the temporal correlations between the neighboring signal samples and is termed as correlated DBP (CDBP) [27, 28]. This approach will improve the accuracy in computing the nonlinear phase shift and there by enhances the nonlinear mitigation performance. Since the optimization of nonlinear coefficient plays a vital role on the performance of DBP, we proposed an amplitude

of 3.2, we investigated the combined performance of SSNL and EKF-CPANE for mitigating the

The numerical model of PM-16-QAM coherent transmission system over DM link [59] is depicted in Figure 6. Here, a fully compensated periodical DM link with several spans has been considered. Each span consists of 80 km of SSMF and 17 km of dispersion compensating fiber (DCF). The SSMF has the following parameters: α = 0.2 dB/km, D = 17 dB/nm-km, γ = 1.2/W-km. The parameters of DCF are given by: α = 0.5 dB/km, D = 80 dB/nm-km and γ = 5/W-km. In this study, the input power to DCF was set to half of the input power to SSMF. Therefore, the gains of EDFA1 and EDFA2 are adjusted accordingly, to compensate the span losses. The NF of both the EDFAs are set to 4 dB. As described earlier, after coherent detection, the signals are further processed by the SSNL and EKF-CPANE algorithms for mitigating fiber nonlinearities. It has been reported in [59], that the combined performance of SSNL and EKF yields an improved tolerance towards nonlinearities of up to 2 dB for a transmission distance of 1200 km and at a

. Further, their combined performance increases the transmission reach by ≈ 250 km

EDFA2 PBS

x N

Applications of Kalman Filters for Coherent Optical Communication Systems

http://dx.doi.org/10.5772/intechopen.71617

LO

SSNL EKF SSNL+EKF Dual Polarization Coherent Receiver

SSNL/CPANE

BER evaluation 219

DCF

200 400 600 800 1000 1200 1400 1600 1800

Figure 7. BER vs. transmission distance curves for PM-16-QAM coherent transmission system over DM link [59].

Transmission distance (km)

≈250 km

at a launch power of 3 dBm and at a BER of 2\*10<sup>2</sup> as depicted in Figure 7.

EDFA1

Figure 6. Simulation model of PM-16-QAM coherent transmission over DM link with DSP module [59].

fiber nonlinearity in DM links [59].

BER of 2\*10<sup>2</sup>

224 Gb/s

PM-16-QAM Tx SSMF

10*−*<sup>4</sup>

10*−*<sup>3</sup>

BER

10*−*<sup>2</sup>

10*−*<sup>1</sup>

Figure 5. (a) Q-factor vs. DBP step size for the combined performance of EKF and DBP after 960 km of SSMF/NZ-DSF transmission at a launch power of 3 dBm [31] (b) Q-factor vs. number of WDM channels after 960 km of SSMF transmission at launch power of 3 dBm [58].

dependent optimization (AO) [58] of the nonlinear coefficient, according to the discrete amplitude levels present in the higher order modulation formats like 16-QAM. The combined performance of AO-CDBP and EKF-CPANE for WDM systems with varying number of channels has been investigated in [58]. Analogous to the single channel systems, the combined performance of AO-CDBP and EKF yields an improved performance also for the WDM case. However, with increasing impact of the cross phase modulation (XPM) as the number of channels increase, the gain obtained from their combined performance starts vanishing which can be observed in Figure 5(b).

#### 3.3. EKF for mitigation of nonlinearities in dispersion managed links

Since the advent of coherent detection and DSP for coherent optical receivers, CD can be effectively compensated by digital equalization in the electric domain ad thereby, eliminating the need for dispersion compensating fibers (DCF). However, nonlinear mitigation in the dispersion managed (DM) links is also vital in order to upgrade existing links. Although, the computational complexity of DBP is quite high, for DM links, the DBP algorithm can be simplified by assuming that the nonlinear behavior repeats itself every span and therefore, the total nonlinearity after N spans of transmission can be approximated to N times the nonlinearity from a single span [70]. This is termed as distance folded DBP [70] and it reduces the complexity by a factor of N assuming the step size of DBP is equal to the span length and the span length is assumed to be constant. Assuming the dispersion is fully compensated in each span, only the nonlinear term in the nonlinear Schrödinger equation (NLSE) can be solved in the time domain avoiding the Fourier and inverse Fourier transformation (FFT/IFFT) pairs which reduces the computational cost of DBP drastically. We call this approach single step nonlinearity mitigation (SSNL). Similar to the unmanaged links as discussed in the earlier section of 3.2, we investigated the combined performance of SSNL and EKF-CPANE for mitigating the fiber nonlinearity in DM links [59].

The numerical model of PM-16-QAM coherent transmission system over DM link [59] is depicted in Figure 6. Here, a fully compensated periodical DM link with several spans has been considered. Each span consists of 80 km of SSMF and 17 km of dispersion compensating fiber (DCF). The SSMF has the following parameters: α = 0.2 dB/km, D = 17 dB/nm-km, γ = 1.2/W-km. The parameters of DCF are given by: α = 0.5 dB/km, D = 80 dB/nm-km and γ = 5/W-km. In this study, the input power to DCF was set to half of the input power to SSMF. Therefore, the gains of EDFA1 and EDFA2 are adjusted accordingly, to compensate the span losses. The NF of both the EDFAs are set to 4 dB. As described earlier, after coherent detection, the signals are further processed by the SSNL and EKF-CPANE algorithms for mitigating fiber nonlinearities. It has been reported in [59], that the combined performance of SSNL and EKF yields an improved tolerance towards nonlinearities of up to 2 dB for a transmission distance of 1200 km and at a BER of 2\*10<sup>2</sup> . Further, their combined performance increases the transmission reach by ≈ 250 km at a launch power of 3 dBm and at a BER of 2\*10<sup>2</sup> as depicted in Figure 7.

Figure 6. Simulation model of PM-16-QAM coherent transmission over DM link with DSP module [59].

dependent optimization (AO) [58] of the nonlinear coefficient, according to the discrete amplitude levels present in the higher order modulation formats like 16-QAM. The combined performance of AO-CDBP and EKF-CPANE for WDM systems with varying number of channels has been investigated in [58]. Analogous to the single channel systems, the combined performance of AO-CDBP and EKF yields an improved performance also for the WDM case. However, with increasing impact of the cross phase modulation (XPM) as the number of channels increase, the gain obtained from their combined performance starts vanishing which

Figure 5. (a) Q-factor vs. DBP step size for the combined performance of EKF and DBP after 960 km of SSMF/NZ-DSF transmission at a launch power of 3 dBm [31] (b) Q-factor vs. number of WDM channels after 960 km of SSMF

8

10

Q-factor (dB) 12

14

2 4 6 8 10 12 14

AO-CDBP EKF AO-CDBP(3)+EKF AO-CDBP(6)+EKF
