A Dynamic Model for a Class of Semi-Autogenous Mill Systems

In this paper, we develop a state space model of semi-autogenous (SAG) mill systems, which is of great significance to the mineral processing control. According to the practical enterprise production condition, the model is composed of two main parts. One is the ore circuit model. Compared with other models of ore circuit, the present model is simple and efficient. Another part is the SAG subsystem model where the state space model is utilized to depict the grinding mechanism. Furthermore, the grinding process of hold-ups are analyzed and presented with the common outputs such as power draw, pulp density and pulp flow. We consider the mass of rocks, granules and water as states and analyze the observability of the SAG subsystem model. For the ore circuit system and the SAG subsystem established respectively, unknown parameters of these models are identified by appropriate parameter estimation methods. Verified by the actual plant data, the result shows that the SAG mill system works well.


I. INTRODUCTION
With the development of mineral resources and the increasing demand for metals, the quantity of copper ore results in a relative request of preparation equipment. Because of the simple process, convenient configuration and low investment, the semi-autogenous (SAG) mill has been adopted increasingly in mines. The prominent feature of SAG mill is that copper rocks can be crushed by the strong impact on each other. As a rule, a little of water is added into the SAG mill to improve efficiency. Since the crushing process is a complex physical process which is related to the resistance strength, hardness, shape, size, humidity, density, particle size distribution and other factors, an applicable model of SAG system is hard to establish. The connections of variables and parameters has provided a remarkable challenge.
The solid system in the SAG mill cylinder is composed of a large number of discrete heterogeneous particles. The discrete element method (DEM) first proposed by Cundall has become a standard tool for processing such particulate matter. Mishra used DEM to analyse the trajectory of individual particles as well as contact force and the distribution The associate editor coordinating the review of this manuscript and approving it for publication was Lei Jiao. of energy [1]. The ball mill was further studied and DEM was adopted to predict the power [2]. Due to the large industrial equipment and complex distribution, DEM takes too much time to simulate the grinding process. In addition, the detailed granularity data is difficult to obtain. Therefore, it is particularly important to propose a macroscopic model analysis.
The dynamic model of SAG mill has attracted many researchers' attention. Apelt et al. made enormous contributions to the model and control for the SAG mill [3]- [6]. Considering the feed ore size distribution, Amestica et al. presented a dynamic model by taking the masses of water, copper rocks and granules as state variables [7], [8]. In [9], the control model was developed with taking the delivery rate of copper rocks, the feed flow rate of water and the rotation speed as control inputs and considering the power, the filling rate and particle size reduction as output variables. However, there is no sensor to measure the particle size in the SAG mill and the influence of rotation speed on the crushing rate is not considered. In [10], a class of SAG mill model was proposed in consideration of the influence of steel ball. Due to lacking of mining instrument and meter to measure the corresponding discharge, this method is not applicable. In [11], a relatively complete dynamic model was set up. The parameters need to be adjusted for different SAG mill models. In power model and dynamic mathematical model, it was not taken into consideration that the volume of charge change with shoulder angle and toe angle in [12], [13]. When parameters need to be estimated, observability analysis of the model is necessary for the identification of the parameters. Whereas, in most of the existing literature, the observability was often neglected.
A variety of parameters in the SAG mill have been subjected to intensive studies. For example, the crushing rate is a crucial importance parameter in the operational process of SAG mill. According to the characteristics of ore, the total energy requirements of three common crushing machines were calculated in [14], which showed that the specific energy of SAG mill grinding is relatively high. The ore crushing model can be described by differential equations [11]. In the process of model simulation, errors accumulating with low accuracy was solved by matrix of differential equation of ore grinding in [15]. For the wet SAG mill, water was added into the cylindrical shell to drive out the copper granules smaller than the grid size. In [16], the retention of ore was expressed by linear fitting equation considering the parameters of grate open area, the location of apertures, mill speed and mill diameter. In [17], the influence of slurry discharge on power was analyzed, based which we know that slurry has little influence on power in the rolling mill. Latchireddi et al. studied the effect of grate's property, load and rotation speed on retention, then solved the complicated flow problem in rotating flow field by modeling the pulp hoist with fluid dynamics [18], [19]. In [20], time delay neural network was utilized to estimate particle sizes, which has a marked impact on fault diagnosis in the SAG grinding machine. In [21], the effect of mineral hardness on particle shape was analyzed by the detailed industrial data. The work index of characterizing the ore discharge rate affects the grade of copper ore, but the relation is difficult to be described by mathematical expression. In [22], the Rosin Rammler distribution was adopted to predict the dynamic function of particle size distribution in the SAG mill. In [23], the global sensitivity analysis is utilized to study the effect of the particle size distribution and the input variable on the working condition. By the analysis method, the particle size connection with other parameters was explored via data mining rather than mechanism analysis. As such, there are some difficulties to control the SAG subsystem by classical algorithms.
From the existing literature about dynamic model of the SAG mill system, the equipment in most studies included only the SAG mill but not the ore circuit. To the best of the authors' knowledge, this paper is the first study in the existing literature to consider the whole production line of copper ore pulp. Moreover, another contribution is confirming that the break rate of copper ore serves as a factor indirectly related to rotation speed of the SAG mill and the relationship is formulated. In this paper, considering the SAG mill system located in Inner Mongolia, we establish the mathematical model for the SAG mill system. First, we presents the grinding circuit and the model of ore feed system in Section II and Section III. For the ore feeding system, a simple and effective model is given by the superposition of ore amount in the belt. We present the movement model of hold-ups within the SAG mill cylinder in Section IV. In Section V, based on mechanism analysis, special internal variables are selected to describe the measurement output for the SAG mill. In Section VI, the observability of SAG mill model is analyzed and unknown parameters are identified. Illustrative examples are presented in Section VII. Finally, we present the discussion and conclusion in Section VIII.

II. CIRCUIT DESCRIPTION
The grinding circuit to be considered consists of an ore bin, heavy-apron feeders, a SAG mill, vibrating screens and a pulp reservoir etc. This research is focused on the primary closedloop grinding circuit as shown in Fig. 1. The fresh copper rocks are transported to the bin by trucks. There are three heavy-apron feeders under the ore bin, which are used to adjust the ore feeding rate by their motor frequencies. The operational mechanism is that the middle heavyapron feeder is normally operated as a fine-tuned equipment and one of the heavy-apron feeders on the side is operated as the coarse adjustment. Ore bins above every heavyapron feeder are installed with level indicators. Furthermore, the operating conditions of the three heavy-apron feeders are as follows • The bilateral heavy-apron feeders are in operation alternatively.
• Each heavy-apron feeder is equipped with frequency detector to monitor the motor frequency.
According to the actual operation of the copper factory, the heavy-apron feeder NO. 1 and NO. 2 are working. The fresh ore is transferred to the ore feed belt NO. 1 from the heavy-apron feeder NO. 1 and NO. 2. Subsequently, the fresh ore is mixed with unqualified ore which is conveyed by the recycle ore feed belt NO. 3. The fresh copper ore and unqualified ore transferred from the ore feed belt NO. 1 and NO. 3 are sent into the SAG mill by the ore feed belt NO. 2. Ultimately, the pulp and unqualified ore are separated by the VOLUME 8, 2020 vibrating screen. The liquidometer is installed on the pump sump.
All available measurements are as follows: • Frequency of the heavy-apron feeder NO.

III. ORE CIRCUIT MODEL
Since the copper rock is filled in the bin and piled up above the heavy plate, movement of the heavy-apron feeder plates drives the ore down. According to the structure feature of motor and electric machine theory, the speed of motor is proportional to the frequency. Taking delays into consideration, the relationship between the ore feed rate of the belt NO. 1 in the position of weighing transducer and motor frequencies of heavy-apron feeders can be given by where b 1 is the feed rate of fresh ore [t/s], f 1 is the motor frequency of the heavy-apron feeder NO. 1 [Hz], f 2 is the frequency of another heavy-apron feeder [Hz], τ 1 and τ 2 are time delays of the heavy-apron feeder NO. 1 and NO. 2, respectively [s], a 1 and a 2 are constants to be identified later.
In order to adjust the charge mass flow rate effectively in the actual production, frequencies of the ore feed belt NO. 1 and NO. 2 are normally fixed. In addition, the ore feed belt NO. 1 is distinct from the recycle ore feed belt NO. 3 in speed and belt width. Hence, the mass flow rate model of ore feed belt NO. 2 can be established in the similar way as follows

IV. SAG SUBSYSTEM MODEL
The internal mechanism model of SAG mill is established based on several commonly used observations. In particular, for the power model of SAG mill, we carry on the detailed analysis to the movement track of the internal material.

A. PULP DISCHARGE
Apelt et al. made a verification that the discharge flow rate with no pump sump is proportional to the square of the internal load [3]. For the SAG mill, the volume of pulp flow can be written as where x s is the mass of the coper granules which can pass through the SAG mill discharge grate and always transform into the pulp mixing with water [t], x w is the mass of water [t], Q is the output volume flow rate of pulp [m 3 /s], η is a constant, ρ s is the average density of copper rocks [t/m 3 ], ρ w is the water density [t/m 3 ].

B. DENSITY
Pulp density is an important physical quantity which reflects grinding characteristics of the copper mine. The concentration of discharge pulp is given by where C is the discharge pulp density [%].

C. POWER DRAW
The power draw P [KW ] includes not only the idle power P 0 [KW ] but also the load power P 1 [KW ] which is mainly affected by the load and the rotation speed of the SAG mill cylinder shell. Hence, the power can be represented by The SAG mill is operated under three conditions theoretically like Fig. 2. Whereas, it is always operated in the second fall mode in practical production. The cross section of the cylinder shows that each layer movement of hold-ups consists of the ore rocks, granules and water as shown in Fig. 3. For the SAG mill model, there are some assumptions as follows • The hold-ups will not drift to other grinding layers during the movement.
• Speeds of the hold-ups are reduced to zero at point B shown in Fig. 3.
• The rotational speed remains the same in a rotational period.
• The energy consumption caused by friction is negligible. Apparently, when the radius equals to r shown in Fig. 3, the maximal drop height H r is given by Combining with the force analysis of hold-ups at point A when the radius equals to r, the expression of centripetal force is expressed as where m is the mass of the hold-ups at point , ω is the angular velocity of the hold-ups [rad/s], v is the tangential velocity [m/s]. By simplifying (7), the relationship between the shoulder angle and the tangential velocity of the hold-ups is obtained by v 2 = gr cos α.
Substituting the common formula of circular motion v = π rn 30 into (8), the relation between the radical position and the charge shoulder angle is given by where n is the rotation speed [r/min]. When the rotational speed is fixed, the shoulder angle depends on the radical position of charge inner surface radius of mill inside liners. The movement locus of hold-ups is analyzed in rectangular coordinates after the the hold-ups leaving point A shown in Fig. 3. Building the two-dimensional reference frame with the origin at point A, we can get trajectory equations of the parabola and the circle such that where x and y are the horizontal and vertical coordinates. From the simultaneous equations (10a) and (10b), the set of coordinate values at point B is given by where X r B is the horizontal coordinate, Y r B is the vertical coordinate.
Applying Newton's Second Law, the lift height with the radial distance equal to r can be governed by Substituting (11b) and (12) into (6), we obtain the maximum of falling height when the radius equals to r as follows Moreover, we can get Therefore the numerical calculation of the minimum radial position of charge inner surface is gained by dX r B dt = 0. By simple manipulation, the result is given by By solving (15) and substituting the result into (9), the minimum radius can be expressed as where R min is the minimum radical position. From the view of cross section shown in Fig. 4, when one of the layers is under the optimum condition, the others may not be under the optimum conditions. Hence, the polycondensation layer is introduced to concentrate masses of all layers VOLUME 8, 2020 to one layer. The radius R of the polycondensation layer is given by where R max is the maximum radius [m].
The maximum falling height of the polycondensation layer is expressed as Since different layers have different rotational periods, we utilize the period of the ploycondensation layer to represent whole. The rotational period t R 0 is the time that it takes to revolve when the radius equals to R in (16). The SAG mill drives the hold-ups to revolve during the rotational period. Hence, the loading power is expressed as where E is the energy provided the hold-ups to revolve by the SAG motor [KJ ], t R 0 is the rotational period [s]. When the cylinder rotates, the impact energy E of the hold-ups for a certain paraboloid is given by where v B is the velocity at point B [m/s], dm is the differential form of the layer mass [t] which can be described by where σ is the line mass [t/m 2 ], ω is the angular speed of the barrel [rad/s], t 2 is the movement time in circular when the radius equals to r [s]. The total mass of mixture is related to the two time periods. One is the movement time in a cycle period t 1 and the other is the duration time in the paraboloid t 2 . As shown in Fig. 5, because t 1 and t 2 are the movement time periods when the radius equals to r, we utilize the time in every layer to express the whole time period which is closely related to the mass in the circular surface and paraboloid. We utilize a ratio to represent the relationship as below where M 1 is the total mass in the circular surface [t], t 1 is the movement time period in the paraboloid [s]. M is the total mass of hold-ups in the cylinder of the SAG mill [t], i.e.
where x r is mass of copper rocks in the SAG mill [t].
The coordinate of B (11b) is substituted into the trigonometric function of β. The relationship between the shoulder angle and toe angle is given by By simplifying (22), the following relationship is obtained By the assumption that the rotational speed remains the same in a rotational period, we can obtain the duration when the particle moves along the circle as follows Substituting the coordinate (11a) into simple kinematics formulas, the paraboloid time is shown as The total time in a rotational period can be given by When the radius r equals to R in (9), the shoulder angle is rewritten as α = arccos Rn 2 π 2 900g .

Substituting (33), (29) and (25) into
Ultimately, the power draw can be written as where λ is a time-varying parameter. Due the computational process is readily comprehensible but too rigmarole, we can facilely solve with the MATLAB program.

V. DYNAMIC MODEL A. DYNAMIC MODEL OF SAG MILL
In this subsection, we consider the case that the coper ore is classified into two types: the rocks and granules. The granules rather than the rocks can pass through the discharge grate of SAG mill. In this case, we can separate them individually. Thus, the masses of rocks, granules and water in the cylindrical shell will be considered as states.

1) COPPER ROCKS
Copper ore is ground by collision with itself and the wearresistant lining. The mass of copper rocks in the SAG mill cylinder x r can be described as a first-order system as follows where u r is the feed rate of copper rocks [t/s], K r is the crushing rate of copper rocks [s −1 ], q n is the unqualified ore [t/s]. The crushing rate of copper rocks K r is affected by the maximum of falling height in the ploycondensation layer H . In (9), the shoulder angle depends on the rotation speed when the radius equals to r. In (17), the maximum value of H can be gained by dH dα = 0 such that dH dα = d 4.5Rsin 2 α cos α dα = 0.
When H reaches the maximum H max , α = 0.304π . Substituting (9) into (17), the relationship between the maximum of falling height in the ploycondensation layer H and the rotate n can be given as 4.5π 2 n 2 R 2 30 4 g 2 − π 4 n 4 R 2 30 6 g 3 .
(36) VOLUME 8, 2020 Hence, whereK r is the maximal crushing rate of copper ore, γ is a time-varying coefficient which is denoted by

2) GRANULES
The granule is the solid particle which can pass through the discharge grate of SAG mill. The granule formed of the ground copper rocks is blended with the water, and eventually forms the pulp. Hence, the change of granules can be expressed as where c s is the discharge rate of granules [t/s].

3) WATER
The content of water varies with input and output flow rates. The dynamic model of water is given by where u w is the input rate of water [t/s], c w is the water discharge rate accompanying with the pulp [t/s]. The outflow rates of water and granules are described in detail [10]. They can be shown as

B. STATE SPACE EXPRESSION
The state space model is used to describe dynamic changes of internal states including the mass of water, copper rocks and granules in the SAG mill. Denote Combining (35), (38) and (39), the state space model can be expressed as By integrating (3), (4) and (34), the output equation is expressed as (41)

VI. PARAMETER IDENTIFICATION A. ORE FEEDING SYSTEM
Recursive least-square method is used to identify the parameters a 1 , a 2 , a 3 and a 4 . Combining (1) and (2), we obtain that The ore feeding system (42) is described in the following compact form The criterion function is chosen as (43) According to the principle of least square method, the recursive formula of the parameter estimation of the multiple input multiple output system is given by where K is the gain of correction, L is the recursive factor and I is the identity matrix.

B. SAG SUBSYSTEM
Since the values of ρ s and ρ w are gained by the experimental measurement, parametersK r and η in the dynamic model can be identified by extend Kalman filter. By discretizing (40) and (41), the discrete system model is obtained by and Observability represents the ability that the system state can be fully reflected by the output. In order to estimate parameters in (45) and (46), we incorporateK r and η into the state and obtain the following augmented system. Denote After system augmentation, the system matrix is rewritten as Similarly, the output matrix can be given by Taylor expansion of the nonlinear function atX k is given by whereF k andH k are the Jacobian matrix to linearize the system at the operation point calculated by It is easy to verify that the linearized system is completely observable. The system matrix and output matrix after discretization are given bỹ The validation has been done by the SAG mill at Inner Mongolia Huogeqi copper mine. Discrete nonlinear system    equations are approximated as linear functions by the polynomial expansion method as follows  where u(k) = 0 u r (k) u w (k) 0 0 T and y(k) =

P(k) Q(k) C(k)
T . By applying observability criterion to linearized system (48), we find that the states are observable in the system (48). We then adopt extended Kalman filter to estimate parameters which is a recursive algorithm based on prior estimations and posterior estimations. It is composed of prediction module and error correction module as shown in Fig. 6.

VII. SIMULATION
With the real data collected from the mine, we identify parameters to ensure the model practicability.

A. ORE FEED SIMULATION
Before the parameter identification, some basic data provided by operators are given in Table 1. In addition, the delays in the ore belt system model (1) and (2) are averaged by multiple measurements. The estimation parameter θ minimizes the criterion function (43) with the optimal value given by

B. SAG MILL SIMULATION
Before estimating parameters by extended Kalman filter method, we used the real data to illustrate the correctness of the SAG mill model.
The initialisation of the observer is given by Table 2. By iteration operation, estimated values of parameters in the SAG mill model are given by In the standard dynamic system of the SAG mill, state variables of mill are masses of copper rocks x r , granules x s and water x w which are described in Fig. 9.   Utilizing the estimated parameters, we continue to simulate with another set of data to predict system outputs. The sample interval T is 1s. Based on obtained results shown in Fig. 10, it can be confirmed that the predicted power P and real production data match well. As is seen from the result obtained by parameters estimation, the model can adapt to the influence of the change of control input from the perspective of qualitative and quantitative analysis.

VIII. CONCLUSION
In this paper, we proposed a complete model for the SAG mill systems including ore feeding system and SAG subsystem. The ore feeding quantity is usually controlled by adjusting the frequencies of asynchronous motors in the the heavy-apron feeder system. The feeding model proposed in this paper has a strong migration adaptability and can meet the control requirement for the engineering designer. In the paper, a standard dynamic model is established for a class of SAG mill from the perspective of mechanism analysis. The content of copper rocks, granules and water are taken as state variables, which conforms to the practical application system. For the SAG subsystem, we introduced our pioneer work on establishing the relationship between the hold-ups and power. In addition, internal states of the SAG mill system are described in detail. The average falling height and other important time-varying parameters are derived by means of the rotational speed. In this way, the relations between states and outputs in the SAG subsystem can be described, so as to reduce complexity of the model as well as the number of estimated parameters and to improve migration adaptability of the model.
The ore feed control system and the SAG subsystem constitute a complete mine grinding process model which has instructive significance in the whole ore dressing operation. In the SAG subsystem, according to the ore particle size, the copper ore is divided into two categories. One is the granule which can pass through the discharge grate. Another is the copper rock which cannot pass through and is crushed in the cylindrical shell. The classified method is cursory. However, there is no sensor to measure the distribution of particle size when the copper ore is crushing in the SAG mill. From mining to slurry outflow, the whole process industry control design is inseparable from a complete model. Thus, the following work is to apply the optimization algorithm and design the controller based on this model. In 2010, he joined Shandong University, where he is currently an Associate Professor with the School of Control Science and Engineering. His current research interests include the optimal control of engineering and the optimization of combined cooling, heating, and power systems.
SHUAI LIU received the B.E. and M.E. degrees in control theory and engineering from Shandong University, China, in 2004 and 2007, respectively, and the Ph.D. degree in electrical and electronic engineering from Nanyang Technological University, Singapore, in 2012. He was a Senior Research Fellow with Berkeley Education Alliance for Research in Singapore (BEARS), from 2011 to 2017. Since 2017, he has been with the School of Control Science and Engineering, Shandong University. His research interests include modeling, cooperative control, distributed consensus, optimal estimation and control, time-delay systems, fault detection, and estimation.