- Email: [email protected]
Visualization of the controller states of an autogenous mill from time series data
Minerals Engineering 56 (2014) 1–9
Contents lists available at ScienceDirect
Minerals Engineering journal homepage: www.elsevier.com/locate/mineng
Visualization of the controller states of an autogenous mill from time series data C. Aldrich a,b,⇑, J.J. Burchell b, J.W. de V. Groenewald c, C. Yzelle b a
Department of Minerals Engineering and Extractive Metallurgy, Western Australian School of Mines, Curtin University of Technology, Perth, WA, Australia Department of Process Engineering, University of Stellenbosch, Private Bag X1, Matieland 7602, Stellenbosch, South Africa c Anglo Platinum Management Services, PO Box 62179, Marshalltown 2107, South Africa b
a r t i c l e
i n f o
Article history: Received 3 June 2013 Accepted 14 October 2013 Available online 14 November 2013 Keywords: Comminution Modeling Nonlinear time series analysis Neural networks
a b s t r a c t The operational variables of an industrial autogenous mill were embedded in a low-dimensional phase space to facilitate visualization of the dynamic behavior of the mill. This was accomplished by use of a multivariate extension of the method of delay coordinates used in nonlinear time series analysis. In this phase space, the controlled states of the mill could be visualized as separate regions or clusters in the phase space. Comparison of the correlation dimension of the state variable of the mill (the load) embedded in phase space suggested that the dynamic behavior of the mill could not be represented by a linear stochastic model (Gaussian or otherwise). The low dimensionality (62) of the correlation dimension further suggested that the mill load depended on a few variables only and that the underlying generative process had a significant deterministic component. In addition, the operational variables could be used as reliable predictors in a neural network model to identify the controlled states of the mill. As a complementary approach to visualization of the operation of the mill, a different neural network model could be used to reconstruct a corrected power load curve by compensating for the effect of varying operating conditions. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Energy consumption in ore milling processes constitutes up to 40% of the energy expenditure in concentrator plants, but is not more than 15–25% efficient, depending on the approach used to calculate the energy required to reduce the size of the ore particles (Fuerstenau and Abouzeid, 2002). As a consequence, the reduction of energy consumption in comminution circuits has drawn considerable interest over an extended period of time. Energy efficient operation of mills is known to depend to a large extent on the behavior of the mill load (Van Nierop and Moys, 2001), but reliably accounting for the effect of load on mill power consumption is not an easy task. Other process variables, such as media shape (Lameck et al., 2006), mill speed (Powell et al., 2009), ore type (Powell et al., 2009) and mill filling (Kiangi and Moys, 2008) also play a significant role that may be difficult to quantify from first principles. This is especially the case with autogenous and semi-autogenous grinding (SAG) mills, which are popular in new plants, owing ⇑ Corresponding author at: Department of Minerals Engineering and Extractive Metallurgy, Western Australian School of Mines, Curtin University of Technology, Perth, WA, Australia. Tel.: +61 8 9266 4349; fax: +61 8 9358 4912. E-mail address: [email protected] (C. Aldrich). 0892-6875/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.mineng.2013.10.018
to their favorable capital and operating costs. SAG mill operation, in particular, is strongly affected by the type and size of the ore being ground, which makes stable control difficult. In these mills, grinding is driven by the mill load (Powell et al., 2009). One approach to mill control is based on the use of power–load curves (Powell et al., 2001, 2009). This approach involves establishing empirical grind curves that reveal at which mill load the circuit’s power draw peak. These curves have a typical concave parabolic shape when plotted with power and load as the ordinate and abscissa respectively. Once these curves are established, they may be used to guide milling operation with respect to achieving a desired throughput or grind quality. In this context, for example, current control strategies are aimed at the maintenance of maximal mill power draw at or near the top of the power load curve (Powell and Mainza, 2006) by manipulation of the mill feed rate and other process variables, such as the coarse to fine ore feed ratio and the water inlet rate. This aim is challenging, since the power load curve is not static for a particular milling circuit. Instead, the curve is affected by other process variables, such as the ore characteristics, and it may therefore be difficult to track the designated operating point with sufficient accuracy. Therefore, accounting for the effect of multiple variables on the controlled states of the mill would be a major advantage.
2
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9
Nomenclature Symbol Meaning a, b, c, d, h, # threshold values used by expert control system (–) e0 scaling parameter (–) K diagonal eigenvalue matrix (–) ti , tj ith and jth point of a time series embedded in phase space (–) x frequency (–) C correlation matrix (–) CN correlation count (–) dc correlation dimension (–) f ðxÞ Fourier transform (–) I(A,B) mutual information of two random variables A and B (–) IðA; BÞ average mutual information of two random variables A and B (–) k embedding lag of phase space (–) kj embedding lag of jth variable in phase space (–) Lnorm normalized load (–) Lavg average load (ton) Lmax maximum load (ton) Lmin minimum load (ton) LROC load rate of change (ton h1) m number of variables (–) M embedding dimension of phase space (–) Mj embedding dimension of j’th variable in phase space (–) N number of samples in a time series (–)
In this paper, analysis of the dynamic behavior of an autogenous mill is considered, based on reconstruction of the mill variables in phase space, which allows time series variables and their history to be considered simultaneously within a state space context. It is shown that most of these states can be predicted from the operational variables of the mill.
2. Autogenous grinding circuit The milling circuit analyzed here consisted of a fully autogenous closed-loop mill with a recycling load that overflowed from a screen unit, as indicated in Fig. 1. The circuit was under expert control, which manipulated the mill’s control variables, namely its water inlet, coarse ore feed rate and fine ore feed rate, to control the output variables of the mill, i.e. the mill load and the power consumption of the mill. Five mill variables were monitored, viz. the fine ore feed rate (x1), coarse ore feed rate (x2), power consumption (x3), mill load (x4), and the inlet water flow rate (x5).
SCREEN
FEED FAG MILL
Fig. 1. Basic flow diagram of grinding circuit.
NE pA pAB Pnorm Pavg Pmax Pmin PROC q P Pq t Tq XE xEj x1 x2 x3 x4 x5 xt Xt yt z1, z2
number of rows in an embedded time series matrix (–) probability density function of random variable A (–) joint probability density function of random variables A and B (–) normalized mill power consumption (–) average mill power consumption (kW) maximum mill power consumption (kW) minimum mill power consumption (kW) mill power rate of change (kW h1) number of eigenvectors in a reduced dimensionality principal component model (–) eigenvector matrix (–) loading matrix with q columns (–) time (s) score matrix with q columns (–) matrix of embedded variables (–) vector of jth embedded variable (–) fine ore feed rate (ton h1) coarse ore feed rate (ton h1) power consumption (kW) mill load (ton) water inlet rate (ton h1) time series observation at time t (–) random variable at time t (–) surrogate data observation at time t (–) discriminant functions (linear combination of variables x1 to x5) (–)
Each individual time series (xi, i = 1, 2, . . . 5), consisting of 32,564 records, was normalized to zero mean and unit variance before further analysis. Typical time series measurements of these variables normalized between 0 and 1 are shown in Fig. 2. Since samples were collected approximately every 10 s, the time series span a period of almost 89 h of continuous operation. Moreover, state labels were logged by the mill controller for each sample collected in the dataset. These logged controller states, represented by the data in Table 1 and Fig. 3, are related to the operating conditions during the period under investigation, and also reflect the collective, heuristic knowledge of the operation of the comminution circuit. Note that not all the states are necessarily displayed, hence the labeling starting at no 3 and skipping some states in between, as not all the mill states occurred over the period of observation. Moreover, normal operating conditions were not defined formally, but could be considered the default state of the mill, i.e. when none of the other conditions were flagged by the controller. In Fig. 3, operational variables are indicated by FO = fine ore, CO = coarse ore, LD = load, PW = power and WI = water inlet, while the different mill controller states are represented by different color1 lines, as indicated in the legend of the figure. As can be seen from the radar chart (Fig. 3), mill states 3 and 6 have virtually identical profiles and would be difficult to distinguish reliably from the operational variables. This is also shown in Table 1, where the deviations of the variables from the normal state are shown in the last column in order from left to right for the fine ore (FO), coarse ore (CO), mill load (LD), power consumption (PW) and water inlet (WI). The mill variables were embedded into a phase space to facilitate analysis of the mill dynamics, as described in more detail in the next section.
1 For interpretation of color in Figs. 3, 5 and 12 the reader is referred to the web version of this article.
3
Fine Ore
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9 1 0.5 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Coarse Ore
Time
2
x 10
4
1 0.5 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Time
2
x 10
4
Power
1 0.5 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Time
2
x 10
4
Load
1 0.5 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Water Inlet
Time
2
x 10
4
1 0.5 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Time
2
x 10
4
Fig. 2. Typical realizations of the scaled mill variables considered in the analysis.
Table 1 Mill states logged by controller and deviations of associated operating conditions from the norm, in sequence from left to right of fine ore (FO), coarse ore (CO), power (PW), mill load (LD) and water inlet (WI). Mill state 3: Normal operating conditions 4: High loading of mill
Description
Deviations of operating variables
Normal operating conditions, the default state of the mill ðL
L
Þ
av g min Normalized load exceeding specified threshold Lnorm ¼ ðLmax Lmin Þ > a
5: Feed disturbance
Combined ore feed rate (x1 + x2) in excess of the average feed rate by a set threshold
6: Feed limited 7: Low mill loading
LROC < 0 and combined ore feed rate x1 + x2 in excess of a set threshold
8: High power consumption
av g min Having a normalized power larger than a specified threshold P norm ¼ ðPmax P min Þ > c
11: Recycle limited 12: Sand loading
ðL
L
Þ
av g min Having a normalized load smaller than a specified threshold Lnorm ¼ ðLmax Lmin Þ > b
ðP
P
Þ
– P ROC 6 d and ðP=LÞROC 6 h and
Circulating Feed Rate Total Feed Rate
P#
2
3. Phase space reconstruction of grinding circuit variables It is generally accepted that the space in which to view the structure of the dynamics of a system is not that of the one dimensional observations, or time series, but rather a space of higher dimensionality (Abarbanel, 1996). This higher dimensional space can be identified by means of a phase space reconstruction, and consists of coordinates composed of the original time series and their time delayed copies. The task of identifying the phase space from the original observations is addressed by the embedding theorem attributed to Takens (1981). Essentially, the reconstruction depends on two parameters, viz. the time delay or lag and the dimensionality of the phase space. More formally, if data matrix X represents N observations on m mill variables, as indicated in Eq. (1), then this matrix of variables is expanded (XE) to contain the variables, as well as a number of lagged copies of each variable, as indicated by Eqs. (2) and (3). For each variable X Ej , the lag kj, as well as the embedding dimension Mj, is determined individually, before concatenation into the aggregated lagged trajectory matrix XE.
X ¼ ½ x1
X E ¼ xE1
x2
xE2
X 11 6 . xm ¼ 6 4 .. X 1N xEm
. . .. ..
xm1
xmN
3 7 7 5
ð1Þ
ð2Þ
where
X Ej ¼ xj;t
x1;tk;j
x1;tðMj 1Þkj
ð3Þ
3.1. Estimation of embedding lag parameters of variables Optimal lags kj, for j = 1, 2, . . . m, were determined by the average mutual information criterion (Gallager, 1968), while the method of false nearest neighbors (Kennel et al., 1992) was used to determine the embedding dimension of the variables. The mutual information (I) is a measure of the amount of information gained with regard to variable bj, drawn from set B, when variable ai, drawn from set A, is measured:
4
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9 Table 2 Embedding of mill variables. Variable
Description
Units
k
M
x1 x2 x3 x4 x5
Mill load Mill power Inlet water Fine ore feed rate Coarse ore feed rate
[ton] [kW] [m3 s1] [ton s1] [ton s1]
1376 1017 1432 1477 1379
4 5 5 5 4
Fig. 3. Radar chart of the mill control states and operating variables summarized in Table 1 scaled to a range of [0, 1].
IðA; BÞ ¼ log2
pAB ðai ; bj Þ pA ðai ÞpB ðbj Þ
ð4Þ
In Eq. (4), pAB (ai,bj) is the joint probability density of measurements of the random variables ai 2 A and bj 2 B. pA(ai) and pB(bj) are the individual or marginal densities of measurements ai 2 A and bj 2 B. The average over all the measurements is known as the average mutual information and can be expressed as:
IðA; BÞ ¼
X p ðai ; bj Þ pAB ðai ; bj Þlog2 AB pðai ÞpB ðbj Þ a ;b i
ð5Þ
j
The ideal time lag was selected to correspond with the sample lag k, at which the average mutual information between Xc and Xtk reached a minimum.
IðX t ; X tk Þ ¼
X
pðxt ; xtk Þlog2
xt ;xtk
pðxt ; xtx Þ pðxt Þpðxtk Þ
Fig. 4. The eigenspectrum of the principal components of the embedded mill variables.
Following this, the eigenvectors P and eigenvalues ^ are calculated for the covariance matrix C using the eigenvalue decomposition.
C¼
1 E
N 1
T
X E XE ¼ P ^ PT
ð8Þ
ð6Þ
3.2. Estimation of embedding dimensions of variables The false nearest neighbor (FNN) algorithm assesses whether or not neighboring points in the phase space are close to one another owing to their dynamic origins, in which case they are considered true neighbors, or whether they appear to be neighbors, owing to the suboptimal unfolding of the phase space attractor as a result of an incorrect embedding dimension. If so, these points are considered to be false neighbors. The ideal embedding dimension is chosen as the dimension at which the percentage of false nearest neighbors reaches a minimum. As indicated in the last two columns of Table 2, the lags (k) and embedding dimensions (M) of the variables were very similar in magnitude.
Finally, the principal component scores constituting the phase space are calculated using principal components Pq (q columns of eigenvector matrix P), based on the reduced dimensionality q which captures significant variance, according any of a number of such criteria, such as the Kaiser criterion (Yeomans and Golder, 1982), scree tests, and some specified level or variance to be captured (Joliffe, 2011).
CT q ¼ X E P q
ð9Þ
As indicated from the eigenspectrum (K) of the principal components of this embedding in Fig. 4, a significant portion of the variation (approximately 61.5%) could be captured by the first three components, resulting in a three-dimensional subspace. This is sufficient for visualization of the phase space. 4. Tracking and controlling the mill dynamics
3.3. Multivariate embedding of mill variables 4.1. Visualization of controlled states The aggregate trajectory matrix XE of the mill variables had a dimensionality of NE = N max(kj) = 22,363 rows x P5 E M ¼ j¼1 M j = 23 columns, which still required further decorrelation to account for the cross-correlation between the different variables. This was accomplished by use of principal component analysis, the principal components of which then constituted the phase space coordinates of this multivariate embedding, i.e. in general for data set X, the covariance matrix C is constructed, i.e.
C¼
1 E
N 1
T
XE XE
ð7Þ
The reconstructed state space of the mill can be augmented by superimposing the controlled states of the mill onto the plot. This is shown in Fig. 5, for a phase space that was reconstructed from the data, after removal of a number of sections where close-to-zero mill power consumption was recorded. The controller states considered here are the states associated with normal operating conditions (NOC), or default state, and states associated with high mill loads, feed disturbances, feed limited conditions, low mill loads, recycle limited conditions, as well as sand loading. In this case, the first three principal component coordinates collectively account
5
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9
Fig. 5. Reconstructed phase space of the autogenous mill with controlled states determined by the mill controller superimposed.
for approximately 73.6% of the variance of the data, indicating a reasonably reliable projection of the data. In Fig. 5, it can be seen that some of the same states occur in different regions of the phase space. The controller state associated with sand loading, for example, occurs as two clusters (red markers) in the phase space, together with several smaller spots in between. The same goes for the low mill load controller state. 4.2. Transitions between controlled states To complement Fig. 5, Table 3 summarizes the time spent in each state over a period of approximately 66 h. The entries in the top left to bottom right diagonal of the table are the samples (at 10 s intervals) associated with each state, while the off-diagonal entries represent the frequencies of transitions between states. For example, from the NOC state, the mill control state moved to a high mill load state on three occasions and to a feed limited state on 57 occasions (1st row of the table). From a high mill load state, it moved to the NOC state only once, and twice to a state associated with a feed disturbance. 4.3. Analysis of mill load dynamics In order to better understand the dynamic behavior of the primary state variable of the mill, i.e. the load, this variable was also embedded into a phase space by itself. The embedding was based on a lag of k = 500 sample intervals, as estimated by use of Eq. (6), as well as an embedding dimension of m = 4, which was the lowest embedding dimension where the number of false nearest neighbors was negligible. The reconstructed attractor is shown in three dimensions in Fig. 6 (scaled data), which is an approximation of the
four-dimensional phase space into which the variable was embedded. Hypotheses about the nature of the underlying model that have generated the mill data can be tested by means of surrogate data analysis (Theiler et al., 1992; Barnard et al., 2001). In this case, the null hypothesis that the load time series data were generated by a Gaussian process undergoing a static transform was tested. The alternative hypothesis would suggest a nonlinear (possibly deterministic) process. The hypotheses can be tested by generating surrogate data that preserve the linear correlational structure of the original time series, as well as its marginal distribution. With the iterative amplitude adjusted Fourier transform (IAAFT) algorithm originally proposed by Schreiber and Schmitz (1996), a two-stage process is applied. It starts with a time series constituting white noise, and in the first stage, the Fourier amplitudes of the time series are replaced by the Fourier amplitudes of the original time series. In the second stage, the derived time series is reordered to the original marginal distribution until convergence of both power spectrum and marginal distribution is attained. Termination of the algorithm is triggered by complete convergence (same reordering in two consecutive steps), or if a maximum number of iterations is reached. More formally, if {xt} is a time series of length N and scaled to have a zero mean, then its Fourier transform is (Eq. (10)): N 1 X f ðxÞ ¼ pffiffiffiffiffiffiffiffiffiffi eixt xt ; 2pN t¼1
for p 6 x 6 p
ð10Þ
qffiffiffiffi For discrete frequencies, xj ¼ 2Np; j ¼ 1; 2 . . . N the transformation can be calculated according to Eq. (11), and the original time series can be recovered by back transformation.
xt ¼
rffiffiffiffiffiffiffi N 2pX ixj t e f ðxj Þ; N j¼1
for t ¼ 1; 2; . . . ; N
ð11Þ
Fourier transform surrogates construct new time series (yt) with the same periodogram as the time series, xt. The Fourier amplitudes are fixed, while the Fourier phases are replaced by uniformly distributed random numbers, urand ðxj Þ 2 ½0; 2p. The AAFT surrogate time series (yt) is realized as
yt ¼
rffiffiffiffiffiffiffi N 2pX ixj t e jf ðxj Þjeiurand ðxj Þ ; for t ¼ 1; 2; . . . ; N N j¼1
ð12Þ
After generation of the AAFT surrogates, the periodogram values are replaced by the periodogram of the original time series and the phases are retained. The amplitudes are subsequently adjusted in the time domain and the two steps are iterated until the
Table 3 Controlled states of the mill with diagonal entries (bold) showing time spent in a particular state and off-diagonal entries showing frequency of state transitions from the particular state in the left hand column to the state in the upper row. Mill states
NOC
High mill load
Feed disturbance
Feed limited
Low mill load
recycle limited
Sand loading
3. NOC 4. High mill load 5. Feed disturbance 6. Feed limited 8. Low mill load 11. Recycle limited 12. Sand loading
16,320 1 4 51 7 3 21
3 669 0 0 0 0 0
4 2 1099 2 0 0 0
57 0 1 3316 0 0 8
0 0 0 2 569 0 9
3 0 0 0 0 303 0
20 0 3 11 4 0 908
Total
16,407
672
1107
3382
880
306
946
6
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9
Fig. 6. Reconstructed attractor of the mill load time series data.
Fig. 7. Correlation dimensions of the mill load time series data (solid curve) and IAAFT surrogate data derived from the mill load time series (broken curves).
periodograms and the amplitude distributions of the original time series and its iterated AAFT surrogate are approximately equal.
When embedded in phase space, the topology of these time series can be represented in a robust way by their correlation dimensions (Grassberger and Procaccia, 1983a,b) and for this reason the correlation dimension is used very widely as a test statistic to compare time series data in phase spaces. In this paper, a version of the analysis proposed by Judd (1994) is used, where the correlation dimension (dc) is plotted for different values of a scale parameter (e0). A plot of the correlation dimensions of the original load time series data (solid curve) and 20 IAAFT surrogate time series (broken curves) derived from the original load time series is shown in Fig. 7. Since the correlation dimensions of the original data and their surrogates, dc(e0), are markedly different, the null hypothesis that the time series was generated by a stochastic process (Gaussian or otherwise) undergoing a static transform can be rejected. It is also interesting to note that the correlation dimension of the load (the solid curve in Fig. 7) is low (less than or equal to 2). Since the dimension of the correlation dimension is an indication of the number of dominant variables present in the evolution of the underlying system dynamics (Boon et al., 2008; Hill et al., 2008), it also suggests that only a few (two) variables would be needed to model the behavior of the mill load. This should not be construed as a generalization for the modeling of load behavior in autogenous mills, as in this case, it is applicable to a comparatively short period of observation of the particular mill run under set conditions of control.
Fig. 8. Prediction of the mill load with a multilayer perceptron neural network.
7
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9 Table 4 Classification of mill states with a multilayer perceptron on validation data. Class
Sand loading
No of samples 61 Correct (%) 54.1
NOC High mill load
Feed Feed Low disturbance limited mill load
Overall
118 72 84.7 94.4
177 94.9
600 86.7
113 91.2
59 81.4
system. Fig. 8 shows the results obtained with of a multilayer perceptron with six hidden layer nodes, trained by the Levenberg– Marquardt algorithm, which was used to predict 500 sample points ahead, based on four inputs, viz. the load at time t, t-500, t-1000 and t-1500. The number of nodes in the hidden layer of the neural network was optimized by use of the Schwarz–Christoffel information criterion (Schwarz, 1978). The neural network was trained on the first 38,128 samples of the time series data, and the prediction shown in Fig. 8 is for the last 5448 sample points.
5. Prediction of mill states from operational variables
Fig. 9. Linear projection of mill states showing maximum separability between the states.
Moreover, the identities of the dominant variables cannot be inferred from the analysis, and it can only serve as a guideline to the construction of models. Such a low dimensionality could also be indicative of a strong deterministic component in the data.
log C N N!1 log e0
dc ¼ lim lim s0!0
C N ðe0 Þ ¼
N 2
ð13Þ
1 X N
I k#i #j < e0 Þ
ð14Þ
06i
Although this would need to be confirmed by more extensive analysis beyond the scope of the present study, the predictability of the load behavior supports the indication of a deterministic
The mill states were predicted from the operational variables in Table 2 by use of a multilayer perceptron neural network configured as a classification model. The network had an input layer with five nodes, corresponding to five inputs x1, x2, . . . x5, as summarized in Table 2, a single hidden layer with nine bipolar sigmoidal activation units, and six output nodes with unipolar sigmoidal activation functions coding for the six mill states. A subset of the data consisting of 4000 samples was randomly divided into a training, test and validation data set, respectively containing 70%, 15% and 15% of the samples. The validation data set was used to assess the accuracy of the model after training with the training and test data sets. The network was trained with the Broyden–Fletcher–Goldfarb–Shanno (BFGS) gradient descent algorithm and a sum of squares error criterion and could predict the state of the mill with 86.7% accuracy on the validation data. A breakdown of the results are given in Table 4, where it can be seen that the model was not able to predict the state labeled Sand Loading meaningfully, but could predict the other states with an accuracy of at least 80%. The mill states are shown in Fig. 9, which shows a linear projection of the variables generated to maximize the separability between the different groupings or states by optimizing the ratio of the between scatter of the samples in the classes to the within scatter of the samples in the classes. From these data, it can be seen that the Sand Loading samples shown by small solid diamonds, are spread across several other groups. Since this class is defined by the rates of change of the power and the power to load ratio (decreasing more than set thresholds), as well as the ratio of feed rates, it is difficult to capture by instantaneous measurement of the variables only (as also suggested by the inability of the neural network model to identify the class).
high loading
1
Scaled Load [-]
NOC
0.8
feed limited
0.6
low loading
0.4 0.2 0
feed disturbance
sand loading
0.5
1
1.5
2
2.5
x 10
4
Time Index NOC
unknown
Scaled Power [-]
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.5
1
1.5
2
2.5
x 10
4
Time Index Fig. 10. Scaled load (top) and power (bottom) measurements with logged mill states representing a period of 89 h of operation.
8
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9
Fig. 11. Power–load curve of the mill with logged controller states.
6. Relationship between mill load and power consumption Fig. 10 shows the load and power measurements as time series data, with the logged controller states superimposed on the data. The relationship between the mill load and power consumption is shown in Fig. 11, again with the logged mill controller states superimposed on the data. As can be seen from Fig. 11, the curve is not well-defined, owing to the effect of the other mill variables not accounted for in the power–load plane.
The effect of disturbances on the power–load curve in Fig. 11 can be eliminated by using neural network models to estimate the curve at fixed (set point) values of the control values. A multilayer perceptron model trained for this purpose managed to account for approximately 71% of the variation in the power consumption, using mill load as a single input only. More than 92% of the variance in the power consumption could be explained by the model, when the water inlet, fine ore and coarse ore feed were also taken into account. More specifically, the model used to generate this smooth power–load curve had six logsigmoidal units in its single hidden layer and a tansigmoidal output unit and was trained with the Broyden–Fletcher–Goldfarb–Shanno (BFGS) gradient descent algorithm and a sum of squares error criterion, as before. A random sample of 70% of the data was used as training data set, while the remainder of the data was used to test and validate the model. Of these additional variables, the water inlet had the most important influence, as shown in Fig. 12. The contributions of the variables were determined by removal of the variable from the model and noting the effect on the model output. The broken red line in Fig. 12 coincides with the 99% significance level for the variable contributions to the model and was generated by use of Monte Carlo simulation of the effects of random input variables in the model, as discussed in more detail by Auret and Aldrich (2011). The solid black curve in Fig. 13 shows the correction when the multilayer perceptron neural network model was used to predict
Fig. 12. Relative importance of variables in power–load curve model.
Fig. 13. Corrected power–load curve generated by neural network model (bottom left) with set point inputs for control variables (solid black line) and neural network model (top left) with 10-sample localized mean values for all inputs (broken red line). The blue dots are the raw data used to validate the models. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9
the power–load curve at the mean values of the control variables (fine ore feed, coarse ore feed and water inlet). Alternatively, the outputs of the model can be averaged over local regions to smooth the curve, as indicated by the broken red line in Fig. 13 with a moving average of 10 samples. However, this approach is not as reliable as that based on the use of a neural network model.
7. Conclusions In this study, the dynamic behavior of an autogenous mill was analyzed by embedding the mill variables in a multivariate phase space where all the variables could be simultaneously considered together with information on the operating states of the mill. In addition, the mill load was also embedded in the phase space by itself. Based on qualitative (visualization) and quantitative (surrogate data) analysis of the phase space embeddings, the following conclusions can be made. The controlled states of the mill could be visualized as separate clusters in low-dimensional (2D or 3D) mappings of the phase space that consisted of multivariate embeddings of the operational variables of the mill. This validates the labels in the mill controller assigned to the different states. These controlled states of the mill could be predicted with a satisfactory degree of accuracy by multilayer perceptron neural networks. Comparison of the correlation dimensions of the load time series data with random time series that had the same power spectrum and marginal distribution as the load data suggest that the load data do not conform to a linear stochastic model. Instead, the data could have a significant deterministic component and were influenced by relatively few (perhaps as few as two) dominant variables. Neural network models were constructed that could predict the power load curve with a high degree of accuracy (R2 > 0.92) and, in principle at least, these models could be used to eliminate the effect of disturbances on the power load curve, which could then be used to assist plant operators to control the mill.
Acknowledgements The authors gratefully acknowledge: Anglo American Platinum’s Advanced Process Control Group for their inputs, as well as making plant data available to the authors.
9
The Anglo American Platinum Centre for Process Monitoring at the University of Stellenbosch in South Africa for use of their computational infrastructure.
References Abarbanel, H.D.I., 1996. Analysis of Observed Chaotic Data, first ed. Springer-Verlag, New York, NY, USA. Auret, L., Aldrich, C., 2011. Empirical comparison of tree ensemble variance importance measures. Chemometrics and Intelligent Laboratory Systems 105, 157–170 (submitted, 30 Jun 10, conditionally accepted, 27 Aug 2010, In press, 22 Dec 2010). Barnard, J.P., Aldrich, C., Gerber, M., 2001. Identification of dynamic pro cess systems with surrogate data methods. AIChE Journal 47 (9), 2064–2075. Boon, M.Y., Henry, B.I., Suttle, C.M., Dain, S., 2008. The correlation dimension: a useful objective measure of the transient visual evoked potential? Journal of Vision 8(1) (6), 1–21. Fuerstenau, D.W., Abouzeid, A.-Z.M., 2002. The energy efficiency of ball milling in comminution. International Journal of Mineral Processing 67, 161–185. Gallager, R.G., 1968. Information Theory and Reliable Communication. John Wiley and Sons, New York, NY, USA. Grassberger, P., Procaccia, I., 1983a. Measuring the strangeness of strange attractors. Physica D: Nonlinear Phenomena 9 (1–2), 189–208. Grassberger, P., Procaccia, I., 1983b. Characterization of strange attractors. Physical Review Letters 50 (5), 346–349. Hill, J., Hossain, F., Sivakumar, B., 2008. Is correlation dimension a reliable proxy for the number of dominant influencing variables for modeling risk of arsenic contamination in groundwater? Stochastic Environmental Research and Risk Assessment 22 (1), 47–55. Joliffe, I.T., 2011. Principal Component Analysis. In: Springer Series in Statistics, second ed. Springer Verlag, New York. Judd, K., 1994. Estimating dimension from small samples. Physica D 71, 421– 429. Kennel, M.B., Brown, R., Abarbanel, H.D.I., 1992. Determining minimum embedding dimension using a geometrical construction. Physical Review A 45, 3403–3411. Kiangi, K.K., Moys, M.H., 2008. Particle filling and size effects on the ball load behavior and power in a dry pilot mill: experimental study. Powder Technology 187 (1), 79–87. Lameck, N.S., Kiangi, K.K., Moys, M.H., 2006. Effect of grinding media shapes on load behavior and mill power in a dry ball mill. Minerals Engineering 19 (13), 1357– 1361. Powell, M.S., Mainza, A.N., 2006. Extended grinding curves are essential to the comparison of milling performance. Minerals Engineering 19 (15), 1487–1494. Powell, M.S., Morrell, S., Latchireddi, S., 2001. Development in the understanding of South African style SAG mills. Minerals Engineering 14 (10), 1143–1153. Powell, M.S., van der Westhuizen, A.P., Mainza, A.N., 2009. Applying grindcurves to mill operation and optimisation. Minerals Engineering 22 (7–8), 625–632. Schreiber, T., Schmitz, A., 1996. Improved Surrogate Data for Nonlinearity Tests. Physical Review Letters 77, 635–638. Schwarz, G.E., 1978. Estimating the dimension of a model. Annals of Statistics 6 (2), 461–464. Takens, F., 1981. Detecting strange attractors in turbulence. Lecture Notes in Mathematics 898, 311–366. Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., 1992. Testing for nonlinearity in time series: the method of surrogate data. Physica D 58, 77–94. Van Nierop, M.A., Moys, M.H., 2001. Exploration of mill power modeled as a function of load behavior. Minerals Engineering 14 (10), 1267–1276. Yeomans, K.A., Golder, P.A., 1982. The Guttman–Kaiser criterion as a predictor of the number of common factors. Journal of the Royal Statistical Society. Series D (The Statistician) 31 (3), 221–229.
Contents lists available at ScienceDirect
Minerals Engineering journal homepage: www.elsevier.com/locate/mineng
Visualization of the controller states of an autogenous mill from time series data C. Aldrich a,b,⇑, J.J. Burchell b, J.W. de V. Groenewald c, C. Yzelle b a
Department of Minerals Engineering and Extractive Metallurgy, Western Australian School of Mines, Curtin University of Technology, Perth, WA, Australia Department of Process Engineering, University of Stellenbosch, Private Bag X1, Matieland 7602, Stellenbosch, South Africa c Anglo Platinum Management Services, PO Box 62179, Marshalltown 2107, South Africa b
a r t i c l e
i n f o
Article history: Received 3 June 2013 Accepted 14 October 2013 Available online 14 November 2013 Keywords: Comminution Modeling Nonlinear time series analysis Neural networks
a b s t r a c t The operational variables of an industrial autogenous mill were embedded in a low-dimensional phase space to facilitate visualization of the dynamic behavior of the mill. This was accomplished by use of a multivariate extension of the method of delay coordinates used in nonlinear time series analysis. In this phase space, the controlled states of the mill could be visualized as separate regions or clusters in the phase space. Comparison of the correlation dimension of the state variable of the mill (the load) embedded in phase space suggested that the dynamic behavior of the mill could not be represented by a linear stochastic model (Gaussian or otherwise). The low dimensionality (62) of the correlation dimension further suggested that the mill load depended on a few variables only and that the underlying generative process had a significant deterministic component. In addition, the operational variables could be used as reliable predictors in a neural network model to identify the controlled states of the mill. As a complementary approach to visualization of the operation of the mill, a different neural network model could be used to reconstruct a corrected power load curve by compensating for the effect of varying operating conditions. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Energy consumption in ore milling processes constitutes up to 40% of the energy expenditure in concentrator plants, but is not more than 15–25% efficient, depending on the approach used to calculate the energy required to reduce the size of the ore particles (Fuerstenau and Abouzeid, 2002). As a consequence, the reduction of energy consumption in comminution circuits has drawn considerable interest over an extended period of time. Energy efficient operation of mills is known to depend to a large extent on the behavior of the mill load (Van Nierop and Moys, 2001), but reliably accounting for the effect of load on mill power consumption is not an easy task. Other process variables, such as media shape (Lameck et al., 2006), mill speed (Powell et al., 2009), ore type (Powell et al., 2009) and mill filling (Kiangi and Moys, 2008) also play a significant role that may be difficult to quantify from first principles. This is especially the case with autogenous and semi-autogenous grinding (SAG) mills, which are popular in new plants, owing ⇑ Corresponding author at: Department of Minerals Engineering and Extractive Metallurgy, Western Australian School of Mines, Curtin University of Technology, Perth, WA, Australia. Tel.: +61 8 9266 4349; fax: +61 8 9358 4912. E-mail address: [email protected] (C. Aldrich). 0892-6875/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.mineng.2013.10.018
to their favorable capital and operating costs. SAG mill operation, in particular, is strongly affected by the type and size of the ore being ground, which makes stable control difficult. In these mills, grinding is driven by the mill load (Powell et al., 2009). One approach to mill control is based on the use of power–load curves (Powell et al., 2001, 2009). This approach involves establishing empirical grind curves that reveal at which mill load the circuit’s power draw peak. These curves have a typical concave parabolic shape when plotted with power and load as the ordinate and abscissa respectively. Once these curves are established, they may be used to guide milling operation with respect to achieving a desired throughput or grind quality. In this context, for example, current control strategies are aimed at the maintenance of maximal mill power draw at or near the top of the power load curve (Powell and Mainza, 2006) by manipulation of the mill feed rate and other process variables, such as the coarse to fine ore feed ratio and the water inlet rate. This aim is challenging, since the power load curve is not static for a particular milling circuit. Instead, the curve is affected by other process variables, such as the ore characteristics, and it may therefore be difficult to track the designated operating point with sufficient accuracy. Therefore, accounting for the effect of multiple variables on the controlled states of the mill would be a major advantage.
2
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9
Nomenclature Symbol Meaning a, b, c, d, h, # threshold values used by expert control system (–) e0 scaling parameter (–) K diagonal eigenvalue matrix (–) ti , tj ith and jth point of a time series embedded in phase space (–) x frequency (–) C correlation matrix (–) CN correlation count (–) dc correlation dimension (–) f ðxÞ Fourier transform (–) I(A,B) mutual information of two random variables A and B (–) IðA; BÞ average mutual information of two random variables A and B (–) k embedding lag of phase space (–) kj embedding lag of jth variable in phase space (–) Lnorm normalized load (–) Lavg average load (ton) Lmax maximum load (ton) Lmin minimum load (ton) LROC load rate of change (ton h1) m number of variables (–) M embedding dimension of phase space (–) Mj embedding dimension of j’th variable in phase space (–) N number of samples in a time series (–)
In this paper, analysis of the dynamic behavior of an autogenous mill is considered, based on reconstruction of the mill variables in phase space, which allows time series variables and their history to be considered simultaneously within a state space context. It is shown that most of these states can be predicted from the operational variables of the mill.
2. Autogenous grinding circuit The milling circuit analyzed here consisted of a fully autogenous closed-loop mill with a recycling load that overflowed from a screen unit, as indicated in Fig. 1. The circuit was under expert control, which manipulated the mill’s control variables, namely its water inlet, coarse ore feed rate and fine ore feed rate, to control the output variables of the mill, i.e. the mill load and the power consumption of the mill. Five mill variables were monitored, viz. the fine ore feed rate (x1), coarse ore feed rate (x2), power consumption (x3), mill load (x4), and the inlet water flow rate (x5).
SCREEN
FEED FAG MILL
Fig. 1. Basic flow diagram of grinding circuit.
NE pA pAB Pnorm Pavg Pmax Pmin PROC q P Pq t Tq XE xEj x1 x2 x3 x4 x5 xt Xt yt z1, z2
number of rows in an embedded time series matrix (–) probability density function of random variable A (–) joint probability density function of random variables A and B (–) normalized mill power consumption (–) average mill power consumption (kW) maximum mill power consumption (kW) minimum mill power consumption (kW) mill power rate of change (kW h1) number of eigenvectors in a reduced dimensionality principal component model (–) eigenvector matrix (–) loading matrix with q columns (–) time (s) score matrix with q columns (–) matrix of embedded variables (–) vector of jth embedded variable (–) fine ore feed rate (ton h1) coarse ore feed rate (ton h1) power consumption (kW) mill load (ton) water inlet rate (ton h1) time series observation at time t (–) random variable at time t (–) surrogate data observation at time t (–) discriminant functions (linear combination of variables x1 to x5) (–)
Each individual time series (xi, i = 1, 2, . . . 5), consisting of 32,564 records, was normalized to zero mean and unit variance before further analysis. Typical time series measurements of these variables normalized between 0 and 1 are shown in Fig. 2. Since samples were collected approximately every 10 s, the time series span a period of almost 89 h of continuous operation. Moreover, state labels were logged by the mill controller for each sample collected in the dataset. These logged controller states, represented by the data in Table 1 and Fig. 3, are related to the operating conditions during the period under investigation, and also reflect the collective, heuristic knowledge of the operation of the comminution circuit. Note that not all the states are necessarily displayed, hence the labeling starting at no 3 and skipping some states in between, as not all the mill states occurred over the period of observation. Moreover, normal operating conditions were not defined formally, but could be considered the default state of the mill, i.e. when none of the other conditions were flagged by the controller. In Fig. 3, operational variables are indicated by FO = fine ore, CO = coarse ore, LD = load, PW = power and WI = water inlet, while the different mill controller states are represented by different color1 lines, as indicated in the legend of the figure. As can be seen from the radar chart (Fig. 3), mill states 3 and 6 have virtually identical profiles and would be difficult to distinguish reliably from the operational variables. This is also shown in Table 1, where the deviations of the variables from the normal state are shown in the last column in order from left to right for the fine ore (FO), coarse ore (CO), mill load (LD), power consumption (PW) and water inlet (WI). The mill variables were embedded into a phase space to facilitate analysis of the mill dynamics, as described in more detail in the next section.
1 For interpretation of color in Figs. 3, 5 and 12 the reader is referred to the web version of this article.
3
Fine Ore
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9 1 0.5 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Coarse Ore
Time
2
x 10
4
1 0.5 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Time
2
x 10
4
Power
1 0.5 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Time
2
x 10
4
Load
1 0.5 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Water Inlet
Time
2
x 10
4
1 0.5 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Time
2
x 10
4
Fig. 2. Typical realizations of the scaled mill variables considered in the analysis.
Table 1 Mill states logged by controller and deviations of associated operating conditions from the norm, in sequence from left to right of fine ore (FO), coarse ore (CO), power (PW), mill load (LD) and water inlet (WI). Mill state 3: Normal operating conditions 4: High loading of mill
Description
Deviations of operating variables
Normal operating conditions, the default state of the mill ðL
L
Þ
av g min Normalized load exceeding specified threshold Lnorm ¼ ðLmax Lmin Þ > a
5: Feed disturbance
Combined ore feed rate (x1 + x2) in excess of the average feed rate by a set threshold
6: Feed limited 7: Low mill loading
LROC < 0 and combined ore feed rate x1 + x2 in excess of a set threshold
8: High power consumption
av g min Having a normalized power larger than a specified threshold P norm ¼ ðPmax P min Þ > c
11: Recycle limited 12: Sand loading
ðL
L
Þ
av g min Having a normalized load smaller than a specified threshold Lnorm ¼ ðLmax Lmin Þ > b
ðP
P
Þ
– P ROC 6 d and ðP=LÞROC 6 h and
Circulating Feed Rate Total Feed Rate
P#
2
3. Phase space reconstruction of grinding circuit variables It is generally accepted that the space in which to view the structure of the dynamics of a system is not that of the one dimensional observations, or time series, but rather a space of higher dimensionality (Abarbanel, 1996). This higher dimensional space can be identified by means of a phase space reconstruction, and consists of coordinates composed of the original time series and their time delayed copies. The task of identifying the phase space from the original observations is addressed by the embedding theorem attributed to Takens (1981). Essentially, the reconstruction depends on two parameters, viz. the time delay or lag and the dimensionality of the phase space. More formally, if data matrix X represents N observations on m mill variables, as indicated in Eq. (1), then this matrix of variables is expanded (XE) to contain the variables, as well as a number of lagged copies of each variable, as indicated by Eqs. (2) and (3). For each variable X Ej , the lag kj, as well as the embedding dimension Mj, is determined individually, before concatenation into the aggregated lagged trajectory matrix XE.
X ¼ ½ x1
X E ¼ xE1
x2
xE2
X 11 6 . xm ¼ 6 4 .. X 1N xEm
. . .. ..
xm1
xmN
3 7 7 5
ð1Þ
ð2Þ
where
X Ej ¼ xj;t
x1;tk;j
x1;tðMj 1Þkj
ð3Þ
3.1. Estimation of embedding lag parameters of variables Optimal lags kj, for j = 1, 2, . . . m, were determined by the average mutual information criterion (Gallager, 1968), while the method of false nearest neighbors (Kennel et al., 1992) was used to determine the embedding dimension of the variables. The mutual information (I) is a measure of the amount of information gained with regard to variable bj, drawn from set B, when variable ai, drawn from set A, is measured:
4
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9 Table 2 Embedding of mill variables. Variable
Description
Units
k
M
x1 x2 x3 x4 x5
Mill load Mill power Inlet water Fine ore feed rate Coarse ore feed rate
[ton] [kW] [m3 s1] [ton s1] [ton s1]
1376 1017 1432 1477 1379
4 5 5 5 4
Fig. 3. Radar chart of the mill control states and operating variables summarized in Table 1 scaled to a range of [0, 1].
IðA; BÞ ¼ log2
pAB ðai ; bj Þ pA ðai ÞpB ðbj Þ
ð4Þ
In Eq. (4), pAB (ai,bj) is the joint probability density of measurements of the random variables ai 2 A and bj 2 B. pA(ai) and pB(bj) are the individual or marginal densities of measurements ai 2 A and bj 2 B. The average over all the measurements is known as the average mutual information and can be expressed as:
IðA; BÞ ¼
X p ðai ; bj Þ pAB ðai ; bj Þlog2 AB pðai ÞpB ðbj Þ a ;b i
ð5Þ
j
The ideal time lag was selected to correspond with the sample lag k, at which the average mutual information between Xc and Xtk reached a minimum.
IðX t ; X tk Þ ¼
X
pðxt ; xtk Þlog2
xt ;xtk
pðxt ; xtx Þ pðxt Þpðxtk Þ
Fig. 4. The eigenspectrum of the principal components of the embedded mill variables.
Following this, the eigenvectors P and eigenvalues ^ are calculated for the covariance matrix C using the eigenvalue decomposition.
C¼
1 E
N 1
T
X E XE ¼ P ^ PT
ð8Þ
ð6Þ
3.2. Estimation of embedding dimensions of variables The false nearest neighbor (FNN) algorithm assesses whether or not neighboring points in the phase space are close to one another owing to their dynamic origins, in which case they are considered true neighbors, or whether they appear to be neighbors, owing to the suboptimal unfolding of the phase space attractor as a result of an incorrect embedding dimension. If so, these points are considered to be false neighbors. The ideal embedding dimension is chosen as the dimension at which the percentage of false nearest neighbors reaches a minimum. As indicated in the last two columns of Table 2, the lags (k) and embedding dimensions (M) of the variables were very similar in magnitude.
Finally, the principal component scores constituting the phase space are calculated using principal components Pq (q columns of eigenvector matrix P), based on the reduced dimensionality q which captures significant variance, according any of a number of such criteria, such as the Kaiser criterion (Yeomans and Golder, 1982), scree tests, and some specified level or variance to be captured (Joliffe, 2011).
CT q ¼ X E P q
ð9Þ
As indicated from the eigenspectrum (K) of the principal components of this embedding in Fig. 4, a significant portion of the variation (approximately 61.5%) could be captured by the first three components, resulting in a three-dimensional subspace. This is sufficient for visualization of the phase space. 4. Tracking and controlling the mill dynamics
3.3. Multivariate embedding of mill variables 4.1. Visualization of controlled states The aggregate trajectory matrix XE of the mill variables had a dimensionality of NE = N max(kj) = 22,363 rows x P5 E M ¼ j¼1 M j = 23 columns, which still required further decorrelation to account for the cross-correlation between the different variables. This was accomplished by use of principal component analysis, the principal components of which then constituted the phase space coordinates of this multivariate embedding, i.e. in general for data set X, the covariance matrix C is constructed, i.e.
C¼
1 E
N 1
T
XE XE
ð7Þ
The reconstructed state space of the mill can be augmented by superimposing the controlled states of the mill onto the plot. This is shown in Fig. 5, for a phase space that was reconstructed from the data, after removal of a number of sections where close-to-zero mill power consumption was recorded. The controller states considered here are the states associated with normal operating conditions (NOC), or default state, and states associated with high mill loads, feed disturbances, feed limited conditions, low mill loads, recycle limited conditions, as well as sand loading. In this case, the first three principal component coordinates collectively account
5
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9
Fig. 5. Reconstructed phase space of the autogenous mill with controlled states determined by the mill controller superimposed.
for approximately 73.6% of the variance of the data, indicating a reasonably reliable projection of the data. In Fig. 5, it can be seen that some of the same states occur in different regions of the phase space. The controller state associated with sand loading, for example, occurs as two clusters (red markers) in the phase space, together with several smaller spots in between. The same goes for the low mill load controller state. 4.2. Transitions between controlled states To complement Fig. 5, Table 3 summarizes the time spent in each state over a period of approximately 66 h. The entries in the top left to bottom right diagonal of the table are the samples (at 10 s intervals) associated with each state, while the off-diagonal entries represent the frequencies of transitions between states. For example, from the NOC state, the mill control state moved to a high mill load state on three occasions and to a feed limited state on 57 occasions (1st row of the table). From a high mill load state, it moved to the NOC state only once, and twice to a state associated with a feed disturbance. 4.3. Analysis of mill load dynamics In order to better understand the dynamic behavior of the primary state variable of the mill, i.e. the load, this variable was also embedded into a phase space by itself. The embedding was based on a lag of k = 500 sample intervals, as estimated by use of Eq. (6), as well as an embedding dimension of m = 4, which was the lowest embedding dimension where the number of false nearest neighbors was negligible. The reconstructed attractor is shown in three dimensions in Fig. 6 (scaled data), which is an approximation of the
four-dimensional phase space into which the variable was embedded. Hypotheses about the nature of the underlying model that have generated the mill data can be tested by means of surrogate data analysis (Theiler et al., 1992; Barnard et al., 2001). In this case, the null hypothesis that the load time series data were generated by a Gaussian process undergoing a static transform was tested. The alternative hypothesis would suggest a nonlinear (possibly deterministic) process. The hypotheses can be tested by generating surrogate data that preserve the linear correlational structure of the original time series, as well as its marginal distribution. With the iterative amplitude adjusted Fourier transform (IAAFT) algorithm originally proposed by Schreiber and Schmitz (1996), a two-stage process is applied. It starts with a time series constituting white noise, and in the first stage, the Fourier amplitudes of the time series are replaced by the Fourier amplitudes of the original time series. In the second stage, the derived time series is reordered to the original marginal distribution until convergence of both power spectrum and marginal distribution is attained. Termination of the algorithm is triggered by complete convergence (same reordering in two consecutive steps), or if a maximum number of iterations is reached. More formally, if {xt} is a time series of length N and scaled to have a zero mean, then its Fourier transform is (Eq. (10)): N 1 X f ðxÞ ¼ pffiffiffiffiffiffiffiffiffiffi eixt xt ; 2pN t¼1
for p 6 x 6 p
ð10Þ
qffiffiffiffi For discrete frequencies, xj ¼ 2Np; j ¼ 1; 2 . . . N the transformation can be calculated according to Eq. (11), and the original time series can be recovered by back transformation.
xt ¼
rffiffiffiffiffiffiffi N 2pX ixj t e f ðxj Þ; N j¼1
for t ¼ 1; 2; . . . ; N
ð11Þ
Fourier transform surrogates construct new time series (yt) with the same periodogram as the time series, xt. The Fourier amplitudes are fixed, while the Fourier phases are replaced by uniformly distributed random numbers, urand ðxj Þ 2 ½0; 2p. The AAFT surrogate time series (yt) is realized as
yt ¼
rffiffiffiffiffiffiffi N 2pX ixj t e jf ðxj Þjeiurand ðxj Þ ; for t ¼ 1; 2; . . . ; N N j¼1
ð12Þ
After generation of the AAFT surrogates, the periodogram values are replaced by the periodogram of the original time series and the phases are retained. The amplitudes are subsequently adjusted in the time domain and the two steps are iterated until the
Table 3 Controlled states of the mill with diagonal entries (bold) showing time spent in a particular state and off-diagonal entries showing frequency of state transitions from the particular state in the left hand column to the state in the upper row. Mill states
NOC
High mill load
Feed disturbance
Feed limited
Low mill load
recycle limited
Sand loading
3. NOC 4. High mill load 5. Feed disturbance 6. Feed limited 8. Low mill load 11. Recycle limited 12. Sand loading
16,320 1 4 51 7 3 21
3 669 0 0 0 0 0
4 2 1099 2 0 0 0
57 0 1 3316 0 0 8
0 0 0 2 569 0 9
3 0 0 0 0 303 0
20 0 3 11 4 0 908
Total
16,407
672
1107
3382
880
306
946
6
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9
Fig. 6. Reconstructed attractor of the mill load time series data.
Fig. 7. Correlation dimensions of the mill load time series data (solid curve) and IAAFT surrogate data derived from the mill load time series (broken curves).
periodograms and the amplitude distributions of the original time series and its iterated AAFT surrogate are approximately equal.
When embedded in phase space, the topology of these time series can be represented in a robust way by their correlation dimensions (Grassberger and Procaccia, 1983a,b) and for this reason the correlation dimension is used very widely as a test statistic to compare time series data in phase spaces. In this paper, a version of the analysis proposed by Judd (1994) is used, where the correlation dimension (dc) is plotted for different values of a scale parameter (e0). A plot of the correlation dimensions of the original load time series data (solid curve) and 20 IAAFT surrogate time series (broken curves) derived from the original load time series is shown in Fig. 7. Since the correlation dimensions of the original data and their surrogates, dc(e0), are markedly different, the null hypothesis that the time series was generated by a stochastic process (Gaussian or otherwise) undergoing a static transform can be rejected. It is also interesting to note that the correlation dimension of the load (the solid curve in Fig. 7) is low (less than or equal to 2). Since the dimension of the correlation dimension is an indication of the number of dominant variables present in the evolution of the underlying system dynamics (Boon et al., 2008; Hill et al., 2008), it also suggests that only a few (two) variables would be needed to model the behavior of the mill load. This should not be construed as a generalization for the modeling of load behavior in autogenous mills, as in this case, it is applicable to a comparatively short period of observation of the particular mill run under set conditions of control.
Fig. 8. Prediction of the mill load with a multilayer perceptron neural network.
7
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9 Table 4 Classification of mill states with a multilayer perceptron on validation data. Class
Sand loading
No of samples 61 Correct (%) 54.1
NOC High mill load
Feed Feed Low disturbance limited mill load
Overall
118 72 84.7 94.4
177 94.9
600 86.7
113 91.2
59 81.4
system. Fig. 8 shows the results obtained with of a multilayer perceptron with six hidden layer nodes, trained by the Levenberg– Marquardt algorithm, which was used to predict 500 sample points ahead, based on four inputs, viz. the load at time t, t-500, t-1000 and t-1500. The number of nodes in the hidden layer of the neural network was optimized by use of the Schwarz–Christoffel information criterion (Schwarz, 1978). The neural network was trained on the first 38,128 samples of the time series data, and the prediction shown in Fig. 8 is for the last 5448 sample points.
5. Prediction of mill states from operational variables
Fig. 9. Linear projection of mill states showing maximum separability between the states.
Moreover, the identities of the dominant variables cannot be inferred from the analysis, and it can only serve as a guideline to the construction of models. Such a low dimensionality could also be indicative of a strong deterministic component in the data.
log C N N!1 log e0
dc ¼ lim lim s0!0
C N ðe0 Þ ¼
N 2
ð13Þ
1 X N
I k#i #j < e0 Þ
ð14Þ
06i
Although this would need to be confirmed by more extensive analysis beyond the scope of the present study, the predictability of the load behavior supports the indication of a deterministic
The mill states were predicted from the operational variables in Table 2 by use of a multilayer perceptron neural network configured as a classification model. The network had an input layer with five nodes, corresponding to five inputs x1, x2, . . . x5, as summarized in Table 2, a single hidden layer with nine bipolar sigmoidal activation units, and six output nodes with unipolar sigmoidal activation functions coding for the six mill states. A subset of the data consisting of 4000 samples was randomly divided into a training, test and validation data set, respectively containing 70%, 15% and 15% of the samples. The validation data set was used to assess the accuracy of the model after training with the training and test data sets. The network was trained with the Broyden–Fletcher–Goldfarb–Shanno (BFGS) gradient descent algorithm and a sum of squares error criterion and could predict the state of the mill with 86.7% accuracy on the validation data. A breakdown of the results are given in Table 4, where it can be seen that the model was not able to predict the state labeled Sand Loading meaningfully, but could predict the other states with an accuracy of at least 80%. The mill states are shown in Fig. 9, which shows a linear projection of the variables generated to maximize the separability between the different groupings or states by optimizing the ratio of the between scatter of the samples in the classes to the within scatter of the samples in the classes. From these data, it can be seen that the Sand Loading samples shown by small solid diamonds, are spread across several other groups. Since this class is defined by the rates of change of the power and the power to load ratio (decreasing more than set thresholds), as well as the ratio of feed rates, it is difficult to capture by instantaneous measurement of the variables only (as also suggested by the inability of the neural network model to identify the class).
high loading
1
Scaled Load [-]
NOC
0.8
feed limited
0.6
low loading
0.4 0.2 0
feed disturbance
sand loading
0.5
1
1.5
2
2.5
x 10
4
Time Index NOC
unknown
Scaled Power [-]
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.5
1
1.5
2
2.5
x 10
4
Time Index Fig. 10. Scaled load (top) and power (bottom) measurements with logged mill states representing a period of 89 h of operation.
8
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9
Fig. 11. Power–load curve of the mill with logged controller states.
6. Relationship between mill load and power consumption Fig. 10 shows the load and power measurements as time series data, with the logged controller states superimposed on the data. The relationship between the mill load and power consumption is shown in Fig. 11, again with the logged mill controller states superimposed on the data. As can be seen from Fig. 11, the curve is not well-defined, owing to the effect of the other mill variables not accounted for in the power–load plane.
The effect of disturbances on the power–load curve in Fig. 11 can be eliminated by using neural network models to estimate the curve at fixed (set point) values of the control values. A multilayer perceptron model trained for this purpose managed to account for approximately 71% of the variation in the power consumption, using mill load as a single input only. More than 92% of the variance in the power consumption could be explained by the model, when the water inlet, fine ore and coarse ore feed were also taken into account. More specifically, the model used to generate this smooth power–load curve had six logsigmoidal units in its single hidden layer and a tansigmoidal output unit and was trained with the Broyden–Fletcher–Goldfarb–Shanno (BFGS) gradient descent algorithm and a sum of squares error criterion, as before. A random sample of 70% of the data was used as training data set, while the remainder of the data was used to test and validate the model. Of these additional variables, the water inlet had the most important influence, as shown in Fig. 12. The contributions of the variables were determined by removal of the variable from the model and noting the effect on the model output. The broken red line in Fig. 12 coincides with the 99% significance level for the variable contributions to the model and was generated by use of Monte Carlo simulation of the effects of random input variables in the model, as discussed in more detail by Auret and Aldrich (2011). The solid black curve in Fig. 13 shows the correction when the multilayer perceptron neural network model was used to predict
Fig. 12. Relative importance of variables in power–load curve model.
Fig. 13. Corrected power–load curve generated by neural network model (bottom left) with set point inputs for control variables (solid black line) and neural network model (top left) with 10-sample localized mean values for all inputs (broken red line). The blue dots are the raw data used to validate the models. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
C. Aldrich et al. / Minerals Engineering 56 (2014) 1–9
the power–load curve at the mean values of the control variables (fine ore feed, coarse ore feed and water inlet). Alternatively, the outputs of the model can be averaged over local regions to smooth the curve, as indicated by the broken red line in Fig. 13 with a moving average of 10 samples. However, this approach is not as reliable as that based on the use of a neural network model.
7. Conclusions In this study, the dynamic behavior of an autogenous mill was analyzed by embedding the mill variables in a multivariate phase space where all the variables could be simultaneously considered together with information on the operating states of the mill. In addition, the mill load was also embedded in the phase space by itself. Based on qualitative (visualization) and quantitative (surrogate data) analysis of the phase space embeddings, the following conclusions can be made. The controlled states of the mill could be visualized as separate clusters in low-dimensional (2D or 3D) mappings of the phase space that consisted of multivariate embeddings of the operational variables of the mill. This validates the labels in the mill controller assigned to the different states. These controlled states of the mill could be predicted with a satisfactory degree of accuracy by multilayer perceptron neural networks. Comparison of the correlation dimensions of the load time series data with random time series that had the same power spectrum and marginal distribution as the load data suggest that the load data do not conform to a linear stochastic model. Instead, the data could have a significant deterministic component and were influenced by relatively few (perhaps as few as two) dominant variables. Neural network models were constructed that could predict the power load curve with a high degree of accuracy (R2 > 0.92) and, in principle at least, these models could be used to eliminate the effect of disturbances on the power load curve, which could then be used to assist plant operators to control the mill.
Acknowledgements The authors gratefully acknowledge: Anglo American Platinum’s Advanced Process Control Group for their inputs, as well as making plant data available to the authors.
9
The Anglo American Platinum Centre for Process Monitoring at the University of Stellenbosch in South Africa for use of their computational infrastructure.
References Abarbanel, H.D.I., 1996. Analysis of Observed Chaotic Data, first ed. Springer-Verlag, New York, NY, USA. Auret, L., Aldrich, C., 2011. Empirical comparison of tree ensemble variance importance measures. Chemometrics and Intelligent Laboratory Systems 105, 157–170 (submitted, 30 Jun 10, conditionally accepted, 27 Aug 2010, In press, 22 Dec 2010). Barnard, J.P., Aldrich, C., Gerber, M., 2001. Identification of dynamic pro cess systems with surrogate data methods. AIChE Journal 47 (9), 2064–2075. Boon, M.Y., Henry, B.I., Suttle, C.M., Dain, S., 2008. The correlation dimension: a useful objective measure of the transient visual evoked potential? Journal of Vision 8(1) (6), 1–21. Fuerstenau, D.W., Abouzeid, A.-Z.M., 2002. The energy efficiency of ball milling in comminution. International Journal of Mineral Processing 67, 161–185. Gallager, R.G., 1968. Information Theory and Reliable Communication. John Wiley and Sons, New York, NY, USA. Grassberger, P., Procaccia, I., 1983a. Measuring the strangeness of strange attractors. Physica D: Nonlinear Phenomena 9 (1–2), 189–208. Grassberger, P., Procaccia, I., 1983b. Characterization of strange attractors. Physical Review Letters 50 (5), 346–349. Hill, J., Hossain, F., Sivakumar, B., 2008. Is correlation dimension a reliable proxy for the number of dominant influencing variables for modeling risk of arsenic contamination in groundwater? Stochastic Environmental Research and Risk Assessment 22 (1), 47–55. Joliffe, I.T., 2011. Principal Component Analysis. In: Springer Series in Statistics, second ed. Springer Verlag, New York. Judd, K., 1994. Estimating dimension from small samples. Physica D 71, 421– 429. Kennel, M.B., Brown, R., Abarbanel, H.D.I., 1992. Determining minimum embedding dimension using a geometrical construction. Physical Review A 45, 3403–3411. Kiangi, K.K., Moys, M.H., 2008. Particle filling and size effects on the ball load behavior and power in a dry pilot mill: experimental study. Powder Technology 187 (1), 79–87. Lameck, N.S., Kiangi, K.K., Moys, M.H., 2006. Effect of grinding media shapes on load behavior and mill power in a dry ball mill. Minerals Engineering 19 (13), 1357– 1361. Powell, M.S., Mainza, A.N., 2006. Extended grinding curves are essential to the comparison of milling performance. Minerals Engineering 19 (15), 1487–1494. Powell, M.S., Morrell, S., Latchireddi, S., 2001. Development in the understanding of South African style SAG mills. Minerals Engineering 14 (10), 1143–1153. Powell, M.S., van der Westhuizen, A.P., Mainza, A.N., 2009. Applying grindcurves to mill operation and optimisation. Minerals Engineering 22 (7–8), 625–632. Schreiber, T., Schmitz, A., 1996. Improved Surrogate Data for Nonlinearity Tests. Physical Review Letters 77, 635–638. Schwarz, G.E., 1978. Estimating the dimension of a model. Annals of Statistics 6 (2), 461–464. Takens, F., 1981. Detecting strange attractors in turbulence. Lecture Notes in Mathematics 898, 311–366. Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., 1992. Testing for nonlinearity in time series: the method of surrogate data. Physica D 58, 77–94. Van Nierop, M.A., Moys, M.H., 2001. Exploration of mill power modeled as a function of load behavior. Minerals Engineering 14 (10), 1267–1276. Yeomans, K.A., Golder, P.A., 1982. The Guttman–Kaiser criterion as a predictor of the number of common factors. Journal of the Royal Statistical Society. Series D (The Statistician) 31 (3), 221–229.