The Maximum Correntropy Criterion-Based Identification for Fractional-Order Systems under Stable Distribution Noises

Abstract

This paper studies the identification for fractional-order systems (FOSs) under stable distribution noises. First, the generalized operational matrix of block pulse functions is used to convert the identified system into an algebraic one. Then, the conventional least mean square (LMS) criterion is replaced by the maximum correntropy criterion (MCC) to restrain the effect of noises, and a MCC-based algorithm is designed to perform the identification. To verify the superiority of the proposed method, the identification accuracy is examined when the noise follows different types of stable distributions. In addition, the impact of parameters of stable distribution on identification accuracy is discussed. It is shown that when the impulse of noise increases, the identification error becomes larger, but the proposed algorithm is always superior to its LMS counterpart. Moreover, the location parameter of stable distribution noise has a significant impact on the identification accuracy.

1. Introduction

Fractional-order calculus (FOC), which originated more than 300 years ago, had remained for a long time as an abstract mathematical concept for the lack of acceptable geometric or physical explanations. It was not until the 1970s that FOC was introduced and applied in fields such as electrochemistry, rheology, energy consumption prediction, and epidemiology [1,2,3,4].

Fractional-order systems or models can be regarded as a generalization of integer-order systems or models based on the concept of FOC. Being nonlocal and historically dependent, FOSs can describe the anomalous behavior of dynamical systems more accurately than integer-order systems. The research related to FOSs can be mainly divided into three categories: control, synchronization, and identification [5,6,7,8]. Identification, in which the measured input and output data are used to establish an appropriate mathematical model under certain optimization criteria, provides detailed information for the FOS before control and synchronization. Methods for FOS identification can be classified into frequency domain methods and time domain methods. A number of frequency domain methods have been developed for FOS identification. Malti R et al. combined interval constraint techniques and constraint propagation techniques to identify FOSs based on frequency domain uncertain but bounded data in [9]. Valerio et al. expended the Levy method to identify non-commensurable FOSs from the frequency response in [10]. After analyzing the frequency characteristics of the fractional-order differentiator, Lin et al. proposed an unbiased parameter estimation method for FOSs in a noisy environment using the auxiliary variable [11]. In [12], the recursive least square identification algorithm was extended to identify the commensurate FOS in the frequency domain. In time domain, the idea of most identification methods for FOSs are learnt from existing integer-order identification methods. Equation-error- and output-error-based methods are the two main methods. In [13], a refined instrumental variable method for continuous-time systems was extended to identify FOS. In [14], the recursive least square with state-variable filters and the prediction-error method are proposed to identify continuous-time FOSs. In [15], a metaheuristic optimization-based method is used to identify FOSs. In addition, the orthogonal basis-based method was widely used in the identification of FOSs. Its main idea is to use orthogonal-basis functions to approximate the fractional-order integral operation and then obtain the corresponding operational matrix. With operational matrices, the underlying system can be transformed into an algebraic one, which simplifies the identification process [16,17,18].

Notably, all the above-mentioned methods adopt mean square error (MSE) as the objective function because of its low computational complexity [19,20]. Based on MSE, the LMS criterion has been developed for various applications, which performs well only when the measurement is without or with Gaussian white noise. In practice, however, the Gaussian assumption does not always hold [21]. Noises with more peak or impulse, which follow the stable distribution, are usually encountered during the data measurement process. This will cause an LMS-based algorithm to lose its performance in the identification process. Therefore, it is necessary to introduce a new identification criterion for FOS identification when the data are contaminated by stable distribution noises. Ref. [21] adopted an approximate absolute error function as the cost function and proposed a stochastic gradient method. When the independent variable fluctuates sharply, the differentiation of this cost function approximates to $\pm 1$, which can restrain the influence of impulse noise. Moreover, correntropy has better robustness and has been used in signal processing and machine learning [22,23]. Over the past few years, a series of adaptive filters based on the MCC have been developed for constructing robust adaptive filters [24,25]. In addition, lots of efforts have been reported to improve the performance of MCC-based algorithms [26,27,28]. To the best of our knowledge, research on the application of the MCC in FOS identification has been very limited.

Motivated by the above discussion, this paper aims to design a robust method for identifying FOSs when the output signal is disturbed by stable distribution noises. The generalized operational matrix of block pulse functions (BPFs) is adopted to discretize the FOS to solve the difficulty involved in directly finding the FOC of signals. Then, correntropy is chosen as the objective function over the conventional MSE to restrain the effect of noise on the identification process; correspondingly, an MCC-based stochastic gradient ascent (MCC-SGA) method is designed to identify FOSs. To verify the superiority of the proposed algorithm, an identification algorithm based on MSE is introduced, and the identification accuracy is examined when noise follows different types of stable distribution. Additionally, the effect of parameters in the identification algorithm and stable distribution on identification accuracy are discussed. The contributions of this paper are threefold: (1) the identification problem of FOSs under stable distribution noises are studied; (2) it is verified both theoretically and experimentally that the MCC-SGA algorithm can deliver more excellent performance than the LMS-SGD algorithm; (3) the impact of parameters in the identification algorithm and the distribution of the noise on identification accuracy are analyzed.

The remainder of the paper is organized as follows. In Section 2, a brief mathematical background of FOC, the BPF-based operational matrix, and the stable distribution noise are given. The identification of FOS based on the MCC is proposed in Section 3, and the impact of kernel width on estimation accuracy is discussed in detail. In Section 4, simulation examples are given. Finally, the conclusion remarks are presented in Section 5.

2. Preliminaries

2.1. Definition of Fractional-Order Integral and Derivative

FOC is the generalization of classical integration and differentiation to the non-integer order. There are several definitions for FOC. Due to the unity of definitions, Riemann–Liouville (R-L) definitions [29] are used in this paper, which are given as follows.

Definition 1.

The R-L fractional integral operator of order ν is defined as

$$\begin{array}{c}{I}^{\nu }f\left(t\right)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\Gamma \left(\nu \right)}{\int }_0^t{(t-s)}^{\nu -1}f\left(s\right)\mathrm{d}s=\frac{1}{\Gamma \left(\nu \right)}{t}^{\nu -1}\ast f\left(t\right),}\hfill & \nu >0,\\ f\left(t\right),\hfill & \nu =0,\end{array}\right.\hfill \end{array}$$

(1)

where $(0,t)$ is the integral interval, $\Gamma (·)$ denotes the Gamma function, and ∗ is the convolution operator.

The R-L fractional-order integral is a linear operator, namely

$$I^{\nu }(\lambda f\left(t\right)+\mu g\left(t\right))=\lambda I^{\nu }f\left(t\right)+\mu I^{\nu }g\left(t\right),$$

(2)

where $\lambda$ and $\mu$ are constants.

Definition 2.

The R-L fractional-order derivative is defined as

$$\begin{array}{cc}\hfill D^{\nu }f\left(t\right)& =D^n\left[I^{n-\nu }f\left(t\right)\right]\hfill \\ \hfill & =\frac{1}{\Gamma (n-\nu )}{\left(\frac{\mathrm{d}}{\mathrm{d}t}\right)}^n{\int }_0^t\frac{f\left(s\right)}{{(t-s)}^{\nu +1-n}}\mathrm{d}s,\hfill \end{array}$$

(3)

where ν is a positive number and n is the smallest integer greater than ν.

2.2. Generalized Operational Matrices of Block Pulse Functions

A set of BPFs in the semi-open interval $[0,T)$ can be defined as follows:

$${\psi }_i\left(t\right)=\left\{\begin{array}{cc}1,\hfill & ih\le t<(i+1)h,\\ 0,\hfill & \mathrm{otherwise},\end{array}\right.i=0,\dots ,M-1,$$

(4)

where $h=\frac{T}{M}$ is the sampling period and M is the number of BPFs presented in the set.

According to [30], an arbitrary absolutely integrable function defined on the interval $[0,T)$ can be represented by a linear combination of the BPFs as

$$f\left(t\right)\simeq \sum _{i=0}^{M-1}{f}_i{\psi }_i\left(t\right)={\mathbf{f}}^{\mathtt{T}}{\Psi }_M\left(t\right),$$

(5)

where the superscript $\mathtt{T}$ denotes transposition, ${\Psi }_M\left(t\right)={[{\psi }_0\left(t\right),{\psi }_1\left(t\right),\cdots ,{\psi }_{M-1}\left(t\right)]}^{\mathtt{T}}$ is the block pulse vector of order M, and $\mathbf{f}={[f_0,f_1,\dots ,f_{M-1}]}^{\mathtt{T}}$ is the coefficient vector defined as

$$f_i=\frac{1}{h}{\int }_{ih}^{(i+1)h}f\left(t\right)\mathrm{d}t\approx f\left(ih\right).$$

(6)

Remark 1.

The approximation errors arisen from Equations (5) and (6) approach zero, with the value of M tending to be infinite [31].

Based on Equation (1), the R-L fractional-order integration of the block pulse vector can be obtained as

$$I^{\nu }{\Psi }_M\left(t\right)=\frac{1}{\Gamma \left(\nu \right)}{t}^{\nu -1}·{\Psi }_M\left(t\right),$$

(7)

which can be written in a matrix form, i.e.,

$$I^{\nu }{\Psi }_M\left(t\right)=P^{\nu }{\Psi }_M\left(t\right).$$

(8)

Here $P^{\nu }$ is the operational matrix of R-L integration, and given as

$$P^{\nu }=\frac{h^{\nu }}{\Gamma (\nu +2)}\left(\begin{array}{cccc}{\xi }_1^{\nu }& {\xi }_2^{\nu }& \cdots & {\xi }_M^{\nu }\\ 0& {\xi }_1^{\nu }& \cdots & {\xi }_{M-1}^{\nu }\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\xi }_1^{\nu }\end{array}\right),$$

(9)

${\xi }_1^{\nu }=1,{\xi }_l^{\nu }=l^{\nu +1}-2{(l-1)}^{\nu +1}+{(l-2)}^{\nu +1},l=2,\dots ,M$. Combining with Equation (5) and using the operational matrix, the fractional integration of function $f\left(t\right)$ can be written as

$$I^{\nu }f\left(t\right)={\mathbf{f}}^{\mathtt{T}}{P}^{\nu }{\Psi }_M\left(t\right).$$

(10)

After some matrix operations, Equation (10) can be rewritten as

$$I^{\nu }f\left(t\right)=\frac{h^{\nu }}{\Gamma (\nu +2)}\sum _{i=0}^{M-1}{q}_i{\psi }_i\left(t\right),$$

(11)

where

$$q_i=\sum _{p=0}^i{f}_p{\xi }_{i-p+1}^{\nu },$$

and $f_p$ can be calculated by using Equation (6).

In this way, the fractional integration of a function can be converted into an algebraic operation, and the complexity of directly calculating the fractional integral of the input and output signals is avoided.

Remark 2.

Most of the currently reported methods utilize the definition of the Grünwald–Letnikov (G-L) derivative to calculate FOC, which poses a huge computational burden and amplifies the effect of noise [32]. In view of this, here, the operational matrices of BPFs are adopted to discretize the FOS.

2.3. Stable Distribution for Impulse Noise

In order to model the impulse noise in real environment, stable distribution is studied in this section. Stable distribution is a four-parameter family of distributions and is usually denoted by $S(\alpha ,\beta ,\gamma ,\mu )$. The parameter $\alpha$ is the characteristic exponent satisfying $0<\alpha \le 2$, which controls the thickness of tails in the distribution. $\mu \in (-\infty ,+\infty )$ is the location parameter, which represents the mean value of the distribution. $\beta \in [-1,1]$ is the symmetry parameter, which specifies the skewness of the distribution. The parameter $\gamma >0$ is the dispersion parameter, which acts as the variance and determines the spread of the distribution around its location parameter. The family of stable distribution includes three special cases:

• Gaussian distribution: When $\alpha =2$, the stable distribution turns to be a Gaussian distribution, and the skewness $\beta$ loses its effect. The probability density function of Gaussian distribution is written as:

  $$f\left(x\right)=\frac{1}{\pi \gamma [1+{\left(\frac{x-\mu }{\gamma }\right)}^2]}.$$

  (12)
• Cauchy distribution: When $\alpha =1$ and $\beta =0$, the stable distribution turns to be a Cauchy distribution. The probability density function of Cauchy distribution is written as:

  $$f\left(x\right)=\frac{1}{\sqrt{8{\gamma }^2\pi }{e}^{-\frac{{(x-\mu )}^2}{8{\gamma }^2}}}.$$

  (13)
• Lévy distribution: When $\alpha =\frac{1}{2}$ and $\beta =1$, the stable distribution reduces to be a Lévy distribution, and the probability density function of Lévy distribution is written as:

  $$f\left(x\right)=\sqrt{\frac{\gamma }{2\pi }}\frac{e^{-\frac{\gamma }{2(x-\mu )}}}{{(x-\mu )}^{3/2}}.$$

  (14)

Except for the aforementioned three cases, the probability density function of stable distribution $S(\alpha ,\beta ,\gamma ,\mu )$ cannot be expressed in a closed form, but its characteristic function is given by

$$\varphi (t;\alpha ,\beta ,\gamma ,\mu )=\left\{\begin{array}{cc}exp\left[it\mu -{\mid \gamma t\mid }^{\alpha }\left(1-i\beta \mathrm{sgn}\left(t\right)tan\frac{\pi \alpha }{2}\right)\right],\hfill & \alpha \ne 1,\\ exp\left[it\mu -\mid \gamma t\mid \left(1+i\beta \frac{2}{\pi }\mathrm{sgn}\left(t\right)ln\mid t\mid \right)\right],& \alpha =1,\end{array}\right.$$

(15)

where

$$\mathrm{sgn}\left(t\right)=\left\{\begin{array}{cc}1,\hfill & t>0,\\ 0,\hfill & t=0,\\ -1,\hfill & t<0.\end{array}\right.$$

Based on the characteristic function, the probability density function of stable distribution $S(\alpha ,\beta ,\gamma ,\mu )$ can be given as

$$p\left(x\right)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{-{\left|\omega \right|}^{\alpha }}{e}^{-i\omega x}\mathtt{d}\omega ,$$

(16)

and the probability density functions of the stable distribution with different parameter values are plotted in Figure 1. Figure 1a shows that as $\alpha$ decreases, the tail of the probability density function becomes heavier, which means that the probability of noises with large magnitude occurring increases. Figure 1b demonstrates that when $\beta >0$, the curve is skewed to the right and the right tail is thicker, resulting in most of the noise being positive; when $\beta <0$, the curve is skewed to the left and the left tail is thicker, which results in the noise being mostly negative; when $\beta =0$, the distribution is symmetric around the location parameter. Figure 1c shows that as $\gamma$ increases, the distribution becomes more spread, which will generate noise with high fluctuations. Figure 1d exhibits that the location of the distribution is decided by the value of $\mu$. As $\left|\mu \right|$ increases, the mean of the noise deviates farther from 0.

3. Parameter Identification of FOS Based on MCC

3.1. Problem Formulation

Consider the linear FOS described by the following differential equation:

$$\sum _{i=0}^n{a}_i{D}^{{\beta }_i}x\left(t\right)=bu\left(t\right),$$

(17)

where $u\left(t\right)$ and $x\left(t\right)$ are the input and output signals. The derivative orders ${\beta }_i,(i=0,\dots ,n)$ are known, sorted in descending order. $a_i$ and b are the unknown parameters to be identified. Under zero initial condition, applying the fractional integral of order ${\beta }_n$ to both sides of Equation (17) yields

$$\sum _{i=0}^n{a}_i{I}^{{\beta }_n-{\beta }_i}x\left(t\right)=b{I}^{{\beta }_n}u\left(t\right).$$

(18)

According to Equation (5), $u\left(t\right)$ and $x\left(t\right)$ can be expanded onto the BPFs as

$$u\left(t\right)={\mathbf{U}}^{\mathtt{T}}{\Psi }_M\left(t\right),x\left(t\right)={\mathbf{X}}^{\mathtt{T}}{\Psi }_M\left(t\right),$$

(19)

where $\mathbf{U}={[u_0,u_1,\dots ,u_{M-1}]}^{\mathtt{T}}$ and $\mathbf{X}={[x_0,x_1,\dots ,x_{M-1}]}^{\mathtt{T}}$. Combined with Equation (10), we have

$$\sum _{i=0}^n{a}_i{\mathbf{X}}^{\mathtt{T}}{P}^{{\beta }_n-{\beta }_i}{\Psi }_M\left(t\right)=b{\mathbf{U}}^{\mathtt{T}}{P}^{{\beta }_n}{\Psi }_M\left(t\right).$$

(20)

For the ease of operation, the data are collected at sampling time $t=(k-1)h,k=1,\cdots ,M$, the values of variables ${\mathbf{X}}^{\mathtt{T}}{P}^{{\gamma }_n-{\gamma }_i}{\Psi }_M{\left(t\right)|}_{t=(k-1)h}$ and ${\mathbf{U}}^{\mathtt{T}}{P}^{{\beta }_n}{\Psi }_M{\left(t\right)|}_{t=(k-1)h}$ are defined as $X_i\left(k\right)$ and $U_n\left(k\right)$, respectively. Taking use of Equation (11), $X_i\left(k\right)$ and $U_n\left(k\right)$ can be given as

$$X_i\left(k\right)=\frac{h^{{\beta }_n-{\beta }_i}}{\Gamma ({\beta }_n-{\beta }_i+2)}\sum _{p=1}^k{x}_p{\xi }_{k-p+1}^{{\beta }_n-{\beta }_i}$$

(21)

and

$$U_n\left(k\right)=\frac{h^{{\beta }_n}}{\Gamma ({\beta }_n+2)}\sum _{p=1}^k{u}_p{\xi }_{k-p+1}^{{\beta }_n}.$$

(22)

Then, Equation (18) can be discretized as follows:

$$\sum _{i=0}^n{a}_i{X}_i\left(k\right)=b{U}_n\left(k\right).$$

(23)

Without loss of generality, we let $a_n=1$; then, Equation (23) can be rewritten as

$$X_n\left(k\right)={\varphi }^{\mathtt{T}}\left(k\right)\theta ,$$

(24)

where

$$\varphi \left(k\right)={[-X_0\left(k\right),\dots ,-X_{n-1}\left(k\right),U_n\left(k\right)]}^{\mathtt{T}}$$

and

$$\begin{gathered} \theta ={[a_0,\dots ,a_{n-1},b]}^{\mathtt{T}}, \\ k=1,\dots ,M. \end{gathered}$$

In practice, the output $x\left(t\right)$ is unobtainable, and the measured output $y\left(t\right)$ can be presented as

$$y\left(t\right)=x\left(t\right)+v\left(t\right),$$

(25)

where $v\left(t\right)$ denotes the measurement noise. Then, Equation (18) can be converted to

$$\sum _{i=0}^n{a}_i{I}^{{\beta }_n-{\beta }_i}y\left(t\right)=b{I}^{{\beta }_n}u\left(t\right)+\sum _{i=0}^n{a}_i{I}^{{\beta }_n-{\beta }_i}v\left(t\right).$$

(26)

Similar to Equation (19), $y\left(t\right)$ and $v\left(t\right)$ can be written as

$$y\left(t\right)=Y^{\mathtt{T}}{\Psi }_M\left(t\right),v\left(t\right)=V^{\mathtt{T}}{\Psi }_M\left(t\right),$$

(27)

where $Y={[y_0,y_1,\dots ,y_{M-1}]}^{\mathtt{T}}$ and $V={[v_0,v_1,\dots ,v_{M-1}]}^{\mathtt{T}}$. Combined with Equation (10), we have

$$\sum _{i=0}^n{a}_i{Y}^{\mathtt{T}}{P}^{{\beta }_n-{\beta }_i}{\Psi }_M\left(t\right)=b{U}^{\mathtt{T}}{P}^{{\beta }_n}{\Psi }_M\left(t\right)+\sum _{i=0}^n{a}_i{V}^{\mathtt{T}}{P}^{{\beta }_n-{\beta }_i}{\Psi }_M\left(t\right).$$

(28)

We represent the values of variables $Y^{\mathtt{T}}{P}^{{\beta }_n-{\beta }_i}{\Psi }_M{\left(t\right)|}_{t=(k-1)h}$, $V^{\mathtt{T}}{P}^{{\beta }_n-{\beta }_i}{\Psi }_M{\left(t\right)|}_{t=(k-1)h}$ as $Y_i\left(k\right)$ and $V_i\left(k\right)$, which can be obtained by the way similar to Equations (21) and (22). Subsequently, Equation (28) can be expressed in the following discrete form:

$$Y_n\left(k\right)={\varphi }_c^{\mathtt{T}}\left(k\right)\theta +v_f\left(k\right),$$

(29)

where

$${\varphi }_c\left(k\right)={[-Y_0\left(k\right),\dots ,-Y_{n-1}\left(k\right),U_n\left(k\right)]}^{\mathtt{T}}$$

(30)

and

$$v_f\left(k\right)=\sum _{i=0}^n{a}_i{V}_i\left(k\right).$$

(31)

The goal of identifying FOS (17) is to find a desired $\widehat{\theta }$ that minimizes the objective function

$$J\left(\widehat{\theta }\right)=J\left[e\left(k\right)\right]=J[Y_n\left(k\right)-{\varphi }_c^{\mathtt{T}}\left(k\right)\widehat{\theta }].$$

(32)

Presently, the objective function $J(·)$ is usually chosen as MSE, which cannot achieve satisfactory performance in non-Gaussian situations, especially when the data are disturbed by impulse noises. The main reason is that the second-order statistics contained in MSE will enlarge the contribution of noises with strong impulse. Therefore, searching for an alternative objective function against large outliers or impulse noises is the main research of this paper.

3.2. Correntropy

Correntropy is a nonlinear measure of the similarity between two random variables X and Y, which is defined as [23]

$$V(X,Y)=\mathrm{E}\left[\kappa (X,Y)\right],$$

(33)

where $\kappa (·)$ is a positive definite kernel, and $\mathrm{E}[·]$ is the expectation operator. A widely used kernel in correntropy is the Gaussian kernel

$$\kappa (x,y)=G_{\sigma }\left(e\right)=exp\left(-\frac{e^2}{2{\sigma }^2}\right),$$

(34)

where $e=x-y$ and $\sigma >0$ denotes the kernel width. Gaussian kernel functions with different kernel widths and their differentiations are shown in Figure 2.

As we can see, correntropy with the Gaussian kernel is a measure of how similar two variables are within a small neighborhood determined by the kernel width $\sigma$. The differentiation of correntropy contains an error term related to noise, which approximates to 0. The smaller the value of $\sigma$, the wider the neighbourhood with a differentiation value of 0, resulting in less sensitivity of the MCC-based algorithm to noise.

3.3. Parameter Identification Based on Maximum Correntropy Criterion

Based on the above discussion, the objective function in Equation (32) is chosen as the correntropy with Gaussian kernel, which can be given as follows:

$$J\left(\widehat{\theta }\right)=\mathrm{E}\left[exp\left(-\frac{e^2\left(k\right)}{2{\sigma }^2}\right)\right].$$

(35)

To achieve the goal of parameter identification, we maximize the instantaneous value of correntropy criterion (35) at every iteration k. Combined with the stochastic gradient method, the MCC-SGA algorithm is expressed as

$$\widehat{\theta }\left(k\right)=\widehat{\theta }(k-1)+{\eta }_1\mathrm{grad}\left(J\left(\widehat{\theta }\right)\right){|}_{\widehat{\theta }(k-1)},$$

(36)

where $\mathrm{grad}(·)$ presents gradient, and can be calculated as

$$\mathrm{grad}\left(J\left(\widehat{\theta }\right)\right)=\frac{1}{{\sigma }^2}exp\left(-\frac{e^2\left(k\right)}{2{\sigma }^2}\right)e\left(k\right){\varphi }_c\left(k\right).$$

(37)

Using Equation (32) and substituting Equation (37) into Equation (36) yields

$$\begin{array}{cc}\hfill \widehat{\theta }\left(k\right)& =\widehat{\theta }(k-1)+\frac{{\eta }_1}{{\sigma }^2}exp\left(-\frac{{\left[Y_n\left(k\right)-{\varphi }_c^{\mathtt{T}}\left(k\right)\widehat{\theta }(k-1)\right]}^2}{2{\sigma }^2}\right)\left[Y_n\left(k\right)-{\varphi }_c^{\mathtt{T}}\left(k\right)\widehat{\theta }(k-1)\right]{\varphi }_c\left(k\right)\hfill \\ \hfill & =\widehat{\theta }(k-1)+\eta P_k\left[Y_n\left(k\right)-{\varphi }_c^{\mathtt{T}}\left(k\right)\widehat{\theta }(k-1)\right]{\varphi }_c\left(k\right),\hfill \end{array}$$

(38)

where $\eta =\frac{{\eta }_1}{{\sigma }^2}$ is the step size of updating formula, $P_k=exp\left(-\frac{{\left[Y_n\left(k\right)-{\varphi }_c^{\mathtt{T}}\left(k\right)\widehat{\theta }(k-1)\right]}^2}{2{\sigma }^2}\right)$ satisfies $0<P_k\le 1$. The presence of term $P_k$ makes the MCC-SGA algorithm have following properties:

1.

In a certain noisy environment, when $\sigma \to \infty$, $P_k\to 1$, then the MCC-SGA algorithm degenerates into

$$\widehat{\theta }\left(k\right)=\widehat{\theta }(k-1)+\eta \left[Y_n\left(k\right)-{\varphi }_c^{\mathtt{T}}\left(k\right)\widehat{\theta }(k-1)\right]{\varphi }_c\left(k\right),$$

(39)

which is the LMS-based stochastic gradient descent (LMS-SGD) algorithm.

2.

The MCC-SGA algorithm can be regarded as the LMS-SGD algorithm with a variable step size $\eta \left(k\right)=\eta P_k$, satisfying $\eta \left(k\right)\le \eta$. This leads to the convergence of the MCC-SGA algorithm being slightly slower than that of the LMS-SGD algorithm when using the same step size. But it also makes the MCC-SGA algorithm be more robust to change of step size.

3.

For a given $\sigma$, if the prediction error $e\left(k\right)$ is far from 0, the term $P_{k}e\left(k\right)$ will approximate to 0, then the MCC-SGA algorithm stops updating. Thus, the MCC-SGA algorithm given in Equation (38) has good robustness to the noise. And the smaller the value of $\sigma$, the higher robustness the MCC-SGA algorithm has.

After some simple calculations, Equation (38) can be rewritten as

$$\widehat{\theta }\left(k\right)=\left(I-\eta P_k{\varphi }_c\left(k\right){\varphi }_c^{\mathtt{T}}\left(k\right)\right)\widehat{\theta }(k-1)+\eta P_k{\varphi }_c\left(k\right)Y_n\left(k\right),$$

(40)

which can be viewed as a discrete-time system with the state $\widehat{\theta }\left(k\right)$. To guarantee the convergence of parameter $\widehat{\theta }\left(k\right)$, all the eigenvalues of symmetric matrices $I-\eta P_k{\varphi }_c\left(k\right){\varphi }_c^{\mathtt{T}}\left(k\right)$ should be within the unit circle [33]. Thus, the choice of $\eta$ should satisfy

$$0<\eta P_k\le \frac{2}{{\lambda }_{max}\left[{\varphi }_c\left(k\right){\varphi }_c^{\mathtt{T}}\left(k\right)\right]},$$

(41)

where ${\lambda }_{max}\left[{\varphi }_c\left(k\right){\varphi }_c^{\mathtt{T}}\left(k\right)\right]$ is the maximum eigenvalue of matrix ${\varphi }_c\left(k\right){\varphi }_c^{\mathtt{T}}\left(k\right)$.

Remark 3.

The data are collected at regular time-instants in the identification interval $[0,T]$. The number of sampling points M is selected to be as large as possible to reduce the approximate error generated by Equation (5), and also to ensure that the estimated parameters can reach a steady state.

3.4. Analysis of Estimation Bias

In this subsection, the bias of parameter estimation is deduced, and the impact of kernel width on estimation accuracy is studied. Before discussing, let us make an assumption that $\varphi \left(k\right)$ is not related to $\widehat{\theta }\left(k\right)$ and denote $\mathrm{E}\left[\varphi \left(k\right){\varphi }^{\mathtt{T}}\left(k\right)\right]=\Omega$.

First, substituting Equation (29) into Equation (38) yields

$$\widehat{\theta }\left(k\right)=\widehat{\theta }(k-1)+\eta P_k{\varphi }_c\left(k\right){\varphi }_c^{\mathtt{T}}\left(k\right)\theta -\eta P_k{\varphi }_c\left(k\right){\varphi }_c^{\mathtt{T}}\left(k\right)\widehat{\theta }(k-1)+\eta P_k{\varphi }_c\left(k\right)v_f\left(k\right).$$

(42)

Taking the expectation on both sides of Equation (42), we can obtain

$$\begin{array}{cc}\hfill \mathrm{E}\left[\widehat{\theta }\left(k\right)\right]& =\mathrm{E}\left[\widehat{\theta }(k-1)\right]+\eta P_k\mathrm{E}\left[{\varphi }_c\left(k\right){\varphi }_c^{\mathtt{T}}\left(k\right)\right]\theta \hfill \\ \hfill & -\eta P_k\mathrm{E}\left[{\varphi }_c\left(k\right){\varphi }_c^{\mathtt{T}}\left(k\right)\right]\mathrm{E}\left[\widehat{\theta }(k-1)\right]+\eta P_k\mathrm{E}\left[{\varphi }_c\left(k\right)v_f\left(k\right)\right].\hfill \end{array}$$

(43)

Based on the above assumption, Equation (43) can be transformed into

$$\begin{array}{cc}\hfill \mathrm{E}\left[\widehat{\theta }\left(k\right)\right]& =\mathrm{E}\left[\widehat{\theta }(k-1)\right]+\eta P_k\Omega \theta -\eta P_k\Omega \mathrm{E}\left[\widehat{\theta }(k-1)\right]+\eta P_k\left(G+{\Sigma }_s\right)\theta \hfill \\ \hfill & -\eta P_k\left(G+{\Sigma }_s\right)\mathrm{E}\left[\widehat{\theta }(k-1)\right]+\eta P_k\Lambda ,\hfill \end{array}$$

(44)

where $G=\mathrm{E}\left[\varphi \left(k\right)s^{\mathtt{T}}\left(k\right)\right]+\mathrm{E}\left[{\varphi }^{\mathtt{T}}\left(k\right)s\left(k\right)\right]$, ${\Sigma }_s=\mathrm{E}\left[s\left(k\right)s^{\mathtt{T}}\left(k\right)\right]$ and $\Lambda =\mathrm{E}\left[{\varphi }_c\left(k\right)v_f\left(k\right)\right]$. Then the deviation of MCC-SGA algorithm for each iteration is

$$b_{MCC}\left(k\right)=\eta P_k\left(G+{\Sigma }_s\right)\left(\theta -\mathrm{E}\left[\widehat{\theta }(k-1)\right]\right)+\eta P_k\Lambda .$$

(45)

Taking the absolute value of Equation (45), we have

$$\begin{array}{cc}\hfill |b_{MCC}\left(k\right)|& =\eta P_k\left|\left(G+{\Sigma }_s\right)\left(\theta -\mathrm{E}\left[\widehat{\theta }(k-1)\right]\right)+\Lambda \right|\hfill \\ \hfill & \le \eta \left|\left(G+{\Sigma }_s\right)\left(\theta -\mathrm{E}\left[\widehat{\theta }(k-1)\right]\right)+\Lambda \right|\hfill \\ \hfill & =\left|b_{MSE}\left(k\right)\right|,\hfill \end{array}$$

(46)

where $b_{MSE}\left(k\right)$ is the deviation of the LMS-SGD algorithm for each iteration.

Remark 4.

It can be seen from Equation (46) that the absolute deviation of the MCC-SGA algorithm for each iteration is no more than that of the LMS-SGD algorithm. Moreover, for a certain noisy environment, if σ decreases, $P_k$ decreases and then $\left|b_{MCC}\left(k\right)\right|$ decreases. This means that the kernel width is useful in reducing the estimation bias. In practice, if the data size is large, a small kernel size should be used to achieve the identification results with high accuracy. However, if the data size is small, the kernel size has to be chosen as a compromise between estimation accuracy and convergence speed.

4. Experiments and Results Analysis

In this section, two numerical examples are given. To illustrate the superiority of the proposed method, the LMS-SGD algorithm is utilized to identify the FOS. The relative error $\delta =\Vert \theta -\widehat{\theta }\Vert /\Vert \theta \Vert$ is adopted to objectively evaluate identification accuracy, where $\widehat{\theta }$ is the estimated parameter value and $\theta$ is the true one. Due to the stochastic noise, a Monte Carlo test is carried out, that is to say, both identification algorithms independently run 30 times. In each run, the algorithm terminates when the pre-specified maximum iteration number is reached.

4.1. Example 1

Consider the FOS [21]

$$D^{0.8}x\left(t\right)+a_{0}x\left(t\right)=bu\left(t\right)$$

(47)

with $a_0=1$ and $b=3$. The identification interval is $[0,16]$. A multi-sine function

$$u\left(t\right)=\sum _{k=1}^{25}sin\left(0.1\pi kt-\frac{k(k-1)\pi }{25}\right)$$

(48)

is chosen as the input signal to excite the system. The initial conditions of the FOS are assumed to be $x\left(t\right)=0,t\le 0$. The sampling period is set as $h=0.001$, then the number of BPFs to be used is $M=T/h=$ 16,000, which is also the maximum number of iterations.

Case 1. The effect of parameters in algorithms on the identification accuracy.

First, we set the step size of the MCC-SGA algorithm and the LMS-SGD algorithm setting as $\eta =0.003$, and the kernel width of the correntropy function as $\sigma =\{1.5,15\}$, respectively. The relative error curves of estimation results obtained by the MCC-SGA algorithm and LMS-SGD algorithm are plotted in Figure 3. It can be seen that as expected, for a large $\sigma$, the MCC-SGA algorithm performs the same as the LMS-SGD algorithm. For a small $\sigma$, the MCC-SGA algorithm can obtain more accurate identification results with a bit slow convergence speed.

Meanwhile, the impact of the step size of these two algorithms is studied. The kernel width of the correntropy function is set as $\sigma =1.5$, and the step size of the MCC-SGA algorithm and LMS-SGD algorithm is set as $\eta =\{0.003,0.05,0.07,0.1\}$.

Table 1. Relative errors of estimations obtained by algorithms with different values of $\eta$ in Experiment 1.

	$\mathit{\eta }\mathbf{=}\mathbf{0.003}$	$\mathit{\eta }\mathbf{=}\mathbf{0.005}$	$\mathit{\eta }\mathbf{=}\mathbf{0.007}$	$\mathit{\eta }\mathbf{=}\mathbf{0.01}$
MCC-SGA	$0.49\times {10}^{-2}$	$2.24\times {10}^{-2}$	$3.60\times {10}^{-2}$	$4.82\times {10}^{-2}$
LMS-SGD	$0.60\times {10}^{-2}$	$2.55\times {10}^{-2}$	$3.95\times {10}^{-2}$	∞

lists the relative errors of the identification results after 16,000 iterations. The results show that compared with the LMS-SGD algorithm, the MCC-SGA algorithm has more stable performance to the change of $\eta$, which makes it have a wider choice of step size and broadens its practical application range. Based on the above discussions, the parameters of identification algorithm are set as $\sigma =1.5$ and $\eta =0.003$ in the following experiments.

Case 2. Output signal is disturbed by different types of stable distribution noises.

In this case, FOS (47) is identified in an environment perturbed by Gaussian, Cauchy and Lévy noises, respectively. In the case of the Gaussian noise $v\left(t\right)\sim S(2,0,0.5,0)$; in the case of the Cauchy noise $v\left(t\right)\sim S(1,0,0.1,0)$; in the case of the Lévy noise $v\left(t\right)\sim S(0.5,1,{10}^{-5},0)$. The figures of these three types of noises are shown in Figure 4. It can be seen that as $\alpha$ decreases, the impulsivity of the noise becomes stronger.

The MCC-SGA algorithm and LMS-SGD algorithm are both used to identify the FOS. The evolution processes of the estimated parameters are shown in Figure 5, Figure 6 and Figure 7. The identification results after 16,000 iterations for the above two methods are given in

Table 2. Estimation results under Gaussian, Cauchy, Lévy noises after 30 runs.

Noise	Parameter	True Value	MCC-SGA		LMS-SGD
			avg	std	avg	std
Gaussian	$a_0$	1	$0.9897$	$0.0215$	$0.9891$	$0.0263$
	b	3	$2.9803$	$0.0263$	$2.9793$	$0.0389$
	$\delta$	-	$0.70\times {10}^{-2}$	-	$0.77\times {10}^{-2}$	-
Cauchy	$a_0$	1	$1.0391$	$0.3125$	$0.8418$	$1.2630$
	b	3	$2.9316$	$0.3400$	$3.0740$	$2.0630$
	$\delta$	-	$2.49\times {10}^{-2}$	-	$5.52\times {10}^{-2}$	-
Lévy	$a_0$	1	$0.9649$	$0.5288$	∞	∞
	b	3	$3.0950$	$0.5582$	∞	∞
	$\delta$	-	$3.07\times {10}^{-2}$	-	∞	-

.

The identification results show the following: (1) When $\alpha =2$, the LMS-SGD algorithm achieves the same identification accuracy as the MCC-SGA algorithm. The reason behind that is that the LMS-based identification method has good performance in a Gaussian noise environment. (2) When $\alpha =1$, the identification results achieved by the LMS-SGD algorithm are lower in accuracy. Moreover, when $\alpha =0.5$, the LMS-SGD algorithm cannot obtain identification results. (3) When the noise switches to Cauchy and Lévy noises, the identification accuracy of these two algorithms becomes worse, but the MCC-SGA algorithm can always achieve satisfactory identifications. Therefore, it can be concluded that the performance of the LMS-SGD algorithm deteriorates seriously when a large impulse noise exists. On the contrary, the MCC-SGA algorithm strains the effect of impulse noise, enhances the robustness, and improves the estimation accuracy.

Case 3. The effect of parameters in stable distribution noises on the identification accuracy.

In this case, parameters of the stable distribution noises are set as $\alpha =\{1,2\}$, $\beta =\{-1,0,1\}$, $\gamma =\{{10}^{-2},{10}^{-3}\}$, and $\mu =\{-0.1,0,0.1\}$. Both the MCC-SGA algorithm and LMS-SGD algorithm are used to perform the identification. The relative errors of identification results after 16,000 iterations are listed in

Table 3. Relative errors of estimations (avg ± std) obtained by MCC-SGA under stable distribution noises with different parameter settings after 30 runs.

$\mathit{\gamma }={10}^{-2}$, $\mathit{\mu }=-0.1$				$\mathit{\gamma }={10}^{-3}$, $\mathit{\mu }=-0.1$
	$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$		$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$
$\mathit{\alpha }=\mathbf{1}$	$0.106\pm 3.20\times {10}^{-2}$	$0.085\pm 1.21\times {10}^{-2}$	$0.091\pm 9.90\times {10}^{-2}$	$\mathit{\alpha }=\mathbf{1}$	$0.086\pm 6.60\times {10}^{-2}$	$0.082\pm 3.12\times {10}^{-2}$	$0.082\pm 1.81\times {10}^{-2}$
$\mathit{\alpha }=\mathbf{2}$	$0.083\pm 5.24\times {10}^{-4}$	$0.083\pm 5.14\times {10}^{-4}$	$0.083\pm 6.13\times {10}^{-4}$	$\mathit{\alpha }=\mathbf{2}$	$0.083\pm 5.85\times {10}^{-5}$	$0.083\pm 5.14\times {10}^{-5}$	$0.083\pm 1.59\times {10}^{-5}$
$\mathit{\gamma }={\mathbf{10}}^{-\mathbf{2}}$, $\mathit{\mu }=\mathbf{0}$				$\mathit{\gamma }={\mathbf{10}}^{-\mathbf{3}}$, $\mathit{\mu }=\mathbf{0}$
	$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$		$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$
$\mathit{\alpha }=\mathbf{1}$	$0.027\pm 9.29\times {10}^{-2}$	$0.012\pm 5.66\times {10}^{-2}$	$0.030\pm 6.53\times {10}^{-2}$	$\mathit{\alpha }=\mathbf{1}$	$0.003\pm 1.42\times {10}^{-2}$	$0.005\pm 1.33\times {10}^{-2}$	$0.006\pm 9.10\times {10}^{-3}$
$\mathit{\alpha }=\mathbf{2}$	$0.005\pm 5.77\times {10}^{-4}$	$0.005\pm 5.77\times {10}^{-4}$	$0.005\pm 5.96\times {10}^{-4}$	$\mathit{\alpha }=\mathbf{2}$	$0.005\pm 5.77\times {10}^{-5}$	$0.005\pm 5.77\times {10}^{-5}$	$0.005\pm 5.77\times {10}^{-5}$
$\mathit{\gamma }={\mathbf{10}}^{-\mathbf{2}}$, $\mathit{\mu }=\mathbf{0}.\mathbf{1}$				$\mathit{\gamma }={\mathbf{10}}^{-\mathbf{3}}$, $\mathit{\mu }=\mathbf{0}.\mathbf{1}$
	$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$		$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$
$\mathit{\alpha }=\mathbf{1}$	$0.087\pm 5.87\times {10}^{-2}$	$0.103\pm 2.90\times {10}^{-2}$	$0.122\pm 1.55\times {10}^{-2}$	$\mathit{\alpha }=\mathbf{1}$	$0.105\pm 4.70\times {10}^{-3}$	$0.103\pm 7.90\times {10}^{-3}$	$0.107\pm 1.27\times {10}^{-2}$
$\mathit{\alpha }=\mathbf{2}$	$0.102\pm 6.22\times {10}^{-4}$	$0.102\pm 6.20\times {10}^{-4}$	$0.102\pm 6.22\times {10}^{-4}$	$\mathit{\alpha }=\mathbf{2}$	$0.102\pm 6.23\times {10}^{-5}$	$0.102\pm 6.21\times {10}^{-5}$	$0.093\pm 1.34\times {10}^{-4}$

and

Table 4. Relative errors of estimations (avg ± std) obtained by LMS-SGD under stable distribution noises with different parameter settings after 30 runs.

$\mathit{\gamma }={10}^{-2}$, $\mathit{\mu }=-0.1$				$\mathit{\gamma }={10}^{-3}$, $\mathit{\mu }=-0.1$
	$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$		$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$
$\mathit{\alpha }=\mathbf{1}$	$0.148\pm 1.50\times {10}^{-1}$	$0.150\pm 1.91\times {10}^{-1}$	$0.102\pm 1.46\times {10}^{-1}$	$\mathit{\alpha }=\mathbf{1}$	$0.128\pm 1.22\times {10}^{-1}$	$0.112\pm 3.38\times {10}^{-2}$	$0.114\pm 4.47\times {10}^{-2}$
$\mathit{\alpha }=\mathbf{2}$	$0.112\pm 6.03\times {10}^{-4}$	$0.112\pm 6.19\times {10}^{-4}$	$0.112\pm 6.68\times {10}^{-4}$	$\mathit{\alpha }=\mathbf{2}$	$0.112\pm 6.02\times {10}^{-5}$	$0.112\pm 6.19\times {10}^{-5}$	$0.112\pm 2.22\times {10}^{-5}$
$\mathit{\gamma }={\mathbf{10}}^{-\mathbf{2}}$, $\mathit{\mu }=\mathbf{0}$				$\mathit{\gamma }={\mathbf{10}}^{-\mathbf{3}}$, $\mathit{\mu }=\mathbf{0}$
	$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$		$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$
$\mathit{\alpha }=\mathbf{1}$	$0.054\pm 2.48\times {10}^{-1}$	$0.057\pm 1.25\times {10}^{-1}$	$0.050\pm 5.24\times {10}^{-2}$	$\mathit{\alpha }=\mathbf{1}$	$0.009\pm 2.19\times {10}^{-2}$	$0.006\pm 1.34\times {10}^{-2}$	$0.005\pm 1.26\times {10}^{-2}$
$\mathit{\alpha }=\mathbf{2}$	$0.006\pm 6.66\times {10}^{-4}$	$0.006\pm 6.66\times {10}^{-4}$	$0.006\pm 6.45\times {10}^{-4}$	$\mathit{\alpha }=\mathbf{2}$	$0.006\pm 6.67\times {10}^{-5}$	$0.006\pm 6.67\times {10}^{-5}$	$0.006\pm 6.67\times {10}^{-5}$
$\mathit{\gamma }={\mathbf{10}}^{-\mathbf{2}}$, $\mathit{\mu }=\mathbf{0}.\mathbf{1}$				$\mathit{\gamma }={\mathbf{10}}^{-\mathbf{3}}$, $\mathit{\mu }=\mathbf{0}.\mathbf{1}$
	$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$		$\mathit{\beta }=-\mathbf{1}$	$\mathit{\beta }=\mathbf{0}$	$\mathit{\beta }=\mathbf{1}$
$\mathit{\alpha }=\mathbf{1}$	$0.119\pm 8.71\times {10}^{-2}$	$0.109\pm 7.24\times {10}^{-2}$	$1.56\times {10}^{-1}\pm 0.041$	$\mathit{\alpha }=\mathbf{1}$	$0.111\pm 6.10\times {10}^{-3}$	$0.109\pm 1.38\times {10}^{-2}$	$0.149\pm 1.91\times {10}^{-2}$
$\mathit{\alpha }=\mathbf{2}$	$0.108\pm 7.06\times {10}^{-4}$	$0.108\pm 6.82\times {10}^{-4}$	$0.108\pm 7.06\times {10}^{-4}$	$\mathit{\alpha }=\mathbf{2}$	$0.108\pm 7.06\times {10}^{-5}$	$0.108\pm 6.83\times {10}^{-5}$	$0.108\pm 6.88\times {10}^{-5}$

. These results demonstrate that under various stable distribution noises, the MCC-SGA algorithm can obtain more accurate estimations than its LMS-SGD counterpart. To better illustrate the effect of each parameter in noise, we plot 3D bar graphs in Figure 8 and Figure 9 to show the results listed in Table 3 and Table 4. It can be noted that the location parameter $\mu$ plays a significant role in the identification accuracy. A larger $\left|\mu \right|$ means the measurement output deviates farther from its true value, which results in poor estimations. When location parameter $\mu =0$, the exponent parameter $\alpha$ has an impact on identification, while when $\mu \ne 0$, $\alpha$ has almost no impact on identification. The dispersion parameter $\gamma$ has an minor impact on the identification accuracy, but it affects the stability of identification results. In addition, the value of symmetric parameter $\beta$ has a small impact on the estimation accuracy. In this experiment, the most accurate estimations can be achieved with $\mu =0,\alpha =2,\gamma ={10}^{-3}$. It can be concluded that when the stable distribution noise has a lighter tail, smaller absolute location parameter, and narrower dispersion, accurate identifications are achieved.

4.2. Example 2

To further validate the effectiveness of the proposed method, the FOS

$$D^{2.2}x\left(t\right)+a_1{D}^{1.1}x\left(t\right)+a_{0}x\left(t\right)=bu\left(t\right)$$

(49)

is considered. Here, $a_1=2$, $a_0=5$ and $b=4$. The identification interval is $[0,16]$. A multi-sine function,

$$u\left(t\right)=\sum _{k=1}^{200}sin\left(0.01\pi kt-\frac{k(200-k)\pi }{100}\right)$$

(50)

is chosen as the input signal. The initial conditions of the FOS are assumed to be $x\left(t\right)=0,t\le 0$. The sampling period is set as $h=0.0016$, then the number of BPFs involved in the proposed method and the maximum number of iterations are $M=T/h=$ 10,000.

Simulations without noise are omitted here to avoid repetition, and FOS (49) is identified when the output is disturbed by general impulse noises. In this case, $v\left(t\right)\sim S(1.5,0,0.2,0)$ and $v\left(t\right)\sim S(1.2,0,0.1,0)$, respectively. The MCC-SGA algorithm and LMS-SGD algorithm are both used to identify FOS. The kernel width is set as $\sigma =1.5$; the step size is set as $\eta =0.003$. The evolution processes of the estimated parameters are shown in Figure 10 and Figure 11. The identification results after $10,000$ iterations are given in

Table 5. Estimation results under general impulse noises after 30 runs.

Noise	Parameter	True Value	MCC-SGA		LMS-SGD
			avg	std	avg	std
$S(1.5,0,0.2,0)$	$a_1$	2	$1.9161$	$1.2632$	$1.7403$	$1.9338$
	$a_0$	5	$4.8356$	$1.1608$	$4.4576$	$2.4246$
	b	4	$4.1103$	$0.9572$	$3.5723$	$1.9572$
	$\delta$	-	$3.21\times {10}^{-2}$	-	$1.10\times {10}^{-1}$	-
$S(1.2,0,0.1,0)$	$a_1$	2	$1.8838$	$1.2222$	$0.8418$	$2.5206$
	$a_0$	5	$4.9404$	$1.0354$	$2.9793$	$2.3874$
	b	4	$4.2532$	$0.9320$	$3.0740$	$1.9482$
	$\delta$	-	$4.25\times {10}^{-2}$	-	$1.96\times {10}^{-1}$	-

.

It can be seen that when the location parameter $\mu =0$, the exponent parameter $\alpha$ has an impact on the identification accuracy. The identification error becomes larger with the decrease in $\alpha$. However, as expected, the proposed algorithm can always obtain more satisfactory results than its LMS counterpart.

5. Conclusions

In this paper, the identification problem of FOS is studied. By using block pulse operational matrix, the identified FOS is converted into an algebraic one. The correntropy is utilized as the objective function, and an iterative formula based on the stochastic gradient ascent method is designed to perform the identification. By studying the impact of parameters in the proposed algorithm on identification, appropriate kernel width and step size are recommended. To demonstrate the superiority of the proposed method, an LMS-based identification method is introduced, and these two algorithms are used to identify FOSs under noises of different types of stable distribution. Moreover, the impact of parameters of stable distribution on identification accuracy is examined. Simulation results demonstrate the following: (1) For a large kernel width, the MCC-SGA algorithm has the same performance as the LMS-SGD algorithm; for a small kernel width, the MCC-SGA algorithm can obtain more accurate estimations with a slightly slower convergence speed. Moreover, the MCC-SGA algorithm is more robust to the change in step size. (2) When the impulsivity of noise increases, the identification error becomes larger, but the proposed method is always superior to its LMS counterpart in identification accuracy and robustness. (3) The location of the noise distribution plays a significant role in the identification accuracy.

It can be believed that the proposed MCC-SGA algorithm provides a new approach to identify the FOS. In the future, FOS identification is still a challenging problem. The modification of the stochastic gradient algorithm to improve the convergence rate of proposed method and the identification of the fractional orders are the future research goals.