Generalized Correntropy Criterion-Based Performance Assessment for Non-Gaussian Stochastic Systems

Control loop performance assessment (CPA) is essential in the operation of industrial systems. In this paper, the shortcomings of existing performance assessment methods and indicators are summarized firstly, and a novel evaluation method based on generalized correntropy criterion (GCC) is proposed to evaluate the performance of non-Gaussian stochastic systems. This criterion could characterize the statistical properties of non-Gaussian random variables more fully, so it can be directly used as the assessment index. When the expected output of the given system is unknown, generalized correntropy is used to describe the similarity of two random variables in the joint space neighborhood controlled and take it as the criterion function of the identification algorithms. To estimate the performance benchmark more quickly and accurately, a hybrid-EDA (H-EDA) combined with the idea of “wading across the stream algorithm” is proposed to obtain the system parameters and disturbance noise PDF. Through the simulation of a single loop feedback control system under different noise disturbances, the effectiveness of the improved algorithm and new indexes are verified.


Introduction
With the continuous improvement of automation in the industrial process, the production process is becoming more and more complex. Control loops play the most important role in automation systems. Product quality, operation safety, material consumption and energy consumption are closely related to the performance of control system. Excellent control loop performance will ensure the effectiveness of the control systems, thus ensuring product quality and reducing product cost. The significance of performance assessment is to realize, restore and maintain the optimal performance of the control loops at all times. Therefore, it is of great practical value to evaluate the system performance quickly and accurately in industrial operation. At present, control loop performance assessment (CPA) has become an essential technology to ensure the smooth progress of industrial production [1].
Noise is unavoidable during the performance evaluation [2]. Most of the current performance evaluation methods assume that the noise obeys the Gaussian distribution. The research focuses on the first or second-order moments of the target output (i.e., the mean and variance). Harris [3] proposed a CPA index based on minimum variance control (MVC). This method is regarded as a milestone in the performance assessment research, so we named the index after the researcher, Harris. Later researchers have made great progress on this fundamental. The MVC-based performance evaluation index has been applied to various types of control systems, and this index has been improved for different situations [4][5][6]. At present, methods of system performance assessment based on MVC are very mature and has a remarkable application effect when the noise disturbance conforms generalized maximum correntropy criterion. Wang [18] applied generalized correntropy to the design of sparse Gaussian Hermite orthogonal filters. Qiu [19] replaced the mean square error loss with generalized correntropy loss (GCL) in the original unscented Kalman filter (UKF) framework, which combines the strength of GCL developed in robust information theoretic learning to solve the non-Gaussian interference and the strength of the UKF in handing strong model nonlinearities. Therefore, it is very promising to study the performance assessment of non-Gaussian stochastic systems by replacing entropy criterion with generalized correntropy.
To estimate the performance benchmark better, this paper integrates the idea of wading across the stream algorithm (WSA) into the traditional EDA to improve the local convergence ability of the algorithm. In the iterative process of the algorithm, the center of some excellent individuals is selected to cross with the best solution to fully retain the information of the best individuals. Generalized correntropy is used as the fitness value of H-EDA to identify the parameters of the controlled auto regression and moving average (CARMA) model of the controlled system, which accelerate the algorithm and improves the accuracy. The effectiveness of the hybrid algorithm is verified through the simulation of a single loop feedback control system under different noise disturbances.
In the next section, the CPA index based on MEC is reviewed, and the generalized correntropy index is proposed to overcome the shortcomings of the previous entropy indexes. The third part introduces the performance evaluation process based on the H-EDA in detail. In the fourth section, a detailed simulation of a single input single output (SISO) system is carried out, and the final conclusion will be given in the last part.

Feedback Control Loop
A SISO feedback control system is considered in this paper for convenience of research, the detailed description is given in Figure 1. where, r represents the system set value; u represents the output of the controller; v represents stationary independent identically distributed unmeasurable noise.
The system in Figure 1 could be described as a CARMA model as follows, where y(k) , u(k) and v(k) are the output, input and unknown noise disturbance of CARMA model process respectively, where n a , n b , n c and τ are the structural parameters of the model; n a , n b , and n c are the order of A(z −1 ), B(z −1 ) and C(z −1 ) respectively, τ is the system delay. The parameters a na , b nb , and c nc can be estimated by the model identification, provided the structural parameters are obtained.

Minimum Entropy Index
Assuming that the input of the given system is 0, the output is, where, Gp is the transfer function without delay. By solving the Diophantine equation, the disturbance transfer function G v could be further decomposed, where, F(z −1 is the impulse response coefficient of the disturbance transfer function G v and R(z −1 is the remaining transfer function satisfying the identity (4), where, The control objective of the minimum entropy controller can be achieved by minimizing the entropy of the output variables when analyzing the linear non-Gaussian system [21,22]. The feedback-invariant only depend on the characteristics of the disturbance in the process rather than the function of the process model or controller. The second term on the right of Equation (5) is the feedback-varying, which means the controller G v can determine the entropy value of the process output. Since both sides of Equation (6) are independent, therefore, Just like the MVC method, the equal sign in Equation (7) holds only when L = 0, then the output entropy reaches the minimum, that is the benchmark entropy, The CPA based on MEC is to compare the actual system output entropy H(y t ) with the output entropy H min (y t under MEC. To sum up, the CPA index based on MEC can be expressed as, This indicator is similar to MVC, which is always between 0 and 1. Generally, the closer it is to 1, the closer the system to the ideal case, indicating that the system performs better; conversely, the closer it is to 0, the worse the system performance is, even including unstable control.

Prevenient Entropy Index and Generalized Correntropy
Entropy is a measure to describe the uncertainty of random variables, even for random variables without mean or variance. The size of entropy is only determined by the shape of the distribution rather than its location. There are many ways to describe entropy, such as Shannon entropy (SE), rational entropy (RE), Renyi entropy and delta entropy. In the researches based on the minimum entropy criterion, SE is widely used. For linear Gaussian systems, SE is equivalent to variance. For a continuous random variable x, Shannon entropy can be expressed as, where γ(x) is its probability density function (PDF). SE is one of the CPA criteria based on MEC system in previous researches [10]. However, it does not meet the "consistency" principle, that is, the results must be the same when entropy could be calculated in many different ways. These characteristics show that SE cannot be used as a new evaluation index. Fortunately, the rational entropy (RE) criterion is proposed by Zhou [10], which has most of the properties of Shannon entropy and satisfies the "consistency" principle. The definition of rational entropy is, where R is the domain of random variable x, and γ(x) is its PDF. RE is a suitable benchmark for CPA of linear non-Gaussian systems. By comparing the output rational entropy under MEC process with actual, a CPA index will be obtained. Of course, RE is not perfect. The PDF of variable x is essential to calculate the output RE. If the PDF is known, it is more convenient to use the RE criterion to calculate the CPA index. However, when the PDF of the variable is unknown, this method is very complicated in the actual calculation process. For instance, when the system parameters are estimated by the EDA, RE as the fitness value, needs to be calculated for several times in each iteration (usually 1000 s), and the operation will be repeated for many times in the identification process. In the actual systems, the parameters and data will be more complex, this method will waste a lot of time. Therefore, this method lacks practical engineering significance.
The mixture correntropy index is adopted by Zhang [13] to solve the problem of translation invariance. However, the Gaussian kernel is applied to the kernel function of mixture correntropy in [13], which is not comprehensive in describing the statistical properties of non-Gaussian random variables [16]. As mentioned before, generalized correntropy is a more appropriate indicator.
Correntropy is a measure of local similarity, which is directly related to the similarity of two random variables in the joint space neighborhood controlled by kernel width. Given two random variables X and Y, their correntropy is defined as [16], where E[·] is the expectation value, κ σ is the kernel function with the width of σ(σ > 0) and F XY (x, y) denotes the joint distribution function of (X, Y). Generally, Gaussian kernel is chosen as the kernel function of correntropy, where e = x − y and λ = 1/2σ 2 is the kernel parameter. In this paper, a well-known GGD function will be used as the kernel function in correntropy [23], where Γ(·) means the gamma function, α > 0 represents shape parameter; β > 0 is bandwidth parameter, λ = 1/ β α is the kernel parameter and γ α, β = α/(2 βΓ(1/a)) is the normalization constant.
The GGD function is used as the kernel function in the correntropy, where, to distinguish the correntropy with Gaussian kernel, we define it as generalized correntropy (GC). Clearly, it is an extension of correntropy.
If the joint distribution f (x, y) of X and Y is known, the correntropy could be calculated by the following formula, Generally, the joint distribution of X and Y is hard to known. Similar to correntropy, generalized correntropy is also calculated by sample estimation, Next, some properties of generalized correntropy are summarized as follows to explain the internal conditions for it to be used as a performance index: Property 3: The generalized correntropy contains the absolute higher moments of n! E |X − Y| αn . These characteristics show that generalized correntropy can be applied to CPA. When the expected output of the given system is known, the target PDF is expected to be tracked by the controlled output PDF as much as possible, that is, the expected output R k and the actual output Y k are as consistent as possible, e k = R k − Y k is close to 0 infinitely. Ideally, the target PDF is tracked by the output completely, the performance of the control system reaches the optimum, and the generalized correntropy reaches the maximum value γ α, β . Therefore, the new CPA index based on generalized correntropy can be defined as, where: closer the index is to 1, the better the control system performance is; otherwise, the system performance is poor and needs to be improved. When the expected output of the system is unknown, it is essential to identify the associated systems. In this case, generalized correntropy could be combined with the EDA to generate a new identification algorithm.
In data analysis of regression and classification, a measure called the correntropic loss (C-loss) is always used instead of correntropy [24]. The GCL between X and Y could be defined as, Obviously, minimizing GCL corresponds to maximizing generalized correntropy.
is a sample taken from p XY , the estimated value of GCL is, GCL could be applied as the fitness value of the H-EDA to search for the optimal parameters and estimate the noise PDF of the given system; the specific algorithm will be given in Section 3.

Improvement and Application of EDA in CPA
In essence, the problem of system parameter identification is an alternative problem of high-dimensional parameter space optimization, so the system parameters can be obtained by the methods of high-dimensional parameters optimization. To obtain the noise PDFs, we have to know the order, time delay and the parameters of the given system.
For the estimation of delay, the simplest and most commonly used method is to analyze the correlation between the input u t and output y t signals [25], which can be expressed as,τ When τ = τ d , R y, u (τ) reaches the peak, the corresponding time is the estimated value of delay.
The Akaike information criterion [26] is applied to get the order of the system in this paper. For the given CARMA model, AIC criterion is as follows, n c is the order of the noise model, AIC(n) = Llgσ 2 e + a(n a + n b + n c )

EDA and Improved
Estimation of distribution algorithm (EDA) is a randomized search heuristic algorithm, which is used to create a probabilistic model of solution space. This model is updated iteratively based on the quality of the solutions sampled by the model [27]. In recent years, EDA as an optimization method has received great attention, and applied to scheduling, project, machine learning and identification problems successfully [27][28][29][30]. Wang [28] proposed a hybrid algorithm named estimation of particle swarm distribution algorithm (EPSDA), which combines PSO (the local search method) with the EDA (the global search method) to improve the efficiency of solving nonlinear bilevel programming problems (NBLP). To overcome the instability of model updating in the iterative process, a novel EDA based on the classic compact genetic algorithm (cGA) is proposed by [27] to optimize the benchmark function ONEMAX. Liang [29] developed a new variant of Gaussian EDA (GEDA) to solve the problem of premature convergence in traditional GEDA. In the process of CPA, EDA is mainly devoted to estimate the parameters of the given systems. Through appropriate improvements, the given system can be identified faster and more accurately, so as to obtain a more accurate benchmark.

Parameter Identification Based on Hybrid-EDA
As the previous work shows, the traditional EDA still has some shortcomings in the optimization process, such as a large amount of calculation, poor accuracy, etc. Therefore, based on the traditional EDA, the space for parameter identification is first obtained through preliminary estimation. To improve the local convergence ability, incorporating WSA [31] ideas into the EDA. In the iterative process, the centers of some excellent individuals are selected to cross-operate with the best individuals, the information of outstanding individuals could be retained more fully.

Acquisition of Parameter Identification Space
The initialization of parameter space is of positive significance to improve the identification accuracy and shorten the time required, so it is necessary to estimate the model parameters preliminarily. The recursive extended least squares (RELS) algorithm is selected in this paper.
The above CARMA model can be expressed as, where, Since v(k) is unmeasurable, it can only be replaced by its estimated valuev(k), where, The purpose is to obtain the parameter vectorθ when the objective function J(θ) reaches the minimum value. J(θ) is defined as follows, Based on the estimation of Gaussian system parameters, the initial estimates of varianceσ v and parametersθ could be obtained, then the parameter identification space Ω W could be taken asθ ± 3σ v .

The Idea of Wading Across Stream Algorithm
Wading across stream algorithm (WSA) [31] is inspired by the idea of "crossing the river with stones". This algorithm examines the "shore" carefully to select an initial starting point firstly, several solutions are randomly searched out in the neighborhood near the starting point to get the best one as the iterative result. Then take this result as the starting point and continue search for several solutions in the nearby neighborhood, select the optimal solution as the third iterations result, and so on until the termination condition is satisfied. The WSA is similar to the simulated annealing algorithm but the effect is better, and the algorithm idea is simple with good local convergence ability, so we integrate this algorithm idea into the EDA to improve its computational efficiency.
To calculate the fitness value of these R solutions, select the optimal individualθ * as the initial starting point according to the fitness value. For ease of description, denote the starting point is asθ * = [x 1 * , x 2 * , · · ·, x n * ], where n = na + nb + nc + 1.
(2) Search strategy Due to the complexity of parameter identification in high-dimensional space, the design of an effective local search strategy could improve the quality of the solution significantly. According to the idea of "crossing the river with stones", when you touch a "stone", you must take that "stone" as the starting point to explore other "stones" around it. Similarly, search for m individuals in the neighborhood radius L k around the starting point θ * to form a new parameters space B, B = [θ 1 ;θ 2 ; · · ·;θ m ] =     â 11 · · ·â 1nab 10 · · ·b 1nbĉ 11 · · ·ĉ 1nĉ a 21 · · ·â 2nab 20 · · ·b 2nbĉ 21 · · ·ĉ 2nc . .
where x ij = x j * +L k · r ij , r ij represents uniformly distributed random numbers in the interval [-1, 1]. Generally, the initial value L 0 of neighborhood radius is taken as one tenth of the whole parameters space, that is L 0 = (0.1 ∼ 0.3)σ v . To accelerate the later convergence of the algorithm, reduce the neighborhood radius gradually with the iterative process. The simplest method is to set L k = αL k−1 , α is the shrinkage coefficient, generally [0.90, 0.99]. The latest optimal individualθ * can be obtained by calculating the fitness values of these m solutions.
Calculate the center point of D, Then, the center point is crossed with the best solution, and the optimal individualθ * is modified as follows,θ where: a is a random number between 0 to 1.
Through the crossover operation, the information of excellent individuals could be retained to the maximum extent, which can avoid falling into the local optimum.

Selection of Fitness Value
From Equation (23), for the CARMA model in this paper, the error at time k could be defined as, For the estimated error sequence e = [e 1 , e 2 , . . . e L ], the GCL could be defined as, The optimal model parameters and the relationship between output and disturbance could be obtained by minimizing GCL. When the GCL reaches the minimum, the corresponding parameters are also optimal, namely, where: Ω W is the parameter identification space.

Algorithm Summary
The following Algorithm 1 is the program steps of the H-EDA identification algorithm.
Algorithm 1 Program steps for the identification algorithm.
1. Describe the system by CARMA model; estimate the system delay τ by analyzing the correlation between u(t) and y(t); determine the model order (n a , n b , n c ). After the coefficientF and noisev t of the feedback-invariant are obtained, the benchmark entropy and output entropy could be obtained by Equations (8) and (11), and the CPA index could be obtained by Equation (9),

Algorithm Validation and Sensitivity Analysis of Initial Parameters
In this section, a model example is given to illustrate the initial parameters of improved algorithm, which further analyses the H-EDA and traditional EDA, consider the following system, For the above example, the parameter vector to be identified is θ= [−1 .7, 0.7, 1.5, 0 .9]. Most of the initial parameters are based on the previous works [12][13][14][31][32][33]. Generally, in the EDA, the number of individuals in the initial population is N = 1000; the optimal number of individuals for each iteration is m = 200; and the maximum number of iterations is n max = 120. In WSA, the number of excellent individuals for crossover operation is N * = 30, R = 800.
In this section, we focus on the sensitivity of shrinkage coefficient α and the initial value of neighborhood radius L 0 . As mentioned before, the initial neighborhood radius is generally L 0 = (0.1 ∼ 0.3)σ v , and the shrinkage coefficient α is between 0.90 to 0.99. To find a proper set of parameters as much as possible, we have carried out a large number of experiments and some representative results are selected and attached in Table 1. In addition, we analyze the significance of a to the algorithm, a is a random number between [0, 1]. If a == 1, it means that the crossover operation is not carried out. A novel parameter vector is defined as P = [L 0 , α,a]. The test results show that the H-EDA is obviously superior to the traditional EDA. The global search capability of the algorithm is mainly determined by the neighborhood radius L k . Too small a value of L k prompts the algorithm to get trapped in local optima, while too large of a value of L k slows down the algorithm; α can accelerate the later convergence of the algorithm, but too small of an α leads to rapid convergence and bad accuracy; the main function of a is to reduce iterations and prevent local optimality. In short, the setting of initial parameters must rely on a large number of experiments and the summary of previous works. For this test case, the optimal initial parameter vector can be selected as P = [0.15σ v , 0.95, rand].

Simulation and Verification
To verify the effectiveness of the identification algorithm and the rationality of the performance assessment index, the following system is considered and the numerical simulation is carried out with Gaussian and non-Gaussian noise signals.
The transfer function of the controller is as follows, It is easy to know that the parameter vector of the system is θ = [−1, 1, −1, −0.2], system delay τ = 2, feedback-invariants F = [1, 0.8], and the controller gain in the simulation is K = 1.2.
For more precise description, the follow-up content is divided into the following parts: in Section 4.1, the effect of the improved algorithm and the improved benchmark entropy on the performance evaluation are explained respectively when the expected output is unknown, and their effectiveness are verified respectively. In Section 4.2, the tracking index is used to evaluate the system performance when the expected output of the system is known.

Simulation When Expected Output Distribution Is Unknown
It can be concluded from Section 3 that the system parameter identification and noise PDF estimation should be carried out firstly when the expected output is unknown, then the benchmark value could be obtained according to the feedback invariant coefficient estimationF and noise estimationv, and the performance evaluation index can be obtained finally. The specific process is shown in Figure 2. According to the analysis ideas in Section 3.3, the parameters of the algorithms in this section are set as follows: In the EDA, the number of individuals in the initial population is N = 1000, the optimal number of individuals for each iteration is m = 200, and the maximum number of iterations is n max = 120. Correspondingly, the parameters of H-EDA are set as, N * = 30, L 0 = 0.15σ v , α = 0.95, R = 800, the rests are consistent with the EDA. As the comparison algorithm, the parameters of the improved PSO (IPSO), particles N = 50, acceleration factor c 1 = 2, c 2 = 2, V max = 0.5.
The termination condition is J (l) − J (l−1) ≤ 0.0001 or reaches the specified number of iterations n max , J (l) is the l-generation GCL.
The noises distribution are subject to normal N(0, 0.255), Beta B (2,9), exponential E(0.5) and bimodal, where (1) When the fitness value is fixed as the GCL, compare the H-EDA with the traditional EDA and the IPSO.
Combine the above three algorithms with GCL, respectively. The probability distributions of the four kinds of noise are shown in Figure 3. It can be determined that the other three kinds of noise are non-Gaussian except for normal distribution. The results of the system parameter identification are shown in Table 2. The results in Table 2 show that the H-EDA could estimate the parameters of CARMA model under different noises when the fitness value is determined, and the identification accuracy and speed are better. Moreover, the fitness value of the algorithm is the generalized correntropy loss function, which proves the effectiveness of the novel criterion in the identification process. It is evident that the noises distribution estimated by the H-EDA in Figures 4 and 5 are very close to the real values. IPSO, as the comparison, is inferior to the H-EDA and cannot estimate the noise of bimodal distribution. Therefore, the superiority of the H-EDA is proved by the simulation results. According to the noise PDFv t and the feedback invariant coefficientF, we can get the performance assessment indexes, as shown in Table 3, where η ME represents the theoretical benchmark index;η MJ−EDA is the index obtained by identification results of the traditional EDA based on GCL;η MJ−H−EDA corresponds to the H-EDA index. Obviously, η MJ−H−EDA is closer to the theoretical value, which proves the effectiveness of the hybridalgorithm again.   (2) When the identification algorithm is fixed as the H-EDA, compare the traditional correntropy (TC) with the generalized correntropy criterion (GCL).
In the previous section, we have proved the effectiveness of generalized correntropy criterion in the performance assessment process. Next, we compare the superiority of GCL with TC based on the H-EDA. The parameters of algorithms, identified objects and noise PDFs are consistent with (1). The simulation results are illustrated in Table 4 and Figure 6.  The results show that when the disturbance obeys Gaussian, non-Gaussian and bimodal distribution, the identification results based on GCL are more accurate, and the time required to obtain the evaluation index is shorter.

Simulation When the Expected Output Distribution Is Known
Assumed that the known expected output distribution is Gaussian with mean value of 0 and variance of 1, as described in Section 2, the output PDF is expected to track the target PDF as much as possible.
The system under Gaussian and non-Gaussian noise disturbance are simulated respectively when the expected output distribution of the system is known. The output PDF tracks the input PDF as shown in Figure 7. The noise in (a) obeys normal distribution, while the noise in (b) obeys the non-Gaussian distribution. The generalized correntropy performance assessment index could be calculated by Formula (18) and shown in Table 5.  In Table 5, η GC is the generalized correntropy index and η C is the traditional correntropy index. Assessment results of generalized correntropy are consistent with the traditional correntropy index when the controller gain K is different, which also proves the effectiveness of the new index.

Conclusions
A new CPA method for linear non-Gaussian systems based on generalized correntropy is proposed by analyzing the limitations of existing performance evaluation methods and indicators. The simple tracking index based on generalized correntropy is given directly when the expected output is known. When the expected output is unknown, the generalized correntropy loss function is used as the fitness value of the H-EDA to estimate the system parameters and noise distribution, so as to calculate the system performance benchmark more quickly and accurately, and then the performance assessment results can be obtained. The single variable method is used to simulate the SISO control system under different noise disturbances, and the effectiveness of the new index and H-EDA is verified respectively. The results show that generalized correntropy index is superior to the traditional correntropy in evaluating the control loop performance of linear systems with non-Gaussian stochastic disturbances, and the H-EDA could also identify the system parameters more quickly and accurately. Then, the follow-up research still focuses on the selection of the new performance assessment index and the improvement of identification algorithm.