Moving Horizon Optimal Estimation for ... - ACS Publications

Subscriber access provided by READING UNIV

Process Systems Engineering

Moving Horizon Optimal Estimation for Temperature Distribution of FCCU Riser Reactor mengying LI, Yi Zheng, and Shaoyuan Li Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b01541 • Publication Date (Web): 16 Jul 2018 Downloaded from http://pubs.acs.org on July 21, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Moving Horizon Optimal Estimation for Temperature Distribution of FCCU Riser Reactor Mengying LI, Yi ZHENG, and Shaoyuan LI∗ Department of Automation, Shanghai Jiao Tong University, and Key Laboratory of System Control and Information Processing, Ministry of Education, Shanghai 200240, China E-mail: [email protected]

Abstract The temperature distribution of riser reactor plays an important role in fluid catalytic cracking unit (FCCU) since it directly affects the quantity and quality of the product. However, there are only limited thermocouples installed on the FCCU riser. In order to monitor the axial temperature distribution of the riser, a neighborhood optimization based distributed moving horizon estimation (N-DMHE) is proposed. The riser reactor is divided into several subsystems according to the location of thermocouples, and the temperature distribution of every subsystem is estimated by an individual small-scale moving horizon estimation (MHE). These MHEs are solved serially in distributed framework. Consequently, the computational time is reduced to satisfy the requirement of online implementation. The performance index for each subsystem includes cost functions not only by itself but also by its neighbors, which is called neighborhood optimization. The neighborhood optimization is a cooperative strategy so that the global performance can be efficiently improved. The proposed method is applied to a FCCU in JiuJiang Petrochemical Co. Ltd. The experiment results show the effectiveness of the proposed method.

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1 Introduction Fluid catalytic cracking unit (FCCU) is an industrial process that converts high molecularweight hydrocarbons to lower molecular-weight products. 1 The violent catalytic cracking reaction is completed in short-contact-time riser reactors where the FCCU catalyst is pneumatically conveyed by the hydrocarbon vapor from the bottom to the top of a vertical lifeline. 1,2 The temperature distribution along riser reactor, which varies from an inlet temperature of roughly 860–880K down to an outlet temperature of roughly 780–800K, is crucial for the cracking rates and the quality of the products. However, there are only a limited number of thermocouples installed on the FCCU riser and it is impractical to place a great many thermocouples along the reactor. In order to real-time monitor the status and characteristics of the entire reaction process, it is necessary to estimate the axial temperature distribution of the riser. Moreover, the estimated temperature distribution can be used to predict the product yield more accurate and be used in the advanced temperature control. Therefore, to develop an algorithm which is able to precisely and rapidly estimate the temperature distribution is important. To design an algorithm for estimating the dynamic temperature distribution using limited measurements is a typical estimator design problem which has received significant attention in large-scale complex chemical processes. 3 Such techniques usually include the Extended Kalman filter (EKF), Extend Luenberger estimator, high gain estimator, moving horizon estimation(MHE), 4 etc. Among them, the MHE, which reformulates the estimation problem as a quadratic problem using a moving fixed-size window, at each time instant, has been one of the most popular method for state estimation and has been successfully applied to various linear, 5,6 nonlinear 7,8 and hybrid 9,10 systems in industrial processes. It not only inherits the good optimization performance of the classical Kalman filter, but also has the merit of explicitly handling noise and state constraints. Since FCCU riser temperature distribution estimation problem has the characteristics of nonlinearity, high dimension and many constraints on noise and state, 11 MHE is chosen in this work to 2


Page 2 of 37

Page 3 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


monitor the temperature distribution. Considering that the online implementation of centralized estimation is impractical for the large-scale system due to its excessive online computational demand, the Distributed realization of MHE, Distributed MHE (DMHE), is appeared in literatures for improving computational speed, flexibility and error torllerance. A DMHE algorithm has been developed for discrete-time linear constrained systems and the convergence properties of DMHE are proven. 12 In ref 13 , DMHE is applied to a set of decoupled sensors where when processing and transmitting large amounts of data, the excellent computational capabilities of DMHE can coordinate sensors activity well. Ref 14 further discusses that DMHE algorithm can reduce the computational load of the overall solution. In ref 15,16, the performance of DMHE under bounded system disturbances and measurement noise are discussed. Using a distributed algorithm can reduce the computational load of the overall solution. However, it cannot avoid that the distributed realization of MHE cannot get an optimization performance as good as that of a centralized one. To further improve the estimation performance of the entire system is still a problem that has been considered in DMHE field. This stimulates developing a DMHE and its coordinate strategy for real-time estimating FCCU riser temperature with an improved performance and a rapid computational speed. Movtivated by above considerations, a neighborhood optimization strategy for DMHE is proposed in this work, which makes a trade-off between the computational speed and the estimation performance. To increase the computational speed, the riser reactor is divided into several subsystems and the temperature distribution of each subsystem is estimated by an individual MHE. To improve the global performance of whole estimator, each subsystem-based MHE considers not only the performance of the local subsystem, but also those of its neighbors. 17 These subsystem-based MHEs are solved serially. Finally, an experiment that applied the algorithm to a FCCU in JiuJiang Petrochemical Co. Ltd., Jiangxi, China, is conducted.

3



The contents are organized as follows. Section 2 presents a five-lumped model of FCCU riser reactor and existing difficulties. A N-DMHE strategy is proposed in Sections 3, which includes the subsystem modeling, the formulation of N-DMHE and the potential expansions of the algorithm. Section 4 discusses the simulated case studies. The experimental results are presented in Section 5 followed by a brief conclusion in Section 6.

2 System Modeling and Existing Difficulties 2.1

Description of the System

Riser Reactor Regenerator

Zone2

Catalyst

Zone1 Riser Pipe

Feed oil

Figure 1: Schematic diagram of FCCU. A FCCU riser reactor is illustrated in Fig.1, the riser pipe consists of two reaction zones which are made of different heat resistant layers. Two zones are serially connected. The heavy feed oil that mixed with hot regenerated catalyst enter the riser with an initial temperature of 860–880K, then reacts endothermic with the temperature down to an outlet temperature of 780–800K. In the riser, the first zone is characterized by high temperature and short contact time. The second zone is featured with larger radius, lower reaction temperature and longer retention time. As a result, the main products are obtained and the coke is deposited on the surface of catalyst which would enter the regenerator to 4


Page 4 of 37

Page 5 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


burn the carbon then return to the riser again. The reaction process is continuous and the catalyst-to-oil ratio is the main input of the process. A FCCU in JiuJiang Petrochemical Co. Ltd is studied in this paper. In the FCCU riser reactor, the first and second reaction zone are 10.2 m and 8m in height, 1m and 2.8m in diameter. There are three thermocouples on the first zone in the position of 2.35m, 3.85m, 10.2m. And there are two thermocouples on the second reaction zone, in the position of 3.8m, 8m. The temperature distribution along the riser is important since it directly affects the cracking rates and the quality of the products. If the reaction temperature increases by 10K, the reaction rate increases by 10-20%, and the corresponding products also change significantly. Therefore, the precise estimation of temperature distribution is critical and meaningful. For the temperature estimation of the riser reactor, a mathematical model is developed and the assessment of the N-DMHE has been performed.

2.2

Dynamic Model

FCCU process is featured by complicated hydrodynamics and complex kinetics. 18,19 The lumped kinetics model is often adopted for the modelling of FCCU process in literatures. Among them, the three-lump kinetics model described in ref 20 is the fundamental model of FCCU riser reactor. This model is developed and extended to four, five, six, ten lump models by other researchers. 21 By increasing the number of lumps, reaction schemes of greater complexity are considered. 22,23 For the sake of ensuring accuracy and reducing computational complexity at the same time, this paper adopts a five-lump dynamical reaction model according to ref 24. Since the dynamic response of component concentrate is much faster than that of temperature, 25 the variation of components over time are neglected, and only the axial distribution of components is considered. On the contrary, the dynamic response of temperature is relatively slow, so both the variation over time and the axial distribution of temperature are considered. 5



Page 6 of 37

Taking the element of height of the riser for material balance calculation and heat balance calculation, the dynamic model of the first reaction zone is given by: dY1 = ST1 · b · Pra1 · φ(tc1 ) · Θ · K · Y1 dX1 1 ∂T ∂T 1 Λ1 + =− b · Pra1 · φ(tc1 ) · Θ · KT · Y1 ∂t ST1 1 + Γ1 ∂X1 1 + Γ1

(1.a) (1.b)

where Y1 = [y A1 , y D1 , y N1 , yG1 , C1] T Θ = diag[exp(−

EA E E ), exp(− D ), exp(− N ), 0, 0] RT RT RT

−(k AD +k AN +k AG +k AC ) 0 0 −(k DN +k DG +k DC ) 0 k AD  k AN k DN −k NG  k AG k DG k NG k AC k DC 0

 K=

 −∆H KT = 

AR ·( k AD + k AN + k AG + k AC )

−∆HDR ·(k DN +k DG +k DC ) −∆HNR ·k NG 0 0

0 0 0 0 0



0 0 0 0 0

T 

where y A1 , y D1 , y N1 , yG1 represent the unconverted rate, diesel yield, gasoline yield and gas yield based on total feed oil, in the fisrt reaction zone of the riser. C1 represents the coke content on the catalyst, quality rate. T represents the reaction temperature. X1 = x1 /L1 , where x1 is the axial position, L1 is the length of first riser reactor, thus X1 stands for the dimensionless position; ST1 = (Ωra1 · ρoil · L1 )/Foil , where Ωra1 is the cross-sectional area of the first reactor, ρoil is the density of feed oil and Foil is the total flow rate; b is catalyst-to-oil ratio; Pra1 is average pressure of the first riser reactor; φ(tc1 ) is inactivation rate of catalyst; Λ1 , Γ1 are heat capacity correction factors of the first reaction zone. E A , ED , EN are activation energy; k AD , k AN , k AG , k AC , k DN , k DG , k DC , k NG are kinetics front factors; Θ represents the activation energy matrix. K and KT are kinetic constant

6


Page 7 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


matrices. At the entrance of riser reactor, the unconverted rate is assumed as 1. The diesel yield, gasoline yield, gas yield and catalyst carbon content are assumed as 0. The inlet temperature Tin is detected by the thermocouple at the entrance. The boundary condition of the first zone is given as: Y1 | X1=0 = [1, 0, 0, 0, 0] T T | X1=0 = Tin

The dynamic model for the second reaction zone can be derived by similarly by changing the script “1” in Equation (1.a) (1.b) to “2”. Set Y2 = [y A2 , y D2 , y N2 , yG2 , C2] T be the product concentrate of the second reactor ; L2 and Ωra2 are the length and the cross-sectional area. The pressure of the second reaction zone is Pra2 ; φ(tc2 ) is inactivation rate of catalyst; Heat capacity correction factors are Λ2 , Γ2 . The boundary conditions are given by: Y2 | X2=0 = [y A1 | X1=1 , y D1 | X1=1 , y N1 | X1=1 , yG1 | X1=1 , C1| X1=1 ] T T | X2=0 = T | X1=1

2.3

Existing Difficulties

The objective of this work is to real-time estimate the FCCU riser temperature distribution according to the limited measurements of thermocouples installed on the riser. However, FCCU is a complicated system with distributed parameters. It is not easy to estimate the temperature by using a nonlinear partial differential equation (PDE) model (1). Finite element method is usually used as a reliable approach to transfer the PDE model to lumping parameter model, which leads to a high dimension of the transferred state 7



space model. It makes the state estimation very time consuming. Moreover, the transferred model is usually in the form of differential-algebraic equations (DAE), which adds the difficulty of designing estimator. In addition, the constraints on noise and state also make the estimation problem more difficult. To address the above problems, a distributed MHE with neighborhood optimization is deduced and the lumping subsystem model would be first introduced in the next section.

3

Moving Horizon Optimal Estimation Design

This section presents the novel N-DMHE algorithms applied for the temperature estimation of FCCU riser reactor. Since the online implementation of this large-scale nonlinear system’s state estimator is time consuming, a distributed realization of MHE is proposed to satisfy the requirements on computing speed. The riser reactor is divided into several sections (each section represents a subsystem) according to the position of thermocouples. As shown in Fig.2, each subsystem contains the part between two thermocouples and includes the upper thermocouple whose measurement is taken as the output of the subsystem. Particularly, each subsystem has only one thermocouple. For each subsystem, the lumping subsystem model is first introduced for the building of the estimators. Then an individual MHE estimator is built for the temperature distribution of each subsystem. The DAE model in each MHE is solved by using the estimated state sequence calculated in previous time instant, then these MHEs are solved serially and a distributed MHE (DMHE) estimator is formed. In addition, a neighborhood optimization strategy that allows the estimators to use more measurements, is proposed to improve the performance of the entire estimator. The specification is organised as follows.

8


Page 8 of 37

Page 9 of 37

Position Nns

S=ns 2.8m

8m

2 1 Nms

1m

2 1 N2

10.2m

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


... S=ms ... S=2

2 1

N1

S=1

2 1

Feed

Temperature

Figure 2: Distributed state estimator scheme.

3.1

Lumping Subsystem Model

To transfer the PDE model (1) to a state-space model (lumping parameter model), the finite difference method 26,27 is employed. According to this method, each subsystem could be divided into several smaller segments in axial direction. If the segments are small enough, each small segment can be approximated as a temperature point. Therefore, each subsystem consists of a series of temperature points that can be augmented as state variables, thus a state estimator could be built to estimate the temperature distribution of each subsystem. Suppose there are ms thermocouples installed on the first reaction zone and (ns − ms ) thermocouples on the second zone. Then the first reaction zone is divided into ms subsystems and the second reaction zone is divided into (ns − ms ) subsystems. The whole system has ns subsystems. As shown in Fig.2, the sth (s = 1, · · · , ns ) subsystem is divided into Ns segments and the Ns th segment represent the thermocouple. Further considering that the measurements are usually provided digitally with a sampling time ∆t in refinery factory, the discrete-time version of the subsystem is derived. Hence, the dynamic model of the sth subsystem in the first reaction zone is expressed

9



Page 10 of 37

by: ST1 s b · Pra1 · φ(tc1 ) · K · Θ · Y1,i −1 Ns Ns ∆t Ns ∆t s Tis (k + 1) = (1 − ) · Tis (k) + · T s (k) + b · f (Y1,i −1 ) ST1 (1 + Γ1 ) ST1 (1 + Γ1 ) i−1

s s Y1,i − Y1,i −1 =

(2.a) (2.b)

where A1 D1 N1 G1 s = [ys,i , ys,i , ys,i , ys,i , C1s,i ] T Y1,i

EA E E ), exp(− Ds ), exp(− Ns ), 0, 0] s RTi−1 RTi−1 RTi−1

Θ = diag[exp(− f (ζ ) = −

Λ1 ∆t · Pra1 · φ(tc1 ) · KT · Θ · ζ 1 + Γ1

j

where s = 1, 2, · · · , ms . i = 1, 2, · · · , Ns is the segments of the sth subsystem. ys,i ( j = A1, D1, N1, G1) represents the product concentration of the i segments in sth subsystem. The dynamic equation for the second reaction zone can be derived by changing the script “1” in Equation (2.a) (2.b) to “2”. To build the estimator of each subsystem, define algebratic equation (2.a) as g, define each segment’s temperature Tis (i = 1, 2, · · · , Ns ) as state variables, define the product yield s (or Y s ) (i = 1, 2, · · · , N ) as algebraic variables, let catalyst-to-oil and the coke content Y1,i s 2,i

ratio b and inlet temperature Tin as the manipulated variable u, let the measurements as the output vector ys . Model (2) can be further expressed by:

xs (k + 1) = As · xs (k ) + As,s−1 · xs−1 (k ) + Bs (k, ws (k )) · u(k ) + µs (k ) 0 = g( xs (k), ws (k ))

(3.a) (3.b)

ys (k ) = Cs · xs (k ) + ηs (k )

(3.c)

where

10


Page 11 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


s T xs = [ T1s , T2s , · · · , TN ] , s = 1, · · · , ns ; s s s ws = [Y1,1 , · · · , Y1,N ]T , s = 1, · · · , ms , s s s (ws = [Y2,1 , · · · , Y2,N ]T , s = ms + 1, · · · , ns ); s

Tin ] T .

u = [b

Here, measurement errors in the measured outputs at time instant k are expressed by ηs (k ). A simple additive state noise term with µs (k ) accounts for unknown disturbances on the system states. As , Bs (k, ws (k)), Cs are the system matrices. As,s−1 is the coupling matrices between

(s − 1)th and sth subsystem. These matrices are given by: 

1− NSs ∆t

  As =  

T1 Ns ∆t ST1

0 ···

0

0

1− NSs ∆t 0 ···

0

0

0

.. .

.. .

0

0

T1

.. . . . . 0 ···

Ns ∆t ST1



  ..  . 

1− NSs ∆t T1

Ns × Ns

A1,0 = 0 

As,s−1

0  0  = .  ..   0

0 ··· 0 0 ··· .. . . . .

0 .. .



Ns ∆t ST1 

0 .. .

      

0 ··· 0 0 Ns × Ns    f (ws−1,Ns (k)) β(k )    f (ws,1 (k )) 0    Bs (k, ws (k )) =   ..   .     f (ws,Ns −1 (k)) 0

11


(4)

Ns ×2


Page 12 of 37

Ns ∆t β(k) = 0, β(k ) = (s = 1, ms + 1) ST1 Cs = 0 0 · · · 0 1 where s = 1, 2, · · · , ns and when ms < s ≤ ns , the matrices can be derived by replacing the script “1” to “2”. It should be noted that Bs (k, ws (k)) in equation (3.a) is related to ws (k) in algebraic equation (3.b) , so model (3) is implicit. Models that compose of both differential and algebraic equations are called DAE models. A simplified solution to the above implicit DAE model is to solve ws (k ) by algebraic equation (3.b) first since the estimated state xs (k ) is known at each time instant, then substitute the solved ws (k) into differential equation (3.a). In the next section, the formulation of N-DMHE will be described based on this lumping parameter model.

3.2

N-DMHE Formulation

The proposed N-DMHE method, in view of its fundamental principle, is based on the local MHE problem. The local MHE minimization problem is first introduced in this section. Then, the formulation of the N-DMHE will be deduced. MHE is a receding optimization method with a fixed information window. For a given moving horizon size M, at time instant k = M, M + 1, · · · , the local MHE optimization problem is given by: k −1

min

x s ( k − M ),{ µ s }

Js (k) =

min

∑

x s ( k − M ),{ µ s } t = k − M

kys (t) − Cs xs (t)k2Rs +kµs (t)k2Qs +θs (k)

(5)

where Rs and Qs are weighting matrices. The MHE only considers the information in the fixed horizon M, so it needs to determine the arrival cost θs (k). A common way to calculate the arrival cost is to apply the extended kalman filter. 28 12


Page 13 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In distributed framework, each subsystem optimize the local optimization function Js (k ) subject to the local dynamic model. The DMHE algorithm can greatly reduce the computational load. However, the optimal solution of the local optimization problem collectively is not equal to the whole system’s optimal solution. Thus the estimation performance of DMHE algorithm needs to be improved. Therefore, in this paper, a neighborhood optimization based DMHE (N-DMHE) is proposed to improve the global estimation performance and reduce the computational load as well. In the estimation of temperature distribution, since the subsystem’s output contains useful information of its upstream neighbor ( subsystem near the entrance), this information can be used to improve the estimation of the upstream neighbors’ states. To this end, define a set Ps = {i1 , · · · , is0 } ⊂ {1, · · · , ns } as the neighborhood of subsystem s. In the neighborhood optimization, neighborhood is defined by the strength of the coupling relationship between subsystems. In the riser, the adjacent subsystems have a strong coupling relationship, so subsystems with adjacent positions are generally selected as neighboring subsystems. The elements in Ps are the selected neighbor of sth subsystem and s ∈ Ps . The new performance index for each subsystem can be improved by

J s (k) =

∑

Ji (k )

(6)

i ∈ Ps

Define the new vectors x s = [ xiT1 , · · · , xiT0 ] T , ws = [wiT1 , · · · , wiT0 ] T , µs = [µiT1 , · · · , µiT0 ] T , ys = s

s

[yiT1 , · · · , yiT0 ]T , η s = [ηiT1 , · · · , ηiT0 ]T , i1 , · · · , is0 ∈ Ps . Define the new system matrices s

s

13


s




Page 14 of 37



0 ··· 0   A i1    Ai ,i Ai  · · · 0  2 1  2 As =  . , ..  .. ..  ..  . . .     0 · · · Ais0 ,is0 −1 Ais0 Bs = Block − diag( Bi1 , · · · , Bis0 ), G s = Block − diag( Gi1 , · · · , Gis0 ), C s = Block − diag(Ci1 , · · · , Cis0 ) The local dynamic model of each subsystem’s neighborhood can be expressed as the following collective form according to equation (3).

x s (k + 1) = As · x s (k) + Bs (k, ws (k)) · u(k) + G s · µs (k ) 0 = g( x s (k ), ws (k ))

(7.a) (7.b)

ys (k ) =C s · x s (k ) + η s (k )

(7.c)

The state prediction of each subsystem’s neighborhood, in the moving horizon l = 1, 2, · · · , M, is able to be calculated by: l

k − M + l −1

x s (k − M + l |k ) = As · x s (k − M |k ) +

∑

As

k − M − j + l −1

Bs ( j, ws ( j))u( j) + G s µs ( j)

j=k− M

(8) It should be noticed that solving Bs ( j, ws ( j)) in equation (8) must solve the ws ( j) first, then plug ws ( j) into equation (4) to get Bs ( j, ws ( j)). There are two ways to calculate ws ( j): • Use the estimated states xˆ s ( j)( j = k − M, · · · , k − 1). Since these states have already been estimated in previous time instant, ws ( j) can be solved quickly through equation

14


Page 15 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(7.b). The calculation equation is given by ws ( j) = arg{ g( xˆ s ( j), ws ( j)) = 0} j = k − M, · · · , k − 1

(9)

• Use the iterative values x s ( j|k) in predictive model (8). ws ( j) can be solved step by step as equation (10) shows.

ws ( j) = arg{ g( x s ( j|k ), ws ( j|k )) = 0} j = k − M, · · · , k − 1

(10)

Since the second approach is iterative, the computation amount is greatly increased. From the perspective of computation, the first approach is adopted in this paper to reduce the complexity and amount of computation. Define the collective weighting matrices Rs = Block − diag( Ri1 , · · · , Ris0 ), Qs = Block − diag( Qi1 , · · · , Qis0 ). The N-DMHE optimization problem is able to be expressed as Φ=

min

−1 x s (k − M ),{µs (t)}kt= k− M

k −1

=

min

(

∑

x s ( k − M ),{ µ s } t = k − M

J s (k)

kys (t) − C s x s (t)k2R +kµs (t)k2Q ) + θ s (k) s

s

(11)

s.t. Equation(8), (9); x L ≤ x s (t) ≤ xU , t = k − M, · · · , k − 1; µ L ≤ µs (t) ≤ µU , t = k − M, · · · , k − 1. The term weighted by symmetric positive definite matrices Rs and Qs are measurement errors and state disturbance. It not only contains the information of the subsystem but also 15



Page 16 of 37

its neighbor. The states should satisfy the constraints of a lowest temperature x L and a highest temperature x H . The disturbance should satisfy the constraints of µ L and µ H . The term θ s (k ) in (11) is the collective arrival cost that is given by θ s (k ) = k x s (k − M) − xs∗ (k − M )k2Πs (k− M)

(12)

where xs∗ (k − M ) denotes the optimal estimate at time k − M given all of the measurements ys (k ) from time 0 to k − M − 1 and Πs (k − M ) is computed from the Kalman filter covariance update as follows T

T

Πs (k ) = − As Πs (k − 1)C s ( Rs + C s Πs (k − 1)C s )−1 C s Πs (k − 1) As T

+ As Πs (k − 1) As + G s Qs G s

T

T

(13)

In the following we sketch the steps that have to be carried out in practice, in order to apply the proposed N-DMHE algorithm. Algorithm 1 N-DMHE algorithm: 1. Initialization: Each subsystem initials the state xs (0), ws (0), subsystem matrices As , As,s−1 , Bs (0, ws (0)), Cs , and the collective states and matrices x s (0), ws (0), As , As,s−1 , Bs (0, ws (0)), C s . Let s = 1 in Step 2,3,4; 2. Calculation: At time instant k, based on the known estimated state x s (k − M), · · · , x s (k − 1), the subsystem calculates ws (k − M ), · · · ws (k − 1) through equation(7.b), then plug it into equation(4) to get Bs (k − M, ws (k − M)), · · · , Bs (k − 1, ws (k − 1)); 3. Optimization: Each subsystem solves the minimization problem (11) based on local measurements ys (k − M), · · · , ys (k − 1) and neighborhood measurements to get the optimal xˆ s (t − M), {µˆ s }; The state sequence xˆ s (k − M|k ), · · · , xˆ s (k |k) stemming from xˆ s (k − M) and {µˆ s } could be derived from the predictive model (8). The estimation xˆ s (k − M + 16


Page 17 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1|k), · · · , xˆ s (k|k ) could be derived correspondingly. Furthermore, each subsystem updates the current estimation xˆ s (k ) = xˆ s (k |k) and xˆ s (k ) = xˆ s (k |k); 4. Communication: In distributed framework, each subsystem transmits the state sequence xˆ s (k − M + 1|k ), · · · , xˆ s (k |k) to the downstream subsystem s+1 through the coupling matrices As+1,s ; Let s = s + 1, go back to Step 2-4 until s=ns ; If s=ns , go to Step 5; 5. Receding horizon: Moving horizon to the next sampling time, that is k → k + 1, go to Step 1, and repeat the above steps. Through the above steps, the on-line estimation of FCCU riser reactor is realized by several small-scaled cooperative MHEs. The computational load is greatly reduced and the neighborhood strategy in a distributed structure can improve the entire performance. The major contribution of this paper is designing a N-DMHE algorithm and applying it to FCCU riser temperature estimation, the theoretical proof of the convergence of N-DMHE is out of the range of this paper.

3.3

Some Discussions

The estimator design problem for FCCU process should be able to address challenges including the large scale, strong nonlinearity and hard constraints, as well as implicit DAE. By applying MHE, which outperforms the EKF in the presence of constraints, N-DMHE is able to deal with these complications with high performance. For the distributed form, the neighborhood optimization strategy enables the usage of the neighboring measurements and performance to coordinate subsystems. Real-time estimation of the FCCU riser temperature distribution can further provide the following advantages. • The estimated temperature distribution is real-time corrected by the measurements, thus the estimated temperature is more reliable and believable than that obtained by method without feedback. 17



• The estimated temperature distribution could be used to predict the product yeild such as diesel yield and gasoline yield, which helps obtaining products with better quality. • The estimated state can be further used in a model based control algorithm, e.g., model predictive control (MPC) to control the process. 29,30 To validate the proposed strategy, both the numerical simulation and experimental results are presented in section 4 and 5.

4 Case Studies A joint simulation using both Matlab and gProms is conducted to verify the performance of the proposed N-DMHE in estimating the temperature distribution of FCCU riser reactor, where gProms is used to simulate FCCU plant in JiuJiang Petrochemical Co. Ltd and matlab is used to run the N-DMHE algorithm. The matlab program accesses the data generated by gProms at each sampling time. Table 1: Industrial Operating Data of the Riser Process variable ρoil Foil Pra Λ1 Λ2 Γ1 Γ2

plant data 6.8 168 360 0.09 0.03 136 24

Unit kg/m3 t/h kPa (kgK )/kJ (kgK )/kJ -

There are three thermocouples installed at 2.35m, 3.85m, 10.2m in the first reaction zone and two thermocouples installed at 3.8m, 8m in the second zone. When performing the finite difference processing on the riser, the difference equations should meet accuracy requirements. So we set the segments of the first and second zone as 40 and 30. According to the position of thermocouples, the segments of five subsystems are set as 18


Page 18 of 37

Page 19 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


N1 = 10, N2 = 6, N3 = 24, N4 = 15, N5 = 15. The state of first reaction zone subjects to a lowest temperature of 770K and a highest temperature of 870K. The state of second zone subjects to a lowest temperature of 760K and a highest temperature of 810K. The industrial operating data of the riser are shown in Table 1. The kinetic parameters of the FCCU reaction are shown in Table 2. The kinetic constants and activation energy are identified under steady-state operation by a modified LM method based on the industrial data in two month, which is not within the scope of this paper. Table 2: Kinetic Constants and Activation Energy Item kg/(kg · Pa · s) Item kJ/mol k AD 1074926.3 k AN 471168.5 EA 71211.9 k AG 1574081.7 k AC 195180.5 ED 72710.9 k DN 215453.9 k DG 246014.3 EN 125346.5 k DC 102892.6 k NG 106786.8 Let sth subsystem’s neighborhood Ps = {s, s + 1}. Let the sampling period ∆t = 2.2s and the optimization horizon length M = 150. The weighting matrices of the objective function Rs = diag(0.5, 0.5..., 0.5), Qs = diag(1, 1..., 1), and Πs (0) = diag(1, 1..., 1). Since the disturbances mainly come from the inlet temperature Tin and the catalyst-tooil ratio b, step signals of Tin and b are respectively sent into the system to validate the performance of the proposed N-DMHE estimator. In comparison, the simulation using DMHE is conducted too. Three cases of simulations are carried out: Case 1: Tin varying, no measurement error; Case 2: b varying, no measurement error; Case 3: b varying and the measurement disturbance existing in outputs. In these cases, different initial state vectors are given to process models and estimators to validate the convergence of the proposed method. Meanwhile, N-DMHE and DMHE use 19



870 860

Temperature(K)

Temperature(K)

850

850 840 830 820

840 830 820 810

0

100

200

300

400

0

100

200

(a)

400

300

400

805

Temperature(K)

Temperature(K)

300

(b)

830 820 810 800 790

800 795 790 785 780

0

100

200

300

0

400

100

200

(d)

(c) 795

Temperature(K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 37

N-DMHE DMHE Model

790

785

780 0

100

200 (e)

300

400

Sampling instant : k

Figure 3: Comparison of temperatures estimated by gProms model, N-DMHE and DMHE in Case 1. the same initial state vectors for comparison. Two estimators are serially calculated. The simulation results are shown in Fig.3–5, where subfigures ( a)–(e) in these figures represent the 1th –5th subsystems respectively. In the subfigures ( a)–(e), each curve describes the state of a temperature point. These curves together form the temperature distribution of the whole system. Fig.3 shows the estimated states using N-DMHE and DMHE in the case of a step changing in the inlet temperature Tin . The estimated states of five subsystems are con20


Page 21 of 37

870 860

Temperature(K)

Temperature(K)

850

850 840 830 820

840 830 820 810

0

100

200

300

400

0

100

(a)

300

400

300

400

805

Temperature(K)

Temperature(K)

200

(b)

830 820 810 800 790

800 795 790 785 780

0

100

200

300

400

0

100

(c)

200

(d)

795

N-DMHE DMHE

Temperature(K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


790

Model

785 780 775 0

100

200 (e)

300

400


Figure 4: Comparison of temperatures estimated by gProms model, N-DMHE and DMHE in Case 2. vergent to the states of plant model. The convergence time of five subsystems are 50s, 78s, 160s, 220s and 280s, respectively, which are gradually increasing since subsystems are serially connected. It is worth noting that each estimated state using N-DMHE goes to a small neighborhood of system state very fast. The curves of states estimated by N-DMHE have the same trend as curves calculated by system model, while the curves resulting from DMHE do not have these trend before convergence. The mean square error

21



870 860

Temperature(K)

Temperature(K)

850

850 840 830 820

840 830 820 810

0

100

200

300

400

0

100

(a)

300

400

300

400

805

820

Temperature(K)

Temperature(K)

200

(b)

830

810 800

800 795 790 785

790 0

100

200

300

780

400

0

100

(c)

200

(d)

795 N-DMHE

Temperature(K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 37

DMHE

790

Model

785

780

0

100

200 (e)

300

400


Figure 5: Comparison of temperatures estimated by gProms model, N-DMHE and DMHE in Case 3. (MSE) between the state estimated by N-DMHE and the state calculated by system model is significantly smaller than that of DMHE. As shown in table 3, in Case 1, the MSE of N-DMHE is 2.4, whereas that of DMHE is 7.1. This result illustrates that N-DMHE has a better estimation performance than DMHE. Furthermore, the N-DMHE algorithm takes 1.2s for each optimization step on average and DMHE algorithm takes 1.0s, while the centralized MHE algorithm takes 8.1s. It demonstrates that N-DMHE could obviously increase the computational speed. 22


Page 23 of 37

840 830

Temperature(K)

Temperature(K)

820

820 810 800 790

810 800 790 780

0

50

100

150

200

250

0

300

50

100

810

200

250

300

200

250

300

790

800

Temperature(K)

Temperature(K)

150

(b)

(a)

790 780

780

770

770 0

50

100

150

200

250

760

300

0

(c)

50

100

150

(d) N-DMHE DMHE

Temperature(K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


770

Model

765

760

0

50

100

150 (e)

200

250

300


Figure 6: Temperatures estimated by gProms model, N-DMHE and DMHE with constraints. Fig.4 and 5 show the estimated states using N-DMHE and DMHE in the cases that a step changing occurs in the catalyst-to-oil ratio b. Similarly, the estimated states of all subsystems converge to the states of plant model and the convergence time is increasing one by one. The MSE between N-DMHE and the state calculated by system model are smaller than that of DMHE, as shown in table 3. In case 3, a Gaussian random noise is added to the thermocouple signals with zero mean to investigate the designed estimator’s characteristic of against noisy. The simulation result shows that the proposed estimator is of good anti-noise performance. 23



840 830

Temperature(K)

Temperature(K)

820

820 810 800 790

810 800 790 780

0

50

100

150

200

250

0

300

50

100

(a)

200

250

300

200

250

300

790

800

Temperature(K)

Temperature(K)

150

(b)

810

790 780

780

770

770 760 0

50

100

150

200

250

300

0

50

100

(c)

150

(d) N-DMHE DMHE

Temperature(K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 37

770

Model

765

760

0

50

100

150 (e)

200

250

300


Figure 7: Temperatures estimated by gProms model, N-DMHE and DMHE without constraints. Table 3: Mean Square Error of N-DMHE and DMHE Algorithm Case 1 Case 2 Case 3 N–DMHE 2.4 2.5 2.7 DMHE 7.1 6.4 6.8 In addition, we have tests where N-DMHE and DMHE are stochastically given the same initial value and we calculate the MSE of two estimators under 100 different initial values. The statistical results are shown in the table 4. Table 4 shows the ratio of the MSE of N-DMHE and the MSE of DMHE. The simulation results show that N-DMHE has smaller

24


Page 25 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


MSE than DMHE and performs better. Table 4: The ratio between MSEs of two estimators with same stochastic initial value Ratio of MSE Maximum Minimal Average N − DMHE 0.62 0.33 0.49 DMHE Furthermore, we have a test to demonstrate MHE’s capability to handle constraints. Fig.6 and 7 show the comparison between estimated states with and without constraints. As shown in Fig.6, the constraints work. the states of 3th subsystem are limited by the lowest temperature 770K and the states of 5th subsystem are limited by the lowest temperature 760K. Fig.7 shows the corresponding states without constraints. The MSE between the constrained estimation of N-DMHE and the state calculated by system model is 6.4, which is smaller than the 6.8 without constraints. Similarly, the MSE between the constrained estimation of DMHE and the state calculated by system model is smaller than that without constraints. The simulation result shows that the MSE is smaller with constrains than without. From these simulation results, it can be seen that the estimations of both estimators converge to the results of plant model even in the presence of disturbances. N-DMHE and DMHE able to incorporate constraints in the formulation. The MSE of N-DMHE is smaller compared with DMHE. To summarize, the proposed N-DMHE has better performance than DMHE in estimating the temperature distribution of the FCCU riser reactor.

5 Experimental Results To further validate the proposed method, the proposed method is implemented to the FCCU in JiuJiang Petrochemical Co. Ltd. As shown in Fig.8, the estimator is running on the server which is equipped with 2.66GHz 8 cores CPU, 16G memory and high speed network. The estimator accesses the real-time database of the FCCU plant in Jiujiang petrochemical company every sampling 25



FCCU

Client Terminal IP21 Database

IIS&Server

Figure 8: The structure of platform used in the experiment. time to obtain the real-time industrial plant measurements. The sampling period is 2.2s. The measurements include inlet temperature, catalyst circulation, feed oil flowrate and the correspondingly measurements of five thermocouples. The catalyst circulation and feed oil flowrate are used to calculate the catalyst-to-oil ratio which is one of the input of system model. The N-DMHE estimator estimates the temperature distribution of FCCU reactor at each time instant and real-time updates the results stored in server. The results are able to be displayed and monitored by any client terminals through the IIS service in the server.

5.1

Model Validation

To test the effectiveness of the designed subsystem state-space model (3), the process data from Jiujiang Petrochemical Com. Ltd. are utilized. The results obtained by process model present a good fittness to industrial measurement data, as shown in Fig.9. The five subfigures (a)–(e) in Fig.9 show the comparison of the outputs of model and the measurements of the five thermocouple-points. The solid lines represent the measurements of five thermocouples and blue dotted lines represent the outputs of system model. Except for some noisy points, the error between the output of process model and the measured value is less than 0.8K, which is acceptable. The mean square errors between the outputs of process model and the measurements at five thermocouple-points are 0.28, 0.31, 0.31, 0.36 and 0.34, respectively. It should be noticed that the model outputs are smoother than

26


Page 26 of 37

Page 27 of 37

the measurements, which is reasonable, because of the unknown system disturbances and

817

796

826.5

816.5

795.5

826 825.5 825 824.5

Temperature(K)

827

Temperature(K)

Temperature(K)

measurement noises in practice.

816 815.5 815 814.5

824 100

200

300

400

100

200

(a)

300

400

(b)

786

782

Temperature(K)

782.5

785.5 785 784.5

794.5 794

793 0

786.5

795

793.5

814 0

Temperature(K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0

100

200

300

400

(c)

Model Measurements

781.5 781 780.5

784

780 0

100

200

300

400

0

100

200 (e)

(d)

300

400


Figure 9: Comparison of the outputs of model and the measurements.

5.2

Performance of N-DMHE

In the experiment, the input of the estimator are catalyst-to-oil ratio and inlet temperature. The catalyst-to-oil ratio is calculated by dividing the catalyst circulation by the feed oil flowrate. At each time instant, the N-DMHE estimator accesses catalyst circulation, feed oil flowrate, inlet temperature and the corresponding measurements of five thermocouples, and online estimates the temperature distribution based on these real-time data. The weighting matrices Rs = diag(5, 5..., 5), Qs = diag(1, 1..., 1), and Πs (0) = diag(1, 1..., 1). The estimation result timed from the start of the estimator and last for 600 sampling moments, is presented. Fig 10–11 show the input data of the estimator, and Fig 12–13 show 27



the output. Fig.10 shows the actual catalyst circulation, feed oil flowrate, and the calculated catalyst-to-oil ratio which fluctuates from 5.8–6.4. Fig.11 describes the fluctuation of the

catalyst(t/h)

inlet temperature, the range of which is between 858K and 864K. 1100

1000

900 0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

feed oil(t/h)

180 170 160

catalyst-to-oil ratio

150 7

6

5


Figure 10: The catalyst circulation, feed oil amount and input: catalyst-to-oil ratio.

865 864 863

Inlet Temperature(K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 37

862 861 860 859 858 857 856 855 0

100

200

300

400

500

600


Figure 11: The input: inlet temperature. The estimation result is shown in Fig.12. The red surface represents the estimated temperature distribution of riser reactor, and the blue surface represents the measurements. 28


Page 29 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Since there are limited thermocouples, the blue surface is formed by straight connecting the measured temperature. It can be seen that the estimated temperature distribution converges to the measurements. And the temperature curves along the riser are very smooth in both zone 1 and 2. At the junction of reaction zone 1 and 2 (segments 40), the slope of the surface is obviously changing because these two zones have different heat-resistant layers.

Figure 12: The temperature distribution estimated by N-DMHE and the measurements of the five thermocouples. For the temperature estimation at five thermocouple-points, the comparison between temperature estimated by N-DMHE and measurements is shown in Fig.13. The subfigures (a)–(e) show the temperature at five thermocouple-points, respectively. The estimated states are convergent to the measurements, and the convergence speed is fast. After convergence, the mean square errors between the estimations and the measurements are 0.24, 0.25, 0.27, 0.33 and 0.31, respectively. These results show that the proposed method gets good results. Additionally, the proposed N-DMHE can be further used to predict product yield which are difficult to detect in practice. Fig.14 shows the satisfactory estimates of the product concentration obtained by using the estimated temperature state. In Fig.14, the horizontal axis represents the axial direction of riser reactor and vertical axis represents the percent concentration. Solid and dotted lines represent the concentration of two reaction 29



835 830 825

800

Temperature(K)

820

Temperature(K)

Temperature(K)

840

815

810

795

790

820 815

805 0

200

400

600

785 0

200

(a)

400

600

0

200

(b)

800

400

600

(c)

790 N-DMHE

795

Temperature(K)

Temperature(K)

790 785

785

Measurements

780

780 775 775 0

200

400

600

0

200

400

600

(e)

(d)


Figure 13: The estimated temperature of N-DMHE and the measurements at five thermocouple-point. 1 yA1

0.9

yA2 yN1

0.8

Percent concentration(%)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 37

yN2 yD1

0.7

yD2

0.6

yG1 G2

y C1 C2

0.5 0.4 0.3 0.2 0.1 0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x

Figure 14: The estimated product concentration of two reaction zones and the outlet product concentration obtained from plant. zones respectively. The dots at the end represent the outlet product concentration obtained from plant. To illustrate the accuracy of the product prediction, the diesel is used as

30


Page 31 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


example. The average outlet diesel concentration obtained from plant is 21.2%, and the average outlet diesel concentration estimated by N-DMHE is 20.8%. The percentage error of 1.8% illustrates an accurate estimation. In a conclusion, the N-DMHE is an effective method for temperature estimation of FCCU process with improved performance and fast computational speed.

6 Conclusions A neighborhood optimization based distributed moving horizon estimation is proposed in this paper for real-time monitoring the axial temperature distribution of the riser reactor. By using this method, the whole system is divided into several subsystems and the temperature distribution is serially estimated by subsystem-based small-scale MHEs so that the computational load is obviously reduced. In addition, a neighborhood optimization strategy that allows each subsystem to consider the performance of itself and its neighbor together to improve the global estimation performance, is employed. In the numerical example, we investigate the performance of the N-DMHE and DMHE algorithm. The simulation results show both approach converge well and the proposed approach is better. Finally, an experiment which applies the proposed N-DMHE to a FCCU in JiuJiang Petrochemical Co. Ltd., is conducted. The experiment result shows the efficiency of the proposed method.

Acknowledgement The authors are grateful for the support from National Nature Science Foundation (NSFC) of China (61590924,61673273).

31



References (1) Theologos, K. N.; Lygeros, A. I.; Markatos, N. C. Feedstock atomization effects on FCC riser reactors selectivity. Chemical Engineering Science 1999, 54, 5617–5625. (2) Chen, G.-Q.; Luo, Z.-H. New insights into intraparticle transfer, particle kinetics, and gas–solid two-phase flow in polydisperse fluid catalytic cracking riser reactors under reaction conditions using multi-scale modeling. Chemical Engineering Science 2014, 109, 38–52. (3) Dochain, D.; Couenne, F.; Jallut, C. Enthalpy based modelling and design of asymptotic observers for chemical reactors. International Journal of Control 2009, 82, 1389–1403. (4) Mohd Ali, J.; Ha Hoang, N.; Hussain, M. A.; Dochain, D. Review and classification of recent observers applied in chemical process systems. Computers and Chemical Engineering 2015, 76, 27–41. (5) Rao, C. V.; Rawlings, J. B.; Lee, J. H. Constrained linear state estimation a moving horizon approach. Automatica 2001, 37, 1619 – 1628. (6) Alessandri, A.; Awawdeh, M. Moving-horizon estimation with guaranteed robustness for discrete-time linear systems and measurements subject to outliers. Automatica 2016, 67, 85 – 93. (7) Rao, C. V.; Rawlings, J. B.; Mayne, D. Q. Constrained state estimation for nonlinear discrete-time systems: stability and moving horizon approximations. IEEE Transactions on Automatic Control 2003, 48, 246–258. (8) MÃijller, M. A. Nonlinear moving horizon estimation in the presence of bounded disturbances. Automatica 2017, 79, 306 – 314. (9) Bemporad, A.; Mignone, D.; Morari, M. Moving horizon estimation for hybrid systems

32


Page 32 of 37

Page 33 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and fault detection. Proceedings of the 1999 American Control Conference (Cat. No. 99CH36251). 1999; pp 2471–2475 vol.4. (10) Ferrari-Trecate, G.; Mignone, D.; Morari, M. Moving horizon estimation for hybrid systems. IEEE Transactions on Automatic Control 2002, 47, 1663–1676. (11) Haseltine, E. L.; Rawlings, J. B. Critical Evaluation of Extended Kalman Filtering and Moving-Horizon Estimation. Industrial & Engineering Chemistry Research 2005, 44, 2451–2460. (12) Farina, M.; Ferrari-Trecate, G.; Scattolini, R. Distributed moving horizon estimation for sensor networks, contract number INFSO-ICT-223854. IFAC Proceedings Volumes 2009, 42, 126 – 131, 1st IFAC Workshop on Estimation and Control of Networked Systems. (13) Farina, M.; Ferrari-Trecate, G.; Scattolini, R. Distributed moving horizon estimation for linear constrained systems. IEEE Transactions on Automatic Control 2010, 55, 2462–2475. (14) Chen, T.; Zhou, D.; Tran, T.; Kastner, C.; Ling, K.; Tseng, K.; Maciejowski, J. M. Distributed moving horizon estimation for power systems. Power & Energy Society General Meeting, 2015 IEEE. 2015; pp 1–5. (15) Zhang, J.; Liu, J. Distributed moving horizon state estimation for nonlinear systems with bounded uncertainties. Journal of Process Control 2013, 23, 1281–1295. (16) Zeng, J.; Liu, J. Distributed moving horizon estimation subject to communication delays and losses. 2015 American Control Conference (ACC). 2015; pp 5533–5538. (17) Zhang, Y.; Li, S. Networked model predictive control based on neighbourhood optimization for serially connected large-scale processes. Journal of Process Control 2007, 17, 37 – 50.

33



(18) McFarlane, R. C.; Reineman, R. C.; Bartee, J. F.; Georgakis, C. Dynamic simulator for a model IV fluid catalytic cracking unit. Computers & chemical engineering 1993, 17, 275–300. (19) Sadeghbeigi, R. Fluid catalytic cracking handbook: An expert guide to the practical operation, design, and optimization of FCC units; Elsevier, 2012. (20) Weekman, V. W.; Nace, D. M. Kinetics of catalytic cracking selectivity in fixed, moving, and fluid bed reactors. AIChE Journal 1970, 16, 397–404. (21) Dasila, P. K.; Choudhury, I. R.; Singh, S.; Rajagopal, S.; Chopra, S. J.; Saraf, D. N. Simulation of an Industrial Fluid Catalytic Cracking Riser Reactor Using a Novel 10Lump Kinetic Model and Some Parametric Sensitivity Studies. Industrial & Engineering Chemistry Research 2014, 53, 19660–19670. (22) Corella, J.; Frances, E. Analysis of the riser reactor of a fluid cracking unit: model based on kinetics of cracking and deactivation from laboratory tests; ACS Publications, 1991. (23) Arbel, A.; Huang, Z.; Rinard, I. H.; Shinnar, R.; Sapre, A. V. Dynamic and control of fluidized catalytic crackers. 1. Modeling of the current generation of FCC’s. Industrial & engineering chemistry research 1995, 34, 1228–1243. (24) Luo, X.; Yuan, P.; Lin, S. Dynamic Model of Fluid Catalytic Cracking Unit I. Reactor section. ACTA PETROLEI SINICA PETROLEUM PROCESSING SECTION 1998, 14, 34–40. (25) Feng, X.; Yeye, W.; Xionglin, L. Soft sensor for inputs and parameters using nonlinear singular state observer in chemical processes. Chinese Journal of Chemical Engineering 2013, 21, 1038–1047. (26) Mitchell, A. R.; Griffiths, D. F. The finite difference method in partial differential equations; John Wiley, 1980. 34


Page 34 of 37

Page 35 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(27) Smith, G. D. Numerical solution of partial differential equations: finite difference methods; Oxford university press, 1985. (28) Qu, C. C.; Hahn, J. Computation of arrival cost for moving horizon estimation via unscented Kalman filtering. Journal of Process Control 2009, 19, 358–363. (29) Zheng, Y.; Li, S.; Wang, X. Distributed model predictive control for plant-wide hotrolled strip laminar cooling process. Journal of Process Control 2009, 19, 1427–1437. (30) Zheng, Y.; Li, S.; Li, N. Distributed model predictive control over network information exchange for large-scale systems. Control engineering practice 2011, 19, 757–769.

35



Graphical TOC Entry

36


Page 36 of 37

Page 37 of 37


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Moving Horizon Optimal Estimation for ... - ACS Publications

Recommend Documents