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

Authors: Li Zhang; Chengjin Zhang; Qingyang Xu; Chaoyang Wang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper is taking actual control demand of rotary kiln as background and builds a calcining belt state space model using PO-Moesp 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 system modeling. 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.


---

## Body

## 1. 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–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–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–13]. In comparison with the conventional methods, it has obvious advantage. This is because the model can be directly got from the input/output (I/O) 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–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 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–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 Section2, 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.


## 2. 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.Figure 1
Technological process chart of rotary kiln.The main calcinations process parameters are listed in Table1.Table 1
Process parameters of rotary kiln.
Parameter
The parameter of kiln
Raw material granularity
10~40 mm
Preheater temperature
900~1100°C
Kiln head temperature
800~950°C
Kiln tail temperature
1100 ± 20°C
Calcining zone temperature
1260~1340°C
Second air temperature
750~900°C
Kiln head cover temperature
800~950°C
Exhaust gas temperature
250°C
Lime temperature after cooling
≤150°C
Kiln head pressure range
−10~−20 pa
Gas pressure
≥6 KPa
Calcinations process pressure
Negative pressure 100–150 Pa
Pulse jet cleaning pressure
0.35~0.5 MPa
3920 dust removel pressure difference
≤1.5 kPa
Kiln body inclination
3~5%
Kiln speed range
0.2~1.28 r/min
CaO content
>88%
Activity degree
>350 mL
Exhaust gas residual carbon
<=1%
Exhaust gas oxygen content
0.7%~2.5%
Dust collector air volume
225000 m3/h
Material level in preheater
5~8.8 m
Material level in cooler
≤0.55 m
production
600 T~650 T
SiO2
<1.5% or <2%
MgO
5%
Power consumption
65 kw·h/t~60 kw·h/t
Heat consumption
1630~1650 kJ/kg


## 3. Problem Formulation and Assumptions

A discrete-time model in state space form is described by the equations:(1)x(k+1)=Ax(k)+Bu(k-d)+w(k),y(k)=Cx(k)+Du(k-d)+v(k),
where input uk∈Rm and output yk∈Rl, m and l are input and output dimensions, respectively; xk∈Rn is the system state at time k, and n is the system order to be identified; process noise vector wk∈Rn and measurement noise vector vk∈Rl are zero mean value white noise series.We introduce the following hypotheses (van Overschee and de Moor, 1995) [29].(a)
The system is asymptotically stable.(b)
The system is observable and reachable.(c)
All the process noisew(k) and measurement noise v(k) are statistically independent of the input u(k), that is mean: E-{u(i)w(j)}=E-{u(i)v(j)}=0 where E-{·}=lim⁡N→∞(1/N)∑E(·).The goal of subspace methods consist in the online estimating ordern, delay d, and system matrices A, B, C, and D.


### 3.1. Matrices Construct

Suppose at timek the past input and output data samples are available; then construct past/future output vectors and the Hankel output matrices as follows:(2a)Up=(u0u1⋯uj-1u1u2⋯uj⋮⋮⋱⋮ui-1ui⋯ui+j-2),⋮⋮(2b)22Yf=(yiyi+1⋯yi+j-1yi+1yi+2⋯yi+j⋮⋮⋱⋮y2i-1y2i⋯y2i+j-2),where ui=[ui1,ui2,…,uil], yi=[yi1,yi2,…,yim], p represent the past information, and f represent the future information, i≥n, j≫max⁡(il,im).The state matrixXp and Xf is defined as follows: Xp=[x0x1⋯xj-1], Xf=[xixi+1⋯xj+i-1].The augmented observation matrixΓi(3)Γi=[CCA⋯CAi]T.Block-triangular matrix Toeplitz matrixHi and His(4)Hi=[D0⋯0CBD⋱⋮⋮⋮⋱0Ci-2BCi-1B⋯D],His=[00⋯0C0⋱⋮⋮⋮⋱0Ci-2Ci-1⋯0].The augmented observation inverse matrix of{A,B} is Δi=[Ai-1BAi-2B⋯B]. And the augmented observation inverse matrix of {A,K} is Δis=[Ai-1KAi-2K⋯K].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:(5a)Xf=AiXp+ΔiUp+ΔisEf,(5b)Yp=ΓiXp+HiUp+HisEp,(5c)Yf=ΓiXf+HiUf+HisEf.The key step of Moesp method is to estimate the augmented observability matrix through the projection futureI/O data onto past I/O data. The augmented observability matrix and state space vector estimation can be got through SVD decomposition.


### 3.2. 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 Γi and estimate X^i of the state sequence Xi can be retrieved. Then the order is got from SVD decomposition. In the second step, the system matrices A, B, C, D and Q, S, R are determined through least square method. The concrete Moesp algorithm is as follows.(1)
Project theYf row space into the orthogonal complement of the Uf row space:
(6)Yf/Uf⊥=ΓiXf/Uf⊥+HidUf/Uf⊥+HisEf/Uf⊥.Since it is assumed that the noise is uncorrelated with the inputs, soEf/Uf⊥=Ef and Uf/Uf⊥=0; therefore Yf/Uf⊥=ΓiXf/Uf⊥+HisEf.(2)
Select the weighting matricesW1 and W2:
(7)W1YfUf⊥W2=W1ΓiXfUf⊥W2+W1HisEfW2.The weighting matrices can be chosen appropriately according to different subspace methods, including N4SID, MOESP, CVA, basic-4SID, and IV-4SID [31–34].Then, we can get(8)oi=W1YfUf⊥W2=W1ΓiXfUf⊥W2.(3)
Carry SVD decomposition:(9)oi=(U1U2)(S1000)(V1TV2T).And then take the number of nonzero eigenvalue as the system orderrank⁡(oi)=n.(4)
The augmented observability matrixΓi=W1-1U1S11/2 or X^i=Xi/Uf⊥W2 is derived from the third step.(5)
Extract estimateA, B, C, D from Γi or X^i.Remark. By reference to [30], the weighting matrices W1 and W2 should satisfy the following three conditions:(1)
rank⁡(W1·Γi)=rank⁡Γi;(2)
rank⁡(Xi/Uf⊥·W2)=rank⁡Xi;(3)
W1·(HisMf+Nf)·W2=0.The first two conditions guarantee that the rank-n property of ΓiXi is preserved after projection onto Uf⊥ and weighting by W1 and W2. The third condition expresses that W2 should be uncorrelated with the noise sequences w(k) and v(k). By choosing the appropriate weighting matrices W1 and W2, all subspace algorithms for LTI systems can be interpreted in the above framework, including N4SID, MOESP, CVA, Basic-4SID, and IV-4SID.


## 3.1. Matrices Construct

Suppose at timek the past input and output data samples are available; then construct past/future output vectors and the Hankel output matrices as follows:(2a)Up=(u0u1⋯uj-1u1u2⋯uj⋮⋮⋱⋮ui-1ui⋯ui+j-2),⋮⋮(2b)22Yf=(yiyi+1⋯yi+j-1yi+1yi+2⋯yi+j⋮⋮⋱⋮y2i-1y2i⋯y2i+j-2),where ui=[ui1,ui2,…,uil], yi=[yi1,yi2,…,yim], p represent the past information, and f represent the future information, i≥n, j≫max⁡(il,im).The state matrixXp and Xf is defined as follows: Xp=[x0x1⋯xj-1], Xf=[xixi+1⋯xj+i-1].The augmented observation matrixΓi(3)Γi=[CCA⋯CAi]T.Block-triangular matrix Toeplitz matrixHi and His(4)Hi=[D0⋯0CBD⋱⋮⋮⋮⋱0Ci-2BCi-1B⋯D],His=[00⋯0C0⋱⋮⋮⋮⋱0Ci-2Ci-1⋯0].The augmented observation inverse matrix of{A,B} is Δi=[Ai-1BAi-2B⋯B]. And the augmented observation inverse matrix of {A,K} is Δis=[Ai-1KAi-2K⋯K].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:(5a)Xf=AiXp+ΔiUp+ΔisEf,(5b)Yp=ΓiXp+HiUp+HisEp,(5c)Yf=ΓiXf+HiUf+HisEf.The key step of Moesp method is to estimate the augmented observability matrix through the projection futureI/O data onto past I/O data. The augmented observability matrix and state space vector estimation can be got through SVD decomposition.


## 3.2. 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 Γi and estimate X^i of the state sequence Xi can be retrieved. Then the order is got from SVD decomposition. In the second step, the system matrices A, B, C, D and Q, S, R are determined through least square method. The concrete Moesp algorithm is as follows.(1)
Project theYf row space into the orthogonal complement of the Uf row space:
(6)Yf/Uf⊥=ΓiXf/Uf⊥+HidUf/Uf⊥+HisEf/Uf⊥.Since it is assumed that the noise is uncorrelated with the inputs, soEf/Uf⊥=Ef and Uf/Uf⊥=0; therefore Yf/Uf⊥=ΓiXf/Uf⊥+HisEf.(2)
Select the weighting matricesW1 and W2:
(7)W1YfUf⊥W2=W1ΓiXfUf⊥W2+W1HisEfW2.The weighting matrices can be chosen appropriately according to different subspace methods, including N4SID, MOESP, CVA, basic-4SID, and IV-4SID [31–34].Then, we can get(8)oi=W1YfUf⊥W2=W1ΓiXfUf⊥W2.(3)
Carry SVD decomposition:(9)oi=(U1U2)(S1000)(V1TV2T).And then take the number of nonzero eigenvalue as the system orderrank⁡(oi)=n.(4)
The augmented observability matrixΓi=W1-1U1S11/2 or X^i=Xi/Uf⊥W2 is derived from the third step.(5)
Extract estimateA, B, C, D from Γi or X^i.Remark. By reference to [30], the weighting matrices W1 and W2 should satisfy the following three conditions:(1)
rank⁡(W1·Γi)=rank⁡Γi;(2)
rank⁡(Xi/Uf⊥·W2)=rank⁡Xi;(3)
W1·(HisMf+Nf)·W2=0.The first two conditions guarantee that the rank-n property of ΓiXi is preserved after projection onto Uf⊥ and weighting by W1 and W2. The third condition expresses that W2 should be uncorrelated with the noise sequences w(k) and v(k). By choosing the appropriate weighting matrices W1 and W2, all subspace algorithms for LTI systems can be interpreted in the above framework, including N4SID, MOESP, CVA, Basic-4SID, and IV-4SID.


## 4. 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:(10)JMRSE(n)=1L∑k=1L∑j=1ny(yk(j)-y^k(j))2∑j=1nyyk(j)2,
where yk(j)-y^k(j) is the model prediction error and L 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 n, for example, n∈[0⋯20], the order of the model will be the one which makes the following AIC index minimum:
(11)AICn(n)=N(m(1+ln⁡⁡2π)+ln⁡|Σn|)+2δnMn,
where
(12)Σn=1N∑i=1Ne(k)e(k)T,e(k)=yn(k)-y^n(k),Mn=2nm+m(m+1)2+nl+ml,δn=NN-((Mn/m)+((m+1)/2)).For calculating theAIC(n) criterion, we first suppose the upper bound nmax⁡ of the system order and then calculate the AICN(1),AICN(2),…,AICN(nmax⁡) sequence; the appropriate system order n is the one which decrees the AIC index obviously, and the order should be as small as possible. The MRSE and AIC index JMRSE(n) 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 order-delay double parameters error criterion, which identifies the two key structural parameters at the same time. That means that the index J(n) is changed into the J(n,d) form.For each given individuald, a state space model can be identified using the Moesp algorithm described in Section 2. Then, its model error can be deserved as Σn,d=(1/N)∑i=1Ne(k)e(k)T, then the AIC(n,d) with respect to individual (n,d) as
(13)AICN(n,d)=N(m(1+ln⁡⁡2π)+ln⁡⁡|Σn,d|)+2δnMn.Also, the MRSE criterionJMSE2 has the similar form:
(14)JMRSE(n,d)=1L∑k=1L∑j=1ny(yn,d,k(j)-y^n,d,k(j))2∑j=1nyyn,d,k(j)2.Other performance indexJ(·) such as SVC IVC, NIC criteria mentioned in [25] can also be modified as this method.The original performance index just identifies the ordern. Suppose n∈[1,nmax⁡]; then J(1),J(2),J(3),…,J(nmax⁡) are calculated, respectively; then the inflection point n* is the best order and the corresponding system matrices A*, B*, C*, D* are the best model. After improvement, we add the delay as other optimal parameters. So the system order n∈[1,nmax⁡] and input delay d∈[1,dmax⁡] are all embedded in J(n,d). Then calculate J(1,1),…,J(nmax⁡,1),J(1,2),…,J(nmax⁡,2), J(1,dmax⁡),…,J(nmax⁡,dmax⁡), respectively. By searching the minimum point n*,d* of surface J(n,d), the best order, delay, and system matrices can all be got. Thus, taking the delay as another parameter in the modelling methods, it can effectively avoid high order results in the delay system.


### 4.2. 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 changed in the performance criterion J(n,d) 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.Figure 2
Sliding time window sketch map.We use sliding window to change delay, which is shown in Figure3. Suppose the output data length is L; select the input data region u(k-dmax⁡),…,u(k),…,u(L). Then the input data can be moved according to different delay from d to dmax⁡.Figure 3
Modeling data sliding sketch map.


### 4.3. The Algorithm Description

The detailed procedure of subspace identification based on ODC algorithm can be expressed as Algorithm1.Algorithm 1
Initializationn=2, d=1, the modelling data after pretreatment: u(k-dmax⁡),…,u(k),…,u(L) and y(k),…,y(L)
(1)  forn=2 to n=nmax⁡
(2)    ford=1 to d=dmax⁡
(3)    Rolling the modelling input datau based on the hypothesis delay d, get the data set u(k-d),…,u(k),…,u(L-d),
y(k),…,y(L);
(4)    Construct input and output Hankel matricesUp, Yp, Uf, Yf;
(5)    CalculateA, B, C, D by Moesp method in Section 3 based on Hankel matrices;
(6)    Substitute theA, B, C, D and d into the formula (1);
(7)    Calculate and store the model error and performance indexJ(n,d);
(8)end for
(9)  end for
(10) Search the inflection point of theJ(n,d) surface.Then the best combinationn*, d* and the optimum matching model parameters are all obtained.


## 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:(10)JMRSE(n)=1L∑k=1L∑j=1ny(yk(j)-y^k(j))2∑j=1nyyk(j)2,
where yk(j)-y^k(j) is the model prediction error and L 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 n, for example, n∈[0⋯20], the order of the model will be the one which makes the following AIC index minimum:
(11)AICn(n)=N(m(1+ln⁡⁡2π)+ln⁡|Σn|)+2δnMn,
where
(12)Σn=1N∑i=1Ne(k)e(k)T,e(k)=yn(k)-y^n(k),Mn=2nm+m(m+1)2+nl+ml,δn=NN-((Mn/m)+((m+1)/2)).For calculating theAIC(n) criterion, we first suppose the upper bound nmax⁡ of the system order and then calculate the AICN(1),AICN(2),…,AICN(nmax⁡) sequence; the appropriate system order n is the one which decrees the AIC index obviously, and the order should be as small as possible. The MRSE and AIC index JMRSE(n) 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 order-delay double parameters error criterion, which identifies the two key structural parameters at the same time. That means that the index J(n) is changed into the J(n,d) form.For each given individuald, a state space model can be identified using the Moesp algorithm described in Section 2. Then, its model error can be deserved as Σn,d=(1/N)∑i=1Ne(k)e(k)T, then the AIC(n,d) with respect to individual (n,d) as
(13)AICN(n,d)=N(m(1+ln⁡⁡2π)+ln⁡⁡|Σn,d|)+2δnMn.Also, the MRSE criterionJMSE2 has the similar form:
(14)JMRSE(n,d)=1L∑k=1L∑j=1ny(yn,d,k(j)-y^n,d,k(j))2∑j=1nyyn,d,k(j)2.Other performance indexJ(·) such as SVC IVC, NIC criteria mentioned in [25] can also be modified as this method.The original performance index just identifies the ordern. Suppose n∈[1,nmax⁡]; then J(1),J(2),J(3),…,J(nmax⁡) are calculated, respectively; then the inflection point n* is the best order and the corresponding system matrices A*, B*, C*, D* are the best model. After improvement, we add the delay as other optimal parameters. So the system order n∈[1,nmax⁡] and input delay d∈[1,dmax⁡] are all embedded in J(n,d). Then calculate J(1,1),…,J(nmax⁡,1),J(1,2),…,J(nmax⁡,2), J(1,dmax⁡),…,J(nmax⁡,dmax⁡), respectively. By searching the minimum point n*,d* of surface J(n,d), the best order, delay, and system matrices can all be got. Thus, taking the delay as another parameter in the modelling methods, it can effectively avoid high order results in the delay system.


## 4.2. 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 changed in the performance criterion J(n,d) 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.Figure 2
Sliding time window sketch map.We use sliding window to change delay, which is shown in Figure3. Suppose the output data length is L; select the input data region u(k-dmax⁡),…,u(k),…,u(L). Then the input data can be moved according to different delay from d to dmax⁡.Figure 3
Modeling data sliding sketch map.


## 4.3. The Algorithm Description

The detailed procedure of subspace identification based on ODC algorithm can be expressed as Algorithm1.Algorithm 1
Initializationn=2, d=1, the modelling data after pretreatment: u(k-dmax⁡),…,u(k),…,u(L) and y(k),…,y(L)
(1)  forn=2 to n=nmax⁡
(2)    ford=1 to d=dmax⁡
(3)    Rolling the modelling input datau based on the hypothesis delay d, get the data set u(k-d),…,u(k),…,u(L-d),
y(k),…,y(L);
(4)    Construct input and output Hankel matricesUp, Yp, Uf, Yf;
(5)    CalculateA, B, C, D by Moesp method in Section 3 based on Hankel matrices;
(6)    Substitute theA, B, C, D and d into the formula (1);
(7)    Calculate and store the model error and performance indexJ(n,d);
(8)end for
(9)  end for
(10) Search the inflection point of theJ(n,d) surface.Then the best combinationn*, d* and the optimum matching model parameters are all obtained.


## 5. 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.


### 5.1. Example  1: MIMO Process with Input Delay

The first model (15) is a MIMO system:
(15)x(k+1)=[0.6030.60300-0.6030.6030000-0.603-0.603000.603-0.603]x(k)+[1.1650-0.69650.62681.69610.07510.05910.35161.7971]u(k-5),y(k)=[0.2641-1.44621.24600.57740.8717-0.7012-0.6390-0.3600]x(k)+[-0.1356-1.2704-1.34930.9846]u(k-5).The time delay isd=5 and the sample is L=1000. We use the zero mean white noise as the input u1 and u2 sequence to excite the system. The input and output curves are shown in Figure 4.PO-Moesp subspace system identification input/output signal.
(a)
Input signalu1
(b)
Input signalu2
(c)
Output signaly1
(d)
Output signaly2Firstly, we use conventional subspace model selection method, which is dependant on the number of eigenvalues. Figure5(a) is the situation of delay = 0 and Figure 5(b) is the situation of delay = 5. According to this strategy, the order of model increases from 4 to 13 when the system has delay.The number of eigenvalues generated whend = 0 and d = 5, respectively, (a) d = 0 and (b) d = 5.
(a)
d
=
0 situation
(b)
d
=
5 situationAlso the MSE (mean squared error), MRSE (mean relative squared error), and AIC criteria are all tested as shown in Figures6 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.Figure 6
Obtain the input delay system order by original error criterion methods.Figure 7
Obtain the input delay system order by AIC index.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. TheJMSE1 of output y1 is shown in Figure 8 and the JMSE2 of output y2 is shown in Figure 9. Three axes are the order n, delay d, and J(·), 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 n=4,d=5. This is the same with the actual system model.Figure 8
The corresponding mean square error surfaceJMSE1 generated by ODC.Figure 9
The corresponding mean square error surfaceJMSE2 generated by ODC.Figure 10
The corresponding AIC criterion surface generated by ODC.The corresponding identified model matrices(16)A^=[0.50170.60470.20240.1514-0.62560.51640.3929-0.08280.22280.3763-0.4884-0.60900.1883-0.09360.6069-0.5297],B^=[-0.03861.8602-1.00280.0809-0.5982-1.0787-0.0583-0.4779],C^=[-1.42761.0033-1.11610.3224-1.1801-0.57790.4116-0.4398],D^=[-0.1356-1.2704-1.34930.9846].For verifying the model, the model output error is drawn ase1 and e2 in Figure 11.Modeling error curvese1 (a) and e2 (b).
(a)
e
1
(b)
e
2


### 5.2. 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 isTs=1s; then use Moesp method modelling the kiln based on the input/output data after preprocessing. Here, three main practical problems are solved.


#### 5.2.1. 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 aboutd1, d2, and n should be carried for solving J(n,d1,d2).In order to obtain the inflection point information more directly, we identify the ordern firstly. According to the industrial field situation, choose the possibility maximum values which are d1max⁡=150, d2max⁡=150. At first, travel all the possible d1, d2 and compute the smallest JMSE1, JMSE2, and JMRSE corresponding to the different order n. 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 n=5, so the order is got. Then specify n=5; the triple loop is reduced to double loop, which is to compute JMSE1 and JMSE2 corresponding to d1 and d2 between [0,150].Figure 12
The corresponding minimumJ(n,d1,d2) curves of each order when traversal d1, d2.


#### 5.2.2. The Second Problem

From Figure13 we notice another problem that the outputs y1 and y2 generate different inflection point.Figure 13
J
MSE
1 and JMSE2 surface generated by ODC when n=5.In order to solve this problem, the criterion should choose MRSE and AIC, which take into accounty1 and y2 both together. Then get Figures 14 and 15.Figure 14
AIC surface generated by ODC whenn=5.Figure 15
J
MRSE surface generated by ODC when n=5.It can be conducted thatn=5,d1=30,d2=100.


#### 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 time-varying process, changing the sampling time from 1 s to 10 s will not affect the model accuracy.From Table2 it can be seen that when set Ts=10s, the order and inflexion point is still n=5. We can also change the delay to d1=30/10=3, d2=100/10=10. The JMRSE surface can be got as in Figure 16; the results show that d1=3,d2=10. This is in agreement with the analysis before.Table 2
WhenTs = 10 s the order, the minmum JMRSE and the delay.
n
1
2
3
4
5
6
7
8
d
1
4
3
3
3
3
3
3
3
d
2
5
10
10
10
10
10
10
10
J
MRSE
698.1251
59.5572
49.9145
38.4619
23.9931
23.9931
23.9931
23.9931Figure 16
J
MRSE surface generated by ODC when n=5 and Ts=10s.The corresponding rotary kiln calcining zone temperature model is(17)x(k+1)=Ax(k)+Bu(k-d1),y(k)=Cx(k)+Du(k-d2)
and the order n=5; delay d1=3,d2=10;
(18)A=[0.99360.0015-0.0062-0.0007-0.00010.00630.99060.0170-0.0017-0.0001-0.00090.00260.9793-0.01470.02090.0000-0.00010.00180.9985-0.0069-0.00000.0001-0.00150.00620.9848],B=10-3×[-0.0908-0.0129-0.16500.09570.0680-0.0329-0.07340.0191-0.0093-0.0012],C=[-5.79952.9327-0.5841-0.01570.0260-6.7753-2.3024-0.0247-0.0596-0.0063],D=10-15×[0.12050.2699-0.10660.0143].Compare the measured value and the predicted value generated by the model in Figures17 and 18.Figure 17
Comparison of calcination temperature measured curve and model predictive curve.Figure 18
Comparison of kiln tail temperature measured curve and model predictive curve.


## 5.1. Example  1: MIMO Process with Input Delay

The first model (15) is a MIMO system:
(15)x(k+1)=[0.6030.60300-0.6030.6030000-0.603-0.603000.603-0.603]x(k)+[1.1650-0.69650.62681.69610.07510.05910.35161.7971]u(k-5),y(k)=[0.2641-1.44621.24600.57740.8717-0.7012-0.6390-0.3600]x(k)+[-0.1356-1.2704-1.34930.9846]u(k-5).The time delay isd=5 and the sample is L=1000. We use the zero mean white noise as the input u1 and u2 sequence to excite the system. The input and output curves are shown in Figure 4.PO-Moesp subspace system identification input/output signal.
(a)
Input signalu1
(b)
Input signalu2
(c)
Output signaly1
(d)
Output signaly2Firstly, we use conventional subspace model selection method, which is dependant on the number of eigenvalues. Figure5(a) is the situation of delay = 0 and Figure 5(b) is the situation of delay = 5. According to this strategy, the order of model increases from 4 to 13 when the system has delay.The number of eigenvalues generated whend = 0 and d = 5, respectively, (a) d = 0 and (b) d = 5.
(a)
d
=
0 situation
(b)
d
=
5 situationAlso the MSE (mean squared error), MRSE (mean relative squared error), and AIC criteria are all tested as shown in Figures6 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.Figure 6
Obtain the input delay system order by original error criterion methods.Figure 7
Obtain the input delay system order by AIC index.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. TheJMSE1 of output y1 is shown in Figure 8 and the JMSE2 of output y2 is shown in Figure 9. Three axes are the order n, delay d, and J(·), 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 n=4,d=5. This is the same with the actual system model.Figure 8
The corresponding mean square error surfaceJMSE1 generated by ODC.Figure 9
The corresponding mean square error surfaceJMSE2 generated by ODC.Figure 10
The corresponding AIC criterion surface generated by ODC.The corresponding identified model matrices(16)A^=[0.50170.60470.20240.1514-0.62560.51640.3929-0.08280.22280.3763-0.4884-0.60900.1883-0.09360.6069-0.5297],B^=[-0.03861.8602-1.00280.0809-0.5982-1.0787-0.0583-0.4779],C^=[-1.42761.0033-1.11610.3224-1.1801-0.57790.4116-0.4398],D^=[-0.1356-1.2704-1.34930.9846].For verifying the model, the model output error is drawn ase1 and e2 in Figure 11.Modeling error curvese1 (a) and e2 (b).
(a)
e
1
(b)
e
2


## 5.2. 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 isTs=1s; then use Moesp method modelling the kiln based on the input/output data after preprocessing. Here, three main practical problems are solved.


### 5.2.1. 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 aboutd1, d2, and n should be carried for solving J(n,d1,d2).In order to obtain the inflection point information more directly, we identify the ordern firstly. According to the industrial field situation, choose the possibility maximum values which are d1max⁡=150, d2max⁡=150. At first, travel all the possible d1, d2 and compute the smallest JMSE1, JMSE2, and JMRSE corresponding to the different order n. 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 n=5, so the order is got. Then specify n=5; the triple loop is reduced to double loop, which is to compute JMSE1 and JMSE2 corresponding to d1 and d2 between [0,150].Figure 12
The corresponding minimumJ(n,d1,d2) curves of each order when traversal d1, d2.


### 5.2.2. The Second Problem

From Figure13 we notice another problem that the outputs y1 and y2 generate different inflection point.Figure 13
J
MSE
1 and JMSE2 surface generated by ODC when n=5.In order to solve this problem, the criterion should choose MRSE and AIC, which take into accounty1 and y2 both together. Then get Figures 14 and 15.Figure 14
AIC surface generated by ODC whenn=5.Figure 15
J
MRSE surface generated by ODC when n=5.It can be conducted thatn=5,d1=30,d2=100.


### 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 time-varying process, changing the sampling time from 1 s to 10 s will not affect the model accuracy.From Table2 it can be seen that when set Ts=10s, the order and inflexion point is still n=5. We can also change the delay to d1=30/10=3, d2=100/10=10. The JMRSE surface can be got as in Figure 16; the results show that d1=3,d2=10. This is in agreement with the analysis before.Table 2
WhenTs = 10 s the order, the minmum JMRSE and the delay.
n
1
2
3
4
5
6
7
8
d
1
4
3
3
3
3
3
3
3
d
2
5
10
10
10
10
10
10
10
J
MRSE
698.1251
59.5572
49.9145
38.4619
23.9931
23.9931
23.9931
23.9931Figure 16
J
MRSE surface generated by ODC when n=5 and Ts=10s.The corresponding rotary kiln calcining zone temperature model is(17)x(k+1)=Ax(k)+Bu(k-d1),y(k)=Cx(k)+Du(k-d2)
and the order n=5; delay d1=3,d2=10;
(18)A=[0.99360.0015-0.0062-0.0007-0.00010.00630.99060.0170-0.0017-0.0001-0.00090.00260.9793-0.01470.02090.0000-0.00010.00180.9985-0.0069-0.00000.0001-0.00150.00620.9848],B=10-3×[-0.0908-0.0129-0.16500.09570.0680-0.0329-0.07340.0191-0.0093-0.0012],C=[-5.79952.9327-0.5841-0.01570.0260-6.7753-2.3024-0.0247-0.0596-0.0063],D=10-15×[0.12050.2699-0.10660.0143].Compare the measured value and the predicted value generated by the model in Figures17 and 18.Figure 17
Comparison of calcination temperature measured curve and model predictive curve.Figure 18
Comparison of kiln tail temperature measured curve and model predictive curve.


## 5.2.1. 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 aboutd1, d2, and n should be carried for solving J(n,d1,d2).In order to obtain the inflection point information more directly, we identify the ordern firstly. According to the industrial field situation, choose the possibility maximum values which are d1max⁡=150, d2max⁡=150. At first, travel all the possible d1, d2 and compute the smallest JMSE1, JMSE2, and JMRSE corresponding to the different order n. 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 n=5, so the order is got. Then specify n=5; the triple loop is reduced to double loop, which is to compute JMSE1 and JMSE2 corresponding to d1 and d2 between [0,150].Figure 12
The corresponding minimumJ(n,d1,d2) curves of each order when traversal d1, d2.


## 5.2.2. The Second Problem

From Figure13 we notice another problem that the outputs y1 and y2 generate different inflection point.Figure 13
J
MSE
1 and JMSE2 surface generated by ODC when n=5.In order to solve this problem, the criterion should choose MRSE and AIC, which take into accounty1 and y2 both together. Then get Figures 14 and 15.Figure 14
AIC surface generated by ODC whenn=5.Figure 15
J
MRSE surface generated by ODC when n=5.It can be conducted thatn=5,d1=30,d2=100.


## 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 time-varying process, changing the sampling time from 1 s to 10 s will not affect the model accuracy.From Table2 it can be seen that when set Ts=10s, the order and inflexion point is still n=5. We can also change the delay to d1=30/10=3, d2=100/10=10. The JMRSE surface can be got as in Figure 16; the results show that d1=3,d2=10. This is in agreement with the analysis before.Table 2
WhenTs = 10 s the order, the minmum JMRSE and the delay.
n
1
2
3
4
5
6
7
8
d
1
4
3
3
3
3
3
3
3
d
2
5
10
10
10
10
10
10
10
J
MRSE
698.1251
59.5572
49.9145
38.4619
23.9931
23.9931
23.9931
23.9931Figure 16
J
MRSE surface generated by ODC when n=5 and Ts=10s.The corresponding rotary kiln calcining zone temperature model is(17)x(k+1)=Ax(k)+Bu(k-d1),y(k)=Cx(k)+Du(k-d2)
and the order n=5; delay d1=3,d2=10;
(18)A=[0.99360.0015-0.0062-0.0007-0.00010.00630.99060.0170-0.0017-0.0001-0.00090.00260.9793-0.01470.02090.0000-0.00010.00180.9985-0.0069-0.00000.0001-0.00150.00620.9848],B=10-3×[-0.0908-0.0129-0.16500.09570.0680-0.0329-0.07340.0191-0.0093-0.0012],C=[-5.79952.9327-0.5841-0.01570.0260-6.7753-2.3024-0.0247-0.0596-0.0063],D=10-15×[0.12050.2699-0.10660.0143].Compare the measured value and the predicted value generated by the model in Figures17 and 18.Figure 17
Comparison of calcination temperature measured curve and model predictive curve.Figure 18
Comparison of kiln tail temperature measured curve and model predictive curve.


## 6. 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.



---

*Source: 816831-2014-10-23.xml*