Load State Identification Method for Ball Mills Based on Improved EWT, Multiscale Fuzzy Entropy and AEPSO_PNN Classification

Abstract

To overcome the difficulty of accurately determining the load state of a wet ball mill during the grinding process, a method of mill load identification based on improved empirical wavelet transform (EWT), multiscale fuzzy entropy (MFE), and adaptive evolution particle swarm optimization probabilistic neural network (AEPSO_PNN) classification is proposed. First, the concept of a sliding frequency window is introduced based on EWT, and the adaptive frequency window EWT algorithm, which is used to decompose the vibration signals recorded under different load states to obtain the intrinsic mode components, is proposed. Second, a correlation coefficient threshold is used to select the sensitive mode components that characterize the state of the original signal for signal reconstruction. Finally, the MFE of the reconstructed signal is used as the characteristic vector to characterize the load state of the mill, and the partial mean value of MFE is calculated. The results show that the mean value of MFE under different load states varies. To further identify the load state, a characteristic mill load vector is constructed from the MFE values of the reconstructed signal under different load conditions and is used as the input of the AEPSO_PNN model, which then outputs the predicted ball mill load state. Thus, a load state identification model is established. The feasibility of the method is verified based on grinding experiments. The results show that the overall recognition rate of the proposed method is as high as 97.3%. Compared with the back propagation (BP) neural network, Bayes discriminant method, and PNN classification, AEPSO_PNN classification increases the overall recognition rate by 8%, 5.3%, and 3.3%, respectively, which indicates that this method can be used to accurately identify the different load states of a ball mill.

1. Introduction

As the main type of mechanical equipment used for ore grinding, ball mills are widely used in the beneficiation process in mining operations [1]. It is imperative but challenging to develop effective modeling, monitoring, and control techniques for complex industrial systems [2,3,4]. Due to their complexity, it is difficult to investigate the internal charge dynamics of ball mills. Energy consumption is obviously related to rotational speed and mill load, and scholars have examined the influence of rotational speed on the energy consumption of mills and achieved good results [5]. For the mill load, it is important to be able to quickly and accurately identify the internal load of a ball mill to ensure that the mill is operating under the best possible working conditions, not only to reduce energy consumption during mineral processing, but also to ensure high grinding efficiency and output [6,7]. Therefore, a method of increasing the load recognition rate for ball mills would have great application value for improving the stability and economic benefits of the grinding process, and efforts to develop such methods have attracted the attention of many scholars at home and abroad [8,9]. To this end, studies have shown that the vibration signal generated by a ball mill during the grinding process is correlated with the load [10].

The vibration signal of a ball mill is nonlinear and nonstationary. Currently, the most widely used methods for processing such signals include the wavelet packet algorithm, empirical mode decomposition (EMD), variable mode decomposition (VMD), local mean decomposition (LMD), and the complete integrated empirical decomposition algorithm (CEEMDAN) [11,12,13,14]. Liu et al. [15] combined the EMD algorithm with principal component analysis (PCA) to extract the vibration signal from the cylinder of a wet ball mill. The results showed that this method can distinguish among different load states, but that the recognition rate requires improvement. Tang et al. [16,17] reported a method of extracting the vibration signal characteristics of a ball mill based on ensemble empirical mode decomposition (EEMD) and interval partial least squares (iPLS) modeling and extended this method to the study of ball mill sound signals. Although the signal features were successfully extracted, there was residual noise in the intrinsic mode functions (IMFs) after decomposition, and white noise with a different amplitude was added each time. Although the above methods can be used to successfully extract signal features, they face problems related to noise residuals and computational burden. Therefore, the key to mill load identification is to find an effective method of extracting the characteristic information of the vibration signal of the ball mill cylinder. The proposed empirical wavelet transform (EWT) algorithm effectively compensates for the above shortcomings. This algorithm not only suppresses the modal aliasing problem and reduces residual noise, but also improves the completeness of decomposition. In reference [18], the EMD, EEMD, and EWT algorithms were compared and analyzed. The EWT algorithm was found to have the best processing effect. Specifically, the EWT algorithm had a better processing speed and better ability to extract modal component signals than the other algorithms. However, in practical engineering, especially under the harsh working conditions of a ball mill, the Fourier spectrum of the EWT segmentation signal easily encounters interference from background noise and must be further improved. In this paper, the adaptive frequency window is used to improve EWT. Compared with traditional EWT and other signal processing algorithms, the denoising effect is more significant.

In recent years, many nonlinear dynamic methods, such as multiscale entropy (MSE), singular value entropy (SVE), permutation entropy (PE), and fuzzy entropy (FE), have been widely used for fault diagnosis, classification, and recognition because of their good performance in terms of feature extraction [19,20,21]. Miao Y et al. [22] applied SVE to the identification of the optimal frequency band. Zhao L et al. [23] completed the fault diagnosis of a gearbox using PE optimization and modified the modal decomposition algorithm. Chang J L et al. [24] applied MSE for load recognition in machine tools. Liu H et al. [25] reported an example of MSE applied for the fault diagnosis of rolling bearings, but the recognition accuracy required further improvement. To diagnose the problem of rolling bearing faults, Zheng H D et al. [26] adopted the method of multiscale fuzzy entropy (MFE), which effectively overcame the defect in the MSE mutation, and the diagnosis result was improved. Compared with the above methods, MFE has some advantages for feature extraction because of its unique performance and ability to accurately reflect the feature information of the original signal.

As a tool for recognition and classification, an artificial neural network is a model abstracted based on neural network theory that originates from the field of physiology. Such models can be used for arbitrary data clustering and pattern classification and are widely used for tasks such as pattern recognition [27,28,29]. Specifically, a probabilistic neural network (PNN) is an artificial neural network with the advantages of a fast training speed, simple parameter adjustment, and good classification performance [30]. However, the classification effect of a probabilistic neural network is greatly influenced by the smoothing parameter σ, and if the selection of σ is not appropriate, then inaccurate results may be obtained. To solve this problem, an adaptive evolutionary particle swarm optimization (AEPSO) algorithm is proposed in this paper to optimize the smoothing parameters in a probabilistic neural network (PNN) so that the optimized network can identify the load state of a ball mill. In this paper, the AEPSO algorithm is used to improve the PNN clustering method; compared with the traditional PNN clustering method and other clustering methods, it has the advantages of high speed and high accuracy.

Considering the nonstationary and nonlinear characteristics of the vibration signal from the cylinder of a ball mill, a load identification method for ball mills is proposed in this study based on improved EWT, MFE, and AEPSO_PNN classification. First, the vibration signals are decomposed using improved EWT, and the mode components of the reconstructed signals are selected using a correlation coefficient threshold. Then, the load state of the ball mill is determined based on the magnitude of the calculated MFE. Finally, AEPSO_PNN is used for learning and classification to enable the recognition of a different load state.

2. Principles of the Load State Identification Method

2.1. Principles of Improved EWT

2.1.1. Principles of EWT

EWT is a widely used method for the adaptive segmentation of signals [31]. The segmentation principle involves adaptively segmenting the Fourier spectrum by marking maximum points in the frequency domain, and a set of bandpass filters suitable for processing signals is constructed in the frequency domain to extract amplitude modulation and frequency modulation (AM-FM) components from the Fourier spectrum.

The Fourier axis [0, π] is divided into n consecutive parts, that is, ${\mathrm{\Lambda }}_n=[{\omega }_{n-1},{\omega }_n]({\omega }_0=0,{\omega }_n=\pi )$, where ${\omega }_n$ is the boundary point between two parts and the corresponding value is the minimum between the two adjacent maximum values in the Fourier spectrum of the signal. Figure 1 [32] shows the division diagram of the Fourier axis. In the figure, ${\omega }_n$ is defined as the center point of ${\mathrm{\Lambda }}_n$. Then, a transition region with a width of $T_n=2{\tau }_n$ is obtained.

Referring to the wavelet construction method of Littlewood–Paley and Meyer, the empirical wavelet function is constructed. After ${\mathrm{\Lambda }}_n$ is determined, the empirical wavelet is used as a bandpass filter. The formulas of the empirical wavelet function ${\stackrel{\wedge }{\psi }}_n(\omega )$ and the empirical scale function ${\stackrel{\wedge }{\varphi }}_n(\omega )$ are as follows [33]:

$${\stackrel{\wedge }{\psi }}_n(\omega )=\{\begin{array}{l}1,(\left|\omega \right|\le (1-\gamma ){\omega }_n)\\ \cos \left\{\frac{\pi }{2}\beta \left[\frac{1}{2\gamma {\omega }_n}(\left|\omega \right|-(1-\gamma ){\omega }_n)\right]\right\},(1-\gamma ){\omega }_n\le \omega \le (1+\gamma ){\omega }_{\mathrm{n}}\\ 0,(\mathrm{others})\end{array}$$

(1)

$${\stackrel{\wedge }{\varphi }}_n(\omega )=\{\begin{array}{l}1,(\left|\omega \right|\le (1-\gamma ){\omega }_n)\\ \cos \left\{\frac{\pi }{2}\beta \left[\frac{1}{2\gamma {\omega }_n}(\left|\omega \right|-(1-\gamma ){\omega }_n)\right]\right\},(1-\gamma ){\omega }_n\le \omega \le (1+\gamma ){\omega }_{\mathrm{n}}\\ 0,(\mathrm{others})\end{array}$$

(2)

where

$$\begin{array}{l}\beta \left(x\right)=x^4(35-84x+70x^2-20x^3);\\ {\tau }_n=\gamma {\omega }_n\gamma <{\min }_n(\frac{{\omega }_{n+1}-{\omega }_n}{{\omega }_{n+1}+{\omega }_n})\end{array}$$

(3)

After the EWT, the approximation coefficient $W_{f1}(0,t)$ and the detail coefficient $W_{f2}(n,t)$ can be expressed as follows.

$$\begin{array}{ll}{W}_{f1}(0,t)& =<x,{\varphi }_1>={\displaystyle \int x(\tau )}\overline{{\varphi }_1(\tau -t)}\mathrm{d}\tau \\ & =F^{-1}[x(\omega ){\stackrel{⌢}{\varphi }}_1(\omega )]\end{array}$$

(4)

$$\begin{array}{ll}{W}_{f2}(n,t)& =<x,{\psi }_n>={\displaystyle \int x(\tau )}\overline{{\psi }_n(\tau -t)}\mathrm{d}\tau \\ & =F^{-1}[x(\omega ){\stackrel{⌢}{\psi }}_n(\omega )]\end{array}$$

(5)

Then, the functional expression of the reconstructed original signal a is as follows:

$$\begin{array}{l}x(t)=W_{f1}(0,t)\ast {\phi }_1(t)+{\displaystyle \sum _{n=1}^N{W}_{f2}(n,t)\ast {\psi }_n(t)}\\ =F^{-1}\left[\stackrel{\wedge }{W_f}(n,\omega ){\phi }_1(\omega )+{\displaystyle \sum _{n=1}^N\stackrel{\wedge }{W_f}(n,\omega )\stackrel{\wedge }{{\psi }_n}(\omega )}\right]\end{array}$$

(6)

where “$\ast$” is a convolution operation and $\stackrel{\wedge }{W_{f1}}(0,\omega )$ and $\stackrel{\wedge }{W_{f2}}(n,\omega )$ are the Fourier transforms of the approximate coefficient $W_{f1}(0,t)$ and the detail coefficient $W_{f2}(n,t)$, respectively. Finally, the signal f is decomposed into the sum of several single component signals.

$$x(t)={\displaystyle \sum _{k=0}^{N-1}{x}_k(t)}$$

(7)

2.1.2. Principle of the Adaptive Frequency Window EWT Algorithm

The division rules of the spectral boundaries of the traditional EWT algorithm are determined by the frequency domain extreme points, but a ball mill is vulnerable to strong noise, resulting in the disorderly arrangement of frequency domain extreme points. Considering these deficiencies, this paper uses the adaptive frequency window EWT to divide the spectral boundaries, as shown in Figure 2.

In Figure 2, the frequency window is represented as $\left[{\omega }_a,{\omega }_b\right]$, where ${\omega }_a,{\omega }_b$ is the central frequency of the lower cutoff band of the window. The shaded area represents the transition region of the segmented portion of the spectrum with width 2$\mathrm{\tau }$. The range of the d support interval is $\left[0,\pi \right]$. The frequency window can slide freely in the interval, and the width range is adaptively variable.

After the frequency window segmentation is improved, Equations (1) and (5) are modified as follows.

$${\stackrel{\wedge }{\psi }}_n(\omega )=\{\begin{array}{l}1,({\omega }_a+\tau \le \left|\omega \right|\le {\omega }_b-\tau )\\ \cos \left\{\frac{\pi }{2}\beta \left[\frac{1}{2\tau }(\left|\omega \right|-{\omega }_b+\tau )\right]\right\},({\omega }_b-\tau \le \left|\omega \right|\le {\omega }_b+\tau )\\ \sin \left\{\frac{\pi }{2}\beta \left[\frac{1}{2\tau }(\left|\omega \right|-{\omega }_a+\tau )\right]\right\},({\omega }_a-\tau \le \left|\omega \right|\le {\omega }_a+\tau )\\ 0,(\mathrm{others})\end{array}$$

(8)

$$\begin{array}{ll}{W}^{\prime }(t)& =<x,\psi >={\displaystyle \int x(\tau )}\overline{\psi (\tau -t)}\mathrm{d}\tau \\ & =F^{-1}[x(\omega )\stackrel{⌢}{\psi }(\omega )]\end{array}$$

(9)

Additionally, Equation (8) must be satisfied as follows.

$$\{\begin{array}{l}\beta (x)=x^4(35-84x+70x^2-20x^3)\\ \tau =\gamma {\omega }_a\gamma <({\omega }_b-{\omega }_a)/({\omega }_b+{\omega }_a)\end{array}$$

(10)

Therefore, the modal component signal can be reconstructed as follows:

$$x^{*}(t)=W^{\prime }(t)\ast \psi (\omega )=F^{-1}\left[\stackrel{\wedge }{W}(\omega )\stackrel{\wedge }{\psi }(\omega )\right]$$

(11)

where $\stackrel{\wedge }{\psi }(\omega )$ is a Fourier transform of $\psi (\omega )$ and $x^{*}(t)$ is an AM–FM component signal for improving EWT extraction.

2.1.3. Simulation and Comparative Analysis of Improved EWT

To verify the ability of the improved EWT method to extract the feature components of the signal, a simulation with the improved EWT approach is performed, and the results are compared with those of the traditional EWT. The simulation signal x (t) is constructed as follows in Equation (12):

$$\{\begin{array}{c}{x}_1(t)=2t^2\\ x_2(t)=1.1\sin (34\pi t)\\ x_3(t)=\{\begin{array}{c}0.7\cos (56\pi t)\begin{array}{cc}& 0<t<0.5\end{array}\\ 0.8\cos (64\pi t)\begin{array}{cc}& t\ge 0.5\end{array}\end{array}\\ x(t)=x_1(t)+x_2(t)+x_3(t)\end{array}$$

(12)

where x (t) is white noise, the signal-to-noise ratio (SNR) is set to 3, and $t\in [0,1]$. Figure 3 shows the improved EWT and traditional EWT decomposition results for the simulation signal x (t).

In Figure 3, the components f2–f5 correspond to the signals x₃(t)~x₁(t), respectively. Figure 3a shows that the noise contained in the signal is well decomposed by the improved EWT and that the degree of coincidence of each component is close to 90%. The two modes that originally belonged to the same component are decomposed because the two modal components obviously have distinct energy signals and can be regarded as two independent modes. Figure 3b shows that traditional EWT can decompose noise, but the components x₁(t), x₂(t), and x₃(t) are deformed because the traditional EWT segmentation method is too simple. When analyzing local noise or nonstationary signals, some local maxima generated by noise and nonstationary components may appear and erroneously remain in the peak sequence, and some useful maxima may not be kept in the peak sequence, resulting in improper segmentation. The improved EWT uses the adaptive frequency window for spectrum segmentation, which can reduce the effects of noise and nonstationary components and greatly increase the reliability of spectrum segmentation.

This comparative study of simulated signals indicates that the improved EWT method can effectively detect the modal components in power spectra, extract components similar to the original signal components, and suppress modal aliasing. Thus, the decomposition effect of the improved EWT method is better than that of the traditional EWT method.

2.2. Principle of Multiscale Fuzzy Entropy

2.2.1. Principle of Fuzzy Entropy

FE is the probability of identifying a new pattern in a time series when the dimension changes, which reflects the complexity and irregularity of the time series. The larger the probability of the time series, the greater the FE value [34]. During the operation of a ball mill, the change in the load state will cause the characteristics of the vibration signal of a cylinder to change in an obvious manner, and FE can effectively characterize the state characteristics of the signal in each frequency band during the sampling time. Therefore, it is feasible to introduce FE as the characteristic parameter of the vibration signal of a ball mill cylinder. The algorithm steps are as follows.

1. The m-dimensional vector is obtained by processing the time series:

$$\begin{array}{l}{X}_i^m=\{u(i),u(i+1),\dots ,u(i+m+1)\}-u_0(i)\\ u_0(i)=\frac{1}{m}{\displaystyle \sum _{j=0}^{m-1}u(i+j)}i=1,2,\dots ,i+m+1\end{array}$$

(13)

where $X_i^m$ is the result of removing the mean $u_0(i)$ of the time series.

2. Calculate the maximum distance between $X_i^m$ and $X_j^m$:

$$d_{ij}^m=d[X_i^m,X_j^m]=\underset{k\in (0,m-1)}{\max }\{|u(i+k)-u0(i)-(u(j+k)-u0(j))|\}$$

(14)

where $i,j=1,2,\cdot \cdot \cdot ,N-m,i\ne j$.

3. The similarity between $X_i^m$ and $X_j^m$ is defined by a fuzzy function as follows:

$$D_{ij}^m=u(d_{ij}^m,n,r)=e^{-{(d_{ij}^m/t)}^n}$$

(15)

where $u(d_{ij}^m,n,r)$ is an exponential fuzzy membership function and n and r are the boundary gradients and widths of the fuzzy membership functions, respectively.

4. Define the functions as follows:

$${\varphi }^m(n,r)=\frac{1}{N-m}{\displaystyle \sum _{i=1}^{N-m}(\frac{1}{N-m-1}{\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^{N-m}{D}_{ij}^m})}$$

(16)

where $i,j=1,2,\cdot \cdot \cdot ,N-m,i\ne j$.

5. The m+1 vector is constructed based on the above four steps.

$${\varphi }^{m+1}(n,r)=\frac{1}{N-m}{\displaystyle \sum _{i=1}^{N-m}(\frac{1}{N-m-1}{\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^{N-m}{D}_{ij}^{m+1}})}$$

(17)

6. The calculation formula of the FE value can be summarized as follows:

$$FuzzyEn(m,n,r)=\underset{N\to \infty }{\lim }[\ln {\varphi }^m(n,r)-\ln {\varphi }^{m+1}(n,r)]$$

(18)

where $i,j=1,2,\cdot \cdot \cdot ,N-m,i\ne j$.

7. When N is limited, Equation (18) is transformed into the following formula.

$$FuzzyEn(m,n,r,N)=\ln {\varphi }^m(n,r)-\ln {\varphi }^{m+1}(n,r)$$

(19)

2.2.2. Principle of Multiscale Fuzzy Entropy

The characteristic frequency band and complexity of the vibration signal of a cylinder under different load conditions in a ball mill are different at different scales. Considering the FE of the vibration signal at different scales can improve the recognition accuracy, therefore the concept of multiple scales is introduced based on FE. The steps in the MFE algorithm are as follows.

1. Construct a new coarse granularity vector for the original time series $X_i=\{x_1,x_2,\cdot \cdot \cdot ,x_n\}$ as follows:

$$y_j(\tau )=\frac{1}{\tau }{\displaystyle \sum _{i=(j-1)\tau +1}^{j\tau }{x}_i}1\le j\le \frac{N}{\tau }$$

(20)

where $\tau =1,2,\cdot \cdot \cdot ,n$ represents the scale factor. When $\tau =1$, the coarse-grain vector is the original sequence. For a given $\tau$, the original sequence is divided into coarse granularity vectors of length $N/\tau$, and Figure 4 shows the coarse granularity process for $\tau =2$ and $\tau =3$.

2. The FuzzyEn of each coarse-grained sequence is determined by the standard deviation of the original sequence. The FuzzyEn value can be expressed as a function of the scale factor in MFE analysis.

2.2.3. Parameter Selection for MFE

According to the definition of MFE, the calculation of MFE is related to the embedding dimension m, similarity tolerance r, exponential function gradient n, and data length N. The selection rules are as follows.

• A large embedding dimension m produces more information when the time series is dynamically reconstructed, and the data sequence $N={10}^m\sim {30}^m$; thus, m is set to 2.
• The similarity tolerance r represents the width of the boundary of the exponential function. If r is too large, then a large amount of statistical information will be lost, and if r is too small, then the sensitivity to noise will be high. r is usually set from 0.1 SD to 0.25 SD (SD denotes the standard deviation of the original time series). Considering the working characteristics of the ball mill, r is set to 0.15 SD.
• n is a weighting factor in the calculation of FE vector similarity. A large n will result in a large gradient, but an overly small n will lead to the loss of detail. To obtain as much detailed information as possible, a small integer is usually used, and n is set to 2 in this case.
• To obtain an accurate MFE calculation result, the data length N should be greater than $100{\tau }_{\max }$. In addition, the maximum scale factor ${\tau }_{\max }$ should also be considered when calculating the MFE, and the value of ${\tau }_{\max }$ is usually between 10 and 20; thus, a = 20 is used in this study.

2.3. Principle of the AEPSO_PNN

2.3.1. PNN Principle

A PNN is a type of radial basis network that was first proposed by Dr. D.F. Speeht in 1989. The PNN is a supervised network classifier based on the Bayes minimum risk criterion [35]. As a feed-forward network, a PNN has the advantages of a fast training speed and simple parameter adjustment. Currently, PNNs are widely used in pattern classification [36]. Compared with other network classifiers, a PNN can not only guarantee real-time performance, but also produce classification and recognition results that are minimally influenced by complex parameter settings.

The signal sample vector can be represented as $X=\left[x_1,x_2,\dots ,x_i,\dots ,x_n\right]$ with states $Y=\left[y_1,y_2,\dots ,y_i,\dots ,y_n\right]$. Then, the prior probability, posterior probability, and class-specific probability density functions for each state can be represented by $P(y_i)$, $P(y_i/X)$, and $P(X/y_i)$, respectively. For a given identification target, $P(y_i)$ is a known parameter, and $P(X/y_i)$ can be estimated using the Parzen function. The corresponding formula is as follows:

$$P(X/y_i)=\frac{{\displaystyle \sum _{j=1}^{N_i}\exp (-\frac{{\Vert X-x_{ij}\Vert }^2}{2{\sigma }^2})}}{N_i{(2\pi )}^{\frac{d}{2}}{\sigma }^d}$$

(21)

where N_i is the total number of samples of the ith class, d is the dimensionality of the feature vector, x_ij is the jth sample of the ith class, and σ is the width of the Parzen function window, that is, the smoothing parameter.

The following formula is obtained from probabilistic and statistical theory.

$$P(y_i/X)=P(X/y_i)P(y_i)/P(X)$$

(22)

If the possibility of misjudgment is not considered, the Bayes rule can be expressed as follows.

$$\vee j\ne i=1,2,3,\cdot \cdot \cdot ,m,ifP(y_i/X)>P(y_j/X),X\in y_j$$

(23)

However, because misjudgment can readily occur in real-world situations, it is necessary to introduce the risk coefficient λ_ij, yielding the following risk function R for the decision conditions.

$$R(y_i/X)={\displaystyle \sum _{j=1}^m{\lambda }_{ij}}P(y_j/X)$$

(24)

In summary, the Bayes minimum risk criterion can be expressed as follows.

$$ifR(y_i/X)>R(y_j/X),X\in y_j$$

(25)

In this paper, the minimum risk criterion is used as the basis of the feed-forward network that serves as the mill load state recognition model. By setting reasonable smoothing parameters, the network is trained on a set of sample feature vectors to estimate the probability densities of three distinct load states and enable the recognition of the mill load state.

2.3.2. Principle of AEPSO

The optimization speed and position updating formulas of the traditional particle swarm optimization algorithm [37] are as follows:

$$V_p{}^{k+1}=\omega V_p^k+c_1{r}_1(W_p-X_p^k)+c_2{r}_2(W_g-X_p^k)$$

(26)

$$X_p{}^{k+1}=X_p^k+V_p^{k+1}$$

(27)

where k is the number of iterations; $\omega$ is the inertial weight of the particle; $c_1,c_2$ are the learning factors of the particle, of which the former is the individual factor and the latter is the global factor; and $r_1,r_2$ are random numbers in the interval [0, 1], which make the particles independent and diverse.

To address the nonlinear problem of ball mill load identification, the AEPSO algorithm introduces a nonlinear adaptive time-varying inertial weight.

$${\omega }_t={\omega }_{start}-({\omega }_{start}-{\omega }_{end})\times \exp (-\frac{1}{1+2t/t_{\max }})$$

(28)

For the learning factors $c_1,c_2$ of particles, the traditional particle colony algorithm usually sets $c_1=c_2=2$, but this approach ignores the phase difference of the algorithm during training. The AEPSO algorithm adopts the strategy of managing the learning factor in segment, and the formula is as follows.

$$\{\begin{array}{c}{c}_1=2.5,c_2=1.5t<t_{\max /2}\\ c_1=1.5,c_2=2.5t\ge t_{\max /2}\end{array}$$

(29)

To enhance the adaptability of particle swarm optimization after iteration, the AEPSO algorithm introduces the local search operator $\beta$ in Equation (13). The revised formula is as follows:

$$X_p{}^{k+1}=X_p^k+\beta \times V_p^{k+1}$$

(30)

where $\beta =rand()[rand()+0.5]$ and $rand()$ is a random number in [0, 1].

2.3.3. Optimization of the PNN by AEPSO

The smoothing parameter σ in the PNN has a considerable influence on the training effect. The improper selection of the σ value makes it easy to misjudge the recognition of the mill load state. Therefore, this paper uses an AEPSO algorithm to optimize the smoothing parameters of the PNN so that the optimized network can effectively identify the state of the mill load. The specific steps in the algorithm are as follows.

• The parameters of the PSO algorithm are initialized, the smoothing parameters σ of the PNN are used as the population particles, the number of iterations is set to 500, and a set of data (σ) is randomly generated as an initial parameter vector.
• The training set samples are input, and the fitness function is used to calculate the fitness value. Then, the optimal individual fitness value and the global optimal fitness value of the group are traversed by comparing each particle (σ). Finally, the particles are adjusted.
• After calculating each particle in the population, the termination condition is determined to be satisfied or not. If not, the state is updated according to the speed and position updating formula; then, the algorithm returns to step 2. Otherwise, the algorithm iterates until termination and outputs the search results.
• The PNN model trained by the optimal parameter combination (σ) is used to classify the test sample set and output the target category.

The network structure of the AEPSO_PNN includes four parts: the input layer, the mode layer, the summation layer, and the output layer, as shown in Figure 4.

As Figure 5 shows, the training step of the load state identification model of a ball mill based on the AEPSO_PNN is as follows.

• The input layer multiplies the received feature vector of the training sample by the weighting coefficient Wj and transmits the result to the mode layer for training. The number of neurons in this layer is the dimension of the feature vector.
• The mode layer first uses the exponential function gj as the activation function. Then, the probability density of each neuron is determined, and finally, the result is transmitted to the summation layer.
• The probability density is the weighted average of the summation, and the resulting estimated probability density is transmitted to the output layer.
• Based on the Bayes minimum risk criterion, the output layer selects the category with the largest posterior probability as the final classification result of the sample.

3. Design of the Load State Identification Method for a Ball Mill

Based on the research on the EWT algorithm, MFE theory, and the PNN clustering algorithm combined with the characteristics of ball mill vibration signals, a feature extraction algorithm for vibration signals is proposed based on modified EWT, MFE, and AEPSO_PNN classification. The specific steps in the algorithm are as follows.

1. Decompose the recorded vibration signal from the cylinder of the ball mill via the adaptive frequency window EWT algorithm to obtain ${\mathrm{AM}-\mathrm{FM}}_i(i=1,2,\cdots ,n)$.

2. Calculate the correlation coefficients for all ${\mathrm{AM}-\mathrm{FM}}_i$ components and the original signal in accordance with Equation (31), and select the sensitive $\mathrm{AM}-\mathrm{FM}$ components based on the threshold given in Equation (32).

$${\mathrm{\rho }}_{xy}=\frac{{\displaystyle \sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y}})}{\sqrt{{\displaystyle \sum _{i=1}^N{(x_i-\overline{x})}^2}}\sqrt{{\displaystyle \sum _{i=1}^N{(y_i-\overline{y})}^2}}}$$

(31)

The correlation coefficient threshold is calculated as

$${\mu }_h=\frac{\max ({\mu }_i)}{10\times \max ({\mu }_i)-3}$$

(32)

where ${\mu }_h$ is the threshold, ${\mu }_i$ is the correlation coefficient between the ith AM-FM component and the original signal, and max is the maximum correlation coefficient value. Each AM-FM component for which the value of the correlation coefficient with the original signal is greater than the threshold ${\mu }_h$ is retained as a sensitive AM-FM component. Each AM-FM component for which the correlation coefficient is smaller than the threshold ${\mu }_h$ is removed as a spurious component.

3. The sensitive AM-FM components are used to obtain the reconstructed vibration signals of different loads in the ball mill.

4. The MFE of the reconstructed vibration signal is calculated, and the result is used as the characteristic vector for the load classification of the ball mill.

5. The characteristic vector matrix is used as the input of AEPSO_PNN, and the load state is used as the output. Then, the load state of the mill is identified.

Thus, the overall flow of the ball mill load identification model that is proposed in this paper based on the modified EWT, MFE, and AEPSO_PNN classification methods can be summarized as shown in Figure 6.

4. Experimental Analysis of Mill Load State Recognition

4.1. Data Collection

To verify the method proposed in this paper, a grinding experiment was performed using a Φ305 × 305 mm Bond index experimental ball mill. The experimental device is shown in Figure 7. The material used in the experiment was tungsten ore from a mine in Jiangxi, China, with a Protodyakonov scale of hardness of 14–18, a density of 1.8 t/m³, and five grades of particle sizes: 1–3 mm, 3–6 mm, 6–9 mm, 9–11 mm, and >11 mm. The experimental parameters considered were the fill rate, powder-to-ball ratio, and grinding concentration. The vibration signal acquisition system of the mill cylinder consisted of a DH5922N dynamic data acquisition instrument and a DH131 acceleration sensor, which were used to collect the signals of various load parameters under three different load conditions. According to the literature, the mill load was divided into the following states: the underloaded state, corresponding to a fill rate of 10–20%; the normal load state, corresponding to a fill rate of 20–40%; and the overloaded state, corresponding to a fill rate of 40–60% [38].

4.2. Decomposition of the Cylinder Vibration Signal

First, we present the typical working conditions corresponding to the three load conditions considered in this analysis: working condition 1 (a fill rate of 10%, a powder-to-ball ratio of 0.4, and a grinding concentration of 0.5), working condition 2 (a fill rate of 30%, a powder-to-ball ratio of 0.6, and a grinding concentration of 0.5), and working condition 3 (a fill rate of 50%, a powder-to-ball ratio of 0.8, and a grinding concentration of 0.5). The waveforms of the cylinder vibration signals recorded under these three working conditions are shown in Figure 8.

As Figure 8 shows, there is a large amount of noise in the vibration signal from the mill cylinder in all three load states, which makes it difficult to effectively extract feature information. To extract the characteristics of the vibration signal of the cylinder, the original signal must be preprocessed. The preprocessing steps are as follows.

1. The improved EWT algorithm was used to adaptively decompose the original signals under three typical working conditions. Then, 10 AM-FM components were obtained.

2. The correlation coefficients between the AM-FM components and the original cylinder vibration signal were calculated using Equation (29), and the threshold values were then calculated in accordance with Equation (30), yielding the following results: 0.21437 for working condition 1, 0.23872 for working condition 2, and 0.19905 for working condition 3. The correlation coefficient values and the threshold values of the vibration signals from the cylinder body under the three working conditions are shown in Figure 9.

As shown in Figure 9, the correlation coefficients between the AM-FM1, AM-FM2, and AM-FM5 components and the original signal for working condition 1 were greater than the threshold value of 0.21437. Thus, these components were identified as sensitive AM-FM components that characterize the vibration signal of the cylinder. For working condition 2, the AM-FM2, AM-FM4, and AM-FM6 components, with correlation coefficients greater than the threshold value of 0.23872, were selected as the sensitive mode components. For working condition 3, the AM-FM1, AM-FM3, and AM-FM4 components, with correlation coefficients greater than the threshold value of 0.23872, were selected as the sensitive mode components. All AM-FM components with correlation coefficients smaller than the corresponding threshold were removed.

3. The selected sensitive modal components are reconstructed, and the results are shown in Figure 10.

Based on a comparison of Figure 8 and Figure 10, the trend of the reconstructed signal waveform is basically the same as that of the original signal. Compared with the original signal, the impact profile of the reconstructed signal curve is obviously distinct, but it preserves the characteristic information of the original signal while effectively denoising the signal. To further quantitatively highlight the preprocessing effect in this paper, the EMD algorithm, EWT algorithm, and improved EWT algorithm are used to decompose the original signals of the three working conditions, and the sensitive components are reconstructed by the correlation coefficient method. Additionally, the SNR is introduced into the comparative analysis before and after processing to qualitatively analyze the comparison results, and the results are shown in

Table 1. Signal-to-noise ratio (SNR) before and after signal processing.

Working Conditions	The Original Signal (SNR/db)	Reconstructed Signals of Three Algorithms (SNR/db)
		EMD	EWT	Improved EWT
1	7.91	13.97	17.22	21.23
2	9.58	15.35	18.94	22.36
3	7.02	14.61	19.07	24.54

.

As Table 1 shows, compared with the original signal, the SNR of the reconstructed signal processed by the improved EWT algorithm increases by 13.32 dB, 12.78 dB, and 17.52 dB under three typical working conditions, which indicates that the noise is considerably reduced after applying the improved EWT algorithm. Compared with those of the EMD algorithm and the EWT algorithm, the SNR of the reconstructed signal processed by the improved EWT algorithm increases the most. Therefore, the preprocessing effect of the improved EWT algorithm is best.

4.3. Decomposition of the Cylinder Vibration Signal

The FE of the reconstructed signal is calculated, and five groups of samples are assessed for each type of ball mill load state. The average value of the FE of the three-state data is calculated, as shown in

Table 2. Fuzzy entropy values of three types of load state vibration signals.

Sample	Underloaded	Normal Load	Overloaded
1	1.19	1.01	0.45
2	1.31	0.88	0.59
3	1.03	0.92	0.45
4	1.42	0.73	0.38
5	1.30	1.11	0.57
Mean	1.25	0.93	0.48

.

Table 2 shows that the FE value of the vibration signal varies by load state and that the FE value of the vibration signal under the same load state fluctuates back and forth near the average value. By comparing the FE values of three different load vibration signals, the FE values of the underloaded state are found to be relatively large, which is due to the relatively small amount of steel ball and mineral material in the cylinder under this condition. Additionally, the collision frequency between the mineral material and the steel ball in grinding production increases with the movement of the cylinder body to a certain height under the action of friction, and the collision frequency with other steel balls, minerals, and the cylinder walls is high in the process of falling. Energy is mainly consumed in the collisions between the steel ball and the tube wall and between the steel ball and other steel balls; thus, the vibration signal is complex, and the signal is highly random. However, the FE value under overloaded conditions is relatively small because there are more steel balls and minerals in the cylinder under these conditions, causing the steel ball and minerals to undergo peristalsis in the grinding process. In this case, the randomness of the signal is small. Under a normal load, energy is mainly used for grinding the quantity of minerals, and so the complexity of generating a signal is relatively moderate. For underloaded conditions and a normal load state, the sample entropy values are similar, and individual overlap occurs, which results in a discriminating effect. Therefore, MSE is introduced into the analysis of the mill vibration signal. The MFE of the reconstructed signal that can characterize the characteristic information of the vibration signal under three different load conditions is calculated. To highlight the superiority of the feature extraction method used in this paper, four combination methods (EWT-MSE, EWT-MFE, improved EWT-MSE, and improved EWT-MFE) are used to analyze the vibration signals of the cylinder of the ball mill under three load conditions. The mean value and standard deviation curve of the three states (20 samples per group) are shown in Figure 11. The parameter selection process of the algorithm is as described above.

In Figure 11, it is evident that the order of the mean value of FE of the vibration signal of the ball mill cylinder under three working conditions displays the following order: underloaded > normal load > overloaded. Specifically, as the ball mill load increases, the amplitude of each component of the vibration signal in the spectrum obviously increases, which leads to a decrease in entropy. Although the variation trend of the FE of the cylinder vibration signal with the scale factor is the same in different load states, the fluctuation range of the entropy value varies, which indicates that FE can be effectively used to identify the load state. By comparing the four graphs, we see that there are obvious fluctuations and interval intersections between the EWT-MSE method and the EWT-MFE method. Although the entropy curve of the improved EWT-MSE method is smooth and the three states are distinguished to a certain extent, there are still overlap and intersection issues at small scales, which may lead to judgment errors. However, the MSE curves of the three load states obtained with the improved EWT-MFE method are smooth, and the fluctuation intervals have obvious limits. This finding indicates that the improved EWT-MFE method can effectively distinguish among the three load states of the ball mill.

4.4. Training and Testing

To verify the effectiveness of the proposed load identification model for a ball mill, 3 × 100 samples were randomly selected from each of the three classes of vibration signals, including 150 as training samples and 150 as test samples. The selected samples were first decomposed via the improved EWT method. The sensitive mode component signal with load state information was screened by the correlation coefficient method and reconstructed. Then, the MFE of the reconstructed signal was normalized as the input of the load state identification model of the ball mill based on AEPSO_PNN, and the load state of the ball mill was output. To highlight the superiority of AEPSO_PNN classification and identification, three clustering methods, namely, PNN classification, back propagation (BP) neural network, and Bayes identification methods, were trained and tested with the abovementioned samples. Then, the identification results were compared with the AEPSO_PNN identification results. For simplicity of description, the underloaded, normal load, and overloaded conditions are indicated by working condition numbers 1, 2 and 3, respectively. The identification effects of various classification methods are shown in Figure 12 and

Table 3. Singular value entropy results for the three working conditions. BP: back propagation. probabilistic neural network (PNN): probabilistic neural network.

Classification Method	Correct Identifications	Load Recognition Accuracy
BP neural network	134	89.3%
Bayes identification method	138	92.0%
P NN classification	141	94.0%
AEPSO_PNN classification	146	97.3%

.

Figure 12 and Table 3 show that the predicted load state of the AEPSO_PNN model of ball mill load state recognition is largely consistent with the real state. Only four samples are misdiagnosed, and the overall recognition accuracy is 97.3%. Specifically, the recognition accuracy of AEPSO_PNN classification for three different load states is 96%, 98%, and 98%, all of which are high recognition levels. The BP neural network, Bayes discriminant method, and PNN classification can also achieve effective load identification. The highest accuracies of these methods are 89.3%, 92.0%, and 94.0%. Compared with the back propagation (BP) neural network, Bayes discriminant method, and PNN classification, AEPSO_PNN classification increases the overall recognition rate by 8%, 5.3%, and 3.3%. The results show that the mill load identification method based on the improved EWT-MFE method and AEPSO_PNN classification is effective, and the identification effect is excellent. Thus, this method provides a new approach for ball mill load identification.

5. Conclusions

By combining the improved EWT algorithm, MFE feature extraction, and AEPSO_PNN clustering, a load identification model of a ball mill is constructed. The main contributions to this work are as follows:

(1)

The strong background noise, nonlinearity, and nonstationarity of the vibration signal of a ball mill cylinder hinder the recognition accuracy. The improved EWT algorithm proposed in this paper can effectively denoise the original signal and retain the feature information.

(2)

The MFE algorithm has obvious advantages in terms of feature extraction. Notably, the MFE difference between underloaded, normal load, and overloaded conditions is large, and the proposed method can distinguish among the load states of the mill.

(3)

The AEPSO_PNN classifier is introduced into the load recognition model of the ball mill to improve the recognition effect. Compared with the BP neural network, the Bayes discriminant method, and PNN classification, AEPSO_PNN classification provides a better recognition effect and the highest load recognition accuracy.

(4)

The effectiveness of the method is verified based on a grinding experiment performed with a Bond work index ball mill in the laboratory.

In future research, the algorithm, structure, and parameter setting process of the proposed model will be optimized and improved to enhance the ability of the model to identify the ball mill load state.