A rank correlation coefficient based particle filter to estimate parameters in non-linear models

Abstract

Particle filtering algorithm has found an increasingly wide utilization in many fields at present, especially in non-linear and non-Gaussian situations. Because of the particle degeneracy limitation, various resampling methods have been researched. The estimation process of particle filtering algorithm is a series of weighted calculation processes, which can be regarded as weighted data fusion. This article proposed an improved particle filtering algorithm combining with different rank correlation coefficients to overcome the shortcomings of degeneracy. By simulating iteration operation in MATLAB, it discovers that the proposed algorithm provides better accuracy in comparison with particle filtering, Gaussian sum particle filter, and Gaussian mixture sigma-point particle filter in Gaussian mixture noise. A practical seven-dimensional harmonic model is also implemented in the simulation. After comparing the performances of different algorithms, we found that the proposed method had more accuracy than the widely used extended Kalman filtering algorithm.

Introduction

Since Gordon et al. ¹ research in 1993, particle filtering (PF) algorithm also known as sequential Monte Carlo (SMC) method has become a recent technique to perform filtering and smoothing for non-linear and non-Gaussian systems. This methodology has been adopted in various fields, including signal processing, navigation, target tracking, robotics, image processing, control, wireless communications, and economics.^2,3

PF algorithms have already been implemented to deal with multiple integrals in the fields of statistics and physics as early as the middle 1950s and been draw into fields of automation around the 1970s. Limited by the particle degeneracy and calculation speed then, the method was unvalued. Gordon et al. ¹ implemented sequential importance resampling (SIR) while they proposed the bootstrap filter. It is an appropriate solution of the particle degeneracy phenomenon. Although the SIR algorithm is advantageous, it has limitations. Much research on the improved algorithms from different perspective has been done and made some progress. Doucet and Johansen ² proposed that various particle methods could be reinterpreted as different instances of SMC methods. Pitt and Shephard proposed the auxiliary particle filter (APF) in Pitt and Shephard, ⁴ which introduced instrumental variable to correct the particle weights according to the likelihood. Van der Merwe and Wan ⁵ proposed the unscented particle filter (UPF). The algorithm takes advantage of the unscented transformation (UT) and the unscented Kalman filter (UKF) to achieve the proposal distribution. The sigma-point particle filter (SPPF) uses a sigma-point Kalman filter (SPKF) for proposal distribution generation and is an extension of the original UPF. ⁵ This kind of SPKF also includes its square-root form, such as square-root UKF (SRUKF) and square-root central difference KF (SRCDKF). They also proposed the Gaussian mixture sigma-point particle filter (GMSPPF) in Van der Merwe et al. ⁶ Kotecha and Djuric ⁷ applying the Gaussian distribution instead of the posterior probability density, which is called the Gaussian particle filter (GPF), avoided the resampling step. Similarly, in Kotecha and Djuric, ⁷ they proposed the Gaussian sum particle filter (GSPF) based on the GPF. ⁸

Initially, particle filter was called the bootstrap filter which is a powerful tool for Bayesian state estimation in non-linear systems. The key idea of this algorithm is to construct the posterior density function (pdf) of the state variables by a set of random samples (particles) with associated weights recursively. After decades of research carried out by many scholars, there is variety of PF algorithms. Almost all of them consist of three important operations: particle propagation, weight computation, and resampling. Particle propagation means the generation of particles, and weight computation amount to the generation of particles and assignment of weights, whereas resampling replaces one set of particles and their weights with another set. ⁹ Inspired by the idea of PF algorithm in Zhou et al., ¹⁰ multi-dimensional simulation in Zhou et al., ¹¹ and data fusion approaches in Zhou et al.,^12,13 this article proposed a novel approach for parameter estimation using PF algorithm and correlation coefficients. More details are in the following section.

This paper is organized as follows. In section “Problem statement,” we will present the hidden Markov model (HMM) and a harmonic parameter model in power system. Section “Proposed PF algorithm” shows the definitions and background information of traditional PF algorithm and different correlation coefficients. The new method that we combining the traditional PF algorithm with different correlation coefficients is also illustrated in this section. In section “Simulation performance and harmonic parameter estimation,” we present a numerical simulation and a practical seven-dimensional harmonic model using MATLAB. And finally, discussions and conclusions are stated in section “Conclusion.”

Problem statement

HMM

General state space models can be described in the forms of HMM as follows

$$x_k=f(x_{k-1},v_k)$$

(1)

$$y_k=h(x_k,w_k)$$

(2)

Equation (1) is the process equation and equation (2) is the measurement equation. In the above equation, $k=1,2,\dots$ is a discrete time series; the hidden state vector of the system $x_k$ satisfies $x_k\in {\mathbb{R}}^n$ at time k, and $y_k\in {\mathbb{R}}^m$ is the measurement vector produced by the system; $f(·)$ denotes the state space function or system transition function under the Markovian hypothesis and associated with the state variables iteration; $h(·)$ denotes the system measurement (observation) function associated with the series of measurements $(y_1,y_2,\dots )$ ; $v_k$ is the system process noise and $w_k$ is the measurement noise. It is worth noting that both Gaussian and non-Gaussian noises can be applied here.

We can consider the state of $x_k$ as a vector of stochastic processes. What we are most interested in is to estimate the probability density function (PDF) $p(x_k|y_k)$ of the state at each time step k conditioned on the whole set of measurements or observations until time k, $\{y_1,y_2,\dots ,y_k\}$ .

Harmonic parameter model

As is stated in He et al., ¹⁴ a power signal with harmonics in discrete form can be expressed using the following equation

$$y_k=\sum _{h=1}^H{A}_{k,h}\sin (h{\omega }_{k}k\mathrm{\Delta }t+{\phi }_{k,h})$$

(3)

where $A_{k,h}$ , $h{\omega }_k$ , and ${\phi }_{k,h}$ are the time-varying amplitude, frequency, and phase of the hth component, respectively. $h=1,2,\dots ,H$ are the harmonic orders, and H is the highest order of harmonic, and $h=1$ represents the fundamental harmonic. ${\omega }_k=2\pi f$ is the fundamental angular frequency and f denotes the utility frequency. $\mathrm{\Delta }t$ is sample interval associate with sampling frequency and k is the walk step.

In order to estimate the state vector of the power system at the current time step k, given a group of measurements or observations $Y_k=\{y_0,y_1,\dots ,y_k\}$ achieved at time steps $0,1,\dots ,k$ . The state vector $x_k$ includes important parameters of voltage/current in power system, which can be defined as

$$\begin{array}{c}\multicolumn{1}{c}{x_k={[A_{k,1},{\omega }_{k}k\mathrm{\Delta }t+{\phi }_{k,1},A_{k,2},2{\omega }_{k}k\mathrm{\Delta }t+{\phi }_{k,2},\dots ,A_{k,h},h{\omega }_{k}k\mathrm{\Delta }t+{\phi }_{k,h},\dots ,A_{k,H},H{\omega }_{k}k\mathrm{\Delta }t+{\phi }_{k,H},{\omega }_k]}^T}\end{array}$$

(4)

where the symbols are synonymous as above.

Having the measurement equations (2) and (3), process equation (1) and the state vector (4), we may establish the state transition matrix $F_k$ and measurement matrix $H_k$ as follows

$$F_k=\left[\begin{array}{ccccccccccc}1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 0\hfill \\ \multicolumn{1}{c}{}& 1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \mathrm{\Delta }t\hfill \\ \multicolumn{1}{c}{}& \hfill & 1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 0\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & 1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 2\mathrm{\Delta }t\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \ddots \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & 1\hfill & \hfill & \hfill & \hfill & \hfill & 0\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & \hfill & 1\hfill & \hfill & \hfill & \hfill & h\mathrm{\Delta }t\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \ddots \hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 1\hfill & \hfill & 0\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 1\hfill & H\mathrm{\Delta }t\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 1\hfill \end{array}\right]$$

(5)

$$H_k=[\sin \left(x_{k,2}\right),0,\sin \left(x_{k,4}\right),0,\dots ,\sin \left(x_{k,2h}\right),0,\dots \sin \left(x_{k,2H}\right),0,0]$$

(6)

What we should notice particularly is the dimensions of the above matrix and vectors. Since the dimension of state vector $x_k$ is $(2H+1)\times 1$ , the dimensions of state transition matrix $F_k$ and measurement matrix $H_k$ are $(2H+1)\times (2H+1)$ and $1\times (2H+1)$ , respectively. This would satisfy the HMM equations (1) and (2). The measurement matrix is different from that in He et al. ¹⁴ The reason is that the algorithm in He’s paper is based on extended Kalman filtering (EKF), which was the first-order approximation to the optimum. The method proposed in this article is based on Bayesian framework, which do not need the similar linearization and approximation process.

Since the high-order harmonics are very small in real application, we only choose the first (fundamental), third, and fifth harmonics for simulation model in section “Simulation performance and harmonic parameter estimation.” Thus, the dimension of the filter is 7-by-7.

Proposed PF algorithm

This section is divided into two parts, which introduce the traditional PF algorithm and the proposed PF algorithm, respectively.

Classical PF algorithm

For non-linear and/or non-Gaussian systems, PF algorithms introduced by Gordon, which is also known as SIR, provide a practical and effective framework to perform filtering and smoothing for dynamic systems.

For state estimation problems, what we are actually concerned about is finding an appropriate estimator of $x_k$ . Based on Bayesian statistics and SMC framework, the PF algorithm can be seen as a method of solving the Bayesian estimation problem. The state space variable $x_k$ can be regard as a random variable with certain distribution $p(x_k)$ . In order to find the distribution $p(x_k)$ , we have to calculate the posterior distribution $p(x_k|y_{1:k})$ , as we mentioned above in the previous section. Applying the Bayes rule, we can obtain

$$p(x_k|y_{1:k})=\frac{p(y_k|x_k)p(x_k|y_{1:k-1})}{p(y_k|y_{1:k-1})}$$

(7)

where $p(y_k|x_k)$ is obtained through the measurement equation (2) and the distribution $p(y_k|y_{1:k-1})$ is thought as a normalizing constant. Then, we have the Chapman–Kolmogorov equation

$$p(x_k|y_{1:k-1})=\int p(x_k|x_{k-1})p(x_{k-1}|y_{1:k-1})\mathrm{d}{x}_{k-1}$$

(8)

The probability density $p(x_k|x_{k-1})$ , also known as the transition density, can be obtained with the system model (1). A prediction-update procedure is formed using equations (7) and (8). This provides us a way to obtain the posterior PDF of the state, recursively.

However, the posterior PDF $p(x_k|y_{1:k})$ is usually hidden and cannot be easily obtained. And therefore, as alternative to the posterior PDF, we draw $N_s$ random samples (particles) $\{x_k^{(i)}{\}}_{i=1}^{N_s}$ with associated weights $\{{\omega }_k^{(i)}{\}}_{i=1}^{N_s}$ . Please note that the superscripts of $x_k$ and ${\omega }_k$ represent the ordinals of particles, not power exponentiation operation. As the number of samples (particles) increases, the estimates become close to the functional description of the PDF. Thus, we get the sample-based posterior PDF approximately

$$\hat{p}(x_k|y_{1:k})=\sum _{i=1}^{N_s}{\omega }_k^{\left(i\right)}\delta (x_k-x_k^{\left(i\right)})$$

(9)

where $\delta (·)$ is the Dirac delta function which is defined as

$$\delta \left(x\right)=\{\begin{array}{cc}+\infty ,\hfill & x=0\hfill \\ \multicolumn{1}{c}{0,}& x\ne 0\hfill \end{array}$$

(10)

The idea of sequential importance sampling (SIS) is the basis of PF algorithm. We choose an importance distribution such that it factors similarly to the posterior density

$$q(x_k|y_{1:k})=q(x_k|x_{1:k},y_{1:k})q(x_{k-1}|y_{1:k-1})$$

(11)

making the weights satisfy

$${\omega }_k^{\left(i\right)}\propto \frac{p(x_k^{\left(i\right)}|y_{1:k})}{q(x_k^{\left(i\right)}|y_{1:k})}$$

(12)

by considering equation (7), we have

$$p(x_k|y_{1:k})=\frac{p(y_k|x_k)p(x_k|x_{k-1})}{p(y_k|y_{1:k-1})}p(x_{k-1}|y_{1:k-1})$$

(13)

and finally

$${\omega }_k\propto \frac{p(y_k|x_k)p(x_k|x_{k-1})}{q(x_k|x_{1:k},y_{1:k})}\frac{p(x_{k-1}|y_{1:k-1})}{q(x_{k-1}|y_{1:k-1})}$$

(14)

$$=\frac{p(y_k|x_k)p(x_k|x_{k-1})}{q(x_k|x_{1:k},y_{1:k})}{\omega }_{k-1}$$

(15)

It has been proven that using the SIS method above alone may lead to the particle degeneracy problem. That is to say, after a few iterations, most of the particle weights reduce close to 0, only one weight increases close to 1. This makes it a waste of computing resources and causes the algorithm become invalid. To avoid this phenomenon, the PF algorithm needs another step, referred to as resampling.

Normally, we use effective sample size $N_{\mathrm{eff}}$ to define the degree of particle degeneracy

$$N_{\mathrm{eff}}=\frac{1}{\sum _{i=1}^{N_s}{\left({\omega }_k^{\left(i\right)}\right)}^2}$$

(16)

During each time step, the posterior distribution is resampled: $N_s$ particles are chosen from the approximated posterior density. With this procedure, it is assured that particles with small or negligible weight will be discarded while particles with more weight will multiply. The pseudocodes for complete PF algorithm are summarized in Algorithm 1.

Algorithm 1. The classical particle filter.
1. $\left[{\{x_k^{\left(i\right)},{\omega }_k^{\left(i\right)}\}}_{i=1}^{N_s}\right]=\mathrm{PF}[{\{x_{k-1}^{\left(i\right)},{\omega }_{k-1}^{\left(i\right)}\}}_{i=1}^{N_s},y_k]$
2. for $i=1:N_s$ do
3.  Draw $x_k^{\left(i\right)}~q(x_k|x_{k-1}^{\left(i\right)},y_k)$ ;
4.  Assign the particles with corresponding weights, ${\omega }_k^{\left(i\right)}\propto {\omega }_{k-1}^{\left(i\right)}\frac{p(y_k|x_k^{\left(i\right)})p(x_k^{\left(i\right)}|x_{k-1}^{\left(i\right)})}{q(x_k^{\left(i\right)}|x_{1:k}^{\left(i\right)},y_{1:k})}$
5. end for
6. Calculate total weight: ${\Omega }_s=\sum _{i=1}^{N_s}{\omega }_k^{\left(i\right)}$
7. for $i=1:N_s$ do
8.  Normalize: ${\omega }_k^{\left(i\right)}={\omega }_k^{\left(i\right)}/{\Omega }_s$
9. end for
10. Calculate effective particle set size: $S_{\mathrm{eff}}=\frac{1}{\sum _{i=1}^{N_s}{\left({\omega }_k^{\left(i\right)}\right)}^2}$
11. If $S_{\mathrm{eff}}<S_{\mathrm{thr}}$ then
12.  Implement the resampling procedure;
13.    ${(x_k^{\left(i\right)},{\omega }_k^{\left(i\right)})}_{i=1}^{N_s}=\mathrm{resample}{(x_k^{\left(i\right)},{\omega }_k^{\left(i\right)})}_{i=1}^{N_s}$ ;
14.    ${\hat{x}}_k=\sum _{i=1}^{N_s}{x}_k^{\left(i\right)}{\omega }_k^{\left(i\right)}$ ;
15. end if

Proposed PF algorithm

The traditional PF algorithm introduced above has averted the particle degeneracy phenomenon by implementing resampling process. When implementing resampling, the idea of traditional PF algorithm is to draw samples (particles) from the PDF $p(x_k|y_{1:k})$ theoretically and set the weights to be equal. The proposed algorithm provides a new approach combining several different correlation coefficient with the resampling process to update particle weights.

Unlike other methods used by other researchers, we introduced several non-linear correlation coefficients into the algorithm. The advantage of this method is that it has good estimation accuracy without increasing the computational complexity of the algorithm (compared with the classical PF). In addition, most of the applicable systems of the algorithm are non-linear; therefore, the use of non-linear correlation coefficients can also avoid the inappropriateness of linear correlation coefficient in some non-linear cases.

Correlation coefficients

Correlation coefficient is defined and used quite frequently in statistics or mathematical analysis. It is a measure that determines the degree of two arrays’ variation tendency. There are several types of correlation coefficients that are perhaps the most widely used: Pearson’s linear correlation coefficient (Pearson’s r), Spearman’s rank-based coefficient (Spearman’s $\rho$ ), and Kendall’s concordance coefficient (Kendall’s $\tau$ ). In addition to the three coefficients above, Xu et al. ¹⁵ introduced a new correlation coefficient named order statistics correlation coefficient, which is based on order statistics and rearrangement inequality. Pearson’s r is a measure of the linear associations between two variables. Spearman’s $\rho$ and Kendall’s $\tau$ are two special cases of rank correlation coefficient which is a measure of the ordinal association between rankings of different ordinal variables, where “ranking” is the assignment of the range in order of numeric value. Although Pearson’s r has fast computation speed, its performance in some non-linear cases is unsatisfactory. On the other hand, neither Spearman’s $\rho$ nor Kendall’s $\tau$ can be implemented as fast as Pearson’s r, but even so, their “ranking” property makes them suitable in non-linear situation. Xu’s order statistics correlation coefficient is a moderate measure that combines the advantages of the coefficients above to achieve a reasonable speed while being well performed in both linear and non-linear situation.

As described above, correlation coefficients provide a measure of the strength and direction of the linear relationship between two time series. We will depict the definitions and structures of several different correlation coefficients below.

Define $(x_i,y_i)$ , $i=1,2,\dots ,N$ to be two time series with length N. Rearranging pairwise the two time series with respect to the magnitudes of $x_i$ , we get two new series denoted by $(x_{(i)},y_{[i]})$ , where $x_{(1)}⩽x_{(2)}⩽\cdots ⩽x_{(N)}$ are called the order statistics of x and $y_{[1]},y_{[2]},\dots ,y_{[N]}$ the associated concomitants. Reversing the roles of x and y, we also define the order statistics of y and the corresponding concomitants which are denoted by $y_{(1)},y_{(2)},\dots ,y_{(N)}$ and $x_{[1]},x_{[2]},\dots ,x_{[N]}$ , respectively. Suppose $x_j$ is at the kth position in the sorted series $x_1,x_2,\dots ,x_N$ , the number $1⩽k⩽N$ is termed the rank of $x_j$ and is denoted by $p_j(=k)$ . Similarly, we can get the rank of $y_j$ denoted by $q_j$ . Such operation of obtaining the ranks of all elements in a series is called ranking. Let $(x_i,y_i)$ and $(x_j,y_j)$ with $i=1,2,\dots ,N$ and $j=i+1,i+2,\dots ,N$ be two data-pairs from the original time series. If $p_j-p_i$ and $q_j-q_i$ have the same sign, we say that the two data-pairs are concordant; otherwise, we say that they are discordant. ¹⁶ Let P stands for the number of concordant pairs and Q the number of discordant pairs, it follows that $P+Q=N(N-1)/2$ . Let $\overline{x}$ , $\overline{y}$ , $\overline{p}$ , and $\overline{q}$ be the arithmetic averages of x, y, p, and q, respectively, Pearson’s correlation coefficient $r_P$ , Spearman’s rho $r_S$ , and Kendall’s tau $r_K$ are defined as follows

$$r_P(x,y)=\frac{\sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\sum _{i=1}^N{(x_i-\overline{x})}^2\sum _{i=1}^N{(y_i-\overline{y})}^2}}$$

(17)

$$r_S(x,y)=1-\frac{6\sum _{i=1}^N{(p_i-q_i)}^2}{N^3-N}$$

(18)

$$r_K(x,y)=\frac{2(P-Q)}{N(N-1)}$$

(19)

As proposed by Xu et al., ¹⁵ the order statistics correlation coefficient can be defined as

$$r_X(x,y)=\frac{\sum _{i=1}^N(x_{\left(i\right)}-x_{(N-i+1)})y_{\left[i\right]}}{\sum _{i=1}^N(x_{\left(i\right)}-x_{(N-i+1)})y_{\left(i\right)}}$$

(20)

Correlation coefficients provide a measure of the strength and direction of the changing correlation relationship between two time series. Pearson’s r computes very fast; however, it may lead to a significant margin of error if non-Gaussian is involved in the model. The other three coefficients can be used under non-Gaussian condition with lower computation efficiency. Xu’s order statistic correlation coefficient provides the solution while offering a wider scope and more efficient application of correlation calculation.

Proposed PF algorithm

Zhou et al. ¹⁷ summarized kinds of methods for choosing the importance density directly or compositely and proposed a new way to choose the particles from importance density, which is based on Pearson correlation coefficient. Unlike the algorithm proposed by Zhou et al., the proposed algorithm in this article provides a new approach combining several different correlation coefficients with the resampling process to update particle weights. We also select other correlation coefficients instead of Pearson correlation coefficient due to its linear characteristic.

SIS method causes the degeneracy phenomenon and provides estimates whose variance increases exponentially with the discrete time index k. Resampling techniques are effective methods to solve this problem. Besides, applying the proposed methods in this article can also refresh the particles weights.

When implementing resampling, the idea of traditional PF algorithm is to draw samples (particles) from the PDF $p(x_k|y_{1:k})$ theoretically and set the weights to be equal. The proposed algorithm provides a new approach combining several different correlation coefficients with the resampling process to update particle weights.

During resampling process, particles chosen from the approximated posterior density according to the effective sample size $N_{\mathrm{eff}}$ have their weights normalized indistinguishably. In other words, for particle sets $\{x_k^{(i)},{\omega }_k^{(i)}{\}}_{i=1}^{N_s}$ , the weights $\{{\omega }_k^{(1)},{\omega }_k^{(2)},\cdots ,{\omega }_k^{(N_s)}\}$ associated with the resampled particles $\{x_k^{(1)},x_k^{(2)},\dots ,x_k^{(N_s)}\}$ are all equal to $1/N_s$ , in most cases. The proposed PF algorithm makes adjustments to the resampling process. It introduces the concept of correlation coefficient to recalculate the weights of resampled particles according to a short time series of observations or measurements $\{y_k{\}}_{k=j-L+1}^j$ with length L. When applying the proposed approach to resampling, particle weights with strong correlation relationship are increased, while particle weights with weak correlation relationship are reduced before normalization process, respectively. Therefore, under the premise of retaining the particle diversity, and with reasonable computational complexity, the proposed method can be significantly more selective and specific about the resampled particles.

In HMM mathematic framework, since the independent observations of the true state $\{y_n{\}}_{n=m,m+1\cdots l}$ and the observations of the particles $\{y_n^{★}{\}}_{n=m,m+1\cdots l}$ vary together ( $1⩽m<l⩽n$ ), the method of correlation coefficient could be applied to check the synchronization of both variable series. Because of the complexity of the computational process of the rank correlation coefficient and normalization, we apply the method in some selected iteration which is dominated by the judgment index $S_{\mathrm{thr}}$ . The implementation of the proposed algorithm is as follows:

• Step 1: Particles propagation. For $N_s$ particles, sample $x_k^{(i)}$ is drawn from the prior density $p(x_k|x_{k-1}^{(i)})$ .

• Step 2: Weights calculation. As illustrated in equation (14), particle weights can be recursively computed. Once we have the particle weights ${\omega }_k^{(i)}$ , we can calculate the normalized weights using the following equation

$${\omega }_k^{(i)}=\frac{{\omega }_k^{(i)}}{\sum _{i=1}^{N_s}{\omega }_k^{(i)}}$$

(21)

• Step 3: Initial resampling process. Compute the effective sample size $N_{\mathrm{eff}}$ according to equation (16), set a suitable value for $S_{\mathrm{thr}}$ , and make comparison between the two parameters to decide whether the algorithm carries out the resampling process.

• Step 4: Prepare two observation (measurement) arrays from the same time series.

Let k denote the discrete time index of the iterations, $X_k=\{x_j:j=k-L+1,\dots ,k\}$ and $X_k^{★}=\{x_j^{★}(i):j=k-L+1,\dots ,k;i=1,\dots ,N_s\}$ denote the true state and a set of all particles in recent time from $k-L+1$ to k, respectively, where L is a constant specified in advance and $N_s$ is the particles number. Accordingly, by calculating with the system observation function (equations (22) and (23)), we obtain $Y_k=\{y_j:j=k-L+1,\dots ,k\}$ and $Y_k^{★}=\{y_j^{★}(i):j=k-L+1,\dots ,k;i=1,\dots ,N_s\}$

$$y_j=h(x_j,0)$$

(22)

$$y_j^{★}\left(i\right)=h(x_j^{★}\left(i\right),0)$$

(23)

• Step 5: Calculate the correlation coefficient.

According to the definition of different rank correlation coefficients above, we get the corresponding equations as follows

$$r_S\triangleq 1-\frac{6\sum _{i=j-L+1}^j{(p_i-q_i)}^2}{L(L^2-1)}$$

(24)

$$r_K\stackrel{\mathrm{\Delta }}{=}\frac{2(P-Q)}{L(L-1)}$$

(25)

$$r_X\triangleq \frac{\sum _{i=j-L+1}^j(y_{\left(i\right)}^{★}-y_{(L-i+1)}^{★})y_{\left[i\right]}}{\sum _{i=j-L+1}^j(y_{\left(i\right)}^{★}-y_{(L-i+1)}^{★})y_{\left(i\right)}}$$

(26)

where $r_S$ , $r_K$ , and $r_X$ denote Spearman’s $\rho$ , Kendall’s $\tau$ , and Xu’s order statistics correlation. Let $(y_i,y_i^{★})$ , $i=j-L+1,\dots ,j$ , be two time series of length L. Rearranging pairwise the two time series with respect to the magnitudes of $y_i$ , we get two new series denoted by $(y_{(i)},y_{[i]}^{★})$ , where $y_{(j-L+1)}⩽\cdots ⩽y_{(j)}$ are called the order statistics of $y_i$ and $y_{[j-L+1]},\dots ,y_{[j]}$ the associated concomitants. ⁷ P and Q denote the number of concordant and discordant pairs, respectively. Here, “Concordant” means that $p_b-p_a$ and $q_b-q_a$ have the same sign (i.e. $\{p_b-p_a>0\cap q_b-q_a>0\}$ ), where $p_i$ is the ordinal number of $y_i$ in the sorted time series $y_{(j-L+1)},\dots ,y_{(j)}$ . Analogously, we get $q_a$ and $q_b$ from time series $y_i^{★}$ .

• Step 6: Adjust the parameter of corresponding weights.

When calculated the rank correlation coefficients of the two time series, we take the exponential function with an adjustable parameter $\alpha$ to expand the numerical range of above results, not limited to $-1⩽\mathrm{cc}⩽1$ . For the purpose of simplifying the calculation, we take a threshold $S_{\mathrm{thr}}$ as reference, redistribute the weights using the proposed PF method for the condition of particle degeneration screened by the threshold. The updated particle weights are proportional to the correlation coefficients of the corresponding time series.

The proposed algorithm proceeds as follows. We can notice that the processes before resampling remain essentially unchanged. Hereafter, we abbreviate Spearman correlation coefficient resampling based particle filtering to SCCRPF, Kendall correlation coefficient resampling based particle filtering to KCCRPF, order statistic correlation coefficient resampling based particle filtering to OCCRPF, and collectively call them XCCRPF.

Algorithm 2. The proposed correlation coefficient PF method.
1. $\left[{\{x_k^{\left(i\right)},{\omega }_k^{\left(i\right)}\}}_{i=1}^{N_s}\right]=\mathrm{XCCRPF}[{\{x_{k-1}^{\left(i\right)},{\omega }_{k-1}^{\left(i\right)}\}}_{i=1}^{N_s},y_k]$
2. for $i=1:N_s$ do
3.  Draw $x_k^{\left(i\right)}~\mathrm{Prior}q(x_k|x_{k-1}^{\left(i\right)},y_k)=p(x_k|x_{k-1}^{\left(i\right)})$ ;
4.  Assign the particle a weight, ${\omega }_k^{\left(i\right)}\propto {\omega }_{k-1}^{\left(i\right)}\frac{p(y_k|x_k^{\left(i\right)})p(x_k^{\left(i\right)}|x_{k-1}^{\left(i\right)})}{q(x_k^{\left(i\right)}|x_{k-1}^{\left(i\right)},y_k)}$
5. end for
6. Calculate total weight: ${\Omega }_s=\sum _{i=1}^{N_s}{\omega }_k^{\left(i\right)}$
7. for $i=1:N_s$ do
8.  Normalize: ${\omega }_k^{\left(i\right)}={\omega }_k^{\left(i\right)}/{\Omega }_s$
9. end for
10. Calculate effective particle set size: $S_{\mathrm{eff}}=\frac{1}{\sum _{i=1}^{N_s}{\left({\omega }_k^{\left(i\right)}\right)}^2}$
11. if $S_{\mathrm{eff}}<S_{\mathrm{thr}}$ then
12.  Calculate corresponding correlation coefficient between theoretical and practical system observations $c{c}_k$ ;
13.   $\mathrm{cc}{e}_k=e^{\alpha \times c{c}_k}$ ;
14.   ${\omega }_k^{\left(i\right)}=\mathrm{cc}{e}_k·p(y_k|x_k\left(i\right))$ ;
15.  Normalize: ${\omega }_k^{\left(i\right)}={\omega }_k^{\left(i\right)}/\sum _{i=1}^{N_s}{\omega }_k^{\left(i\right)}$ ;
16. end if

Justification

Smith and Gelfand ¹⁸ proved that if samples have upper boundary, the cumulative distribution function (CDF) of corresponding samples also has upper boundary under the customary bootstrap (particle filters). In this article, the screening procedure affected by the threshold $S_{\mathrm{thr}}$ behaves similarly to the resampling procedure essentially.

Suppose that samples $\{x_k^{★(i)};i=1,2,\dots ,N_s\}$ are available from a continuous density function $G(x)$ and that samples are required from the PDF proportional to $L(x)G(x)$ , where $L(x)$ is a known function. According to the proof process given by Smith and Gelfand, the sample drawn from the discrete distribution $\{x_k^{★(i)};i=1,2,\dots ,N_s\}$ with probability mass $L(x_k^{★(i)})/\sum _{j=1}^{N_s}L(x_k^{★(j)})$ on $x_k^{★(i)}$ , converges to the distribution of required density as $N_s$ tends to infinity.

Density function $G(x)$ is identified with the prior $p(x_k|y_{1:k-1})$ . In the proposed algorithm, $L(x)$ is identified with $\mathrm{cc}{e}_k·p(y_k|x_k)$ (the likelihood) where $\mathrm{cc}{e}_k$ is the above exponential transform of the corresponding correlation coefficient. With the normalized particle weights (equation (18)), samples (particles) drawn from the discrete distribution could also tend in distribution to the required density as $N_s$ tends to infinity.

Simulation performance and harmonic parameter estimation

In this section, firstly we apply a one-dimensional model that has been analyzed in many publications before. Then we apply the proposed method to a given seven-dimensional state space model. The effectiveness of the proposed method is verified by comparing the parameter estimation results of different research methods.

One-dimensional model

The following example of simulation is well known and has been analyzed in many publications before. It is noteworthy that we set the process noise as a Gaussian mixture heavy-tailed distribution instead of a common Gaussian distribution in order to test the effectiveness and reliability of the proposed method more rigorously

$$x_k=\frac{x_{k-1}}{2}+\frac{25x_{k-1}}{1+x_{k-1}^2}+8\cos (1.2(k-1))+v_{k-1}$$

(27)

$$y_k=\frac{x_k^2}{20}+w_k$$

(28)

where $n_k~N(0,1)$ is zero mean Gaussian random variables with standard deviation (SD) of 1. The initial state $x_0=0.1$ , particles number $N_s=50$ , and the adjustable parameter $\alpha =10$ . The simulation step $N=100$ . The process noise distribution $v_k$ is a Gaussian mixture given by

$$p\left(v_k\right)=ϵN(v;0,{\sigma }_{v1}^2)+(1-ϵ)N(v;0,{\sigma }_{v2}^2)$$

(29)

where $ϵ=0.8$ , ${\sigma }_{v1}^2=1$ , and ${\sigma }_{v2}^2=10$ , that can model the heavy-tailed distributions well.

The simulation algorithm was coded using MATLAB, and the program was performed on an Inter Core i5-4570 PC of 3.2 GHz and 4 GB memory. In order to assess the performance of different algorithms, we can use the root-mean-square error (RMSE) as criterion. The RMSE is given by

$$\mathrm{RMSE}=\sqrt{\frac{\sum _{k=1}^N{({\hat{x}}_k-x_k)}^2}{N}}$$

(30)

where ${\hat{x}}_k$ is the estimation of true state variable value $x_k$ , N is the simulation step size, and here, $N=100$ . We repeated the state estimation simulations for 1000 times in MATLAB, in order to get more accurate simulation results. Average results of 1000 random simulation runs are shown in Table 1. AvgT is the average time required by corresponding filtering algorithm for processing 1000 iterations. GSPF and GMSPPF are introduced as contrast to have a comprehensive evaluation of performance of the proposed method.

Table 1.

Simulation results of different algorithms.

Algorithm	RMSE mean	RMSE variance	AvgT (ms)
PF	9.400	0.722	18.752
GSPF	7.114	3.242	73.068
GMSPPF	9.687	4.669	299.248
SCCRPF	5.534	2.126	25.197
KCCRPF	5.511	1.486	237.779
OCCRPF	4.745	0.468	25.210

RMSE: root-mean-square error; PF: particle filtering; GSPF: Gaussian sum particle filter; GMSPPF: Gaussian mixture sigma-point particle filter; SCCRPF: Spearman correlation coefficient resampling based particle filtering; KCCRPF: Kendall correlation coefficient resampling based particle filtering; OCCRPF: order statistic correlation coefficient resampling based particle filtering.

Figure 1 shows the state estimation results of different algorithms in one simulation experiment. Different algorithms are marked in several different colors for the sake of convenience as illustrated in the legend. RMSE values of 10 random simulation runs of different algorithms are plotted in Figure 2. It is shown that the performance of the proposed algorithms is better than the PF, GSPF, and GMSPPF algorithms, while the three proposed algorithms perform similarly. In Table 1, KCCRPF’s extra time consumption is due to the computational complexity of Kendall’s $\tau$ . Even though the time consumption of Kendall’s $\tau$ was higher, the computation accuracy was better. As long as the time consumption is reduced, KCCRPF is still available. With similar simulation time with GMSPPF, the simulation accuracy of KCCRPF is much higher than that of GMSPPF. Furthermore, we are also expected to shorten the calculation time of Kendall’s $\tau$ in the following work. According to the results from this simulation, we may find that SCCRPF and OCCRPF are more accurate than other algorithms, in condition of consuming similar computation time, and OCCRPF has lower RMSE variance.

Figure 1.

State estimation results of different algorithms.

Figure 2.

RMSE of different algorithms.

Seven-dimensional model

As the above-mentioned model established in section “Harmonic parameter model,” we consider establishing a concrete seven-dimensional model to estimate the harmonic parameters in power system. The simulation signal is set as follows

$$\begin{array}{c}\multicolumn{1}{c}{y_k=A_{k,1}\sin ({\omega }_{k}k\mathrm{\Delta }t+{\phi }_{k,1})}\\ \multicolumn{1}{c}{+A_{k,3}\sin (3{\omega }_{k}k\mathrm{\Delta }t+{\phi }_{k,3})}\\ \multicolumn{1}{c}{+A_{k,5}\sin (5{\omega }_{k}k\mathrm{\Delta }t+{\phi }_{k,5})}\end{array}$$

(31)

The state equation and observation equation in matrix are as follows

$$x_{k+1}=\left[\begin{array}{ccccccc}1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 0\hfill \\ \multicolumn{1}{c}{}& 1\hfill & \hfill & \hfill & \hfill & \hfill & 1/f_s\hfill \\ \multicolumn{1}{c}{}& \hfill & 1\hfill & \hfill & \hfill & \hfill & 0\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & 1\hfill & \hfill & \hfill & 3/f_s\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & 1\hfill & \hfill & 0\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & 1\hfill & 5/f_s\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & \hfill & 1\hfill \end{array}\right]x_k+v_k$$

(32)

$$\begin{array}{c}{y}_{k+1}=[\sin \left(x_{k,2}\right),0,\sin \left(x_{k,4}\right),0,\hfill \\ \sin \left(x_{k,6}\right),0,0]x_k+w_k\hfill \end{array}$$

(33)

where $[A_{k,1},A_{k,3},A_{k,5}]=[1,0.3,0.2]$ , ${\omega }_k=2\pi f$ , and $f=50$ is the utility frequency, $\mathrm{\Delta }t=1/f_s$ and $f_s=3200$ is the sampling frequency, $[{\phi }_{k,1},{\phi }_{k,3},{\phi }_{k,5}]=[\pi /3,\pi /4,\pi /6]$ .

The initial state was set to

$$\begin{array}{cc}\multicolumn{1}{c}{x_1=}& [1,2\pi f/f_s+\pi /3,0.3,3\times 2\pi f/f_s\hfill \\ \multicolumn{1}{c}{}& +\pi /4,0.2,5\times 2\pi f/f_s+\pi /6,2\pi f\hfill \end{array}]$$

(34)

In practical harmonic control, odd harmonics appear more frequently than even harmonics. In a balanced three-phase system, even harmonics have been eliminated due to the symmetrical relationship, and only odd harmonics exist. Because of this, we only choose odd harmonics in the simulation. The process noise is $v_k~N(0,0.1)$ , and the measurement noise $v_k$ should be calculated for an analog of the signal-to-noise ratio (SNR) to be 40 dB. SNR is defined as follows

$$\mathrm{SNR}\left(\mathrm{dB}\right)=20{\log }_{10}\frac{A_s}{A_n}$$

(35)

where $A_s$ and $A_n$ denote the amplitude of signal and noise, respectively. For different kinds of state variables (amplitude and phase), we set similar process noise parameters. This is more conducive to our observation of the impact of non-linearity on the estimation process and results.

We introduced the EKF algorithm as a contrast to show the superiority of the proposed algorithm. The EKF applies the Kalman filter to non-linear systems by simply applying the first-order approximation of Taylor series of the non-linear models. We wanted to do this because the relative simplicity and practical efficacy had made the EKF one of the most widely used algorithms in multiple fields.

The simulation runs for 200 steps, and the results are shown in Figure 3. Observations $o_1$ , $o_3$ , and $o_5$ denote the amplitude of three different frequencies, which should be 1, 0.3, and 0.2 theoretically. $o_2$ , $o_4$ , and $o_6$ denote the estimated phase angles, which have been translated into degree and 60, 45, and 30 are the theoretical value. $o_7$ is the utility frequency of $50\mathrm{Hz}$ . The last figure in Figure 3 is the simulated power signal given by equation (31).

Figure 3.

State estimation results of harmonic model.

As can be seen from Figure 3, the EKF algorithm tends to be stable after a period of oscillation, but its estimation value is biased and has a large error. Especially, for the estimation of phase angle, the deviation is large. The PF and OCCRPF algorithms perform better than the EKF algorithm. The estimation results of the amplitudes by PF and OCCRPF algorithms are similar to those by the EKF algorithm, while the fluctuation is smaller in the estimation process. For the estimation of phase angles, the PF and OCCRPF algorithms are obviously better than EKF. For the estimation of the utility frequency, the difference between the three algorithms is not significant, and the stability estimation values of the three algorithms are within fractions of the actual value.

Talking to the estimation effect of one single algorithm, EKF is the worst, especially for the estimation of the phase angle. The reason can be seen from the state transition matrix equation (5) and the measurement matrix equation (6). The estimations of odd terms of the state vector maintain linear variation, while the estimations of even terms of the state vector are originally from the function of sinusoid, which exhibit its non-linear characteristics. It is also well known that EKF can only give biased estimation for non-linear systems. This leads to errors of the estimation results of EKF. There is another reason why the accuracy of performance of the EKF algorithm was poor in the simulation. The issue is the parameter setting of process noise. We set the same noise parameters for the amplitude and phase variables. The phase vector has a non-linear calculation process of trigonometric function, which enlarges the deviation between the estimated and true value, thus further enlarging the phase error. As mentioned previously, there are no such problems for PF and OCCRPF algorithms.

Conclusion

In this article, a new particle filter algorithm combing with different rank correlation coefficients is proposed. Different rank correlation coefficients are introduced to measure the variation tendency between the particles and the state, and therefore obtain the particle weights. Compared with the PF, GSPF, and GMSPPF, the proposed method has better performance of estimation accuracy, and the computation time consumes about the same. Moreover, the computational complexity of the proposed algorithm is lower than GSPF and GMSPPF, which can be reflected by time consumption in Gaussian mixture noise.

A practical harmonic model of power signal is also simulated in this article. The simulation compared the EKF, PF, and OCCRPF algorithms. Because of non-linearity, EKF made biased estimates of state variables, especially for the phase variables. When it comes to PF and OCCRPF, the estimation of phase variables is much better than PF. Error analysis and simulation results show that PF and OCCRPF are superior to EKF for non-linear systems.

The algorithm proposed in this article leads to at least two interesting directions for continued research. First, because the computation accuracy of KCCRPF is well but causes high time consumption, new approach will be analyzed to reduce the time consumption of calculating Kendall’s $\tau$ . Second, as the algorithm proposed is a feasible tool for tracking variables in non-linear and non-Gaussian systems, we might consider using it to estimate other state parameters in different power quality disturbance signals.