Optimizing the control system of cement milling: process modeling and controller tuning based on loop shaping procedures and process simulations

PROCESS SYSTEMS ENGINEERING

Optimizing the control system of cement milling: process modeling and controller tuning based on loop shaping procedures and process simulations

D. C. Tsamatsoulis

Halyps Building Materials S.A., Italcementi Group, Phone: [phone], 17th Klm Nat. Rd. Athens – Korinth, 19300, Aspropyrgos, Greece. E-mail: [e-mail]

ABSTRACT

Based on a dynamical model of the grinding process in closed circuit mills, efficient efforts have been made to optimize PID controllers of cement milling. The process simulation is combined with an autoregressive model of the errors between the actual process values and the computed ones. Long term industrial data have been used to determine the model parameters. The data include grinding of various cement types. The M - Constrained Integral Gain Optimization (MIGO) loop shaping method is utilized to determine PID sets satisfying a certain robustness constraint. The maximum sensitivity is considered as such a criterion. Both dynamical parameters and PID sets constitute the inputs of a detailed simulator which involves all the main process characteristics. The simulation is applied over all the PID sets aiming to find the parameter region that provides the minimum integral of absolute error, which functions as a performance criterion. For each cement type a PID set is selected and put in operation in a closed circuit cement mill. The performance of the regulation is evaluated after a sufficient time period, concluding that the developed design combining criteria of both robustness and performance leads to PID controllers of high efficiency.

Keywords: Dynamics; Cement; Mill; Grinding; Uncertainty; PID; Robustness; Simulation.

INTRODUCTION

In the cement industry, a heavy industry absorbing extremely high energy, the automatic control of the grinding process remains a challenging issue, due to the elevated degree of uncertainties, process non-linearity and frequent change of the set points and the respective model parameters during operation. For productivity and quality reasons, grinding is mostly performed in closed circuits: The ball cement mill (CM) is fed with raw materials. The milled product is fed via a recycle elevator to a dynamic separator. The high fineness stream of the separator constitutes the final circuit product, while the coarse material returns back to the CM to be ground again. A simplified flow sheet, showing the basic components of a closed grinding system, is demonstrated in Figure 1.

In the current automation one of the following is usually selected as a process variable: (1) The power of the recycle elevator. (2) The return flow rate from the separator. (3) An electronic ear in the mill inlet. (4) The mill power or combination of the above. The common field among the majority of designs of controllers is the assumption of a model which describes the process dynamics. Modelling of an industrial milling system is a delicate task due to the multivariable character of the process, the elevated degree of load disturbances, the different cement types ground in the same mill, as well as the incomplete or missing information about some key process characteristics such as clinker hardness, the materials' moisture, grinding media condition, mill holdup, separator efficiency and full particle size of different streams.

Numerous studies describe techniques of mill automation with varying degrees of complexity. Because of the multivariable character and the presence of non – linearities, Model Predictive Control schemes were developed (Efe and Kaynak, 2002, Ramasamy et al., 2005, Prasath et al., 2010). However, it has been mentioned by Astrom and Hagglund (2006) that, in industrial process control, more than 95% of the control loops are of the PID type. As was clearly stated by Astrom (2000), model uncertainty and robustness have been a central theme in the development of the field of automatic control. Therefore, for the model-based control of a CM, not only the values of the dynamic model parameters are necessary, but also their uncertainty.

Austin et al. (1984) developed analytical population models based on the material flows and on the mass balances of the granulometric fractions within the grinding circuit. Tsamatsoulis et al. (1994) applied this model based also on a similar theoretical approach (Tsamatsoulis and Papayanakos, 1994). Van Breusegem et al. (1994) applied a very interesting black box model based on a set of four ordinary differential equations (ODEs') and developed a multivariable linear quadratic controller (MLQC). An analogous approach was followed by Boulvin et al. (1999) and Lapore et al. (2002), leading to the so-called gray-box model. Boulvin et al. (1999) utilized the obtained models to develop MLQC strategies or to implement PI controllers, while Lapore et al. (2002) developed a model predictive control scheme. Tsamatsoulis (2009) developed a black box model between the process and control variables, with some similarities with the one of Van Breusegem et al. (1994). The derived dynamical parameters were used to parameterize PID controllers between the recycle elevator power and the mill feed flow rate (Tsamatsoulis, 2009; Tsamatsoulis and Lungoci, 2010).

The objective of the current study is twofold: (i) to extend the model developed by Tsamatsoulis (2009) by investigating the dynamical parameters not only per CM but per cement type (CEM) as well. Long term data of an industrial mill were used for this purpose. A detailed approach of analyzing the parameter uncertainty is followed and the results are discussed. An attempt is also made to model the errors between the actual process values and the computed ones by the dynamical model. (ii) The results of the first step are the inputs of advanced techniques to parameterize a PID controller in a robust way. The obtained sets of PID parameters satisfying a robustness constraint along with the dynamical results and their uncertainty, constitute the entries of extensive simulators to define the optimum area of PID coefficients as concerns performance. A PID set belonging to the optimum area is selected and placed in operation in the mill investigated. The long term results of the actual performance are evaluated. A similar approach has been presented by Tsamatsoulis (2011) in two consecutive papers analysing a different process: the optimization of PID controllers regulating the quality of cement raw meal. The differences between this previous study and the present one are the following: (i) the method was applied to regulate the chemical modules of the raw meal, while in the current work an attempt is made to regulate the cement mill operation in an optimal way; (ii) the simulation of the raw meal mixing process included the uncertainties of the raw materials, while in this study the uncertainty of the dynamical parameters is directly incorporated in the simulation. (iii) The present investigation is extended to a series of cement products, while the mentioned papers analyze one semifinal product.

PROCESS MODEL

In all cases the dynamics between the process variable and the mill feed flow rate is modelled. The model considers as process variables the power of the recycle elevator or the flow rate of the separator return. The total feedback control loop is demonstrated in Figure 2. The symbols G_c, G_w, G_p and G_f denote the transfer functions of the controller, weight feeders, process and filter, respectively. The feeder output, u₁, constitutes the input of the system, while the filter output, y, is the process variable.

The transfer function between the process output and input contains an integral term as was proven by Tsamatsoulis (2009) after application of a double pulse in the mill feed. Therefore, G_p is given by Equations (1) to (3) in the Laplace domain:

where Q is the mass flow rate, PV represents the power of the recycle elevator or return flow rate, Q₀ is the flow rate corresponding to PV= PV₀ of the steady state, k_v is the gain and T_d is the delay time. The meaning of the coefficient α_x is the following: there are cases where the filter input is not the process value but the ratio between the PV and the maximum one, PV_Max, expressed as a percentage. In these cases α_x=100/ PV_Max. Otherwise, if the controller input is equal to PV, then α_x=1.

To avoid undesirable measurement noise, the signal x is filtered: the transfer function of the filter can be described with one of the two following equations:

First order filter:

Exponentially weighted moving average filter (EWMA):

where T_f symbolizes the first order time constant, T_s is the sampling period of the unfiltered signal x and k is a constant belonging to the interval [0, 1]. The set of the model parameters consists of the gain k_v, delay time T_d, flow rate Q₀ and variable PV₀, corresponding to the system steady state, under the specific operating conditions, such as: (a) cement type, (b) raw materials moisture and grindability, (c) gas flows, (d) pressures, (e) temperatures, (f) condition of the grinding media, (g) separator efficiency etc. The short or long term variance of these conditions generates the model uncertainty as to the parameter values.

The model parameters are estimated using the convolution theorem between the input signal u₁ and the process variable y, expressed by the Equation (6).

where g(t) is the pulse system response. Exclusively routine operational data of the CM are utilized. These data are sampled by applying convenient software. Using the Newton-Raphson non linear technique, the optimum dynamical parameters are computed by minimizing the residual error provided by Equation (7):

where s_res represents the residual error, y_calc is calculated from the model process value, y_exp is the actual one, N is the number of experimental points and k₀ is the number of the independent model parameters. At time I the error between y_calc and y_exp is given by Equation (8). This error is modeled with the autoregressive Equation (9).

To investigate whether this model error is adequate, its regression coefficient was checked and its standard error compared with the residual error of the dynamical model.

IDENTIFICATION OF THE MODEL PARAMETERS

Code in Visual Basic for Application was developed to load and to process industrial data for routine operation of a cement mill, directly extracted from the plant database. In each extraction two days worth of data are loaded, with a sampling period of one minute. Then the software checks for feeder stoppages and finds continuous operational data sets of 250 minutes duration. In parallel it detects the cement type milled during this time interval. Data for CM5 of the Halyps cement plant were utilized, belonging to a long time period from January, 2010, to January, 2011. In this mill three cement types are milled conforming to the European Norm EN 197-1. The following codes are used for the above: Code 1 = CEM A-L 42.5 N, Code 2 = CEM IV B (P-W) 32.5 N, Code 3 = CEM II B-M (P-L) 32.5 N. Some details of the CM circuit are provided in Table 1.

For the selected CM, the power of the recycle elevator constitutes the process variable. Because the process value is a percentage of the maximum power of the elevator, the formula α_x =100/ PV_Max is used. For the given application PV_Max =50 KW and, consequently, α_x = 2. A first order filter is applied with T_f = 180 s. By applying a step increase of the mill feed and sampling the flow rate every second, it is experimentally found that the weight scales follow first order kinetics with a time constant of T_w = 11 s. Therefore, the transfer function G_w is given by Equation (10):

Afterwards the software determines the optimum dynamical parameters for each data set and the corresponding regression coefficient, R. A minimum coefficient, R_Min, is decided and data sets presenting R<R_Min are neglected in subsequent processing. In the current application, an R_Min=0.7 is chosen. In this way a total population of results is obtained, including all the three cement types. According to the cement type in each data set, the software creates three subpopulations. If there are two types for a data set, then this is not attributed to a subpopulation. For the full population and the three subpopulations, the mean value and standard deviation of each dynamical parameter are computed. A minimum number of deviations, N_s, was also employed as concerns k_v and T_d populations. Sets not satisfying the constraints │k_v - Aver_kv│ ≤ N_s∙s_kv or │T_d - Aver_Td│ ≤ N_s∙s_Td are considered as outlying and not taken into account in further processing. A value of N_s equal to 2 is taken in this investigation. The described double screening of the dynamical parameters can be thought of as a procedure to enhance the validity of the results and to reduce the impact of the load disturbances. The processing continues with the calculation of the average and median values, as well as of the standard deviation, which can be considered to be a measure of the parameter uncertainty. The results for the gain and delay time are demonstrated in Table 2. For each dynamical parameter, the coefficient of variation %CV=100×Std. Dev. / Average is also computed. The parameters of the error model are presented in the same table as well. Only the values of A₁, A₂ are shown because A₀ is around zero in all cases.

In Table 2 it can be observed that the average gain per CEM type is different from the one computed for all the sets. The %CV of gain for each CEM type is smaller than that for all types. As a consequence, the calculation of the parameters per cement type contributes to the decrease of the uncertainty. To investigate the adequacy of the autoregressive error model, the average values of s_Err and s_Res are compared and shown in Table 2. Since the ratio of the two errors remains always less than 0.4, it can be concluded that the superposition of the error model to the dynamical one contributes to a noticeable decrease of the estimation error.

Therefore, the incorporation of an autoregressive model for the error enhances the reliability of the dynamical model in case it is used for simulation. As a result, Err(t) shown in Equation (9) is a coloured noise while Std_ModErr is a white noise, with standard deviation s_Err. The error model possesses a remarkable property. The sum of the parameters A₁, A₂ has a significantly lower standard error than each parameter separately; therefore, a correlation exists. The correlation coefficient between A₁ and A₂ is found in all the four cases to be between 0.87 and 0.89. Thus, for simulation purposes the following group of equations could be used in the time domain: the dynamical set (1) to (3) and (6), the filter Equations (4) or (5) and the autoregressive formula (9) where A₁ is a normally distributed variable with average value and standard deviation given in Table 2. Then, the sum A₁+A₂ can be considered as a constant S_A. Consequently, A₂ is a variable that follows the normal distribution and amounts to S_A-A₁.

Processing of the Dynamical Results

The ability of the model to describe the dynamics between the two selected variables was initially examined by constructing the distribution of the regression coefficients of the entire data set population, depicted in Figure 3.

As can be observed, the sets with R ≥ 0.7 constitute a percentage higher than 75% of the total population. Because the regression coefficient is related to the data variance, two functions are generated: (a) the one between R and %CV of the process variable and (b) the one between the percentage of sets with R ≥ 0.7 and the same %CV. The coefficients of variation are grouped in sets. For each set the average R and the percentage of sets with R ≥ 0.7 are computed. The results are shown in Figure 4, where the distribution of %CV of the data set population is also demonstrated. More than 95% of the data sets have %CV ≥ 1.4. After this threshold, the average R of each %CV group always remains higher than 0.75 (Figure 4a). For %CV ≥ 1.4 the percentage of data sets with R ≥ 0.7 is higher than 70% (Figure 4b). It must be elucidated that sets with low %CV present very poor dynamics.

Distributions of Gain and Delay Time

The distributions of the gain and delay time for each cement type are plotted in Figure 5. The gain distributions approach the Gaussian distribution, but the variance is high. The time delay differential distributions deviate noticeably from the normal one. The above is a strong verification of the short and long term variation of the operating conditions and materials properties.

PID CONTROLLER DESIGN

It is crucial in control system design and tuning to assure the stability and performance of the closed loop in the case that strong load disturbances exist or a mismatch between accepted and actual model occurs, i.e. to guarantee robustness. In spite of the huge amount of work presented on tuning of PID controllers, it remains a challenge to parameterize effectively and to operate in the long term existing PID loops installed in grinding systems. A widely applied technique to adjust PID controllers is robust loop shaping (McFarlane and Glover, 1992; Zolotas and Halikias, 1999; Yaniv and Nagurka, 2005). A remarkably effective loop shaping technique was proven to be the called MIGO (M-constrained integral gain optimization) developed by Astrom, Panagopoulos and Hagglund (1998, 2002, 2004). The design approach is to maximize integral gain subject to a constraint on the maximum sensitivity defined by Equation (11). For the specific feedback control loop between recycle elevator power and mill feed designed in Figure 2, the sensitivity is given by Equation (12):

The controller transfer function G_c is provided by Equation (13). The variables k_p, k_i, k_d constitute the proportional, integral and differential gains of the controller correspondingly.

The error, e, is given by the difference y_sp-y, y_sp representing the set point of the process variable. For M_s belonging to the interval [1.3, 2.5] and with a stepwise increase of 0.1, the (k_p, k_i, k_d) sets are computed for each cement type using the average or median dynamical parameters shown in Table 2. The k_p and k_i values as a function of k_d and M_s are shown in Figure 6 for the CEM B-M (P-L) 32.5, using the average dynamical parameters. In this way, for each (M_s, k_d) pair a PID set of parameters is determined.

Impact of the Dynamical Parameter Uncertainty on the Open Loop Properties

The model uncertainty expressed by the respective variance of each parameter has a strong effect on the open loop features. To verify the above, controllers corresponding to [M_s, k_d] = [1.5, 2] are selected for each CEM type and for all the data as well. For these controllers, the M_s are computed for each data set and the cumulative distributions are plotted. The results are depicted in Figure 7.

All the distributions present a noticeable dispersion around the median value. Therefore, the impact of the parameter uncertainty on the maximum sensitivity is very strong, rendering it more necessary for the operation of a robust controller. To identify in more detail this impact in connection with the selected PID parameters on the open loop properties variance, the following procedure was employed:

(i) Design parameters M_s,Des and k_d are selected belonging to the ranges: 1.3 ≤ M_s,Des ≤ 2.5 with a step 0.1 and 0 ≤ k_d ≤ 5.6 with a step 0.4. For each [M_{s, Des}, k_d] pair, the (k_p, k_i, k_d)^T vector is computed and applied over all the dynamical sets.

(ii) For each one of these sets, the maximum sensitivity, gain margin, g_m, and phase margin, φ_m, are calculated. Then for the total population of the three open loop characteristics, the median value X₅₀, and the percentiles X₂₅, X₇₅ are determined, where X = M_s, g_m, φ_m.

(iii) The sharpness index SI(X) of each distribution is defined by Equation (14). SI(X) can be considered to be a measure of the coefficient of variation of each distribution.

(iv) For the main cement types ground in CM5 – CEM codes 2 and 3 – the sharpness indices are calculated for the full range of M_{s, Des} and k_d. Then for each open loop characteristic X and CEM code Y, the ratio RSI(X, Y) is computed according to Equation (15).

The RSI distributions for the three open loop properties are plotted in Figure 8 as two dimensional graphs, in which a significant reduction of the variance of the open loop characteristics is observed. Consequently the parameterization of the controller per cement type results in a significant improvement of the parameter uncertainty.

PID OPTIMIZATION BY PROCESS SIMULATION

The feedback control loop shown in Figure 2 is simulated in order to find the optimum area of the PID sets, determined in the previous section. The Integral of Absolute Error between y and y_sp, as a percentage of the y_sp, is defined as the optimization criterion, given by Equation (16). The sets minimizing %IAE compose the region of the optimum parameters.

The integration is performed in a time interval t₀. The simulator has the ability to select a family of PID vectors among the full sets found by the MIGO technique and simulates the mill operation for each (k_p, k_i, k_d)^T. As concerns the PID implementation, a digital form is introduced, described by Tsamatsoulis (2009). For 50 minutes the mill operates in manual mode with constant feed. Afterwards the mill is transferred to automatic control by applying the selected (k_p, k_i, k_d)^T. The settings are the following: (1) cement type; (2) sampling period, T_s; (3) actuator period, T_a; (4) minimum control variable, Q_Min; (5) maximum control variable, Q_Max; (6) set point of process value, y_sp; (7) mill start up feed, Q_in; (8) time interval of mill operation in automatic mode, T_oper; (9) average or median values and standard deviations of the dynamical parameters; (10) average A₁ coefficient of the model of errors and its standard deviation; (11) number of standard deviations of the parameters (9) and (10); (12) average A₁+A₂ coefficient; (13) average error and standard deviation of the error of the model of errors, Av.s_Err and s s_Err, respectively; (14) process value noise before filtering, err_Noise; (15) time constant of the first order filter, T_f, or constant k of the moving average filter; (16) number of iterations, N_Iter. The simulation can run in two versions: (a) a simplified and approximate one (b) a full one where as much as possible of the uncertainty is taken into account. In the first version the parameters (7), (9), (10) are considered to be constant. Dynamical coefficients are equal to the average or median values. Their standard deviation is neglected. The selection between average or median values depends on the way the PID is parameterized. In this version the noise of the process variable is considered to be negligible and the number of iterations is only one. The full implementation of the simulator tries to express numerically the large range of uncertainties and disturbances encountered in the actual CM operation:

(i) The startup feed Q_in, can take values between two limits [Min, Max]. Then a generator of random numbers selects a value within this interval.

(ii) For each dynamical or error model parameter D {k_v, T_d, Q₀, PV₀, A₁} with average or median value D_M and standard deviation s_D a number of standard deviations N_s is assumed. The software generates a random number, p, in the interval (0, 1). By applying the normal distribution a parameter D is found, satisfying in parallel the condition (17).

Thus, for each iteration, representing a simulated mill automatic operation with duration equal to T_oper, a new set of (k_v, T_d, Q₀, PV₀, A₁) is generated, satisfying the constraints described.

(iii) The average error and the respective standard deviation of the error model are activated. By repeating the technique analyzed in the previous step, a value of s_Err is computed for each simulation.

(iv) The noise of the process output before filtering is introduced. Its value is found experiment-tally from the actual operation, by subtracting the

filter for a certain time interval and by sampling the process variable with high frequency.

(v) For each iteration I the corresponding %IAE(I) is computed. To investigate whether the number of iterations, N_Iter, is adequate to converge IAE, for each I, %IAE(I)_Aver is computed from Equation (18):

(vi) The number of iterations is found using the following convergence criterion: if there is a sufficiently large number of iterations I₀, so that:

Then IAE values converge to the actual average. Otherwise the iterations must be increased and convergence rechecked. A value of α_Conv equal to 0.02 was used.

The settings utilized for the simplified and full implementation of the simulator are shown in Table 3. The simulator is applied using these settings and the complete sets of PID gains.

Initial Simulations

Before the application of the simulator by using the full range of PID sets and parameter uncertainty, some preliminary simulations were performed to investigate the impact of some main factors on the performance and robustness of the control. Initially the simplified simulator was used: The mill starts to operate with Q_in=26 t/h grinding cement of code equal to 2. All the other parameters are those of Table 3. The impact of M_s on the response of the process variable is investigated as follows: For a constant k_d=4, four PID sets were selected so that they correspond to M_s=1. 3, 1.7, 2.1, 2.5.

The evolution of the control and process variables as a function of time for 10 hours of operation in automatic mode is shown in Figure 9, from which the following conclusions can be derived: M_s robustness constraint affects noticeably the closed loop properties as concerns the settling time and overshoot. For low M_s values, the process variable does not settle, even after 10 hours of operation. For the selected value of k_d, as M_s moves from 1.3 to 2.5, the settling time passes through a minimum. Overshoot appears for a certain value of maximum sensitivity, increasing as M_s augments. Thus, for the selected k_d, an M_s=1.7 provides the minimum settling time and it can be considered to be the optimal value. The response of the process variable has also been plotted for another set (M_s, k_d) = (1.9, 1.2). This set generates a response without overshoot and with settling time lower than the set (M_s, k_d) = (1.7, 4), leading to the conclusion that the optimization of the parameters is a two-dimensional problem.

To investigate the shape of the process variable response, when the parameter uncertainty, the model of errors and the noise are incorporated into the simulator, the following simulation was performed: For an initial feed of 26 t/h, the full version was implemented three times with N_iter=1. Grinding of cements with code 2 was considered. All the other parameters were those of Table 3. The controller operates with (k_p, k_i, k_d)^T corresponding to (M_s, k_d) = (1.9, 1.2). The simplified version of the simulator was also applied in order to compare the responses. The results are presented in Figure 10. Depending on the values of the set (k_v, T_d, Q₀, PV₀, A₁), created with a random generator but belonging to the actual population of the dynamical parameters, the %IAE differs from the one calculated using the simplified version. Consequently, the full version of the simulation provides a more realistic picture of the actual variance of the process and of the introduced uncertainties as well.

Optimization of the PID Parameters

For each cement type the simplified and full versions are applied. For every PID set corresponding to a (M_s, k_d) pair the %IAE(N_Iter)_Aver was computed. In the full version, an adequate number of iterations is necessary so that the condition (19) is valid for the last I₀ iterations. This number depends on the N_s level. For N_s=1, a N_Iter ≥ 72 suffices for the convergence of %IAE(N_Iter)_Aver. An example of convergence for CEM code equal to 3 is demonstrated in Figure 11: during the last ~20 iterations, the average error remains in the ±2% zone of the final %IAE(N_Iter)_Aver.

The region of the PID sets giving the optimum %IAE_Aver was computed as follows: (a) the controller providing the minimum %IAE_{Aver, Min} was found; (b) the sets with %IAE_Aver ≤ 1.2 %IAE_{Aver, Min} were determined; (c) this area constitutes the optimum region of the PID coefficients. This procedure was performed for both the simplified and full versions of the simulator. The results for the three CEM codes are shown in Figure 12 and the following conclusions can be extracted.

The optimum area after the application of the simplified version always constitutes a subset of the corresponding area obtained using the full version. The additional area obtained via the full version always contains lower values of M_s in comparison with the common region; thus, it is more robust. This enhancement of the robustness can be explained by the fact that all the main uncertainties have been incorporated in the full implementation. Consequently, if one utilizes

the simplified simulator to speed up the optimization process, one must always select a (k_p, k_i, k_d)^T located in the left border of the optimum region. The combination of a robustness constraint with a performance criterion, maximum sensitivity and integral of absolute error, respectively, leads to the determination of an optimum region of PID coefficients. Sets of this optimum area can be chosen and implemented in the current study of the operation of the grinding system.

APPLICATION OF THE RESULTS TO AN ACTUAL MILLING SYSTEM

From the results demonstrated in Figure 12, three controllers are chosen, one per CEM type. These sets, presented in Table 4, were applied in the CM studied. The gain and phase margins are also depicted in the same table.

These PID parameters per CEM type were implemented continuously for 290 days of CM operation. The performance of the parameterization was evaluated by calculating the daily %IAE_Aver. The data were sampled from the database of the plant with a period of 1 minute. The long-term performance of the techniques applied to tune the controller parameters of CM5 is shown in Figure 13, where three distinct periods are demonstrated. In the first one the PID was tuned by trial and error. In the second one the parameters were adjusted with loop shaping techniques and mainly via the MIGO method, as analyzed by Tsamatsoulis (2009). One PID set was used for all the cement types during the first and second period. In this previous level of development of the implementation of the MIGO technique in the automatic control of mill operation, only the robustness constraint was considered during the design. In the last period the parameters shown in Table 4 were put in operation.

For each period the average %IAE is also plotted, as well as the low and high control limits, LCL and HCL, respectively, calculated using the statistical norm ISO 8258:1991 dealing with Shewhart control charts. The values of the above parameters are shown in Table 5. With the implementation of the loop shaping techniques and mainly of the MIGO method in the second period, a significant decrease of the %IAE was achieved in comparison with the period in which the PID was parameterized with trial and error. In the third period where the controller was tuned per cement type, by utilizing a combination of robustness and performance criteria during the design, a further amelioration of %IAE occurs: Its value reached a level of 1.6%, resulting in a 30% decrease compared with that of the previous period. The high control limit is about equal to the average of the previous period. The current set point of the recycle elevator power is 13.3 kW. The discrimination of the measurement is 0.1 kW, corresponding to 0.75% of the set point. As can be seen, the LCL is very near to this point. Based on this assessment the PID parameters chosen belong to the optimum region.

CONCLUSIONS

Based on a dynamical model of the grinding process in closed circuit mills, efficient efforts have been made to optimize the PID controller regulating the operation. The model involved not only the grinding process, but also the weight feeders and filter dynamics. Simultaneously an autoregressive model of the errors between the actual process values and the computed ones was built. The combination of these two models to simulate the process results in a significant decrease of the standard error. The model parameters and their uncertainty, computed per cement type, feed the M-constrained Integral Gain Optimization technique, to determine PID sets satisfying a certain robustness constraint. The maximum sensitivity was utilized as such criterion.

Both dynamical and PID parameters constitute the inputs of a detailed simulator involving all the main process characteristics. The simulation was applied over all the cement types and (k_p, k_i, k_d)^T parameters aiming to find the PID sets providing the minimum integral of absolute error, which is used as a performance criterion. Two versions of the simulator are implemented: (a) the first one where the parameter uncertainties are neglected and only the average or median values are considered; (b) a second one where all the uncertainties are included and model parameters are chosen using a random number generator. The optimum area of the PID parameters, found with the full version, is wider than the one determined by the simplified version, extended into a more robust region, due to the incorporation of the uncertainties.

For each cement type a PID set is selected and placed in operation in an actual closed circuit cement mill. After a period of 290 operating days the results were assessed using the daily average %IAE and compared with the results of the previous long-term periods: The first one where the PID was tuned by trial and error and the second one where loop shaping techniques were utilized, using only the robustness constraint during the design. For both periods a PID set operated independently of the milled cement type. The developed design combining criteria of both robustness and performance led to PID controllers of high efficiency, leading to an average %IAE equal to 1.6%, while around 98% of the %IAE population was less than 3%.

The same methodology has been successfully applied to other closed circuit cement mills with productivities ranging from 30 t/h to 90 t/h depending on the size of the circuit and the CEM type produced. In cases where the flow rate of the material returned from the separator is measured, then this variable can also be used as a process variable. The corresponding dynamics between feed and return is determined using the same model. The selection of the best process variable depends on the value of the dynamical model regression coefficient. Due to the good efficiency achieved in such applications, the dynamics are expected to describe high capacity circuits processing as much as 150 t/h of fresh feed.

NOMENCLATURE

A₀, A₁, A₂

coefficients of error model

Aver_kv

average of gain population

%∙(10³∙kg)^-1

Aver_Td

average of gain population

60∙s

Err(I)

error between calculated

and actual process values at time I

%

err_Noise

standard error of the noise

%

g(t)

pulse system response

g_m

gain margin

G_c

controller transfer function

G_f

filter transfer function

G_p

process transfer function

G_w

feeders transfer function

%IAE

integral of absolute error

%

k

constant of EWMA filter

k_v

gain

%∙(10³∙kg)^-1

k₀

number of model parameters

k_p

proportional gain of PID controller

k_i

integral gain of PID controller

k_d

derivative gain of PID controller

M_s

maximum sensitivity

N

number of experimental points

N_Iter

number of iterations

N_s

number of standard deviations

PV

elevator power

10³∙W

PV₀

elevator power in steady state

10³∙W

Q

mill inlet

(3.6)^-1∙kg∙s^-1

Q₀

mill inlet in steady state

(3.6)^-1∙kg∙s^-1

R

regression coefficient

R_Min

minimum regression coefficient

RSI

ratio of sharpness indices

S

sensitivity function

SI

sharpness index

s_res

dynamical model residual error

%

s_Err

standard error of the error model

%

s_kv

standard deviation of gain population

%∙(10³∙kg)^-1

s_Td

standard deviation of time delay population

%∙(10³∙kg)^-1

T_d

delay time

60∙s

T_a

actuator period

s

T_f

time constant of the first order filter

s

T_oper

period of simulator operation

3600∙s

T_s

sampling period

s

T_w

feeders time constant

s

x

process output before filtering

%

X₅₀

median value

X₂₅, X₇₅

percentiles 25% and 75%

y

process variable

%

y_sp

set point of the process variable

%

u

controller output

(3.6)^-1∙kg∙s^-1

u₁

process input

(3.6)^-1∙kg∙s^-1

Greek Symbols

α_x

coefficient of transformation

%∙(10³∙W)^-1

α_Conv

convergence error

φ_m

phase margin

Subscripts

calc

calculated value

exp

experimental value

In

initial value

Min

minimum value

Max

maximum value

Abbreviations

CEM

cement type

CM

cement mill

EWMA

exponentially weighted moving average

HCL

high control limit

LCL

low control limit

PID

proportional integral derivative controller

Submitted: December 9, 2012

Revised: May 6, 2013

Accepted: June 2, 2013

• Astrom, K. J., Panagopoulos, H. and Hagglund, T., Design of PI controllers based on non-convex optimization. Automatica, 34, 585-601 (1998).
• Astrom, K. J., Model uncertainty and robust control. Lecture Notes on Iterative Identification and Control Design, Lund University (2000).
• Astrom, K. J. and Hagglund, T., Revisiting the Ziegler-Nichols step response method for PID control. Journal of Process Control, 14, 635-650 (2004).
• Astrom, K. and Hagglund, T., Advanced PID control. Research Triangle Park: Instrumentation, Systems and Automatic Society (2006).
• Austin, L. G., Klimpel, R. R. and Luckie, P. T., Process engineering of size reduction: Ball milling. Society of Mining Engineers, New York (1984).
• Boulvin, M., Renotte, C., Vande Wouver, A., Remy, M., Tarasiewitz, S. and Cesar, P., Modeling, simulation and evaluation of control loops for a cement grinding process. European Journal of Control, 5, 10-18 (1999).
• Efe, M. O. and Kaynak, O., Multivariable nonlinear model reference control of cement mills. 15th Triennial World Congress, Barcelona (2002).
• Hagglund, T. and Astrom, K. J., Revisiting the Ziegler-Nichols tuning rules for PI control. Asian Journal of Control, 4, 364-380 (2002).
• Lepore, R., Vande Wouwer, A. and Remy, M., Modeling and predictive control of cement grinding circuits. 15th Triennial World Congress, Barcelona (2002).
• McFarlane, D. and Glover, K., A loop shaping design procedure using H_∞ synthesis. IEEE Transactions on Automatic Control, 37, 759-769 (1992).
• Panagopoulos, H. Astrom, K. J. and Hagglund, T., Design of PID controllers based on constrained optimization. IEEE Proceedings-Control Theory and Applications, 149, 32-40 (2002).
• Prasath, P., Recke, B., Chidambaram, M. and Jorgensen, J. B., Application of soft constrained MPC to a cement mill circuit. 9th International Symposium on Dynamics and Control of Process Systems, Leuven, Belgium (2010).
• Ramasamy, M., Narayanan, S. S. and Rao, Ch. D. P., Control of ball mill grinding circuit using model predictive control scheme. Journal of Process Control, 15(3), 273-283 (2005).
• Tsamatsoulis, D., Fotopoulos, A., Marinos, I. and Papayannakos, N., Optimization of the cement production process using computer simulations. Applications of Automatic Systems in the Process Industry, Athens (1994).
• Tsamatsoulis, D. and Papayannakos, N., Axial dispersion and holdup in bench scale trickle bed reactor at operating conditions. Chemical Engineering Science, 49, 523-529 (1994).
• Tsamatsoulis, D., Dynamic behavior of closed grinding systems and effective PID parameterization. WSEAS Transactions on Systems and Control, 4, 581-602 (2009).
• Tsamatsoulis, D. and Lungoci, C., Effective optimization of controllers stabilizing closed circuit grinding systems of cement. 12th International Conference on Optimization of Electrical and Electronic Equipment, Brasov, Romania, (2010).
• Tsamatsoulis, D., Effective optimization of the control system for the cement raw meal mixing process: I. PID tuning based on loop shaping. WSEAS Transactions on Systems and Control, 6, 239-253 (2011).
• Tsamatsoulis, D., Effective optimization of the control system for the cement raw meal mixing process: II. Optimizing robust PID controllers using real process simulators. WSEAS Transactions on Systems and Control, 6, 276-288 (2011).
• Van Breusegem, V., Chen, L., Werbrouck, V., Bastin, G. and Wertz, V., Multivariable linear quadratic control of a cement mill: An industrial application. Control Engineering Practice, 2, 605-611 (1994).
• Yaniv, O. and Nagurka, M., Automatic loop shaping of structured controllers satisfying QFT performance. Transactions of the ASME, 125, 472-477 (2005).
• Zolotas, A. C. and Halikias, G. D., Optimal design of PID controllers using the QFT method. IEE Proc. Control Theory Appl., 146, 585-590 (1999).

^*

To whom correspondence should be addressed

Publication Dates

• Publication in this collection
  20 Mar 2014
• Date of issue
  Mar 2014

History

• Received
  09 Dec 2012
• Accepted
  02 June 2013
• Reviewed
  06 May 2013