- Email: [email protected]
Adaptive control of a laboratory grinding circuit
M ~ . ~ Engl~oer#lg, Vol. 2, No. 2, pp. 245-258, 1989 Printed in Great Britain
0892-6875/89 $3.00+0.00 © 1989 Pergamon Preu pic
ADAPTIVE CONTROL OF A LABORATORY GRINDING CIRCUIT
Y. DUB~ $ and D. HODOUIN § t
Dept. d'Ing~nierie, Universit~ du Quebec a Trois-Rivi~res, C.P. 500, Trois-Rivi~res, 0u~bec, Canada G9A 5H7 § Dept. de Mines et Metallurgie, Universite Laval, Ste. Foy, Quebec, Canada GIK 7P4
(Received 5 January 1989)
ABSTRACT A multivarlable adaptive control scheme has been studied and evaluated with the help of a phenomenological simulator which was calibrated on a laboratory grinding circuit. The adaptive control scheme includes a Kalman filter, an adaptive mechanism, and a control law. The adaptive mecb-nism is based on the estimation of a fundamental coefficient of the process since a direct identification of the whole set of parameters is impossible due to the large size of the state-space model. The control law that minimizes a quadratic criterion is obtained by dynamic programming. Simulated results show that the control law exhibits good performances. The control scheme is robust and is able to face IOOZ change on grinding rate parameters, without deterioration of the performances. Keywords Computer modeling, Grinding circuit, Adaptive control, Kalman filtering, Quadratic control. INTRODUCTION The main objective of the implementation of a computer control scheme for a c o n t i n u o u s p r o c e s s is u s u a l l y to m i n i m i z e the v a r i a t i o n s of the product quality while minimizing the variations of the control variables. Such a goal can be a c h i e v e d by f o r m u l a t i n g some p e r f o r m a n c e c r i t e r i o n w h i c h can be mathematlcally optimized when a suitable dynamic model of the process to be controlled is available. The performance of the control strategy will depend on both the quality of the selected control scheme and the accuracy of the process model. F o r the s p e c i f i c case of ore g r i n d i n g circuits, automation should help to maintain operating conditions, such as circulating load r a t i o s and p r o d u c t fineness, close to their optimum values while smoothing the variations of the manlpulated variables such as ore feed rate, water addition rate, pump speed. The definition of a performance criterion should includes these requirements with an adequate weighting of the v a r i o u s a s p e c t s . The f o r m u l a t i o n of a q u a d r a t i c o b j e c t i v e f u n c t i o n is the e a s i e s t way t o h a n d l e this problem, because the theory of optimal control is well developed in this context. As the model of the process is usually based on a linear state-space formulation, this approach leads to LQ control algorithms [1,2]. This type of s t r a t e g y will be tested in the present study. T h e r e are two m a i n t e c h n i q u e s to d e v e l o p a d y n a m i c m o d e l of a g r i n d i n g circuit. The empirical approach consists in the description of the process by a transfer matrix identified by suitable experimental design; then the model is reformulated in terms of state-space variables. It has already been shown that this approach is feasible [3]. The second technique, which will be used here, is based on a phenomenological description of the process [4,5]. The d e v e l o p m e n t of the model requires a good understanding of the mechanisms of g r i n d i n g a n d c l a s s i f i c a t i o n and is m o r e c o m p l e x to m a n i p u l a t e that the e m p i r i c a l a p p r o a c h b e c a u s e it leads to h l g h - o r d e r non l i n e a r equations. H o w e v e r the p h e n o m e n o l o g i c a l m o d e l a l l o w s t h e d e f i n i t i o n of p h y s i c a l parameters in which the prior knowledge of the process can be embedded. This 245
246
Y. Dune and D. Ho~tnN
advantage is used here to design an adaptive mechanism to take into account that the process model changes when the circuit efficiency changes. This gives more robustness to the model and as a consequence to the control algorithm. In the paper a phenomenological state-space model is developed, linearized and c a l i b r a t e d for a l a b o r a t o r y g r i n d i n g c i r c u i t involving a ball mill and a screen for the r e c y c l i n g of the c o a r s e p a r t i c l e s to the m i l l feed. The control state feedback law is then established using a quadratic performance index. The states are observed using a Kalman filter. Both the filtering and control algorithms are made adaptive to cope with grinding efficiency drifts. The m e t h o d is e v a l u a t e d u s i n g a g r i n d i n g c i r c u i t s i m u l a t o r i n v o l v i n g a detailed description of the equipment behaviour and including noise generators [6,7,8].
LABORATORY GRINDING CIRCUIT The laboratory grinding circuit is shown in Figure I. The grinding mill is a grate discharge ball mill. The mill inside diameter is 40 cm and its length 40 cm. The mill is equipped with a c y l i n d r i c a l s c r e e n f i x e d on the m i l l output. The fine p a r t i c l e s p a s s i n g t h r o u g h the s c r e e n (finer t h a n 150 micrometers) are analysed by the Leeds and Northrup Microtrac particle size monitor. The coarse particle flow feeds a constant level sump and is recycled to the input of the mill. The measured sump level value is the input to a sump-pump controller and the speed of the pneumatic pump is controlled by an a u t o m a t i c air v a l v e . The r e c y c l i n g s t r e a m f l o w r a t e is m e a s u r e d by a magnetic flowmeter and a densimeter is also included. The circuit is fed by two automatic conveyors which are equipped with load cells and tachometers to measure the feedrate. Their set points can be selected manually or calculated by the computer. The percentage of solids in the mill is normally selected as 65% and maintained constant by varying the water pump speed. The set point of the water pump is given by the computer, based on the fresh ore feedrate and the recycling stream flowrate measurements.
~
Feeders
~~
~Centrltu pgu,mlp °'"
In Diaphragm pump
[I Fine product~ Fig.1
Centrifugal~ pump
Laboratory grinding circuit
As mentioned in the introduction, the control algorithm will be evaluated on a simulator. B a t c h t e s t s and d y n a m i c m e a s u r e m e n t s were made to adjust the parameters of the simulator [9]. Steady-state experiments give the values of the phenomenological parameters of the mill: the fragment size distribution y (breakage function), the rate of breakage of particles (selection f u n c t i o n )
Adaptive control of a laboratory grinding circuit
247
and the hold-up in the mill [10]. Dynamic measurements are used to verify the reliability of the simulator. Grinding tests have shown that a good agreement is obtained between simulated and measured values [5]. PHENOMENOLOGICAL
PROCESS MODEL
The m a t h e m a t i c a l m o d e l of the l a b o r a t o r y g r i n d i n g c i r c u i t is b a s e d discretized p a r t i c l e sizes and g r i n d i n g p a r a m e t e r s . This a p p r o a c h extensively used [references 11 to 16] and only the necessary parameters simulate the laboratory circuit are presented here.
on is to
Hold-up mill model: To determine the volume of water and solids retained in the mill, the m i l l t r a n s p o r t behaviour is assumed to be modelled by three identical interactive tanks in series. The steady-state model for the hold-up is given by the following equations: V t _ 6LQo
+ Vr
K
where Qo is the volumetric flowrate, V t the total pulp volume, V r the residual pulp volume, and K and L two c o e f f i c i e n t s to be d e t e r m i n e d . The d y n a m i c variations of the pulp volumes in the three tanks are given in [17]. Grinding model: Grinding dynamics are represented by two perfect mixers in series where firstorder breakage occurs [5]. It is assumed that the solids are d i v i d e d into c l a s s e s (index K) which corresponds to grinding properties characterized by the values of the breakage and rate functions. A matrix formulation of the equations describing the process is for the class K: K d(H~ mj) dt v
cK :
dt where
(SH) K g
K Hj m~
(I)
H .K
J J dC~
-- M K m~_1 - K K j-1 Mjmj + (B K- I)
J --
Qj-I v 3
CK ( 9-I
- cK ) ]
HK is the hold-up weight of solids of class K in mixer j; 3 H the total solid weight ~Z,kH~);j M the mass flowrate; m the particle size distribution; B the breakage function matrix (distribution of fragments); C the mass of solids per unit volume; SH the selection function matrix, expressed as the mass of particles of each size class which are broken per unit of time; Qj the volumetric flowrate of the mixer j.
Classification model: The c l a s s i f i c a t i o n of the mill o u t p u t is m a d e by a s c r e e n with openings (150 micrometers). The screen is assumed to have zero-order and the fraction of particles recycled to the mill is given by:
constant dynamics
Mrm r = EMom o w h e r e E is the c l a s s i f i c a t i o n matrix (proportion of feed particles in each size class retained on the screen), M o the mass flowrate at mill output, M r the r e c y c l e d flowrate and m o and m r the particle size distributions at the same locations.
248
Y. DUBE and D. HODOUIN
L I N E A R I Z E D S T A T E - S P A C E MODEL The above d e s c r i b e d model which is used in the simulator is continuous and non-linear. H o w e v e r to d e s i g n a d i g i t a l f i l t e r a n d a d i g i t a l c o n t r o l a l g o r i t h m a l i n e a r d i s c r e t e - t i m e model is required. Two kinds of d i s c r e t e models are c o m m o n l y used, ARMA model and state-space model. Since the statespace approach has been selected, the p h e n o m e n o l o g i c a l m o d e l m u s t be discretized, linearized and w r i t t e n under the u s u a l s t o c h a s t i c s t a t e - s p a c e form [18]:
(2)
Xk+ I = PkXk + QkUk + w k Yk = FkXk + GkUk + Vk
where x is the state-space vector, u the control vector, k the time index, P, Q, F, G m a t r i c e s of coefficients, and w and v random white noises. In order to linearize Eq. (I) it is assumed that the mill behaves as a single perfect mixer and that the mill hold-up has zero-order dynamics. Also the m a t e r i a l properties description is s i m p l i f i e d , a s s u m i n g that a m i x t u r e of ores of various g r i n d a b i l i t i e s (classes K) behaves as a m a t e r i a l that has an average grindability. This leads to the following e q u a t i o n for the overall solids: H dm dt
=
Mfmf - Mfm
+ (B - I)SHm
I
w h e r e Mf is the m i l l m a s s f e e d r a t e and m the m i l l o u t p u t p a r t i c l e s i z e distribution. The final state-space form of the d i s c r e t e - t i m e model for the w h o l e circuit is given in r e f e r e n c e [5]. For the present i l l u s t r a t i o n of the method, the state-space vector x has 22 e l e m e n t s . T h e f i r s t 20 p o s i t i o n s r e p r e s e n t the v a r i a t i o n of t h e m i l l d i s c h a r g e p a r t i c l e - s i z e d i s t r i b u t i o n a r o u n d its s t e a d y - s t a t e v a l u e s . T h e l a s t two p o s i t i o n s r e p r e s e n t the v a r i a t i o n (around their steady-state value) of the mill feedrate and of the m i l l f e e d r a t e d e l a y e d by o n e s a m p l i n g i n t e r v a l . T h e a g r e e m e n t of t h i s d i s c r e t e model with the complete n o n - l i n e a r p h e n o m e n o l o g i c a l simulator and the p i l o t - s c a l e d y n a m i c e x p e r i m e n t a l m e a s u r e m e n t s is d i s c u s s e d in the p r e v i o u s l y cited references. KALMAN F I L T E R I N G The Kalman filter is an a l g o r i t h m that m i n i m i z e s the v a r i a n c e of the filtered v a r i a b l e s and r e c o n s t r u c t the s t a t e - s p a c e v e c t o r which is g e n e r a l l y d i f f e r e n t from the vector of the m e a s u r e d variables. The laboratory circuit m e a s u r e d values are the p a r t i c l e size d i s t r i b u t i o n of the screen undersize flow and the f e e d r a t e to t h e m i l l ( r e c y c l e d p l u s f r e s h o r e f l o w r a t e ) . T h e 20 f i r s t e l e m e n t s of t h e s t a t e v e c t o r r e p r e s e n t the m i l l d i s c h a r g e p a r t i c l e size distribution. Thus the Kalman filter has to r e c o n s t r u c t the mill d i s c h a r g e particle size distribution, b a s e d on the c i r c u i t p r o d u c t p a r t i c l e s i z e d i s t r i b u t i o n and has to e s t i m a t e m i n i m u m v a r i a n c e v a l u e s for t h e c i r c u i t product size distribution. The p a r t i c l e size d i s t r i b u t i o n of the circuit fine product (mD) is related to the mill output p a r t i c l e size d i s t r i b u t i o n (mo) by the following equation: Mf mp = ( I Mp
- E) m
where M D is the mass flowrate of the circuit fine product and I the identity m a t r i x . - In the above e x p r e s s i o n Mp and m are related to the part of the state vector x w h i c h contains the mill output p a r t i c l e size d i s t r i b u t i o n (noted x i for the particle size interval i) as follows: m i = msi + x i Mp = Mf~
(i - Ei)(msi + x i)
w h e r e m s is the s t e a d y - s t a t e v a l u e of m. As v a r i a b l e mp is related to the state v a r i a b l e by:
a consequence
the m e a s u r e d
Adaptive control of a laboratory grinding circuit
(I
mpi =
-
249
E i)
(x i + msi) ~(I - Ei)(msi + x i)
= aix i + b i
(3)
Since a i and b i are f u n c t i o n s of x i, this o b s e r v a t i o n e q u a t i o n is not linear. To o v e r c o m e t h i s p r o b l e m , the x i v a l u e in a i and b i is r e p l a c e d b ~ its e s t i m a t e d value given by Eq. 2 w h e r e x k is taken as the filtered value x k at the p r e v i o u s time. T h e second m e a s u r e d value is the r e c y c l e d f l o w r a t e and its v a r i a t i o n around its s t e a d y - s t a t e value is equal to the state v e c t o r p o s i t i o n 21. The overall o b s e r v a t i o n e q u a t i o n can finally be w r i t t e n as: Yk = FkXk + GkUk The Kalman algorithm, here for convenience: Xk+1
which
can
= (PkXk + QkUk ) + Kk+l[Yk+1
be
found
in many
text-books,
is
summarized
- F k + 1 ( P k X k + QkUk)]
Kk+ I = M k + i F ' k + 1 ( F k + i M k + i F ' k + I + Vk+1) -I Mk+ I = PkZkPk , Z k = (I - K k F k ) M k. w h e r e ~ is the filtered given by: W k = E[WkW'k]
;
s t a t e - s p a c e vector
and W and V are d i a g o n a l m a t r i c e s
Vk = E[VkV'k]
T h e s t a n d a r d d e v i a t i o n of the r e c y c l e d f l o w r a t e is selected as 25% of the current s t e a d y - s t a t e value. The standard d e v i a t i o n of the c i r c u i t p r o d u c t particle size distribution is s e l e c t e d as 15% of the current s t e a d y - s t a t e value. These values are used to c o n s t r u c t the m a t r i x V. T h e m a t r i x W is adjusted empirically to g i v e a s u i t a b l e r e d u c t i o n of the v a r i a n c e by the K a l m a n filter. The m a t r i x W value is selected as: W_- I0-5 i
CONTROL A L G O R I T H M The control a l g o r i t h m is based on a m u l t i v a r i a b l e state f e e d b a c k scheme [I, 19] w h e r e the c o n t r o l law is d e s i g n e d for m i n i m i z i n g a q u a d r a t i c o b j e c t i v e function. T h e m e t h o d is f i r s t p r e s e n t e d in a g e n e r a l c a s e , t h e n it is i l l u s t r a t e d by s e l e c t i n g a specific a p p l i c a t i o n to the g r i n d i n g circuit. The o b j e c t i v e f u n c t i o n is of the f o l l o w i n g form: co
J = Z [U'kRU k + e'kTe k] k--1
(4)
w h e r e ek, the d e v i a t i o n of the process output to a set point v e c t o r Z is: ek = Zk - Yk The m i n i m i z a t i o n of the o b j e c t i v e f u n c t i o n is p e r f o r m e d under the c o n s t r a i n t of the process model given by the f o l l o w i n g linear state-space equation: r
Xk+ I = PkXk + QkUk Yk = FkXk + GkUk T h e q u a d r a t i c o b j e c t i v e function (4) does not g u a r a n t e e a zero s t e a d y - s t a t e error. To avoid this effect n u m e r i c a l integrators are added to the initial s t a t e - s p a c e model. They have the f o l l o w i n g form: Vk÷1
-- Vk + ek = Vk + Zk - Yk
250
Y. DUBE and D. HoDoun~
E x p a n d i n g we obtain: Vk+ I = v k - FkX k - GkU k + Zk
i-) (- °I()(°) u (o}
The final augmented system of equations is:
VK+ I
-F k
vk
-G k
zk
() (, :)() (0:) =
Vk
UK
+
O
(5)
Vk
and more compactly: Xk+1
= AkXk ÷ BkUK + fk
Yk = CkXk + DkUk The new objective function is defined as:
3
=
7
k=1 where:
[e'kTek + U'kRU k]
(z:)()
o) ;
~k =
Tv
Vk
The optimal control law m i n i m i z i n g the new cost function is given by: Uk = - LkXk + ~k w h e r e L and ~ are calculated by a dynamic p r o g r a m m i n g method and given in (I):
Lk+ I = (R + B'kKkB k + DkTDK) -I
[B'kKkA k + D'kTC k]
Uk+1
[D'kTZ k - B' k (gk ÷ Kkfk)]
= (R + B'kKkB k + DkTDk) -I
Kk+ I = ~kKk [A k - BkLk+ I ] + ~kCk gk+1
= ~k(gk + Kk [BkUk+1
(6)
+ f]) - #kZk
and: ~k = A'k - C ' k T D k ( R + D'kTDk )-I B'k ~k = C'T - C'TDk(R + D~TDk )-I D' T k k k The matrix L and the vector U are functions of T and R which must be positive definite. L a r g e T values will lead to large v a r i a t i o n s of the control actions u compared to the output error variations. C o n v e r s e l y if R is large c o m p a r e d to ~, the o u t p u t e r r o r v a r i a t i o n s will be large r e l a t i v e l y to the control a c t i o n s u. M a t r i c e s T and R a r e t h e r e f o r e c h o s e n i t e r a t i v e l y to g i v e s a t i s f a c t o r y closed loop dynamics. To i l l u s t r a t e t h i s type of control strategy an a p p l i c a t i o n to the g r i n d i n g process p r e v i o u s l y d e s c r i b e d will be given in the next sections. Generally the control o b j e c t i v e for a g r i n d i n g system is to m a i n t a i n the product size a n d the m i l l t h r o u g h p u t c o n s t a n t u s i n g t h e t w o a v a i l a b l e manipulated variables: the o r e f e e d r a t e and the s u m p w a t e r a d d i t i o n rate. In the present case the water a d d i t i o n to the circuit does not influence either the screen p e r f o r m a n c e or the mill efficiency. This is due to the fact that the
Adaptive control of a laboratory ~'inding circuit
251
percentage of solids in the mill is maintained constant by a local feed-back loop and that the water is added to the sump only to allow a suitable pumping behaviour. Furthermore the screen performance is almost independent of the operating conditions leading to a naturally stable fine p r o d u c t size. The most sensitive operating variable is here the mill throughput which has to be controlled for a stable operation of both the mill and the screen. The amount of recycled solids is a function of the fresh feed particle size distribution and of the ore grindability, variables which are considered as disturbances, and a function of the fresh ore feed rate which is the manipulated variable. Despite the fact that only the mill throughput is controlled by the ore fresh feed rate, the control strategy is based on the whole set of state variables of the process, including o b v i o u s l y the p a r t i c l e size d i s t r i b u t i o n . The extension to the multivariable case is straightforward since the formalism is exactly the same. For instance one could have added in the criteria another o u t p u t v a r i a b l e s of the m o d e l such as some characteristics of the product particle size, with a suitable weighting factor. ADAPTIVE MECHANISM The filtering and control algorithms are based on the matrices A, B, C and D. When the ore characteristics and the o p e r a t i n g levels v a r y the m o d e l parameters change and the matrices A, B, C and D must be continuously updated to keep track of the parameter variations. In the present case, preliminary simulation runs show that the v a r i a t i o n of the g r i n d i n g rate p a r a m e t e r s m o d i f i e s the process dynamics. As a consequence the optimal controller is detuned when the process model changes. Usual real-time adaptive algorithms, such as r e c u r s i v e least s q u a r e s e s t i m a t i o n of the model matrix elements, cannot be considered here due to the very large number of parameters involved in the matrices A, B, C and D. The method selected is based on a key parameter of the phenomenological model which is estimated using the deviation between the experimental values and the calculated results. It is assumed that the main variations of the mill model can be taken into account by a modification of the selection f u n c t i o n by a multiplicative constant Uk, such that at time k: S k = UkSo The constant u k modifies by the same factor the rate of breakage of all the particles, independently of their size. So it characterizes the grinding mill efficiency and lumps all the effects due to the grlndabillty and partlcle size distribution v a r i a t i o n s of the f r e s h feed as w e l l as the e f f e c t s of the v a r i a t i o n of the o p e r a t i n g c o n d i t i o n s such as the ball load and the pulp theology. This index is recursively calibrated using the difference between the discrete model prediction and the experimental v~lue (here the simulated value) of the mill throughput (respectively M~ and M~). The ek index is given by: P ~k+1
=
Uk
+
8(Mf
-
e Mf)
The speed of the adaptive algorithm is adjusted by the parameter B which empirically selected to give suitable dynamics to the estimation process.
is
The parameter ~ does not explicitly appear in the model matrices A, B, C and D. However it is formally present in the original phenomenological model. As a c o n s e q u e n c e , the m a t r i c e s A, B, C and D h a v e to be rederived from the p h e n o m e n o l o g i c a l m o d e l each time that the p a r a m e t e r ~ is updated by the adaptive mechanism. SIMULATED RESULTS In a first step the adaptation algorithm will be tested. In the absence of disturbances the steady-state value of the simulated mill throughput is 47.5 kg/h for a feedrate of 20 kg/h. At time index k = 0 (the sampling period is 18 seconds) the f r e s h ore g r i n d a b i l i t y is d o u b l e d . This d i s t u r b a n c e is equivalent to a multiplication of the selection function by a factor 2. The fresh ore feedrate is also perturbated according to the scheme in Figure 2.
252
Y, Dime and D. HODOUn~
F u r t h e r m o r e , to compare easily the responses of the discrete model and the adaptive discrete model, the simulation of the measured variable is performed without addition of simulated measurement noises. The comparison between the mill throughput predicted by the discrete model without adaptation mechanism and the mill throughput given by the simulator is given in Figure 2. Figure 3 gives the response of the mill throughput predicted by the adaptive discrete model and by the simulator. A brief comparison of these curves clearly shows the advantage of the adaptive procedure. 52
44 GD
v
36 LLI '~ r",, 2 8 w L~J h 20
'1 12
I
i
I i
'
'
I
I
I
I
52 SIMULATOR DISCRETE MODEL
L--
44
v
g
3s
1"
~ 2s
......f
-~ 20 .J 12 0
Fig.2
I
I
40
80
I
I
120 TIME
I
160 200 INOEX
I
I
240
280
320
Predicted mill throughput without model adaptation
In a second step the performance of the control algorithm depicted in Figure 4 is studied. This algorithm must keep the mill throughput at a constant set point using fresh ore fe~drate as control variable. The p e r f o r m a n c e index given by Eq. (4) combines the squared variations of the mill throughput and of the fresh ore feedrate. The state-space model is obtained by substituting in Eq. (5) the matrices P, Q, F and G corresponding to the state and observation equations of the grinding circuit. 52
SIMULATOR ------
44
DISCRETE
MODEL
I-
~..r ss ~ as "1" I.-
J
J
f
20
12 0
I
I
I
I
40
80
120
160
TIME
Fig.3
I
200
I
240
I
280
300
INDEX
Predicted mill throughput with model adaptation
Adaptive control of a laboratory grinding circuit
grlndiblllty dlitur blnce Fresh ore feed rite
.J
253
miilurement nolle
LABORATORY I Mill
throughput
I I OItINDIHO S M IULATORCIRCUITI - Pr°duc' "eriicl" ~ I slse
KALMAN FILTER Filtered states
I
LINEARIZED
-=,
I
mill
I [
O,T,MA,
I CALCULATIONI-oF B, C I ANo DA'MATRICESI -
Fig.4
DISCRETE
dlltrlbuUon
"
II
1.
I
I/+"
CsIculitldP
throughput
ADAPTIVEI_
MECHANISM
Block diagram of the control strategy
The o b j e c t i v e of the control strategy is to keep the v a r i a t i o n of the m i l l t h r o u g h p u t a r o u n d its s t e a d y - s t a t e value M s as close as p o s s i b l e to a setpoint value. Using the p r e v i o u s l y defined notation the variable: M
-
M s
=
x(21)
÷
u
should track the set point value. T h e o p t i m a l c o n t r o l law, g i v e n in Eq. 6, is now tested w e i g h t i n g factors of the p e r f o r m a n c e c r i t e r i o n selected as:
T =
(2:) ;
by
simulation
for
R = I
0 T h e d i s t u r b a n c e w h i c h is g e n e r a t e d in the simulator is a change of the ore g r i n d a b i l i t y which, for a s t e a d y - s t a t e operation, w o u l d be e q u i v a l e n t to a m u l t i p l i c a t i o n of the s e l e c t i o n f u n c t i o n by a factor 2. F i g u r e 5 shows the r e s p o n s e of the u n c o n t r o l l e d g r i n d i n g c i r c u i t to t h i s g r i n d a b i l i t y d i s t u r b a n c e w h i l e the fresh ore feedrate is m a i n t a i n e d at a value of 20 kg/h. A strong v a r i a t i o n of the mill throughput can be observed. Under control, F i g u r e 6 shows the optimal t r a j e c t o r y of u and the v a r i a t i o n of the mill t h r o u g h p u t around the set point of 47.5 kg/h (steady state value, for a fresh ore feedrate of 20 kg/h, before d i s t u r b i n g the ore grindability). This c a l c u l a t i o n is p e r f o r m e d w i t h o u t any s i m u l a t i o n of m e a s u r e m e n t n o i s e . The c o n t r o l a l g o r i t h m is very e f f i c i e n t in the absence of noise measurement, to m a i n t a i n the mill throughput at the set point value.
50 46 v
42
~: (9
38
jo. o
34
t-.J -'2 3O 26
0
Fig.5
ME 2/2--O
I 16
I 32
J 48
J 64 TIME
I I 80 96 INDEX
112
1 2i 8
1 4i 4
U n c o n t r o l l e d mill t h r o u g h p u t r e s p o n s e to a g r i n d a b i l i t y
increase step
254
Y. Dum~ and D. HODOUIN
4eI ..........
SET
POINT
~44 -
32
24
2O
Fig.6
0
16
Optimal
32
control
I I 12
I
48
64 80 96 TIME INDEX
performances
in a n o i s e - f r e e
I 128
I 144
environment
(R = I)
When noises, with the s t a t i s t i c a l c h a r a c t e r i s t i c s mentioned in t h e K a l m a n filtering section, are simulated and added to the m e a s u r e d variables, the optimal control results of F i g u r e 7 are obtained. Very strong v a r i a t i o n s of the manipulated v a r i a b l e can be o b s e r v e d w h i c h are not a c c e p t a b l e for many p r a c t i c a l reasons. To reduce their a m p l i t u d e a larger value can be given to R. F i g u r e 8 s h o w s the b e h a v i o u r of the control scheme for R = 100 w i t h o u t noise a d d i t i o n and Figure 9 its b e h a v i o u r in the p r e s e n c e of noise. One can s e e f r o m F i g u r e 8 r e s u l t s t h a t the v a r i a t i o n s of the mill t h r o u g h p u t are larger than the v a r i a t i o n s w h e n R w a s g i v e n a v a l u e I. However smaller v a r i a t i o n s of fresh ore feedrate are needed in the p r e s e n c e of noises. Figure 10 and 11 show the responses for R = 2500 and similar remarks can be made. 80
A
70
A
~< 6o
A
A A
A
I-
~.
A
A
A A
A
A
50
AA
A
A
;ET POINT o
A
A
A
~.
A
~
A
A
6
A
/,
A
A
A A
A
A
A
40
A ~A A
A
A A
i•J ,.1
3O
A
i
A
I
20 ~
i
I 0
30
0 0
o
0
I
I
I
00
o
0
o
o
0 0
0
o
oo
o
o0
o O0 000
00 0
0
0
0
00
0
0
0
0
0
0
0
0
0 I( 0
Fig.7
I
o
o
0 0
o w w "
I
40 o
~
I
I 16
Optimal
I 32
I 48
| 64
I 80
TIME
INDEX
I 96
I 112
I 128
I 144
control p e r f o r m a n c e s in the p r e s e n c e m e a s u r e m e n t noises (R = I)
of
Adaptive control of a laboratory grinding circuit
255
I,-
~,..56 [ - ~ - za. 'LR
" _1
SET POINT
40 52 24
~' uJ
16 o bJ LIJ "
8
I 16
0
I 32
I 48
I 64
I 80
TIME Fig.8
O p t i m a l control p e r f o r m a n c e s
I 96
I 112
I 128
I 144
INDEX
in a n o i s e - f r e e e n v i r o n m e n t
(R = 500)
In the p r e s e n c e of simulated m e a s u r e m e n t noises, it is d i f f i c u l t to see the effect of R on the true mill t h r o u g h p u t due to the large standard d e v i a t i o n of the m e a s u r e m e n t noises. To clarify this aspect of the control strategy, the r e s u l t s are p r e s e n t e d in F i g u r e 12 w i t h o u t the m e a s u r e m e n t noise, w h i l e the optimal control law is c a l c u l a t e d using t h e n o i s y m e a s u r e m e n t of the m i l l t h r o u g h p u t a n d p r o d u c t p a r t i c l e size d i s t r i b u t i o n . This is p o s s i b l e o n l y b e c a u s e the v a r i a b l e s are c a l c u l a t e d by a simulator then d i s t u r b e d by a noise generator. A m i n i m u m v a r i a t i o n for the mill t h r o u g h p u t is o b t a i n e d for R =100. For R = I the control a l g o r i t h m is v e r y s e n s i t i v e to the n o i s e a n d p r o d u c e s large v a r i a t i o n s on control variable. For R = 2500 the v a r i a t i o n s of the mill t h r o u g h p u t are slowly corrected.
72 A v
64
-
A
A
A
A
A
I= 0. "I(9
56
o(¢
48
-J _1
40
A
A A
A
A
-r" I-
A
A
A
A
/&
A
.,.
A
A
~&~ A A
A
A
A
:E
~
SET POINT
A
A
A
A
I
32
t~
b'l"3 2 J
0
24 ua'-"
A
A A
A
A
m
~.1 u.
A
51
A
A
A
I
I
I
00000000000000000
I
I
I
000000000000000000000
oO
L_oo,.o
i., o IRJ-o
v
j/6
i
32
I
48
i
64 TIME
Fig.9
I
80
I
96
I
112
I
128
I
144
INDEX
O p t i m a l control p e r f o r m a n c e s in the presence of m e a s u r e m e n t noises (R = 500)
256
Y. DUBEand D. HODOUIN
k-
Q-
56,
f-
~:~'r 40
I
t
I
I
I
I
I
I
I
I
I
I
32
24
""
I,LI
oc
a I.LI bJ h
16
8
0 Fig.10
Optimal
I
I
I
16
32
48
control
64 80 96 T I M E INDEX
performances
in a n o i s e - f r e e
These s i m u l a t i o n tests show that R, and also Q, can p a r a m e t e r s of the control law. The choice of Q and R of t h e n o i s e s a n d on the p r o c e s s d y n a m i c s of the relative c h o i c e of t h e s e s m a t r i c e s m u s t be d o n e d i f f i c u l t when noises with large v a r i a n c e s affect the a p h e n o m e n o l o g i c a l simulator, as the one used here, these parameters.
I
I
I
112
128
144
environment
(R = 2500)
be c o n s i d e r e d as tuning depends on the v a r i a n c e s controlled loop. The iteratively and can be process. In these cases is very useful to a d j u s t
72 A
,-,, I-m o.. "I(.9 O oc -t'I-_.1 _1
A
& A
Z~ Z~
A
A
56 A A
A
A
48
A
A
A
A
A
A
A
A
A
AA
~
Zl
A
A
A
A
A
AA AA
SET POINT
A
40
& A
A
&
A
:E 52 I
"
I
I
I
I
I
I
I
I
32
v
0
W
00000
~0000000000000000000000000000000000 O~
0000000 w -,
'~
16 0 Fig.11
I
I
I
16
32
48
I
I
I
64 80 96 TIME INDEX
I
i
I
112
128
144
O p t i m a l control p e r f o r m a n c e s in the p r e s e n c e m e a s u r e m e n t noises (R = 2500)
of
Adaptive control of a laboratory grinding circuit
257
561 R=I
.-. ¢-
3¢
,,='-
ool I.
56L
:z: t.,9 ::::) 0
48=,.
-,I-
40
I
I
I
I
I
I
I
I
I
I
I
I
R = I00
n_1 "
56-
R=2500
48~ 40 0
I
I
I
16
32
48
64 80 96 T I M E INDEX
I
I
I
112
128
144
Fig.12 Variations of the simulated throughput (clean of measurement noises) under optimal control in a noisy environment, for different R values
CONCLUSIONS This paper evaluates by simulation the feasibility and the performance of a model-based optimal control of grinding processes. Application is performed for a laboratory grinding circuit. The following conclusions can be drawn: A p h e n o m e n o l o g i c a l m o d e l involving the description of the complete particle size distribution can be used as a basis for optimal control, However this model has to be simplified, linearized and discretized to make it s u i t a b l e for the a p p l i c a t i o n of l i n e a r - q u a d r a t i c o p t i m a l control theory, The linearization of the model used in the control low necessitates the adaptation of its parameter in order to account for the mismatch b e t w e e n the model and the circuit when the ore grindability or mill efficiency changes, The adaptive scheme is feasible when a very small number of parameters are recursively estimated; here a single scale factor of the rate of breakage is estimated on-line with the process, The adaptive parameter is used to update the model in the control law as well as in the Kalman filter, The w e i g h t e d c r i t e r i o n used is very useful to tune the controller; the manipulated variable can be weighted accordingly to the amplitude of t h e n o i s e m e a s u r e m e n t to g i v e s a t i s f a c t o r y d y n a m i c s to the controlled grinding circuit, The control scheme is robust and can accommodate strong variations of ore grindability in a highly noisy environment.
258
Y. Dime and D. HODOUIN
The approach which is used is based on a physical model of the grinding mill and as a consequence can be transposed to any other circuit. This opens the door to the design of a general grinding circuit adaptive controller, based on a physical u n d e r s t a n d i n g of g r i n d i n g dynamics, rather than on b l a c k - b o x models. ACKNOWLEDGEMENTS The authors w i s h to thank the N a t u r a l S c i e n c e s Council of Canada for its financial support.
and
Engineering
Research
REFERENCES I. 2. 3.
4. 5. 6. 7. 8. 9. 10.
11. 12. 13.
14. 15. 16. 17. 18. 19.
Foulard, C., Gentil, S., Sandraz, J.P., C o m m a n d e et r e g u l a t i o n par calculateur numerique: de la theorie aux applications, Eyrolles, 1987. Sage, A.P., White, C.C., Optimum Systems Control Prentice-Hall, 1977. Lanthier, R., Hodouin, D., Empirical Models for Multivariate Filtering of Closed G r i n d i n g C i r c u i t s R e a l - T i m e Data. Proceedings of Mathematical Modelling of Materials Processing Operations - Met. Soc. of AIME, 1987, p. 75-90. Rajamani, K., O p t i m a l C o n t r o l of Ball Mill Grinding, Ph.D. Thesis, University of Utah, 1979. Dube, Y., Montage, simulation et algorithme pour la commande automatique du procede de broyage, Ph.D. Thesis, Laval University, 1985. Dube, Y., Lanthier, R., Hodouin, D., Computer-Aided Dynamic Analysis and Control Design for Grinding Circuits, CIM Bulletin, September 1987. Dube, Y., Hodouin, D., D e s i g n of an A d a p t i v e Filter for a Laboratory Grinding Circuit, Proceedings, 5th IFAC Symposium on Automation in Mining Mineral and Metal Processing, Japan, 1986. Hodouin, D., Dube, Y., Lanthier, R., Stochastic Modelling and Simulation of Filtering and Control Strategies for Grinding Circuits, International Journal of Mineral Processing 22 261-274, 1988 Hodouin, D., Berube, M.A., Everell, M.D., Modelling of Twelve Continuous Grinding Experiments on a N e w - B r u n s w i c k S u l f i d e Ore, R e s e a r c h Report, CANMET, Energy, Mines and Resources, Ottawa, 1977. Gupta, V.K., Hodouin, D., Berube, M.A., Everell, M.D., The Estimation of Rate and Breakage Distribution Parameters from Batch Grinding Data for a Complex Pyritic Ore Using a Back-Calculation Method. Powder Technology, 28 (1981) 97-106. Lynch, A.J., The Use of Simulation Models in the Design, Optimization and C o n t r o l of W e t - G r i n d i n g Circuits, A I M E Annual Meeting, San Francisco, California, 1972. Herbst, J.A., Mika, T.S., Rajamani, A., A Comparison of Distributed and L u m p e d P a r a m e t e r Models for O p e n Circuit Grinding Mills, ~th European Symposium of Comminution, Nurnberg, R.P., dechema-monographien, 1975. Garden, R.P., V e r g h e s e , K., T a n k s - i n - s e r i e s T r a n s i e n t Models for the Determination of Model Simulation Parameters in Continuous, Closed-Circuit C o m m i n u t i o n Processes, 4th E u r o p e a n Symposium on Comminution Nurnberg, Dechma-Monographien, 1975. Smith, H.W., Guerin, D., Dynamic Modelling, Simulation and Control of a G r i n d i n g Circuit, 18th A n n u a l C o n f e r e n c e of M e t a l l u r g i s t s , Sudbury, Ontario, Canada, 1979. Austin, L.G., Weller, K.A., Simulation and Scale-up of Wet Ball Milling, XIV International Mineral Processin E Congress, Toronto, Canada, 1982. Herbst, J.A., M o d e l l i n g and Simulation of Crushing Circuit for Design, Control and Optimization, Short Course, Society of Mining Engineers, 1979. Tanaka, T., P r e d i c t e d D y n a m i c R e s p o n s e of H o l d - U p in Continuous Ball Mills, Proc., Australas. Inst. Min. Metall., No. 225, September 1975. Takahashi, Y., Rabins, M.J., Auslander, D.M., Control and Dynamic Systems, Addlson-Wesley, Mass., 1970. Goodwin, G.C., Sin, K.S., A d a p t i v e F i l t e r i n g P r e d i c t i o n and Control, Prentice-Hall, New-Jersey, 1984.
0892-6875/89 $3.00+0.00 © 1989 Pergamon Preu pic
ADAPTIVE CONTROL OF A LABORATORY GRINDING CIRCUIT
Y. DUB~ $ and D. HODOUIN § t
Dept. d'Ing~nierie, Universit~ du Quebec a Trois-Rivi~res, C.P. 500, Trois-Rivi~res, 0u~bec, Canada G9A 5H7 § Dept. de Mines et Metallurgie, Universite Laval, Ste. Foy, Quebec, Canada GIK 7P4
(Received 5 January 1989)
ABSTRACT A multivarlable adaptive control scheme has been studied and evaluated with the help of a phenomenological simulator which was calibrated on a laboratory grinding circuit. The adaptive control scheme includes a Kalman filter, an adaptive mechanism, and a control law. The adaptive mecb-nism is based on the estimation of a fundamental coefficient of the process since a direct identification of the whole set of parameters is impossible due to the large size of the state-space model. The control law that minimizes a quadratic criterion is obtained by dynamic programming. Simulated results show that the control law exhibits good performances. The control scheme is robust and is able to face IOOZ change on grinding rate parameters, without deterioration of the performances. Keywords Computer modeling, Grinding circuit, Adaptive control, Kalman filtering, Quadratic control. INTRODUCTION The main objective of the implementation of a computer control scheme for a c o n t i n u o u s p r o c e s s is u s u a l l y to m i n i m i z e the v a r i a t i o n s of the product quality while minimizing the variations of the control variables. Such a goal can be a c h i e v e d by f o r m u l a t i n g some p e r f o r m a n c e c r i t e r i o n w h i c h can be mathematlcally optimized when a suitable dynamic model of the process to be controlled is available. The performance of the control strategy will depend on both the quality of the selected control scheme and the accuracy of the process model. F o r the s p e c i f i c case of ore g r i n d i n g circuits, automation should help to maintain operating conditions, such as circulating load r a t i o s and p r o d u c t fineness, close to their optimum values while smoothing the variations of the manlpulated variables such as ore feed rate, water addition rate, pump speed. The definition of a performance criterion should includes these requirements with an adequate weighting of the v a r i o u s a s p e c t s . The f o r m u l a t i o n of a q u a d r a t i c o b j e c t i v e f u n c t i o n is the e a s i e s t way t o h a n d l e this problem, because the theory of optimal control is well developed in this context. As the model of the process is usually based on a linear state-space formulation, this approach leads to LQ control algorithms [1,2]. This type of s t r a t e g y will be tested in the present study. T h e r e are two m a i n t e c h n i q u e s to d e v e l o p a d y n a m i c m o d e l of a g r i n d i n g circuit. The empirical approach consists in the description of the process by a transfer matrix identified by suitable experimental design; then the model is reformulated in terms of state-space variables. It has already been shown that this approach is feasible [3]. The second technique, which will be used here, is based on a phenomenological description of the process [4,5]. The d e v e l o p m e n t of the model requires a good understanding of the mechanisms of g r i n d i n g a n d c l a s s i f i c a t i o n and is m o r e c o m p l e x to m a n i p u l a t e that the e m p i r i c a l a p p r o a c h b e c a u s e it leads to h l g h - o r d e r non l i n e a r equations. H o w e v e r the p h e n o m e n o l o g i c a l m o d e l a l l o w s t h e d e f i n i t i o n of p h y s i c a l parameters in which the prior knowledge of the process can be embedded. This 245
246
Y. Dune and D. Ho~tnN
advantage is used here to design an adaptive mechanism to take into account that the process model changes when the circuit efficiency changes. This gives more robustness to the model and as a consequence to the control algorithm. In the paper a phenomenological state-space model is developed, linearized and c a l i b r a t e d for a l a b o r a t o r y g r i n d i n g c i r c u i t involving a ball mill and a screen for the r e c y c l i n g of the c o a r s e p a r t i c l e s to the m i l l feed. The control state feedback law is then established using a quadratic performance index. The states are observed using a Kalman filter. Both the filtering and control algorithms are made adaptive to cope with grinding efficiency drifts. The m e t h o d is e v a l u a t e d u s i n g a g r i n d i n g c i r c u i t s i m u l a t o r i n v o l v i n g a detailed description of the equipment behaviour and including noise generators [6,7,8].
LABORATORY GRINDING CIRCUIT The laboratory grinding circuit is shown in Figure I. The grinding mill is a grate discharge ball mill. The mill inside diameter is 40 cm and its length 40 cm. The mill is equipped with a c y l i n d r i c a l s c r e e n f i x e d on the m i l l output. The fine p a r t i c l e s p a s s i n g t h r o u g h the s c r e e n (finer t h a n 150 micrometers) are analysed by the Leeds and Northrup Microtrac particle size monitor. The coarse particle flow feeds a constant level sump and is recycled to the input of the mill. The measured sump level value is the input to a sump-pump controller and the speed of the pneumatic pump is controlled by an a u t o m a t i c air v a l v e . The r e c y c l i n g s t r e a m f l o w r a t e is m e a s u r e d by a magnetic flowmeter and a densimeter is also included. The circuit is fed by two automatic conveyors which are equipped with load cells and tachometers to measure the feedrate. Their set points can be selected manually or calculated by the computer. The percentage of solids in the mill is normally selected as 65% and maintained constant by varying the water pump speed. The set point of the water pump is given by the computer, based on the fresh ore feedrate and the recycling stream flowrate measurements.
~
Feeders
~~
~Centrltu pgu,mlp °'"
In Diaphragm pump
[I Fine product~ Fig.1
Centrifugal~ pump
Laboratory grinding circuit
As mentioned in the introduction, the control algorithm will be evaluated on a simulator. B a t c h t e s t s and d y n a m i c m e a s u r e m e n t s were made to adjust the parameters of the simulator [9]. Steady-state experiments give the values of the phenomenological parameters of the mill: the fragment size distribution y (breakage function), the rate of breakage of particles (selection f u n c t i o n )
Adaptive control of a laboratory grinding circuit
247
and the hold-up in the mill [10]. Dynamic measurements are used to verify the reliability of the simulator. Grinding tests have shown that a good agreement is obtained between simulated and measured values [5]. PHENOMENOLOGICAL
PROCESS MODEL
The m a t h e m a t i c a l m o d e l of the l a b o r a t o r y g r i n d i n g c i r c u i t is b a s e d discretized p a r t i c l e sizes and g r i n d i n g p a r a m e t e r s . This a p p r o a c h extensively used [references 11 to 16] and only the necessary parameters simulate the laboratory circuit are presented here.
on is to
Hold-up mill model: To determine the volume of water and solids retained in the mill, the m i l l t r a n s p o r t behaviour is assumed to be modelled by three identical interactive tanks in series. The steady-state model for the hold-up is given by the following equations: V t _ 6LQo
+ Vr
K
where Qo is the volumetric flowrate, V t the total pulp volume, V r the residual pulp volume, and K and L two c o e f f i c i e n t s to be d e t e r m i n e d . The d y n a m i c variations of the pulp volumes in the three tanks are given in [17]. Grinding model: Grinding dynamics are represented by two perfect mixers in series where firstorder breakage occurs [5]. It is assumed that the solids are d i v i d e d into c l a s s e s (index K) which corresponds to grinding properties characterized by the values of the breakage and rate functions. A matrix formulation of the equations describing the process is for the class K: K d(H~ mj) dt v
cK :
dt where
(SH) K g
K Hj m~
(I)
H .K
J J dC~
-- M K m~_1 - K K j-1 Mjmj + (B K- I)
J --
Qj-I v 3
CK ( 9-I
- cK ) ]
HK is the hold-up weight of solids of class K in mixer j; 3 H the total solid weight ~Z,kH~);j M the mass flowrate; m the particle size distribution; B the breakage function matrix (distribution of fragments); C the mass of solids per unit volume; SH the selection function matrix, expressed as the mass of particles of each size class which are broken per unit of time; Qj the volumetric flowrate of the mixer j.
Classification model: The c l a s s i f i c a t i o n of the mill o u t p u t is m a d e by a s c r e e n with openings (150 micrometers). The screen is assumed to have zero-order and the fraction of particles recycled to the mill is given by:
constant dynamics
Mrm r = EMom o w h e r e E is the c l a s s i f i c a t i o n matrix (proportion of feed particles in each size class retained on the screen), M o the mass flowrate at mill output, M r the r e c y c l e d flowrate and m o and m r the particle size distributions at the same locations.
248
Y. DUBE and D. HODOUIN
L I N E A R I Z E D S T A T E - S P A C E MODEL The above d e s c r i b e d model which is used in the simulator is continuous and non-linear. H o w e v e r to d e s i g n a d i g i t a l f i l t e r a n d a d i g i t a l c o n t r o l a l g o r i t h m a l i n e a r d i s c r e t e - t i m e model is required. Two kinds of d i s c r e t e models are c o m m o n l y used, ARMA model and state-space model. Since the statespace approach has been selected, the p h e n o m e n o l o g i c a l m o d e l m u s t be discretized, linearized and w r i t t e n under the u s u a l s t o c h a s t i c s t a t e - s p a c e form [18]:
(2)
Xk+ I = PkXk + QkUk + w k Yk = FkXk + GkUk + Vk
where x is the state-space vector, u the control vector, k the time index, P, Q, F, G m a t r i c e s of coefficients, and w and v random white noises. In order to linearize Eq. (I) it is assumed that the mill behaves as a single perfect mixer and that the mill hold-up has zero-order dynamics. Also the m a t e r i a l properties description is s i m p l i f i e d , a s s u m i n g that a m i x t u r e of ores of various g r i n d a b i l i t i e s (classes K) behaves as a m a t e r i a l that has an average grindability. This leads to the following e q u a t i o n for the overall solids: H dm dt
=
Mfmf - Mfm
+ (B - I)SHm
I
w h e r e Mf is the m i l l m a s s f e e d r a t e and m the m i l l o u t p u t p a r t i c l e s i z e distribution. The final state-space form of the d i s c r e t e - t i m e model for the w h o l e circuit is given in r e f e r e n c e [5]. For the present i l l u s t r a t i o n of the method, the state-space vector x has 22 e l e m e n t s . T h e f i r s t 20 p o s i t i o n s r e p r e s e n t the v a r i a t i o n of t h e m i l l d i s c h a r g e p a r t i c l e - s i z e d i s t r i b u t i o n a r o u n d its s t e a d y - s t a t e v a l u e s . T h e l a s t two p o s i t i o n s r e p r e s e n t the v a r i a t i o n (around their steady-state value) of the mill feedrate and of the m i l l f e e d r a t e d e l a y e d by o n e s a m p l i n g i n t e r v a l . T h e a g r e e m e n t of t h i s d i s c r e t e model with the complete n o n - l i n e a r p h e n o m e n o l o g i c a l simulator and the p i l o t - s c a l e d y n a m i c e x p e r i m e n t a l m e a s u r e m e n t s is d i s c u s s e d in the p r e v i o u s l y cited references. KALMAN F I L T E R I N G The Kalman filter is an a l g o r i t h m that m i n i m i z e s the v a r i a n c e of the filtered v a r i a b l e s and r e c o n s t r u c t the s t a t e - s p a c e v e c t o r which is g e n e r a l l y d i f f e r e n t from the vector of the m e a s u r e d variables. The laboratory circuit m e a s u r e d values are the p a r t i c l e size d i s t r i b u t i o n of the screen undersize flow and the f e e d r a t e to t h e m i l l ( r e c y c l e d p l u s f r e s h o r e f l o w r a t e ) . T h e 20 f i r s t e l e m e n t s of t h e s t a t e v e c t o r r e p r e s e n t the m i l l d i s c h a r g e p a r t i c l e size distribution. Thus the Kalman filter has to r e c o n s t r u c t the mill d i s c h a r g e particle size distribution, b a s e d on the c i r c u i t p r o d u c t p a r t i c l e s i z e d i s t r i b u t i o n and has to e s t i m a t e m i n i m u m v a r i a n c e v a l u e s for t h e c i r c u i t product size distribution. The p a r t i c l e size d i s t r i b u t i o n of the circuit fine product (mD) is related to the mill output p a r t i c l e size d i s t r i b u t i o n (mo) by the following equation: Mf mp = ( I Mp
- E) m
where M D is the mass flowrate of the circuit fine product and I the identity m a t r i x . - In the above e x p r e s s i o n Mp and m are related to the part of the state vector x w h i c h contains the mill output p a r t i c l e size d i s t r i b u t i o n (noted x i for the particle size interval i) as follows: m i = msi + x i Mp = Mf~
(i - Ei)(msi + x i)
w h e r e m s is the s t e a d y - s t a t e v a l u e of m. As v a r i a b l e mp is related to the state v a r i a b l e by:
a consequence
the m e a s u r e d
Adaptive control of a laboratory grinding circuit
(I
mpi =
-
249
E i)
(x i + msi) ~(I - Ei)(msi + x i)
= aix i + b i
(3)
Since a i and b i are f u n c t i o n s of x i, this o b s e r v a t i o n e q u a t i o n is not linear. To o v e r c o m e t h i s p r o b l e m , the x i v a l u e in a i and b i is r e p l a c e d b ~ its e s t i m a t e d value given by Eq. 2 w h e r e x k is taken as the filtered value x k at the p r e v i o u s time. T h e second m e a s u r e d value is the r e c y c l e d f l o w r a t e and its v a r i a t i o n around its s t e a d y - s t a t e value is equal to the state v e c t o r p o s i t i o n 21. The overall o b s e r v a t i o n e q u a t i o n can finally be w r i t t e n as: Yk = FkXk + GkUk The Kalman algorithm, here for convenience: Xk+1
which
can
= (PkXk + QkUk ) + Kk+l[Yk+1
be
found
in many
text-books,
is
summarized
- F k + 1 ( P k X k + QkUk)]
Kk+ I = M k + i F ' k + 1 ( F k + i M k + i F ' k + I + Vk+1) -I Mk+ I = PkZkPk , Z k = (I - K k F k ) M k. w h e r e ~ is the filtered given by: W k = E[WkW'k]
;
s t a t e - s p a c e vector
and W and V are d i a g o n a l m a t r i c e s
Vk = E[VkV'k]
T h e s t a n d a r d d e v i a t i o n of the r e c y c l e d f l o w r a t e is selected as 25% of the current s t e a d y - s t a t e value. The standard d e v i a t i o n of the c i r c u i t p r o d u c t particle size distribution is s e l e c t e d as 15% of the current s t e a d y - s t a t e value. These values are used to c o n s t r u c t the m a t r i x V. T h e m a t r i x W is adjusted empirically to g i v e a s u i t a b l e r e d u c t i o n of the v a r i a n c e by the K a l m a n filter. The m a t r i x W value is selected as: W_- I0-5 i
CONTROL A L G O R I T H M The control a l g o r i t h m is based on a m u l t i v a r i a b l e state f e e d b a c k scheme [I, 19] w h e r e the c o n t r o l law is d e s i g n e d for m i n i m i z i n g a q u a d r a t i c o b j e c t i v e function. T h e m e t h o d is f i r s t p r e s e n t e d in a g e n e r a l c a s e , t h e n it is i l l u s t r a t e d by s e l e c t i n g a specific a p p l i c a t i o n to the g r i n d i n g circuit. The o b j e c t i v e f u n c t i o n is of the f o l l o w i n g form: co
J = Z [U'kRU k + e'kTe k] k--1
(4)
w h e r e ek, the d e v i a t i o n of the process output to a set point v e c t o r Z is: ek = Zk - Yk The m i n i m i z a t i o n of the o b j e c t i v e f u n c t i o n is p e r f o r m e d under the c o n s t r a i n t of the process model given by the f o l l o w i n g linear state-space equation: r
Xk+ I = PkXk + QkUk Yk = FkXk + GkUk T h e q u a d r a t i c o b j e c t i v e function (4) does not g u a r a n t e e a zero s t e a d y - s t a t e error. To avoid this effect n u m e r i c a l integrators are added to the initial s t a t e - s p a c e model. They have the f o l l o w i n g form: Vk÷1
-- Vk + ek = Vk + Zk - Yk
250
Y. DUBE and D. HoDoun~
E x p a n d i n g we obtain: Vk+ I = v k - FkX k - GkU k + Zk
i-) (- °I()(°) u (o}
The final augmented system of equations is:
VK+ I
-F k
vk
-G k
zk
() (, :)() (0:) =
Vk
UK
+
O
(5)
Vk
and more compactly: Xk+1
= AkXk ÷ BkUK + fk
Yk = CkXk + DkUk The new objective function is defined as:
3
=
7
k=1 where:
[e'kTek + U'kRU k]
(z:)()
o) ;
~k =
Tv
Vk
The optimal control law m i n i m i z i n g the new cost function is given by: Uk = - LkXk + ~k w h e r e L and ~ are calculated by a dynamic p r o g r a m m i n g method and given in (I):
Lk+ I = (R + B'kKkB k + DkTDK) -I
[B'kKkA k + D'kTC k]
Uk+1
[D'kTZ k - B' k (gk ÷ Kkfk)]
= (R + B'kKkB k + DkTDk) -I
Kk+ I = ~kKk [A k - BkLk+ I ] + ~kCk gk+1
= ~k(gk + Kk [BkUk+1
(6)
+ f]) - #kZk
and: ~k = A'k - C ' k T D k ( R + D'kTDk )-I B'k ~k = C'T - C'TDk(R + D~TDk )-I D' T k k k The matrix L and the vector U are functions of T and R which must be positive definite. L a r g e T values will lead to large v a r i a t i o n s of the control actions u compared to the output error variations. C o n v e r s e l y if R is large c o m p a r e d to ~, the o u t p u t e r r o r v a r i a t i o n s will be large r e l a t i v e l y to the control a c t i o n s u. M a t r i c e s T and R a r e t h e r e f o r e c h o s e n i t e r a t i v e l y to g i v e s a t i s f a c t o r y closed loop dynamics. To i l l u s t r a t e t h i s type of control strategy an a p p l i c a t i o n to the g r i n d i n g process p r e v i o u s l y d e s c r i b e d will be given in the next sections. Generally the control o b j e c t i v e for a g r i n d i n g system is to m a i n t a i n the product size a n d the m i l l t h r o u g h p u t c o n s t a n t u s i n g t h e t w o a v a i l a b l e manipulated variables: the o r e f e e d r a t e and the s u m p w a t e r a d d i t i o n rate. In the present case the water a d d i t i o n to the circuit does not influence either the screen p e r f o r m a n c e or the mill efficiency. This is due to the fact that the
Adaptive control of a laboratory ~'inding circuit
251
percentage of solids in the mill is maintained constant by a local feed-back loop and that the water is added to the sump only to allow a suitable pumping behaviour. Furthermore the screen performance is almost independent of the operating conditions leading to a naturally stable fine p r o d u c t size. The most sensitive operating variable is here the mill throughput which has to be controlled for a stable operation of both the mill and the screen. The amount of recycled solids is a function of the fresh feed particle size distribution and of the ore grindability, variables which are considered as disturbances, and a function of the fresh ore feed rate which is the manipulated variable. Despite the fact that only the mill throughput is controlled by the ore fresh feed rate, the control strategy is based on the whole set of state variables of the process, including o b v i o u s l y the p a r t i c l e size d i s t r i b u t i o n . The extension to the multivariable case is straightforward since the formalism is exactly the same. For instance one could have added in the criteria another o u t p u t v a r i a b l e s of the m o d e l such as some characteristics of the product particle size, with a suitable weighting factor. ADAPTIVE MECHANISM The filtering and control algorithms are based on the matrices A, B, C and D. When the ore characteristics and the o p e r a t i n g levels v a r y the m o d e l parameters change and the matrices A, B, C and D must be continuously updated to keep track of the parameter variations. In the present case, preliminary simulation runs show that the v a r i a t i o n of the g r i n d i n g rate p a r a m e t e r s m o d i f i e s the process dynamics. As a consequence the optimal controller is detuned when the process model changes. Usual real-time adaptive algorithms, such as r e c u r s i v e least s q u a r e s e s t i m a t i o n of the model matrix elements, cannot be considered here due to the very large number of parameters involved in the matrices A, B, C and D. The method selected is based on a key parameter of the phenomenological model which is estimated using the deviation between the experimental values and the calculated results. It is assumed that the main variations of the mill model can be taken into account by a modification of the selection f u n c t i o n by a multiplicative constant Uk, such that at time k: S k = UkSo The constant u k modifies by the same factor the rate of breakage of all the particles, independently of their size. So it characterizes the grinding mill efficiency and lumps all the effects due to the grlndabillty and partlcle size distribution v a r i a t i o n s of the f r e s h feed as w e l l as the e f f e c t s of the v a r i a t i o n of the o p e r a t i n g c o n d i t i o n s such as the ball load and the pulp theology. This index is recursively calibrated using the difference between the discrete model prediction and the experimental v~lue (here the simulated value) of the mill throughput (respectively M~ and M~). The ek index is given by: P ~k+1
=
Uk
+
8(Mf
-
e Mf)
The speed of the adaptive algorithm is adjusted by the parameter B which empirically selected to give suitable dynamics to the estimation process.
is
The parameter ~ does not explicitly appear in the model matrices A, B, C and D. However it is formally present in the original phenomenological model. As a c o n s e q u e n c e , the m a t r i c e s A, B, C and D h a v e to be rederived from the p h e n o m e n o l o g i c a l m o d e l each time that the p a r a m e t e r ~ is updated by the adaptive mechanism. SIMULATED RESULTS In a first step the adaptation algorithm will be tested. In the absence of disturbances the steady-state value of the simulated mill throughput is 47.5 kg/h for a feedrate of 20 kg/h. At time index k = 0 (the sampling period is 18 seconds) the f r e s h ore g r i n d a b i l i t y is d o u b l e d . This d i s t u r b a n c e is equivalent to a multiplication of the selection function by a factor 2. The fresh ore feedrate is also perturbated according to the scheme in Figure 2.
252
Y, Dime and D. HODOUn~
F u r t h e r m o r e , to compare easily the responses of the discrete model and the adaptive discrete model, the simulation of the measured variable is performed without addition of simulated measurement noises. The comparison between the mill throughput predicted by the discrete model without adaptation mechanism and the mill throughput given by the simulator is given in Figure 2. Figure 3 gives the response of the mill throughput predicted by the adaptive discrete model and by the simulator. A brief comparison of these curves clearly shows the advantage of the adaptive procedure. 52
44 GD
v
36 LLI '~ r",, 2 8 w L~J h 20
'1 12
I
i
I i
'
'
I
I
I
I
52 SIMULATOR DISCRETE MODEL
L--
44
v
g
3s
1"
~ 2s
......f
-~ 20 .J 12 0
Fig.2
I
I
40
80
I
I
120 TIME
I
160 200 INOEX
I
I
240
280
320
Predicted mill throughput without model adaptation
In a second step the performance of the control algorithm depicted in Figure 4 is studied. This algorithm must keep the mill throughput at a constant set point using fresh ore fe~drate as control variable. The p e r f o r m a n c e index given by Eq. (4) combines the squared variations of the mill throughput and of the fresh ore feedrate. The state-space model is obtained by substituting in Eq. (5) the matrices P, Q, F and G corresponding to the state and observation equations of the grinding circuit. 52
SIMULATOR ------
44
DISCRETE
MODEL
I-
~..r ss ~ as "1" I.-
J
J
f
20
12 0
I
I
I
I
40
80
120
160
TIME
Fig.3
I
200
I
240
I
280
300
INDEX
Predicted mill throughput with model adaptation
Adaptive control of a laboratory grinding circuit
grlndiblllty dlitur blnce Fresh ore feed rite
.J
253
miilurement nolle
LABORATORY I Mill
throughput
I I OItINDIHO S M IULATORCIRCUITI - Pr°duc' "eriicl" ~ I slse
KALMAN FILTER Filtered states
I
LINEARIZED
-=,
I
mill
I [
O,T,MA,
I CALCULATIONI-oF B, C I ANo DA'MATRICESI -
Fig.4
DISCRETE
dlltrlbuUon
"
II
1.
I
I/+"
CsIculitldP
throughput
ADAPTIVEI_
MECHANISM
Block diagram of the control strategy
The o b j e c t i v e of the control strategy is to keep the v a r i a t i o n of the m i l l t h r o u g h p u t a r o u n d its s t e a d y - s t a t e value M s as close as p o s s i b l e to a setpoint value. Using the p r e v i o u s l y defined notation the variable: M
-
M s
=
x(21)
÷
u
should track the set point value. T h e o p t i m a l c o n t r o l law, g i v e n in Eq. 6, is now tested w e i g h t i n g factors of the p e r f o r m a n c e c r i t e r i o n selected as:
T =
(2:) ;
by
simulation
for
R = I
0 T h e d i s t u r b a n c e w h i c h is g e n e r a t e d in the simulator is a change of the ore g r i n d a b i l i t y which, for a s t e a d y - s t a t e operation, w o u l d be e q u i v a l e n t to a m u l t i p l i c a t i o n of the s e l e c t i o n f u n c t i o n by a factor 2. F i g u r e 5 shows the r e s p o n s e of the u n c o n t r o l l e d g r i n d i n g c i r c u i t to t h i s g r i n d a b i l i t y d i s t u r b a n c e w h i l e the fresh ore feedrate is m a i n t a i n e d at a value of 20 kg/h. A strong v a r i a t i o n of the mill throughput can be observed. Under control, F i g u r e 6 shows the optimal t r a j e c t o r y of u and the v a r i a t i o n of the mill t h r o u g h p u t around the set point of 47.5 kg/h (steady state value, for a fresh ore feedrate of 20 kg/h, before d i s t u r b i n g the ore grindability). This c a l c u l a t i o n is p e r f o r m e d w i t h o u t any s i m u l a t i o n of m e a s u r e m e n t n o i s e . The c o n t r o l a l g o r i t h m is very e f f i c i e n t in the absence of noise measurement, to m a i n t a i n the mill throughput at the set point value.
50 46 v
42
~: (9
38
jo. o
34
t-.J -'2 3O 26
0
Fig.5
ME 2/2--O
I 16
I 32
J 48
J 64 TIME
I I 80 96 INDEX
112
1 2i 8
1 4i 4
U n c o n t r o l l e d mill t h r o u g h p u t r e s p o n s e to a g r i n d a b i l i t y
increase step
254
Y. Dum~ and D. HODOUIN
4eI ..........
SET
POINT
~44 -
32
24
2O
Fig.6
0
16
Optimal
32
control
I I 12
I
48
64 80 96 TIME INDEX
performances
in a n o i s e - f r e e
I 128
I 144
environment
(R = I)
When noises, with the s t a t i s t i c a l c h a r a c t e r i s t i c s mentioned in t h e K a l m a n filtering section, are simulated and added to the m e a s u r e d variables, the optimal control results of F i g u r e 7 are obtained. Very strong v a r i a t i o n s of the manipulated v a r i a b l e can be o b s e r v e d w h i c h are not a c c e p t a b l e for many p r a c t i c a l reasons. To reduce their a m p l i t u d e a larger value can be given to R. F i g u r e 8 s h o w s the b e h a v i o u r of the control scheme for R = 100 w i t h o u t noise a d d i t i o n and Figure 9 its b e h a v i o u r in the p r e s e n c e of noise. One can s e e f r o m F i g u r e 8 r e s u l t s t h a t the v a r i a t i o n s of the mill t h r o u g h p u t are larger than the v a r i a t i o n s w h e n R w a s g i v e n a v a l u e I. However smaller v a r i a t i o n s of fresh ore feedrate are needed in the p r e s e n c e of noises. Figure 10 and 11 show the responses for R = 2500 and similar remarks can be made. 80
A
70
A
~< 6o
A
A A
A
I-
~.
A
A
A A
A
A
50
AA
A
A
;ET POINT o
A
A
A
~.
A
~
A
A
6
A
/,
A
A
A A
A
A
A
40
A ~A A
A
A A
i•J ,.1
3O
A
i
A
I
20 ~
i
I 0
30
0 0
o
0
I
I
I
00
o
0
o
o
0 0
0
o
oo
o
o0
o O0 000
00 0
0
0
0
00
0
0
0
0
0
0
0
0
0 I( 0
Fig.7
I
o
o
0 0
o w w "
I
40 o
~
I
I 16
Optimal
I 32
I 48
| 64
I 80
TIME
INDEX
I 96
I 112
I 128
I 144
control p e r f o r m a n c e s in the p r e s e n c e m e a s u r e m e n t noises (R = I)
of
Adaptive control of a laboratory grinding circuit
255
I,-
~,..56 [ - ~ - za. 'LR
" _1
SET POINT
40 52 24
~' uJ
16 o bJ LIJ "
8
I 16
0
I 32
I 48
I 64
I 80
TIME Fig.8
O p t i m a l control p e r f o r m a n c e s
I 96
I 112
I 128
I 144
INDEX
in a n o i s e - f r e e e n v i r o n m e n t
(R = 500)
In the p r e s e n c e of simulated m e a s u r e m e n t noises, it is d i f f i c u l t to see the effect of R on the true mill t h r o u g h p u t due to the large standard d e v i a t i o n of the m e a s u r e m e n t noises. To clarify this aspect of the control strategy, the r e s u l t s are p r e s e n t e d in F i g u r e 12 w i t h o u t the m e a s u r e m e n t noise, w h i l e the optimal control law is c a l c u l a t e d using t h e n o i s y m e a s u r e m e n t of the m i l l t h r o u g h p u t a n d p r o d u c t p a r t i c l e size d i s t r i b u t i o n . This is p o s s i b l e o n l y b e c a u s e the v a r i a b l e s are c a l c u l a t e d by a simulator then d i s t u r b e d by a noise generator. A m i n i m u m v a r i a t i o n for the mill t h r o u g h p u t is o b t a i n e d for R =100. For R = I the control a l g o r i t h m is v e r y s e n s i t i v e to the n o i s e a n d p r o d u c e s large v a r i a t i o n s on control variable. For R = 2500 the v a r i a t i o n s of the mill t h r o u g h p u t are slowly corrected.
72 A v
64
-
A
A
A
A
A
I= 0. "I(9
56
o(¢
48
-J _1
40
A
A A
A
A
-r" I-
A
A
A
A
/&
A
.,.
A
A
~&~ A A
A
A
A
:E
~
SET POINT
A
A
A
A
I
32
t~
b'l"3 2 J
0
24 ua'-"
A
A A
A
A
m
~.1 u.
A
51
A
A
A
I
I
I
00000000000000000
I
I
I
000000000000000000000
oO
L_oo,.o
i., o IRJ-o
v
j/6
i
32
I
48
i
64 TIME
Fig.9
I
80
I
96
I
112
I
128
I
144
INDEX
O p t i m a l control p e r f o r m a n c e s in the presence of m e a s u r e m e n t noises (R = 500)
256
Y. DUBEand D. HODOUIN
k-
Q-
56,
f-
~:~'r 40
I
t
I
I
I
I
I
I
I
I
I
I
32
24
""
I,LI
oc
a I.LI bJ h
16
8
0 Fig.10
Optimal
I
I
I
16
32
48
control
64 80 96 T I M E INDEX
performances
in a n o i s e - f r e e
These s i m u l a t i o n tests show that R, and also Q, can p a r a m e t e r s of the control law. The choice of Q and R of t h e n o i s e s a n d on the p r o c e s s d y n a m i c s of the relative c h o i c e of t h e s e s m a t r i c e s m u s t be d o n e d i f f i c u l t when noises with large v a r i a n c e s affect the a p h e n o m e n o l o g i c a l simulator, as the one used here, these parameters.
I
I
I
112
128
144
environment
(R = 2500)
be c o n s i d e r e d as tuning depends on the v a r i a n c e s controlled loop. The iteratively and can be process. In these cases is very useful to a d j u s t
72 A
,-,, I-m o.. "I(.9 O oc -t'I-_.1 _1
A
& A
Z~ Z~
A
A
56 A A
A
A
48
A
A
A
A
A
A
A
A
A
AA
~
Zl
A
A
A
A
A
AA AA
SET POINT
A
40
& A
A
&
A
:E 52 I
"
I
I
I
I
I
I
I
I
32
v
0
W
00000
~0000000000000000000000000000000000 O~
0000000 w -,
'~
16 0 Fig.11
I
I
I
16
32
48
I
I
I
64 80 96 TIME INDEX
I
i
I
112
128
144
O p t i m a l control p e r f o r m a n c e s in the p r e s e n c e m e a s u r e m e n t noises (R = 2500)
of
Adaptive control of a laboratory grinding circuit
257
561 R=I
.-. ¢-
3¢
,,='-
ool I.
56L
:z: t.,9 ::::) 0
48=,.
-,I-
40
I
I
I
I
I
I
I
I
I
I
I
I
R = I00
n_1 "
56-
R=2500
48~ 40 0
I
I
I
16
32
48
64 80 96 T I M E INDEX
I
I
I
112
128
144
Fig.12 Variations of the simulated throughput (clean of measurement noises) under optimal control in a noisy environment, for different R values
CONCLUSIONS This paper evaluates by simulation the feasibility and the performance of a model-based optimal control of grinding processes. Application is performed for a laboratory grinding circuit. The following conclusions can be drawn: A p h e n o m e n o l o g i c a l m o d e l involving the description of the complete particle size distribution can be used as a basis for optimal control, However this model has to be simplified, linearized and discretized to make it s u i t a b l e for the a p p l i c a t i o n of l i n e a r - q u a d r a t i c o p t i m a l control theory, The linearization of the model used in the control low necessitates the adaptation of its parameter in order to account for the mismatch b e t w e e n the model and the circuit when the ore grindability or mill efficiency changes, The adaptive scheme is feasible when a very small number of parameters are recursively estimated; here a single scale factor of the rate of breakage is estimated on-line with the process, The adaptive parameter is used to update the model in the control law as well as in the Kalman filter, The w e i g h t e d c r i t e r i o n used is very useful to tune the controller; the manipulated variable can be weighted accordingly to the amplitude of t h e n o i s e m e a s u r e m e n t to g i v e s a t i s f a c t o r y d y n a m i c s to the controlled grinding circuit, The control scheme is robust and can accommodate strong variations of ore grindability in a highly noisy environment.
258
Y. Dime and D. HODOUIN
The approach which is used is based on a physical model of the grinding mill and as a consequence can be transposed to any other circuit. This opens the door to the design of a general grinding circuit adaptive controller, based on a physical u n d e r s t a n d i n g of g r i n d i n g dynamics, rather than on b l a c k - b o x models. ACKNOWLEDGEMENTS The authors w i s h to thank the N a t u r a l S c i e n c e s Council of Canada for its financial support.
and
Engineering
Research
REFERENCES I. 2. 3.
4. 5. 6. 7. 8. 9. 10.
11. 12. 13.
14. 15. 16. 17. 18. 19.
Foulard, C., Gentil, S., Sandraz, J.P., C o m m a n d e et r e g u l a t i o n par calculateur numerique: de la theorie aux applications, Eyrolles, 1987. Sage, A.P., White, C.C., Optimum Systems Control Prentice-Hall, 1977. Lanthier, R., Hodouin, D., Empirical Models for Multivariate Filtering of Closed G r i n d i n g C i r c u i t s R e a l - T i m e Data. Proceedings of Mathematical Modelling of Materials Processing Operations - Met. Soc. of AIME, 1987, p. 75-90. Rajamani, K., O p t i m a l C o n t r o l of Ball Mill Grinding, Ph.D. Thesis, University of Utah, 1979. Dube, Y., Montage, simulation et algorithme pour la commande automatique du procede de broyage, Ph.D. Thesis, Laval University, 1985. Dube, Y., Lanthier, R., Hodouin, D., Computer-Aided Dynamic Analysis and Control Design for Grinding Circuits, CIM Bulletin, September 1987. Dube, Y., Hodouin, D., D e s i g n of an A d a p t i v e Filter for a Laboratory Grinding Circuit, Proceedings, 5th IFAC Symposium on Automation in Mining Mineral and Metal Processing, Japan, 1986. Hodouin, D., Dube, Y., Lanthier, R., Stochastic Modelling and Simulation of Filtering and Control Strategies for Grinding Circuits, International Journal of Mineral Processing 22 261-274, 1988 Hodouin, D., Berube, M.A., Everell, M.D., Modelling of Twelve Continuous Grinding Experiments on a N e w - B r u n s w i c k S u l f i d e Ore, R e s e a r c h Report, CANMET, Energy, Mines and Resources, Ottawa, 1977. Gupta, V.K., Hodouin, D., Berube, M.A., Everell, M.D., The Estimation of Rate and Breakage Distribution Parameters from Batch Grinding Data for a Complex Pyritic Ore Using a Back-Calculation Method. Powder Technology, 28 (1981) 97-106. Lynch, A.J., The Use of Simulation Models in the Design, Optimization and C o n t r o l of W e t - G r i n d i n g Circuits, A I M E Annual Meeting, San Francisco, California, 1972. Herbst, J.A., Mika, T.S., Rajamani, A., A Comparison of Distributed and L u m p e d P a r a m e t e r Models for O p e n Circuit Grinding Mills, ~th European Symposium of Comminution, Nurnberg, R.P., dechema-monographien, 1975. Garden, R.P., V e r g h e s e , K., T a n k s - i n - s e r i e s T r a n s i e n t Models for the Determination of Model Simulation Parameters in Continuous, Closed-Circuit C o m m i n u t i o n Processes, 4th E u r o p e a n Symposium on Comminution Nurnberg, Dechma-Monographien, 1975. Smith, H.W., Guerin, D., Dynamic Modelling, Simulation and Control of a G r i n d i n g Circuit, 18th A n n u a l C o n f e r e n c e of M e t a l l u r g i s t s , Sudbury, Ontario, Canada, 1979. Austin, L.G., Weller, K.A., Simulation and Scale-up of Wet Ball Milling, XIV International Mineral Processin E Congress, Toronto, Canada, 1982. Herbst, J.A., M o d e l l i n g and Simulation of Crushing Circuit for Design, Control and Optimization, Short Course, Society of Mining Engineers, 1979. Tanaka, T., P r e d i c t e d D y n a m i c R e s p o n s e of H o l d - U p in Continuous Ball Mills, Proc., Australas. Inst. Min. Metall., No. 225, September 1975. Takahashi, Y., Rabins, M.J., Auslander, D.M., Control and Dynamic Systems, Addlson-Wesley, Mass., 1970. Goodwin, G.C., Sin, K.S., A d a p t i v e F i l t e r i n g P r e d i c t i o n and Control, Prentice-Hall, New-Jersey, 1984.