Optimizing the Control System of Clinker Cooling: Process Modeling and Controller Tuning

: This paper aims to present efﬁcient efforts to optimize the proportional-integral-differential (PID) controller of clinker cooling in grate coolers, which have a ﬁxed grate and at least two moving ones. The process model contains three transfer functions between the speed of the moving grate and the pressures of the static and moving grates. The developed software achieves the identiﬁcation of the model parameters using industrial data and by implementing non-linear regression methods. The design of the PID controller follows a loop-shaping technique, imposing as a constraint the maximum sensitivity, M s , of the open-loop transfer function and providing a set of PIDs that satisfy a range of M s . A simulator determines the optimal PID sets among those calculated at the design step using the integral of absolute error ( IAE ) as a performance criterion. The combination of a robustness constraint with a performance criterion, M s and IAE respectively, leads to an area of controllers with M s belonging to the range of 1.2 to 1.35. The IAE is between 4.2% and 4.8%, depending on the set-point value. PID sets located near the middle of this area can be chosen and implemented in the cooler’s routine operation.


Introduction
Grate coolers constitute a key component of the clinker's dry production lines, having the rapid cooling of clinker and heat recovery as their main function. The stability degree of the cooler's operation affects both the performance of the entire burning line and the product's quality. A simplified representation of the flows and the basic components of a rotary kiln (RK) is demonstrated in Figure 1a. The raw meal, consisting mainly of ground limestone and clay, is fed to the preheater, where it is heated and partially calcined by the hot gases from the kiln and the pre-calciner in an anti-current flow. The decomposed material enters the kiln, where after the calcination is completed, it is transformed into clinker. The air volumes needed for the calcination and combustion come from the primary, secondary, and tertiary air. The very hot clinker then falls onto the grates of the cooler through the hood of the kiln and is cooled by air streams passing vertically from the grates. The cooler includes a fixed grate and at least two moving ones. Figure 1b shows a simplified graphic of a grate cooler. Due to their oscillating motion, the grates cause the solid material to flow forward and spread. After the final moving grate, clinker is first passed through a crusher and then transferred to the corresponding silo.
During clinker cooling, a heat exchange occurs between the solid material and vertically incoming fresh air. A part of the hot air volume enters the kiln as secondary air, another part is led to the pre-calciner, and the remaining goes towards the de-dusting filter as exhaust air. The effectiveness of clinker cooling is linked to two factors: (a) the recovery of heat from the material entering the cooler and (b) its cooling rate, because rapid cooling helps to increase the quality of the product [1]. Another particularly important goal is to keep the temperatures and flow rates of the secondary and tertiary air as stable as possible because any significant disturbance to them affects kiln operation. The above is possible by controlling and stabilizing the thickness of the clinker bed on the first moving grate, taking into account that this thickness depends on the speed of the grate and the mass flow rate from the kiln. Because the latter cannot be controlled by the cooler, as any change is a disturbance, the only control variable is the speed of the first moving grate. The primary air is also not controlled by the cooler. It is fed to the main burner with the fuel and determines the shape of the flame. Subsequent grates usually move 10-20% faster than the first one. For a very thick layer, cooling air cannot blow through the clinker, the cooling is uneven, and the clinker remains too hot at the outlet of the cooler. In the case of a thin layer, the air passes rapidly through it and the heat recovery is poor, resulting in reduced temperatures of secondary and tertiary air. In addition, there is a risk of damage to the grate's plates if a very hot wave of clinker comes from the kiln, so a sufficient material thickness protects the cooler. In practice, reliable continuous measurement of bed thickness is difficult. For this reason, the pressure under the first moving grate is used, considered as proportional to the bed's layer on the grate, as the process variable. The latter also depends on the flow rate of air entering the grate. In most cases, in industrial practice, the speed of cooling fans remains relatively constant and is regulated only under specific conditions. Thus, the main loop of the cooler's automatic control, in an attempt to attenuate disturbances, is between the speed and the pressure of the first moving grate. In most automation systems in modern cement plants, conventional proportional-integraldifferential (PID) controllers are installed in the subsystem of the clinker cooling, using the previous loop. The difficulty arises in finding optimal or at least satisfactory values of the PID parameters, which adequately regulate the grate pressure. In the process control, more than 95% of the control loops are of PID type [2] (p. 1). Astrom [3] clearly stated that a central theme in automatic control is model uncertainty and robustness. Therefore, the automatic control of the cooler requires a reliable model, which includes control and process variables, as well as variables that describe disturbances and noises entering the system. In terms of dynamic parameters, both their mean values and uncertainty are equally necessary.
Tsamatsoulis [4] developed a simplified dynamic model of clinker cooling. He modeled the response of three different pressures to changes in the speed of the moving grate using time-delayed integral transfer functions and long-term industrial data. The simplified transfer functions indicate that there is some effect of the speed of the moving grate on the pressure of the fixed one. However, the modeling did not take into account that any change in the pressure of the fixed grate causes a disturbance in the moving one. Internal model control (IMC) was applied to design a PID in [5] for stabilizing the grate pressure. The authors state that the use of IMC overcomes several drawbacks of the PID controller. The model used two transfer functions between (1) grate speed and pressure as well as (2) kiln feed flow rate and grate pressure. The output of the second function causes a disturbance in the pressure. The first-order with time delay (FOTD) model was applied to both functions. Yu et al. [6] carried out research on the modeling of grate coolers based on industrial data. Among several transfer functions, they extracted a third-order one, between speed and pressure of the grate. In parallel, they analyzed the dominant factors affecting this pressure. In [4][5][6], the description of the models is provided in the Laplace domain.
Xiaohong et al. [7] developed a predictive control scheme for grate coolers based on a grey model, attempting to build a controller more effective than PID when the fluctuation of clinker from the kiln head is high. Jun et al. [8] implemented a dynamic matrix control (DMC) algorithm to control a grate cooler by creating a prediction model, designing the controller based on this model, and finally correcting the resulting error. According to the authors, the results show that DMC can forecast the changes in grate pressure and allow the operator to adjust the grate speed. Fuzzy logic was also utilized for control purposes in some cases. Wang and Yuan [9] developed a self-tuning fuzzy PID control method based on a Kalman filter to overcome the frequently varying operating conditions and high noise-to-signal ratio of pressure signals. Zhang and Lu [10] implemented two fuzzy controllers simultaneously in two different points of the grate cooler. They implemented the controllers in an industrial cooler and compared the results with those of the previously applied controller, concluding that significant improvement had been achieved. In recent years, several researchers have implemented machine learning techniques to model the cooler's operation. Karikov et al. [11] developed a neural network model as a stage of full automation. An Elman network [12]-a particular case of a multi-layer recurrent network-was trained and tested using industrial data. The authors conclude that the obtained models were suitable for solving the synthesis of control methods. Seraj and Shooredeli [13] deployed soft sensor models based on neural networks, concluding that these models help to study the behavior of a cooler in various conditions and design controllers. Geng et al. [14] implemented a multi-model fusion neural network to predict the cooler's pressure. The authors attempt to improve the accuracy by balancing the prediction advantages of single algorithms through design. Wang et al. [15] applied support vector machine techniques to optimize the grate pressure setting, aiming to ensure the cooling and quality of clinker.
The objectives of the present study are as follows: (i) to extend the simplified model developed by Tsamatsoulis [4] by finding a more precise model and by connecting the dynamics of the moving grate to that of the fixed one. This study attempts modeling of the errors between the actual process values and the calculated ones by the dynamic model. The continuous operation results have been retrieved from the kiln cooler of the Halyps cement plant. (ii) To determine sets of PID gains by applying a loop-shaping technique that takes into account the dynamic parameters and a strong robustness condition, first deployed by Astrom et al. [16] and further developed in [2] (pp. 206-221). (iii) To build a process simulator, which takes as inputs the dynamic parameters and their uncertainty, modeling errors and noise, load disturbances, and PID gains, and determines the range of PID values to achieve optimal performance of control. The author applied a similar approach to controlling various processes in the cement industry: (1) raw meal quality control [17][18][19]; (2) cement mill operation [20,21]; and (3) kiln pre-calciner operation [22]. The author developed all software of the controller optimization in C#.
The structure of the paper is as follows. Section 2 provides the hardware setup and a brief description of the operation. Section 3 describes the process model and the method applied to find the optimal dynamic parameters. Section 4 analyzes the determination of the model parameters and the controller design based on the loop-shaping technique. Section 5 refers to optimizing the controller gains by using a simulator that incorporates process characteristics and disturbances. Finally, Section 6 summarizes the main findings and conclusions of this study.

Hardware Setup and Software Operation
The operating system of the cooler is PCS7, Version 6.1, by Siemens AG (Munich, Germany). The software extracting data from the plant process server is PI, Version 3.2.0.0, by OSIsoft (San Leandro, CA, USA). The PI client acquires the process data with a sampling period of 5 s and stores them on the PI server. The data saved in this server constitute the plant database. The developed software, in C#, extracts from the PI server, using specific queries, specified data of the cooler with a 20 s period, and stores them in a Structured Query Language (SQL) database. Afterwards, it processes the SQL data, attempting to determine the optimal PID sets constituting a PID optimizer. Figure 2 presents a simplified diagram of this process. The main tasks of this optimizer are:

•
Data extraction from PI and storing to SQL; • Data extraction from SQL and determination of the best models; • PID design; • Simulation of the process to find the optimal PID gains.

Process Model
The process model contains three transfer functions between the speed of the moving grate and the pressures of the static and moving grate, named G p1 , G p1g , and G g . The controller transfer function, G c , is also included to realize the feedback control loop presented in Figure 3. All the signals are expressed as percentages. The transfer function between the speed of the grate, Z, and static grate output, X 1 , is given by the N-order Equation (1) in the Laplace domain.
where S is the speed of the moving grate and P 1 represents the pressure of the fixed grate. The corresponding symbols of the steady state of the mentioned variables are S 0 and P 10 , respectively. The dynamic parameters are the gain, k vp , the time constant, T p , and the exponent, N p . A load disturbance, d, and noise, n p , are added to X 1 , providing the input X 1D to the second transfer function, G p1g , as shown in Equation (2).
A first-order dynamic describes the transfer function between static and moving grates, according to Equation (3). The output of this function, X, if different from zero, constitutes a disturbance for the moving grate operation.
where the gain, k vpg , and the time constant, T pg , are the dynamic parameters. The signals X and Z are summed up and inserted into the transfer function, G g , of the moving grate, expressed by the N-order Equation (4).
where P g is the pressure of the grate and P g0 is the corresponding steady-state pressure for a speed equal to S 0 . U is the sum of the signals X and Z. The set of dynamic parameters consists of the gain, k vg , the time constant, T g , and the exponent, N g . The addition of the variable n, representing the noise, to the signal, P g , derives the pressure, P. Equation (5) computes the control error which enters the controller transfer function, G c , described by Equation (6).
The proportional, integral, and derivative gain of the PID controller are given by the set (k P , k I , and k D ).
The sets of dynamic and steady-state parameters are the model parameters whose values need to be determined using actual process data. Process conditions and material characteristics affecting these parameters are:

•
The flow rate and granulometry of clinker dropping from the kiln to the cooler; • The flow rate, the number, and the quality of conventional and alternative fuels; • The possible existence of rings inside the kiln; • The pressures and probable clogging inside the preheater.
The variability of such conditions generates uncertainly as to parameters' values. The parameters of each transfer function are estimated using the convolution theorem between the input signal, u, and the output variable, y, expressed by Equation (7).
where g(t) is the pulse system response. The following operating data were collected: (a) the speed, S, and pressure, P, of the first moving grate and (b) the pressure, P 1 , of the static grate. Appropriate software sampled these data. The optimal model parameters have been estimated by minimizing the residual error, s Res , using the Levenberg-Marquardt algorithm [23].
where M is the number of data in the defined period, k 0 is the number of independent parameters, P(I) and P 1 (I) for I = 1 to M are the actual pressures of the moving and fixed grates, and P Calc (I) and P 1,Calc (I) for I = 1 to M are the calculated pressures for the same period. The regression coefficient, R Mod , is calculated by Equation (9).
where s P1 and s P are the standard deviations of the populations of P 1 and P, respectively. The sampling period corresponding to the time interval between two consecutive indices is T s . At the time I, the errors between calculated and actual values are given by Equation (10) and represent both the signal noise and the modeling error. Autoregressive (AR) Equations (11) and (12) model both errors.
The disturbance, d, is modeled by a rectangular pulse of height, D 0 , and duration, T D .

Determining the Best Models
A module of the software is dedicated to loading and processing data of the cooler, directly extracted from the plant database with a sampling period, T s , of 20 s. The selection of 20 s was made so that the sampling period is significantly shorter than the system time constants. The results should verify this decision. Along with the timestamp, the data of speed, pressure of the moving grate, and pressure of the fixed grate were extracted and then stored in an SQL database. Afterwards, the software downloads periods during which the kiln was operating from the database, removes the short downtime periods, and separates the remaining sets at 150-minute intervals. Table 1 shows the main characteristics of the cooler. Typical data sets, downloaded from the database and separated in 150-minute intervals, are demonstrated in Figure 4.  The following equations transform the process data into signals expressed as percentages: where S a is the speed of the first moving grate in RPM, P 1a the pressure of the fixed grate in mbar, and P a the pressure of the moving grate in mm H 2 O. Additionally, S aMin = 9.5 RPM, S aMax = 25 RPM, S Min = 33%, S Max = 55%, P 1aMax = 100 mbar, and P amax = 800 mm H 2 O. A wide variety of models described by Equations (1), (3) and (4) were studied, selecting an extended range of exponents, N p and N g . N p took all integer values in the interval (1,9), while N g took the integer values in the interval (2,9). Equations (17) and (18) calculate the residual errors for each experimental data series. For every period of 150 min, the software utilizes the results of the first 30 min to build the initial values by convolution, while those of the last 120 min were used to compare experimental and calculated pressures. Therefore, for T s = 20 s, M = 360. Equations (19) and (20) give the corresponding regression coefficients.
For each data set of 150-minute duration, the algorithm finds the optimal dynamic parameters first by minimizing the s Res and then calculates the other two residual errors and the regression coefficients. For some values of N p and N g , the parameters needing determination are (a) the gains, k vp , k vpg , and k vg ; (b) the time constant, T p , T pg , and T g ; and (c) the steady-state values, S 0 , P 10 , and P g0 , resulting in k 0 = 9. The parameters k vp , T p , S 0 , and P 10 are entered into Equation (17), thus k 1 = 4, while k vpg , k vg , T pg , T g , S 0 , and P g0 are included in Equation (18), therefore k 2 = 6. The processed data sets of 150-minute duration amounted to 949, resulting in approximately 450,000 groups of data. The validity of the modeling is assessed by comparing the achieved regression coefficients. The percentages of the data sets with regression coefficients, according to Equations (9), (19) and (20), higher than certain thresholds are presented in Tables 2-4.   The results of these three tables lead to the following conclusions: • As far as the model regression coefficient is concerned, R Mod , for 4 ≤ N g ≤ 7 and 4 ≤ N p ≤ 9, more than 77% of the data sets present a R Mod ≥ 0.7.

•
The regression coefficient, R P , between P, P 1 , and S is higher than 0.8 for more than 67% of the data sets if N g ≥ 4.

•
As to the corresponding regression coefficient between P 1 and S, there exist enough combinations of (N p , N g ) pairs where more than 42% of the total sets have a R P1 ≥ 0.7. • R P1 is generally less than R P due to the load disturbances coming from the kiln side.

•
The modeling verifies the measurable impact of the speed of the moving grate on the pressure of the static grate.
The determination of a suboptimum area in the plane (N p , N g ) is done as follows: in Table 3, suboptimal pairs (N p , N g ) are those where the percent is higher than 90% of the maximum value. Respectively, in Table 4, the pairs (N p , N g ) are defined as those with percentages higher than 95% of the highest value. The section of both regions, shown with the symbol "X" in the cells of Table 5, provides the area of exponents satisfying both suboptimal conditions. After selecting the data sets with a R Mod ≥ 0.7, the mean values and standard deviations of the dynamic and steady-state parameters for (N p , N g ) = (6, 6) and (N p , N g ) = (4, 4) were calculated and are shown in Table 6. The variance in the dynamic parameters is very high, probably causing problems in a robust automatic control strategy. As the exponents of the models increase, the three constants, T p , T pg and T g , decrease. The total time lag depends not only on the time constants but also on the exponent of the dynamics. For each set of exponents, the non-linear regression technique tries to find the time constants so that the adjustment of experimental and calculated values is optimal. The mean time constants are from 4.5 min to 6 min, 13 to 18 times bigger than the sampling period of 20 s. Therefore, the selection of this period is correct.  Figure 5a,b depicts the experimental and calculated results for two models. These figures are part of the user-friendly interface provided by the software. One can observe that the fitting of calculated results to the actual ones is better for the moving grate than for the fixed, regardless of the model. The reason is that the pressure, P 1 , is affected not only by the speed, S, but also by disturbances arriving from the kiln.

PID Controller Design
Guaranteeing robustness is a central theme in modern control theory. During the design and parameterization of the controller, when severe load disturbances or high uncertainty in the model parameters occur, the designer must ensure the stability and efficiency of the loop. A technique for adjusting PID controllers is the robust loopshaping method [24,25]. Such an effective loop-shaping technique is the MIGO (Mconstrained integral gain optimization) method, developed by Astrom, Panagopoulos, and Hagglund [2,3,16,26]. Tsamatsoulis successfully applied this method in other cases of automatic control of cement production processes [17,18,[20][21][22]. The approach is to maximize the integral gain for a certain value of the maximum sensitivity, M s , defined by Equations (21) and (22). The letter "M" in the MIGO method comes from the maximum sensitivity constraint.
where S is the sensitivity as a function of the frequency, ω, and G L the open-loop transfer function, given by Equation (23). This transfer function is calculated by interrupting the signal of pressure, P, before adding it to the set-point, SP, and opening the loop.
Among the models shown in Table 5, the model with N p = 6 and N g = 6 was chosen because it is close to the middle of the region. Then, for M s belonging to the interval of 1.1 to 1.7, and with a stepwise increase of 0.1, the sets (k P , k I , and k D ) are computed. The two independent variables are the set (M s , k D ) and the method that calculates the set (k P , k I ), satisfying the robustness constraint. The algorithm is given in [2] (pp. 206-221) and an analytical implementation in [17].
The Nyquist curve of a system is the locus of the complex number, G L (iω), as ω goes from 0 to ∞ [2] (p. 22). Figure 6 shows the Nyquist curve of the G L function for (k P , k I , and k D ) = (0.6922, 0.0305, and 1), corresponding to M s = 1.3. From this figure, the method becomes obvious: the algorithm for given values of M s , k D finds values of k P , k I , if such values exist, so that the curve G L is tangent to a circle, centered at the point (−1, 0) and with a radius of 1/M s . Figure 7 demonstrates the functions between k P , k I , and M s , k D .

Autoregressive Model
The algorithm checked the adequacy of the autoregressive model for Err P and Err P1 errors for the dynamic models with (N p , N g ) = (6,6) and (N p , N g ) = (4, 4). Table 7 presents (a) the coefficients; (b) average standard errors; (c) regression coefficients, R, of Equations (11) and (12); and (d) the values of s ResP1 , s ResP . Table 7. Autoregressive model parameters. One can conclude the following from these results: (a) the regression coefficients of both autocorrelation models are sufficiently high, implying that the errors s ResP1 and s ResP represent a colored noise; (b) the coefficients are very similar for both dynamic models; and (c) the sum of the two coefficients has a significantly smaller standard deviation than each coefficient separately, meaning that one can consider this sum as constant.

PID Optimization by Simulation
The control variable affects not only the process variable but also a variable that is a process disturbance, making, for this reason, the control loop shown in Figure 3 particularly interesting. In addition, high uncertainty characterizes the dynamics. Therefore, to find the optimal PID coefficients, a simulation is required that considers as many process characteristics as possible. The optimization criterion is the integral of absolute error between P and SP, as a percentage of SP defined by Equation (24). The integration time is an interval, t 0 . IAE is an acceptable criterion of optimization because a sluggish or much faster controller would give higher errors compared to a robust controller, especially in the presence of load disturbances.
The simulator can select the PID sets found by the MIGO technique in Section 4.2 and simulates the cooler for each (k p , k i , and k d ). Equation (25) provides the digital form of PID implementation.
where T s is the sampling and actuator period in seconds, e(I), e(I − 1), and e(I − 2) are the errors, e, given by Equation (5)

Short Description of the Simulator
The simulator uses the dynamic model with (N p , N g ) = (4, 4), which belongs to the optimal models' area and shows the lowest residence time compared to other models due to its lower exponents. According to Astrom [2] (pp. 25, 33), the average residence time of an N-order model with a time constant, T 0 , is N·T 0 . The simulation uses the median values of dynamic parameters as reference values. The steady-state parameters are not independent of each other: when S 0 increases, both P 10 and P g0 decrease. After grouping the steady-state parameters of the model with (N p , N g ) = (4, 4) and a R Mod ≥ 0.7, linear expressions, described by Equations (26)    The objective of the simulator is to find the PID sets, among those shown in Figure 7, which achieve the minimum IAE by performing both set-point tracking and load disturbance attenuation. The simulator starts at a moving grate speed and pressures equal to steadystate values so that the movable grate pressure is 75 mm H 2 O below the set-point. Table  9 shows the data used. The total operating time of the simulator is 360 min. During the first 30 min, the simulator performs the initial convolutions between input signals and frequency responses, and then each new point is calculated from the last 1800/T s points. The cooler runs in manual mode during these first 30 min with constant grate speed and then runs in automatic mode. A random white noise belonging to the interval (−s ResP , s ResP ) is added to the first two points, P, while for the following pressures, P, the noise is added according to Equation (11). The same procedure is followed for the pressure, P 1 , using the residual error, s ResP1 , and Equation (12). Disturbance, D, lasts from T dist,min to T dist,max minutes, and its amplitude is a random number belonging to the interval (−MP 10 /100·P 10 , MP 10 /100·P 10 ). The disturbance enters the system halfway through the operating time. The simulator every T s seconds provides the new speed based on Equation (25). In each of the N I simulations, the parameters change randomly around the reference values, within a range ± ma/100, where ma is the maximum disturbance margin. The value of the ma of each parameter is given in rows 8-11, 14-15, and 20-22 of Table 9. There is also a minimum grate speed, S Low , to prevent the grate from completely stopping when the thickness of the bed is too small. T g , min 7 6 P g0 , % From Table 8  7 S 0 , % From Table 8  8 Margin of disturbance of k vg , % 10 9 Margin of disturbance of T g , % 10 10 Margin of disturbance of P g0 , % 10 11 Margin of disturbance of S 0 , % 10 12 The simulator calculates IAE in each run, then the average IAE NI for all N I runs. Because the count of uncertainties is large and their magnitude is significant, a criterion for IAE convergence is required. If the simulator runs K times and N I times in each run, then the average IAE(K) is: Then if in the next K + 1 iteration satisfies condition (29), IAE is converged to IAE Aver . If not, further iterations are needed.

Initial Simulations
Before implementing the simulator, by using the full range of PID sets and parameter uncertainty, some preliminary simulations were performed to investigate the impact of maximum sensitivity function, M s , on the performance and robustness of the control.
The simulator executed three times with SP a = 400 mm H 2 O, N I = 1, total simulation time equal to 480 min, and the following controller parameters: For the initial simulations, the selection of a run time higher than that of Table 9 is necessary for checking the behavior of the controllers in a better manner in a longer operating time. Figure 8 depicts the controller's action with the different settings in setpoint tracking and attenuating the load disturbances. All controllers stabilize the process, but each has its characteristics. One can observe the following in Figure 8a. The velocity of reaching the set-point is very slow with the low M s controller. In particular, if the load disturbance results in a reduction of the static grate pressure, then it is hard for the controller to approach the set-point within a reasonable time. Such a disturbance occurs when the operator reduces the kiln speed to keep the material in the kiln. The response of the controller with the highest M s , in Figure 8c, has oscillations, and if the load disturbance is in the same direction as the rate of change of the moving grate pressure, the amplitude of the oscillations increases. The controller with M s = 1.3, in Figure 8b, is the more robust one: it reaches the set-point within a reasonable time and reacts with rapidly damped oscillation to an intense load disturbance.

Full Simulation and PID Optimization
The simulation applies all the PID sets calculated in Section 4.2 to the process data and the uncertainty margins shown in Table 9. The required number of iterations to fulfill condition (19) is 5·N I = 1000. Figure 9 shows the values of IAE Aver as a function of M s and k D for SP a = 400 mm H 2 O.
The calculations for determining the PID area, which presents the best IAE Aver , are the following: (a) The PID set providing the minimum IAE Aver,Min was found; (b) The sets with IAE Aver ≤ 1.1·IAE Aver,Min were determined; (c) This area constitutes the optimum region of the PID coefficients.    In conclusion, the combination of a robustness constraint with a performance criterion, maximum sensitivity, and integral of absolute error, respectively, leads to determining an optimal area of PID coefficients. PID sets located near the middle of this area can be chosen and implemented in the routine operation of the cooler under examination.
The study of the distribution of IAE for PIDs belonging to the optimal region was performed by applying the simulation for N I = 1000 and three sets (k P , k I , and k D ) with values shown in Table 10. Figure 11a,b demonstrates the percentiles of IAE for SP a = 400 mm H 2 O and SP a = 375 mm H 2 O, respectively. For SP a = 400 mm H 2 O, the three controllers are almost equivalent: the median IAE is~4%, while 80% of IAEs are less than~5%. On the other hand, for SP a = 375 mm H 2 O, the performance differences of the three controllers are not negligible: the smaller IAEs occur with the slow controller, with M s = 1.2 and k D = 0.5, and the differences become more significant for higher percentiles. Therefore, a weak impact of the set-point value on the optimal PID gains exists.

Conclusions
By implementing dynamic models of clinker cooling, efficient efforts have been made in optimizing the PID controller regulating the process. The modeling utilizes exclusively industrial data of the cooler of the Halyps cement plant. The process model contains three transfer functions between the speed of the moving grate and the pressures of the static and moving grate. The two transfer functions between the speed and the two pressures are expressed by N-order dynamics, while the third between the pressures of fixed and moving grate is expressed by first-order dynamics. The models include dynamic parameters such as gains and time constants, steady-state variables, and exponents of the dynamics. The best models' identification was achieved by implementing non-linear regression methods. The process model uses high-efficiency autoregressive equations to express the signal's noise and model mismatch and a rectangular pulse to describe the load disturbances.
Based on modeling results, the design of the PID controller follows according to a loop-shaping technique, which considers as a constraint the maximum sensitivity, M s , of the open-loop transfer function. The output of this algorithm is a set of PIDs satisfying a range of M s constraints. The control variable affects not only the process variable but also a variable that is a process disturbance, making, for this reason, the control loop particularly interesting. A simulator determines the optimal PID sets, among those calculated at the design step, using the IAE as a performance criterion, where a table of 35 inputs constitutes the simulator's settings. The combination of a robustness constraint with a performance criterion, maximum sensitivity, and integral of absolute error, respectively, leads to controllers with M s belonging to the range of 1.2 to 1.35. The IAE is between 4.2% and 4.8%, depending on the set-point value. These results are a strong argument that the cooler investigated can be effectively regulated by the designed PID. Residual error between pressures P and P Calc , % s ResP1 Residual error between pressures P 1 and P 1,Calc , % StdErr P Standard error of the Equation (11), % StdErr P1 Standard error of the Equation (12), % t 0 Integration time in Equation (24), min t D Duration of the disturbance, d, min T p Time constant of the transfer function G p1 , min T pg Time constant of the transfer function G p1g , min T s Sampling and actuation period, s T g Time constant of the transfer function G g , min U Sum of the signals X and Z, % X Output pressure from the transfer function G p1g , % X 1 Difference between P 1 and P 10 , % Y Difference between P g and P g0 , % X 1D Sum of the signals X 1 , d, and n p , % Z Difference between S and S 0 , %