Non-Gaussian Systems Control Performance Assessment Based on Rational Entropy

Control loop Performance Assessment (CPA) plays an important role in system operations. Stochastic statistical CPA index, such as a minimum variance controller (MVC)-based CPA index, is one of the most widely used CPA indices. In this paper, a new minimum entropy controller (MEC)-based CPA method of linear non-Gaussian systems is proposed. In this method, probability density function (PDF) and rational entropy (RE) are respectively used to describe the characteristics and the uncertainty of random variables. To better estimate the performance benchmark, an improved EDA algorithm, which is used to estimate the system parameters and noise PDF, is given. The effectiveness of the proposed method is illustrated through case studies on an ARMAX system.


Introduction
Poor control loop performance will reduce the effectiveness of the control loop, which may lead to product quality degradation, increase product costs and other issues. There are many factors in the chemical process, which include inadequate parameter tuning and maintained controllers, equipment failure, without or insufficient feedforward compensation, inappropriate control structure design and so on, that can affect control loop performance [1]. The purpose of Control loop Performance Assessment (CPA) is to provide a comprehensive health assessment framework for control loops. Such framework includes assessing, detecting and diagnosing as well as suggesting improvement measures [2]. Harris (1989) proposed a performance indicator based on minimum variance control (MVC) [3], which laid the foundation for the development of CPA research field. Based on the MVC assessment technique, the benchmark value can be estimated from routine operational data as long as the delay is known or estimated, without additional experiments. Later researchers extended the scope of the basic Harris indicator to make the MVC-based performance assessment method suitable for feedforward/feedback control systems [4], unsteady and non-minimum phase systems [5], system for setting changes, time-varying system [6] and MIMO system [7]. Although the MVC-based controller can minimize the variance of the output, the robustness of the control system often fails to meet the requirements. In addition, constraints are often encountered during the actual production process. Therefore, Grimble et al. extended minimum variance (MV) benchmark and proposed the generalized minimum variance(GMV) algorithm to solve these problems. The GMV algorithm uses the sum of the minimized system output variance and the constrained control signal variance as an objective function and adds a dynamic minimum variance weight matrix to it [8]. Grimble (2002) used the

MEC Index
Consider a generic feedback control systems shown in Figure 1, where r t is the set point, u t is the controller output, v t is the unmeasured disturbance. G c , G p and G v denote the transfer functions of the feedback controller, the process and disturbance dynamics, respectively. The set point is set to zero by convenience and the disturbances are assumed to be zero mean.

MEC Index
Consider a generic feedback control systems shown in Figure 1, where is the set point, is the controller output, is the unmeasured disturbance. , and denote the transfer functions of the feedback controller, the process and disturbance dynamics, respectively. The set point is set to zero by convenience and the disturbances are assumed to be zero mean. Let the system under consideration be described by an ARMAX model.
where ( ), ( ), ( ) are the output, the input and the noise of the ARMAX process, respectively. The disturbance transfer function in Figure 1 can be further decomposed as follows by Diophantine equation, where ( −1 ) is the impulse response coefficients of in −1 with order − 1 and ( −1 ) is the remaining transfer function that satisfies the Identity (2).
where = −1 + −̃. The feedback-invariant terms are not functions of the process model or the controller; they depend only on the characteristics of the disturbance acting on the process. The second term is feedback-varying. This means that of the process output entropy (Equation (3)) depends on the structure and parameters of the controller ( ). The entropy of the output variable can reach the minimum value if = 0.
For non-Gaussian variables, unlike Gaussian variables which have the particularity that all distribution information is contained in the first and second moments and the higher moments above the second moment are zero. So, the MVC control that only minimize the second order does not apply to the non-Gaussian systems. Fortunately, the entropy is alternative uncertainty measurement which is more general in representing the system randomness using the probability distribution that all the stochastic information is included. Therefore, all higher order moments including the second one can be optimized using the entropy instead of the mean square error optimization.
For a linear non-Gaussian system, the goal of the minimum entropy controller is to minimize the entropy of the system output variables [25][26][27][28]. Like the conventional MVC, the minimum entropy value will be obtained if and only if = 0, Let the system under consideration be described by an ARMAX model.
where y(t), u(t), v(t) are the output, the input and the noise of the ARMAX process, respectively. The disturbance transfer function G v in Figure 1 can be further decomposed as follows by Diophantine equation, where F q −1 is the impulse response coefficients of G v in q −1 with order τ − 1 and R q −1 is the remaining transfer function that satisfies the Identity (2).
where L = The feedback-invariant terms are not functions of the process model or the controller; they depend only on the characteristics of the disturbance acting on the process. The second term is feedback-varying. This means that of the process output entropy (Equation (3)) depends on the structure and parameters of the controller (G c ). The entropy of the output variable can reach the minimum value if L = 0. For non-Gaussian variables, unlike Gaussian variables which have the particularity that all distribution information is contained in the first and second moments and the higher moments above the second moment are zero. So, the MVC control that only minimize the second order does not apply to the non-Gaussian systems. Fortunately, the entropy is alternative uncertainty measurement which is more general in representing the system randomness using the probability distribution that all the stochastic information is included. Therefore, all higher order moments including the second one can be optimized using the entropy instead of the mean square error optimization.
For a linear non-Gaussian system, the goal of the minimum entropy controller is to minimize the entropy of the system output variables [25][26][27][28]. Like the conventional MVC, the minimum entropy value will be obtained if and only if L = 0, H min (y t ) = H(Fv t ). (4) MEC based assessment compares the actual system-output entropy H(y t ) to the output entropy H min (y t ) as obtained using minimum-entropy controller. And the MEC-based CPA index is represented by where H min (y t ) is the entropy of the output variable with MEC and H(y t ) is the entropy of the output variable with actual controller. This index is similar to the MVC index will be always within the interval [0, 1], where MVC index values close to unity indicate good performance with regard to the theoretically achievable output minimum entropy. "0" means the worst performance, including unstable control. In fact, the relative entropy or Kullback-Leibler (KL) divergence [29,30] can reflect the distance between two probability distributions. It is defined as follows.
The relative entropy can be also used as a performance assessment index if p(·) denotes the probability density function (PDF) of the output variable with MEC and q(·) denotes the PDF of output variable with actual controller.
This index equal to zero indicates that the controller is the minimum entropy controller and when it deviates from zero, it is not a minimum entropy controller. However, the problem of the relative entropy-based assessment index is that the corresponding the relative entropy index is not a convex index. In other words, this index can obtain whether the current controller is the minimum entropy controller but it is difficult to find a suitable threshold to determine the current controller performance is good or bad. In this sense, it is not appropriate to use the KL distance as a control loop performance assessment index.

Rational Entropy
In [25,26], the authors gave a method for calculating the entropy, which is based on the following lemma.
Lemma 2 [25]. For the two random variables X and Y, the entropy or the amount of information is revealed by If X and Y are mutually independent, H(X|Y) = H(X), then H(X, Y) = H(X) + H(Y). However, the conclusion of Equation (6) is wrong because the two lemmas are unsuitable for the considering condition. For the Lemma 1, it is valid only for discrete random variables. But for continuous random variables, we will discuss whether Lemma 1 is established by the following examples.
Suppose the system is as follows, y t = a t + ca t−1 .
According to (6), the entropy is obtained as follows, Then for the same distribution a t , we use the different coefficients c to get the probability distribution as shown in Figure 2. Obviously, the different coefficients lead to the different distribution of ca t . Since the entropy is determined by the shape of the distribution, that is, the coefficients cannot be omitted, H(ca t ) = H(a t )(c = 1).
According to (6), the entropy is obtained as follows, Then for the same distribution , we use the different coefficients to get the probability distribution as shown in Figure 2. Obviously, the different coefficients lead to the different distribution of . Since the entropy is determined by the shape of the distribution, that is, the coefficients cannot be omitted, ( ) ≠ ( )( ≠ 1). As for the Lemma 2, probability theory shows that if the PDFs of and are known and and are mutually independent, the PDF of = + is calculated as follows [28], where (•) , (•) and (•) are the PDFs of , and . The random variable is still a univariate variable but Lemma 2 uses the properties of multivariate random variables. The distribution of the sum of two random variables is not the same as the joint distribution. In other words, ( ) ≠ ( ) + ( ). Due to confusion of concepts, the results of the performance evaluation are not credible.
As aforementioned, the entropy is determined by the shape of the distribution (or the shape of the PDF), then the entropy of the feedback invariant is determined by the shape of the PDF of + 1 −1 + ⋯ + −1 − +1 . In [27], a method of computing the entropy of feedback invariants called the consistent discrete distribution approximation method is introduced, which improves the deficiencies of literature [25]. Although this method ensures the unity and consistency of entropy calculation, different standards are required for different non-Gaussian noise. Therefore, the identification of distribution function is an indispensable step and it will be illustrated in subsequent sections.
In the previous studies, Shannon Entropy is one of performance assessment criteria based on minimum entropy control [28]. It is defined as, It seems to be a new benchmark to describe. But the Shannon entropy of the continuous random variable may be negative or even negative infinite, this means that the SE does not satisfy the "consistency" property. As a result, its uncertainty determines that it cannot be used as a new As for the Lemma 2, probability theory shows that if the PDFs of X and Y are known and X and Y are mutually independent, the PDF of Z = X + Y is calculated as follows [28], where f X (·), f Y (·) and f Z (·) are the PDFs of X, Y and Z. The random variable Z is still a univariate variable but Lemma 2 uses the properties of multivariate random variables. The distribution of the sum of two random variables is not the same as the joint distribution. In other words, H(Z) = H(X) + H(Y). Due to confusion of concepts, the results of the performance evaluation are not credible. As aforementioned, the entropy is determined by the shape of the distribution (or the shape of the PDF), then the entropy of the feedback invariant is determined by the shape of the PDF of v t + n 1 v t−1 + · · · + n τ−1 v t−τ+1 . In [27], a method of computing the entropy of feedback invariants called the consistent discrete distribution approximation method is introduced, which improves the deficiencies of literature [25]. Although this method ensures the unity and consistency of entropy calculation, different standards are required for different non-Gaussian noise. Therefore, the identification of distribution function is an indispensable step and it will be illustrated in subsequent sections.
In the previous studies, Shannon Entropy is one of performance assessment criteria based on minimum entropy control [28]. It is defined as, It seems to be a new benchmark to describe. But the Shannon entropy of the continuous random variable may be negative or even negative infinite, this means that the SE does not satisfy the "consistency" property. As a result, its uncertainty determines that it cannot be used as a new standard. Fortunately, a rational entropy (RE) instead of the SE is proposed by Zhou [24]. This type of entropy exhibits most properties of the Shannon's entropy and, at the same time, satisfies the "consistency" property. In this paper, we use the rational entropy of the process output with MEC and the actual output to calculate performance index. Let x be a random variable in R and γ(x) be its PDF, the rational entropy (RE) is given as [24] Although the expression of RE is similar to that of the relative entropy, RE and relative entropy have different meanings: RE reflects the uncertainty of random variables and relative entropy reflects the distance between two probability distributions.
In fact, [24] gave a CPA index of the output stochastic distribution control (SDC) systems. However, [24] only gave the calculation method of the theoretical benchmark value and did not give the estimation method of this benchmark; On the other hand, the index of [24] was only for the SDC systems, so its theory benchmark was not generic. In other words, a CPA framework that directly uses the characteristic of continuous random variables for general non-Gaussian systems has not yet been established. This paper aims to build a MEC-based CPA index for the general feedback control system with non-Gaussian disturbances.
In the chemical process, its output is generally measurable but it is difficult to accurately obtain the process model due to the lack of complete physicochemical knowledge and random disturbance distribution. It means the actual output entropy can be obtained directly from the collected output samples by Equation (11) but the entropy of MEC process still needs to be estimated through the identification of system parameters and noise PDF estimation. To summarize, the complete algorithm to evaluate the MEC-based index and to assess feedback controls contains the steps described as follows.
(1) Select the time-series-model type. Determine/estimate the system time delay τ and system order.
(2) Identify the closed-loop model from collected output samples.
(3) Estimate the benchmark entropy of process data. (4) Compute the performance index.

System Identification and Noise PDF Estimation
From Section 2, it is easy to know that system parameters and PDF of noise need to be accurately estimated in order to obtain RE-based CPA values.
Let q −1 be the unit backward shift operator and define the three polynomials in q −1 as: Defined A(q), B(q) and C(q) are polynomials in q −1 of order na, nb and nc, respectively. There are many kinds of algorithms to estimate the order [31] and delay [32] of the model. In this paper, we adopt that the order of the model is obtained by the Akaike information criterion [31]. Then, a simple and easy method, called correlation analysis, is applied here to estimateτ.
By means of indirect identification, a recursive extended least squares (RELS) algorithms is adopted then the closed loop non-Gaussian system identification is turned into open-loop Gaussian system identification. We can obtain the initial estimate of the parametersθ and the estimation of noise varianceσ v . Based on the estimation of the parameters of the Gaussian system, we can use theθ ± 3σ v Entropy 2018, 20, 331 7 of 15 as the initialized range of the parameter space. Then the improved EDA algorithm is used to obtain the system parameters and the noise distribution estimation.
The Estimation of distribution algorithm (EDA) is a population evolutionary algorithm based on statistical learning theory. The probability model is used to describe the distribution information of candidate solutions in the search space. A statistical learning method is used to establish a probabilistic model describing the distribution from the perspective of population. Then, a random sampling of the estimated probability distribution model is used to generate some new individuals to replace some of the individuals with poor fitness values in the initial population to form a new generation population. When satisfying the iteration termination condition, the iteration of this algorithm will be terminated and finally the optimal outcomes obtained by employing EDA are the best fitness value of the current population.
Based on the traditional EDA algorithm, the parameters of preliminary estimation and data selection are added to improve the speed of searching and optimization precision. And the system parameter identification problem is transformed into the optimization problem in high dimensional parameter space. In order to make the parameter space covers the real parameters as much as possible, the preliminary estimation is used to determine the range of initialization parameter space.
Then, we adoptθ as seeds. In order to eliminate the seeds with large deviations, adding the screening conditions. Since the mean of the noise is assumed to be zero, the mean of the residuals can be regarded as the screening condition. It can be described as follows.
Finally, we take the error entropy as the fitness value, when the error entropy reaches the minimum while the parameter reaches the optimum. The improved EDA algorithm is summarized as follows.
(1) Preliminary estimation. Rough estimation parameters are obtained by recursive maximum likelihood method. Then the task of initializing the parameters space can be accomplished with a uniform distribution using the initial value of parameter as the value range. (2) Screen seeds. From the first generation randomly select R group seed from the parameter space and calculate the mean value of the residuals. Remove the seeds whose mean value is more than ε and keep the seeds whose mean value is less than ε. If the number of reserved seeds is less than R, re-sampling is performed until the number of seeds retained is not less than R and then the seeds reserved are collected in new parameter spaces Ψ (l) (Q)(Q ≥ R).

Simulation
In order to illustrate the effectiveness of the proposed method, the following system which is the same as that of [25,27] is considered. To demonstrate more clearly, the section is divided as follows: In Section 4.1, parameters estimation and CPA with a fixed controller and unimodal distribution noises is discussed. Then keep all parameters unchanged, a simulation experiment is carried out for the bimodal distribution noises in Section 4.2 and the proposed CPA indices with different gain controllers are given in Section 4.3.

Parameters and CPA Index Estimation with Fixed Controller and Unimodal Distribution Noises
In this simulation, the gain of the controller is chosen as = 1.2. is assumed to follow a unimodal distribution including normal distribution N(0,0.255) , Beta distribution B (2,9) and exponential distribution E(0.5). The corresponding parameter estimation results are as follows

Simulation
In order to illustrate the effectiveness of the proposed method, the following system which is the same as that of [25,27] is considered.
The transfer function of the controller is chosen as G c = To demonstrate more clearly, the section is divided as follows: In Section 4.1, parameters estimation and CPA with a fixed controller and unimodal distribution noises is discussed. Then keep all parameters unchanged, a simulation experiment is carried out for the bimodal distribution noises in Section 4.2 and the proposed CPA indices with different gain controllers are given in Section 4.3.

Parameters and CPA Index Estimation with Fixed Controller and Unimodal Distribution Noises
In this simulation, the gain of the controller is chosen as K = 1.2. v t is assumed to follow a unimodal distribution including normal distribution N(0, 0.255), Beta distribution B (2,9) and exponential distribution E(0.5). The corresponding parameter estimation results are as follows The results show that the parameters of the ARMAX model can be estimated by the improved EDA algorithm under different noises. However, we are more concerned with the noise estimated values and its distribution.
Noise estimated values can be obtained by the improved EDA algorithm. There are many ways to describe the estimation effect, such as histogram and kernel density estimation. Since PDF is used in the MEC-based CPA index, the kernel density estimation method for actual noises and estimated noises will be used to show the estimation effect. The corresponding PDF estimation results of actual noises and estimated noises are shown in Figures 4-6. The results show that the parameters of the ARMAX model can be estimated by the improved EDA algorithm under different noises. However, we are more concerned with the noise estimated values and its distribution.
Noise estimated values can be obtained by the improved EDA algorithm. There are many ways to describe the estimation effect, such as histogram and kernel density estimation. Since PDF is used in the MEC-based CPA index, the kernel density estimation method for actual noises and estimated noises will be used to show the estimation effect. The corresponding PDF estimation results of actual noises and estimated noises are shown in Figures 4-6.   After identification by the improved EDA algorithm, not only the parameters of the system are obtained but also the distribution of noise is estimated. It is clear that the estimated disturbance distribution is close to the true system disturbance distribution. Hence, the simulation results in this case demonstrate the efficiency of the improved EDA algorithm. With this distribution and the feedback-invariant estimation, the MEC-based indices can be easily computed as shown in Table 1.   After identification by the improved EDA algorithm, not only the parameters of the system are obtained but also the distribution of noise is estimated. It is clear that the estimated disturbance distribution is close to the true system disturbance distribution. Hence, the simulation results in this case demonstrate the efficiency of the improved EDA algorithm. With this distribution and the feedback-invariant estimation, the MEC-based indices can be easily computed as shown in Table 1.  After identification by the improved EDA algorithm, not only the parameters of the system are obtained but also the distribution of noise is estimated. It is clear that the estimated disturbance distribution is close to the true system disturbance distribution. Hence, the simulation results in this case demonstrate the efficiency of the improved EDA algorithm. With this distribution and the feedback-invariant estimation, the MEC-based indices can be easily computed as shown in Table 1. Based on Table 1, the CPA index under MVC and MEC maintain at around 0.82 and 0.94, respectively. It indicates that when the system noise obeys the Gaussian distribution or the other unimodal distribution, there is no significant difference between the minimum variance and the minimum entropy criterion.
To illustrate that MEC-based CPA index is more applicable than MVC-based under the circumstances where the variance may fail. Bimodal distribution noises are selected as an example in the next subsection.

Parameters and CPA Index Estimation with Fixed Controller and Bimodal Distribution Noises
Keep all parameters except the distribution of v t unchanged and the distribution of v t is chosen as the following bimodal distribution. Parameter estimation results show that both methods seem to be available as system parameter estimates. However, the noise and its estimated distribution shown in Figure 7 indicate that the RELS method is ineffective in estimating the noise of the bimodal distribution. Based on Table 1, the CPA index under MVC and MEC maintain at around 0.82 and 0.94, respectively. It indicates that when the system noise obeys the Gaussian distribution or the other unimodal distribution, there is no significant difference between the minimum variance and the minimum entropy criterion.
To illustrate that MEC-based CPA index is more applicable than MVC-based under the circumstances where the variance may fail. Bimodal distribution noises are selected as an example in the next subsection.

Parameters and CPA Index Estimation with Fixed Controller and Bimodal Distribution Noises
Keep all parameters except the distribution of unchanged and the distribution of is chosen as the following bimodal distribution. Parameter estimation results show that both methods seem to be available as system parameter estimates. However, the noise and its estimated distribution shown in Figure 7 indicate that the RELS method is ineffective in estimating the noise of the bimodal distribution.   Figure 7 shows that that the proposed method is also effective for the bimodal non-Gaussian system identification. According to the parameter identification results, the CPA value can be calculated as shown in Table 2. Comparing the data in Tables 1 and 2, the estimated system performance using the entropy index is consistent where the noise obeys above distributions. Simulation results show that the performance assessment method based on minimum entropy can effectively evaluate noise obeying random distribution system and reduce the limitations of the existing evaluation methods.
Although in the framework of different noise and the same controller, the proposed MEC-based CPA index can give a consistent conclusion, is there similar consistency for different controllers? We will discuss this problem next.

CPA Index with Different Gain Controllers
In this simulation, the controller gain K will be changed to get the controller performance trends. K = 0.8 is the ME controller. In this subsection, the method of CPA index estimation in [27] is used here for comparison. And two scales (0.05 and 0.01) [27] are selected to the CPA index estimation as shown in Tables 3 and 4, respectively, and the corresponding CPA index estimation value with the proposed method is given in Table 5. Moreover, the CPA index estimation value with the proposed method for other controllers is given in Table 6.  Tables 3-5 provide performance assessment results of the different controllers. When the system controller selects the optimal controller, the system output entropy can get to the minimum value, which means that the system output y t is equivalent to the system's feedback invariant and the corresponding CPA index is very close to 1. Meanwhile, with the increase of K, the CPA values estimated by the two methods the have the same monotonous decreasing trend whether they are unimodal distribution or bimodal distribution. But it is noticeable in Table 3 that there are certain CPA values (the blue ones) which give wrong evaluation conclusions. This is because this large discretization scale (0.05) is not suitable for this type of bimodal distribution. Furthermore, it is easy to see from Tables 5 and 6 that the proposed CPA index can well reflect the performance of the current controller regardless of whether the control gain increases or the control gain decreases. The above simulation results show that the proposed scheme is more suitable for non-Gaussian systems CPA.

Conclusions
In this paper, by analyzing the limitation of MVC-based CPA method, a new MEC-based CPA method is developed for linear non-Gaussian system and the entropy index is obtained by PDF estimation. Furthermore, the improved EDA algorithm is adopted to estimate the parameters and noise distribution of the system and then calculate the entropy of the feedback invariant. The effectiveness of the CPA procedures is tested by many simulations. It can be concluded that the new MEC-based CPA index can be used to assess the control loop performance of linear systems with non-Gaussian stochastic disturbance. However, the efficiency of the improved EDA algorithm is still an open problem. Further research will focus on optimization algorithms or consider other strategies for parameter identification and PDF estimation.