Modelling of Lime Kiln Using Subspace Method with New Order Selection Criterion

This paper is taking actual control demand of rotary kiln as background and builds a calcining belt state space model using POMoesp subspace method. A novel order-delay double parameters error criterion (ODC) is presented to reduce the modeling order. The proposed subspace order identification method takes into account the influence of order and delay on model error criterion simultaneously. For the introduction of the delay factors, the order is reduced dramatically in the systemmodeling. Also, in the data processing part sliding-window method is adopted for stripping delay factor from historical data. For this, the parameters can be changed flexibly. Some practical problems in industrial kiln process modeling are also solved. Finally, it is applied to an industrial kiln case.


Introduction
Rotary kiln is the key equipment in metallurgy, building materials, and many other industries.A lot of rotary kiln control theories and methodologies have been studied, for the reason that kiln belongs to the typical complex process, which has the characteristics of multivariables, time delay, serious coupling, and difficult modelling [1][2][3][4][5].
In the field of rotary kiln modeling, researchers started studying internal status of calcinations process as early as 1960s.The initial researchers established control model based on solid mechanics, heat transfer mechanism, and dry calcined kinetics.Imber and Paschkis calculated the optimum kiln length for heating process in 1961 [6,7].These early approaches are all based upon the rmodynamic equilibrium and analysis of the kiln structure, and the establishment of these equations harshly demands material flow balance.So, it is not a good way to put in use widely.In 1980s dynamic system identification technology was used in rotary kiln modeling, but the identification model can only characterize part states of kiln.So this technology was restricted in application and exploitation.In the 1990s, after the emergence of the new intelligent modeling and control methods, a large number of intelligent modeling ways started to use in building model of kiln successfully [8][9][10].After a long-term development many mature experiences have already been obtained, but there are still many problems such as imprecise model, excess restrictions, and difficulty in promotion.Consequently, building model has become a bottleneck on the development road of kiln process control.
Subspace identification method as an effective identification modelling way for multi-input and multi-output system has drawn much attention recently [11][12][13].In comparison with the conventional methods, it has obvious advantage.This is because the model can be directly got from the input/output (/) data without nonlinear optimization and iteration.Furthermore, the computations are based on robust tools such as QR-factorization and singular value decomposition (SVD), for which numerically reliable algorithms are available [14,15].
Order and delay are the most important construction parameters of industrial process model.There exists an extensive literature for order estimation algorithms of linear, dynamical, state space systems [16][17][18][19][20].Among all of these methods, the most famous one is AIC (Akaike information criterion), which is introduced by Akaike [21].And also, great deal of subsequent researches were done about the properties effects of the penalty term [22,23].And then 2 Mathematical Problems in Engineering Baure introduces another criterion for the Larimore type of procedures, which is similar to AIC [24,25].From above literature, an obvious problem can be found that the delay parameter is always neglected in systems order selection.The information about the delay of any process is valuable for industrial process model.However, there exist only few references dealing with the estimation of the order in delay system [26][27][28].When the system has input delay, many methods attempt to increase the model order as the cost to improve model accuracy.This enhances the system complexity.
In this paper, the calcining belt state space model of rotary kiln is built by using PO-Moesp subspace method.After analyzing the main components technology and calcinations reaction mechanism in detail, the input and output variables of modeling are selected according to industrial field operational procedures.And more, a novel order-selecting method is put forward.The original error criterion is replaced by the new one, which includes the two parameters-the order and the delay.This improvement effectively solves the problem of model order too high when modelling object with time delay, and it also increases the precision of model.
The outline of this paper is as follows: in Section 2, the basic process of rotary kiln is described; in Section 3 the general subspace identification algorithm is introduced and the order selection problem is addressed.Section 4 presents the main results of this paper.The method to identify the order and the delay is deduced in detail.In Section 5, we discuss the performance of the proposed method by means of some simulation studies on an industrial rotary kiln illustration.

Lime Rotary Kiln Process Description
The rotary lime kiln is a long cylinder that has a 3-5 ∘ d slope, and it rotates about its axis at roughly 0.5-1.5 rpm.Lime mud is fed from the elevated cold end and moves down the kiln due to rotation and gravity.The hot end, in which the burner operates, is typically maintained at about 1200 ∘ C by burning fuel [1].Active lime rotary kiln system has many technological processes, such as long kiln technology system, chain grate cooling technology system, and vertical preheater technology system.The process is mainly divided into material system and air flow system.Raw material becomes product after preheating, high temperature calcinations, and cooling.Primary air, secondary air, and coke oven gas constitute the air flow system.The whole technological process is shown in Figure 1.In this paper, the lime kiln under study is 600 t/d, 4.0 m × 60 m.
The main calcinations process parameters are listed in Table 1.

Problem Formulation and Assumptions
A discrete-time model in state space form is described by the equations: where input   ∈   and output   ∈   ,  and  are input and output dimensions, respectively;   ∈   is the system state at time , and  is the system order to be identified; process noise vector   ∈   and measurement noise vector V  ∈   are zero mean value white noise series.
(a) The system is asymptotically stable.
(b) The system is observable and reachable.
The goal of subspace methods consist in the online estimating order , delay , and system matrices , , , and . where represent the past information, and  represent the future information,  ≥ ,  ≫ max(, ).
The state matrix   and   is defined as follows: The augmented observation matrix Γ Block-triangular matrix Toeplitz matrix   and The augmented observation inverse matrix of {, } is Suppose the input and state matrices are all row full rank matrices, and they are row space orthogonal.Then the augmented input/output matrices can be written as follows: The key step of Moesp method is to estimate the augmented observability matrix through the projection future / data onto past / data.The augmented observability matrix and state space vector estimation can be got through SVD decomposition.

General Algorithm.
In this section, we introduce the general subspace method framework which is proffered by Favoreel et al. [30].Subspace identification algorithm consists of two main steps.The first step always performs a weighted projection of the row space of the previously defined data Hankel matrices.From this projection, the augmented observability matrix Γ  and estimate X of the state sequence   can be retrieved.Then the order is got from SVD decomposition.In the second step, the system matrices , , ,  and , ,  are determined through least square method.The concrete Moesp algorithm is as follows.
Then, we can get (3) Carry SVD decomposition: And then take the number of nonzero eigenvalue as the system order rank(  ) = .
(4) The augmented observability matrix or X =   / ⊥   2 is derived from the third step.
Remark.By reference to [30], the weighting matrices  1 and  2 should satisfy the following three conditions: (1) rank( The first two conditions guarantee that the rank- property of Γ    is preserved after projection onto  ⊥  and weighting by  1 and  2 .The third condition expresses that  2 should be uncorrelated with the noise sequences () and V().By choosing the appropriate weighting matrices  1 and  2 , all subspace algorithms for LTI systems can be interpreted in the above framework, including N4SID, MOESP, CVA, Basic-4SID, and IV-4SID.

The Proposed Method
In the classical system identification theory, the actual model structure is usually assumed to be known.However, in practical it is always not clear.Subspace system identification method determines the order of the system by the nonzero eigenvalue of the augmented observability matrix.However, the system nonzero singular values may be very small.This may lead to the wrong system order and large identification error.
4.1.The Order-Delay Double Parameters Error Criterion.The most directly order-selection method is based on the error performance criterion.This idea is to choose the smallest possible order that keeps the error below a certain level.Then the MRSE (mean relative squared error) index is introduced by model error as follows: where   () − ŷ () is the model prediction error and  is the sample number.In [35] use the AIC, which was originally developed by Akaike and then adapted by Larimore for SMI.Given a set of samples, for a sequence of system order , for example,  ∈ [0 ⋅ ⋅ ⋅ 20], the order of the model will be the one which makes the following AIC index minimum: where For calculating the AIC() criterion, we first suppose the upper bound  max of the system order and then calculate the AIC  (1), AIC  (2), . . ., AIC  ( max ) sequence; the appropriate system order  is the one which decrees the AIC index obviously, and the order should be as small as possible.The MRSE and AIC index  MRSE () can be analyzed in the same manner.
However, the performance index based on a single order parameter cannot provide an effective solution to the delay system, which is shown in (1).This led to the problem that the original identification method had to increase the order as cost to improve model accuracy.Here, we introduce an orderdelay double parameters error criterion, which identifies the two key structural parameters at the same time.That means that the index () is changed into the (, ) form.

The Delay Factor
Stripping from Historical Data.The introduction of delay parameters in performance criterion has resulted to a notable problem.The modelling historical data matrices have already included delay information.It is difficult to change  in the performance criterion (, ) artificially.To solve this problem, the sliding-window method is adopted here.Sliding-window principle is shown in Figure 2. When new samples are added to the window, the oldest data inside the window will be discarded.
We use sliding window to change delay, which is shown in Figure 3. Suppose the output data length is ; select the input data region ( −  max ), . . ., (), . . ., ().Then the input data can be moved according to different delay from  to  max .

The Algorithm Description.
The detailed procedure of subspace identification based on ODC algorithm can be expressed as Algorithm 1.
Then the best combination  * ,  * and the optimum matching model parameters are all obtained.

Simulation Results
To demonstrate the superiority of the proposed order selection method in this paper over the conventional method, their performance is evaluated through a numerical example and an industrial illustrate.( The time delay is  = 5 and the sample is  = 1000.We use the zero mean white noise as the input  1 and  2 sequence to excite the system.The input and output curves are shown in Figure 4.
Firstly, we use conventional subspace model selection method, which is dependant on the number of eigenvalues.Figure 5  situation of delay = 5.According to this strategy, the order of model increases from 4 to 13 when the system has delay.
Also the MSE (mean squared error), MRSE (mean relative squared error), and AIC criteria are all tested as shown in Figures 6 and 7.These methods all have the distinct inflection point at 13, which is for the 4-order system with input delay 5. Obviously, for obtaining the ideal model of a delay system, these conventional methods have to increase the order.
Next, the performances of the proposed ODC method in this paper are presented.As the two search parameters are available, so the index performance is shown as a surface.The  MSE1 of output  1 is shown in Figure 8 and the  MSE2 of output  2 is shown in Figure 9. Three axes are the order , delay , and (⋅), respectively.The AIC index surface is shown in Figure 10.The minimum value of these indexes can be got easily; they are also the inflection points.The order and delay results are  = 4,  = 5.This is the same with the actual system model.
The corresponding identified model matrices For verifying the model, the model output error is drawn as  1 and  2 in Figure 11.

Example 2:
The Kiln Industrial Illustration.To demonstrate the superiority of the proposed order selection method in this paper over the conventional method, their performance is evaluated through an industrial illustration.The data come from actual kiln production data of an enterprise.Here, the gas flow and the second air flow are selected as the control input, and the calcination temperature and kiln tail temperature are taken as output variables.The sampling time is   = 1; then use Moesp method modelling the kiln based on the input/output data after preprocessing.Here, three main practical problems are solved.

The First Problem.
The problem is that these two input variables have different delay.They need to be identified, respectively; that is to say, a triple loop about  1 ,  2 , and  should be carried for solving (,  1 ,  2 ).
In order to obtain the inflection point information more directly, we identify the order  firstly.According to the industrial field situation, choose the possibility maximum values which are  1 max = 150,  2 max = 150.At first, travel all the possible  1 ,  2 and compute the smallest  MSE1 ,  MSE2 , and  MRSE corresponding to the different order .These curves are shown in Figure 12.As can be seen from the curves, all of these three error criteria achieve the inflection point at  = 5, so the order is got.Then specify  = 5; the triple loop is reduced to double loop, which is to compute  MSE1 and  MSE2 corresponding to  1 and  2 between [0, 150].

The Second Problem.
From Figure 13 we notice another problem that the outputs  1 and  2 generate different inflection point.
In order to solve this problem, the criterion should choose MRSE and AIC, which take into account  1 and  2 both together.Then get Figures 14 and 15.

5.2.3.
The Third Problem.However, there is still a problem which needs to be resolved.In general, the model with large time delay will increase the difficulty of controller designing.We know that when the sampling frequency is higher than the actual needed frequency, there will be lots of redundant data.And this will raise the model order and the delay.Therefore, the delay can be reduced by properly decreasing sampling frequency.Considering the kiln is a slow timevarying process, changing the sampling time from 1 s to 10 s will not affect the model accuracy.
From Table 2 it can be seen that when set   = 10 s, the order and inflexion point is still  = 5.We can also change the delay to  1 = 30/10 = 3,  2 = 100/10 = 10.The  MRSE surface can be got as in Figure 16; the results show that  1 = 3,  2 = 10.This is in agreement with the analysis before.
The corresponding rotary kiln calcining zone temperature model is Compare the measured value and the predicted value generated by the model in Figures 17 and 18.

Conclusion
In this paper, the calcining belt state space model of rotary kiln is built using PO-Moesp subspace method.And a novel double parameters error performance criterion for the order choosing in subspace modelling is introduced.Since the presented method considering the order and delay simultaneously, it can reduce the model order of the delay system effectively.And also, a strategy for stripping the delay factors from the historical data is also proposed.The algorithm is verified in identifying an industrial lime kiln.In this example, we solve the practical problem of industrial process with multidelay and also reduce the order by adjusting sampling time.Further research could shed more light on the issue of applying the model online.The problem in industrial field is more complex than simulation environment.How to extract problems from industrial practice and guide the direction of modeling research has become the study focus.

Figure 1 :
Figure 1: Technological process chart of rotary kiln.

Figure 6 :
Figure 6: Obtain the input delay system order by original error criterion methods.

4 Figure 7 :Figure 8 :
Figure 7: Obtain the input delay system order by AIC index.

Figure 9 :
Figure 9: The corresponding mean square error surface  MSE2 generated by ODC.

Figure 10 :
Figure 10: The corresponding AIC criterion surface generated by ODC.

Table 1 :
Process parameters of rotary kiln.

Table 2 :
When   = 10 s the order, the minmum  MRSE and the delay.