Nonlinear and Constrained State Estimation Based on the Cubature

Subscriber access provided by UCF Libraries

Article

Nonlinear and constrained state estimation based on the cubature Kalman filter Jafar Zarei, and Ehsan Shokri Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/ie4020843 • Publication Date (Web): 24 Jan 2014 Downloaded from http://pubs.acs.org on January 25, 2014

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 free 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 accessible to all readers and 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.

Industrial & Engineering Chemistry Research 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 38

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

Nonlinear and constrained state estimation based on the cubature Kalman filter Jafar Zarei1, Ehsan Shokri2 Faculty of Electrical and Electronic Engineering, Shiraz University of Technology, Shiraz, Iran. 1

2

[email protected]

[email protected] ([email protected])

Abstract: This paper investigates the use of several nonlinear estimation algorithms such as extended Kalman filter (EKF), unscented Kalman filter (UKF) and cubature Kalman filter (CKF) in the problem of state estimation in chemical processes. Three simulation case studies are considered to evaluate the performance of proposed method. Second case study uses the experimental data to investigate the accuracy of the CKF against the UKF in practical applications. Simulation results confirm the superiority of the CKF to the EKF and UKF. However, all these approaches fail to handle constraint issue in state estimation problems. Subsequently, a modified CKF is introduced to overcome linear constraint in nonlinear estimation problems. The final part of the paper shows simulation results that confirm the effectiveness of the proposed constrained CKF (CCKF). Potential profits that can be achieved while applying the proposed approach in constrained estimation problems is shown compared to conventional moving horizon estimation (MHE) algorithm. Keywords: Process monitoring, Cubature Kalman filter, Constrained cubature Kalman filter, Nonlinear state estimation, Constrained nonlinear state estimation. 1. Introduction The significant aspect of process safety and reliability is diagnosing malfunctions. To guarantee safety operation of any process, it is necessary to diagnose equipment malfunctions, process faults and some special events that may ruin the process as early as possible. Therefore, process monitoring approaches have become very popular in recent decades and play an important role to have a successful process. To achieve this aim, keeping 1 ACS Paragon Plus Environment


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

track of crucial state and parameters of a process and comparing them with upper and lower bounds is a common way for process monitoring. Nevertheless, it may be difficult or impossible to directly measure some states and most parameters of processes. Hence, it is needed to estimate important states and parameters indirectly. Thus, process monitoring works based on state estimation algorithms. In the following, some successful approaches to the state and parameter estimation of processes are reviewed. Kalman filter (KF) is the optimum solution to state and parameter estimation in linear systems.1-3 The theory of nonlinear filters is not comprehensive or complete as much as those for linear systems, and presents suboptimal approaches. Extended Kalman filter (EKF) is a development of KF exploited to state estimation in nonlinear systems.1 In this approach, a first-order linearization of nonlinear models is derived, and then Kalman filter is applied to the linearized model.4 The EKF is used to identify the process and measurement noise covariance for state and parameter estimation.5 In ref 6, the same approach is proposed for the state estimation of time-delay systems. The EKF has been widely applied to batch quality estimation problem.7 However, the EKF contains several drawbacks, for example first-order linearization can cause large errors in calculating the mean and covariance of the state vector; moreover, the system equations may be non-differentiable to derive the Jacobian matrices.8 These concerns may cause divergence of the filter especially in the complex industrial setups.9 The second-order EKF that attempts to reduce the linearization error of the conventional EKF has been proposed to overcome this issue.10,

11

This filter exploits the second order

Taylor expansions of nonlinear model functions. However, the computational efforts increase in the second-order EKF. To speed up the computation, a parallel EKF has been investigated.12 In this filter, the time update and measurement update equations of EKF are decoupled so that the time and measurement updates can be computed on separate processors. Although the mentioned approaches to overcome the flaw of EKF are useful, for non-

2 ACS Paragon Plus Environment

Page 2 of 38

Page 3 of 38

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


differentiable systems the EKF fails to estimate states. Therefore, a sophisticated filter that is free of derivatives and has a better estimation rather than EKF seems necessary. These include Bayesian filter,13 Gaussian sum filter,14 Particle filter,15, 16 Unscented Kalman filter,17, 18

and Cubature Kalman filter.19 Among these methods, a well-known nonlinear filter that overcomes flaws of EKF is

called unscented Kalman filter (UKF).17, 18 The UKF is based on unscented transformation that is a procedure to calculate the statistics of a random variable which undergoes a nonlinear function.20 The UKF can be used as an alternative to the EKF for observing process control applications.21-25 Since, the UKF does not need to calculate Jacobian matrix, it can be applied to nonlinear systems with high degree of nonlinearity. However, to achieve a proper estimation the tuning parameters of UKF should be accurately adjusted. Therefore, to reduce the complexity of adjusting parameter, it is necessary to exploit a nonlinear filter avoiding linearization while it has fixed parameters to be independent of tuning parameters. A nonlinear Kalman filter that is presented lately is called cubature Kalman filter (CKF).19 This filter is used in the state estimation of high-dimensional systems and claims more accuracy and stability than previous ones. The main idea of CKF is using cubature rule to approximate multi-dimensional integrals which are yielded from Bayesian filter under Gaussian assumption. In this filter, cubature rule is applied to spherical-radial integrals, and consequently fixed points, which are called cubature points, with a constant weight are obtained. The CKF is a derivative free filter that propagates cubature points through a nonlinear function to calculate the mean and covariance of state and measurement vectors.19 Because the parameters of CKF are appropriately adjusted, it usually has a better performance compared to the UKF. A comparison between CKF and UKF discussed form different aspects in ref

has been

19

. It is shown that the cubature points are more

accurate and more principled in mathematical terms than the sigma points approach. In ref 26,



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

CKF has been applied as cubature information filter and the square root cubature information filter has been investigated in ref 27. This filter also has been used for different purposes such as mobile robot Monte Carlo localization, spacecraft attitude estimation, and navigation.28-30 In addition, in ref 31 the cubature points have been exploited for sensor fault detection and the superiority of utilizing cubature points than sigma points has been shown. Therefore, in this paper state estimation by CKF in chemical processes will be investigated and it will be compared with conventional nonlinear filters. The performances of the mentioned filters will be studied by root mean square error (RMSE) indicator. However, in practical applications, physical limitations impose some constraints to the states or parameters of the system. For instance, the estimation of concentrations should stay positive, because negative concentration has not physical meaning. Moving horizon estimation (MHE) is a method to state estimation with constraints.32 This approach considers possible constraints on the states and requires the solution of a nonlinear programming problem at each time step. The main disadvantage of this approach is the high computational effort of numerical optimization at each time step. This computation is undesirable in critical applications that require on-line state estimation. Some methods to overcome this problem are investigated in refs

33, 34

. Nonlinear Kalman Filters have less numerical effort than MHE

algorithm. Thus a method considering constraints and a nonlinear Kalman filter simultaneously is required. Motivated by these considerations, in this paper a new approach to the constrained state estimation is proposed. The proposed approach exploits CKF, and is called constrained cubature Kalman filter (CCKF). Although in cubature Kalman filter theory the constraints are not incorporated into estimation problem, the linear constraints can be integrated in the CKF by projecting the unconstrained CKF estimation onto the boundary of the admissible region. The proposed


Page 4 of 38

Page 5 of 38

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


method to handle constraints, exploits the advantages of CKF while linear constraints on the states are considered simultaneously. To this end, first the performance of CKF has been compared with previous well-known nonlinear Kalman filters such as the UKF and the EKF. The first two case studies evaluate the performances of each filter by RMSE indicator. The experimental data is considered for the second case study to evaluate the proposed method in practical applications. Next a modified algorithm of CKF to overcome the constrained problems has been proposed. A constrained case study to compare the CKF and the CCKF has been chosen and simulation results demonstrate the efficiency of the CCKF. The proposed method is also compared with MHE algorithm, and the potential profit of the proposed method is shown. The rest of this paper is organized as follows: in Section 2, the problem of state estimation in nonlinear systems is presented. Section 3 explains the CKF algorithm briefly. The proposed method CCKF to constrained state estimation is elucidated in Section 4. Section 5 illustrates three simulation case studies to evaluate the performance of explained methods in Section 3, and Section 4, and Section 6 concludes the paper. 2. Problem formulation The state space model of a nonlinear dynamic system with additive noises is described by the pair of difference equations in discrete time.

xk = F ( xk −1 , uk −1 ) + wk −1

(1)

zk = H ( xk , uk ) + vk

(2)

Equations (1) and (2) are called process equation, and measurement equation, respectively.13 At discrete time k the state vector of dynamic system is xk ∈ , and the measurement vector is z k ∈ .The known control input uk ∈ could be derived from a controller. w k −1 and vk are zero mean and independent Gaussian noise sequences with 5 ACS Paragon Plus Environment


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 6 of 38

covariance Qk −1 , and Rk−1 , respectively. : → and : → are known functions. To monitor a process, the states or parameters of the system should be known. Nevertheless, it is not always possible or suitable to measure the states or parameters which are required. Therefor the nonlinear filters provide a means for inferring the missing information from the corrupted and noisy measurement. To be able to monitoring state space, all of the states should be available all the time. The objective is to find an estimate xˆ k of xk . To meet this objective, nonlinear Kalman filters use noisy measurement and control input to estimate true states of the dynamic system which is described by equations (1) and (2). 3. Cubature Kalman Filter Cubature Kalman filter is obtained from Bayesian filter under Gaussian assumption and uses the cubature rule to approximate the mean and covariance of a random variable propagated under a nonlinear function. In the Bayesian filter, the posterior density of the state vector can describe completely statistical properties of the state at that time. Under the Gaussian assumption, the multi-dimensional integrals, obtained from Bayesian filter, can be approximated by cubature rule.13 Consider the nonlinear system (1), (2), according to the Bayesian filter, the equations of states, covariance and measurements obtaining form Time update and Measurement update steps are given as follows: Time Update: Compute the prediction of state as: xˆ k k −1 = E F (x k −1 , u k −1 ) + w k −1 Dk −1 

(3)


Page 7 of 38

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


where wk−1 is assumed to be zero mean and uncorrelated noise sequence with measurement, and Dk −1 denotes the history of input and measurement up to k-1 , it is obtained:

xˆ k k −1 = E F(x k −1 , u k −1 ) Dk −1  =

∫ F (x

k −1

, u k −1 ) × p ( x k −1 Dk −1 ) dx k −1

R nx

=

∫

(4)

F ( x k −1 , u k −1 ) × N ( x k −1;xˆ k −1|k −1 , Pk −1|k −1 ) dx k −1

R nx

Compute associated covariance as:

(

)(

)

Pk k −1 = E  x k − xˆ k k −1 x k − xˆ k k −1  × N ( x k −1 ;xˆ k −1|k −1 , Pk −1|k −1 ) dx k −1 − xˆ

T

z k −1  = ∫ F ( x k −1 , u k −1 )F ( x k −1 , u k −1 )  nx R

T

T k |k −1 k |k −1

xˆ

(5)

+ Q k −1

Measurement Update: The prediction of measurement density is given by:

(

p ( z k Dk −1 ) = N z k ; zˆ k k −1 , Pzz , k k −1

)

(6)

Therefore, the measurement, and associated covariance are presented by: zˆ k|k −1 =

∫ H ( x , u ) × N ( x ;xˆ k

k

k

k |k −1

, Pk|k −1 ) dx k

(7)

× N ( x k ;xˆ k |k −1 , Pk |k −1 ) dxk − zˆ k |k −1zˆ Tk |k −1 + R k

(8)

R nx

Pzz , k |k −1 =

∫ H ( x , u )H ( x , u )

T

k

R

k

k

k

nx

and the cross covariance is calculated by:

Pxz ,k |k −1 =

∫ x H (x ,u )

T

k

k

k

× N ( xk ;xˆ k |k −1 , Pk |k −1 ) dxk − xˆ k|k −1zˆ Tk |k −1 + R k

(9)

R nx

When the new measurement zk is received, the Bayesian filter computes state estimation and covariance:

(

xˆ k k = x k k −1 + Wk z k − zˆ k k −1

)

(10)

Pk k = Pk k −1 − Wk Pzz ,k k −1 WkT

(11)

where the Kalman gain given as: 7 ACS Paragon Plus Environment


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

Wk = Pxz ,k k −1Pzz−1,k k −1

Page 8 of 38

(12)

As shown above, equations (4),(5),(7)-(9) may be written in the weighted multidimensional integrals form which integrand is a Nonlinear function × Gaussian density. Hence, the problem reduces to approximate the mentioned multi-dimensional integrals by cubature rules. The cubature rule to approximate n-dimensional Gaussian weighted integral with mean

µ and covariance P is as follows:

∫ F ( x )N ( x; µ , P ) dx ≈

Rn

(

1 m ∑ F µ + Pξi m i =1

)

(13)

T

where m=2n, P = P P , i.e. P is a square factor of the covariance, ξi is the i-th element of set of m cubature points introduced by:

 1   0   0   −1   0   0               0   1   .   0   −1  0   .   .   .   m  .   .   .    ,   ,   ,...,   ,   ,...,    2  .   .   .   .   .   .             .   . . . . .              0 0  −1   0   0   1 

(14)

The heart of CKF is to utilize cubature rules for approximating multi-dimensional integrals. Therefore, the CKF algorithm is summarized below: Step 1. Initialize with: xˆ 0 = E ( x 0 )

(15)

T P0 = E ( x 0 − xˆ 0 )( x 0 − xˆ 0 )   

Step 2. Generate cubature points, i = 1, 2,..., m = 2n : χ i ,k −1|k −1 = Pk −1|k −1ξ i + xˆ k −1|k −1

(16)

Step 3. Compute propagated cubature points, i = 1, 2,..., m = 2n :


Page 9 of 38

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


χ *i ,k |k −1 = F ( χ i ,k −1|k −1 , u k −1 )

(17)

Step 4. Estimate predicted state and predicted error covariance:

xˆ k |k −1 =

1 m * ∑ χ i,k|k −1 m i =1

(18)

Pk |k −1 =

1 m * χ i ,k |k −1χ *i ,Tk |k −1 − xˆ k |k −1xˆ Tk |k −1 + Q k −1 ∑ m i =1

(19)

Step 5. Redraw cubature points: χ i , k |k −1 = Pk |k −1 ξ i + xˆ k |k −1

(20)

Step 6. Estimate predicted measurement and associated covariance: Ζ i ,k |k −1 = H ( χ i ,k |k −1 , u k )

zˆ k |k −1 =

(21)

1 m ∑ Ζi,k|k −1 m i =1

(22)

Pzz ,k |k −1 =

1 m Zi ,k |k −1ΖTi ,k |k −1 − zˆ k|k −1zˆ Tk |k −1 + R k ∑ m i =1

(23)

Pxz ,k |k −1 =

1 m ∑ χ i,k|k −1ΖTi,k|k −1 − xˆ k|k −1zˆ Tk|k −1 m i =1

(24)

Step 7. Compute Kalman gain:

Wk = Pxz ,k|k −1Pzz−1,k|k −1

(25)

Step 8. Update the state estimation and associated covariance: xˆ k |k = xˆ k |k −1 + Wk ( z k − zˆ k |k −1 )

(26)

Pk|k = Pk|k −1 − Wk Pzz ,k|k −1WkT (27)

Step 9. Go to Step 2 for the next sample. In the next section, the CCKF that handle constraints on states is studied. 4. Constrained Cubature Kalman Filter



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 10 of 38

This section conducts a method to deal with the constraints of states in the CKF algorithm. Assume that the state constraint for each state is represented as: xmin ≤ xi ,k ≤ xmax

(28)

where xi ,k denotes the i-th element of state vector xk , xmin and xmax are the lower and upper bounds of xi ,k . For instance, the admissible area of a second order system according to (28) can be represented by a rectangle as shown in Figure 1. The assumed constraint information can be taken into account in the CKF algorithm during the Time Update step. After propagating cubature points, if the unconstrained cubature points which are obtained from (17) are outside of the admissible region, they can be projected onto the boundary of the admissible region and continue the CKF algorithm. For instance, as shown in Figure 1, two propagated cubature points are outside of the admissible area and they have been projected onto boundary lines. The mean and covariance regard to new constrained propagated cubature points are further update in the Time Update and

Measurement Update steps. It is worth noting that, the new propagated cubature points, mean and covariance are included the information of constraints, thus it is able to handle the constraint on states in the CCKF algorithm. Extending the proposed method to higher dimension is straight forward. By this method it is possible to handle linear constraints such as:

Cxk ≤ d

(29)

The propagated cubature points which do not satisfy the inequality (29) can be projected onto the boundary of admissible area. Note that, after generating cubature points in Step 2 and redrawing cubature points in Step 5 during the CKF algorithm, state projecting can be applied. Sometimes constraints are nonlinear instead of linear form in (29) as: 10 ACS Paragon Plus Environment

Page 11 of 38

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


g(xk ) ≤ q

(30)

In these cases Taylor series expansion of the constraint equation around xˆ k −1|k −1 can be performed to obtain:35 g (x k ) ≈ g ( xˆ k −1|k −1 ) + g′ ( xˆ k −1|k −1 )( x k − xˆ k −1|k −1 ) + φ (x k )

(31)

where φ ( x k ) is the remaining terms of Taylor series expansion. By neglecting the second and higher-order terms the constraint equation (30) can be rewritten:36, 37 g′ ( xˆ k −1|k −1 ) x k ≤ q − g ( xˆ k −1|k −1 ) + g ′ ( xˆ k −1|k −1 ) xˆ k −1|k −1

(32)

This equation is equivalent to linear constraint equation in (29) as: C = g′ ( xˆ k −1|k −1 )

d = q − g ( xˆ k −1|k −1 ) + g′ ( xˆ k −1|k −1 ) xˆ k −1|k −1

(33)

Therefore, the mentioned method can be deal with nonlinear constraints after linearization. Note that after linearizing the nonlinear constraints, they can be approximately satisfied rather than exactly. This kind of approximation can be applied when the difference between nonlinear constraint and linearized one is negligible, or where the constraint equation has some uncertainties.35 Consequently, the CCKF algorithm can be written as follows: Step 1. Initialize with: xˆ 0 = E ( x 0 )

(34)

T P0 = E ( x 0 − xˆ 0 )( x 0 − xˆ 0 )   

Step 2. Generating cubature points: χ i ,k −1|k −1 = Pk −1|k −1ξ i + xˆ k −1|k −1

(35)



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 12 of 38

Step 3. After generating cubature points by (35) the points out of admissible area are projected onto boundary to obtain the constrained cubature points by using a projection function as follows: χ Ci ,k −1|k −1 = Pr ( χ i ,k −1|k −1 )

(36)

where Pr(.) denotes the projection function and C is abbreviation of constrained. Step 4. Propagate the constrained cubature points. Transform the constrained cubature points through the process function F (.) : χ *i ,Ck |k −1 = F ( χ Ci ,k −1|k −1 , u k −1 )

(37)

Step 5. Estimate predicted state and predicted error covariance:

xˆ Ck |k −1 =

1 m *C ∑ χ i,k|k −1 m i =1

(38)

PkC|k −1 =

1 m *C χ i , k |k −1χ *i ,Ck |k −1T − xˆ Ck |k −1xˆ Ck |k −1T + Q k −1 ∑ m i =1

(39)

Step 6. Redraw cubature points: χ Ci , k |k −1 = P C k |k −1ξi + xˆ Ck |k −1

(40)

Step 7. Apply checking and projecting points to (40) and after achieving new constrained points continue the ordinary CKF algorithm. If CKF estimates the states out of admissible area in (26) the same method can apply to it. Step 8. Go to Step 2 for the next sample. As shown in equations (36)-(40) and Step 7 the constrained information is interfering to calculate the covariance matrices and estimation of states. In the next section, the CKF and constrained CKF are applied to applications and are compared against previous methods. 5. Case studies


Page 13 of 38

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 this section three different chemical applications are discussed to investigate the performance of the proposed methods. First Case study is an irreversible highly exothermic chemical reaction and has a nonlinear process equation with linear measurement equation. Second one is a pH process with nonlinear measurement equation. Third Case study is a constrained state estimation and uses the proposed method in Section 4. It is attempted to compare proposed methods against previous ones especially UKF algorithm. The tuning parameter κ might have significant influence on UKF performance in some case studies. Nevertheless, the CKF has fixed parameters to generate cubature points. Thus, to make a fair comparison different tuning parameter, κ, is considered. In addition, different process and measurement noise covariance lead to different estimations performances. Therefore, the varying process and measurement noise levels are taken into account. In addition all the results are obtained over 100 Monte Carlo independent runs. A performance indicator that is exploited in this paper is root mean square error (RMSE). The RMSE for variable x is defined as follows:

RMSE ( x ) =

1 N

N

∑ ( x − xˆ ) i

2

(41)

i

i =1

where x is the real data and xˆ is estimated one.

Case study I) A nonlinear chemical process adopted from ref

38

is considered to

investigate the performance of the CKF, UKF and EKF. The behavior of an irreversible highly exothermic chemical reaction A → B which is taken place in a non-adiabatic CSTR can be described by nonlinear ordinary differential equations (ODEs). In addition the reactor wall has sufficient effect on the system dynamics and it should be taken into account. Accordingly, the model of CSTR is described by the following equations:



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 14 of 38

dC A FR = ( C A0 − C A ) − k0 e− Ea / RTR C A dt VR

( hAT )i dTR FR = (T0 − TR ) + dt VR VR ( ρ C p )

(TR − TW ) +

− ( ∆H R )

f

( ρC p ) f

k0 e − Ea / RTR C A (42)

( hAT )i ( hAT )e dTW = (TR − TW ) − (TW − TJ ) dt VW ( ρ C p ) W VW ( ρ C p )W ( hAT )e dTJ FJ = (TJ 0 − TJ ) + dt VJ VJ ( ρ C p )

(TW − TJ ) J

Figure 2 illustrates the schematic of this CSTR. This model can be rewritten in normalized dimensionless form as: dx1 dt dx2 dt dx3 dt dx4 dt

= p1u1 + p1u1u2 − p1u2 x1 − p2 e− p3 /(1+ x2 ) (1 + x1 ) − p1 x1 = p1u3 + p1u2u3 − p1 x2 − p1u2 x2 − p4 x2 + p4 x3 + p5 p2 e − p3 /(1+ x2 ) (1 + x1 )

(43) = p6 x2 − p6 x3 − p7 x3 + p7 x4 = p8u4 + p8u4u5 − p8 x4 − p8u5 x4 + p9 x3 − p9 x4

The state vector is defined as:

 CA − CA,ref TR − TR,ref TW − TW , ref TJ − TJ ,ref  x=  TR,ref TW ,ref TJ ,ref   CA,ref

T

(44)

where CA is the concentration of reactant A. TR ,TW , and TJ are temperatures of reactor, wall and jacket, respectively. The unit of CA is mol/m3, and the unit of TR, TW , and TJ is K. The input vector is:

 CA0 − CA,ref FR − FR, ref T0 − T0,ref TJ 0 − TJ 0,ref FJ − FJ ,ref  u=  FR,ref T0,ref TJ 0,ref FJ ,ref   CA,ref

T

(45)

where CA0 is the feed concentration of reactant A, FR the flow rate. T0 and TJ0 are the feed temperature, and the inlet jacket temperature, respectively, and FJ is the coolant flow rate. The unit of CA0 is mol/m3, the unit of T0 and TJ0 is K, and the unit of FR and FJ is m3/s. The corresponding references values are: 14 ACS Paragon Plus Environment

Page 15 of 38

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


CA0,ref = CA,ref = 3 mol / m3 , FR,ref = 60 ×10−5 m3 / s , FJ ,ref = 15×10−4 m3 / s , FR,ref = 60 ×10−5 m3 / s , TR ,ref = Tw,ref = TJ ,ref = T0,ref = TJ 0,ref = 298 K Available measurements are x1, x2 and x4. The third state x3 which denotes the wall temperature is not measureable, hence the output is assumed to be:

y = Cx with

1 0 0 0  C =  0 1 0 0   0 0 0 1 

(46)

The model parameters used in the present application are summarized in Table 1. In addition, this model has two stable steady states and one unstable equilibrium points, which are presented in Table 2. This case study has been carried out by applying the CKF, the UKF and the EKF to state estimation. The obtained results of filters are compared by RMSE indicator. To have a fair comparison different process and measurement noise levels are taken into account. In addition, in the case of UKF the different selections of tuning parameter κ are considered and we tried to choose a κ that results in higher performance. Note that all of the results obtained over 100 Monte Carlo independent runs. The process equation (43) has been discretized by Rung-Kutta fourth order method with ∆t = 0.1s . The true data is generated by using discrete process equation, the initial state is set as x 0 = [ −0.97540, 0.47245, 0.42690, 0.37934] , the T

process noise covariance and measurement noise covariance are different for each simulation and different choices for κ has been considered. In this case two scenarios are followed. In the first scenario the system works normally and the results of state estimation are given. Second scenario investigates the performance of each filter when a sudden change happens in the parameters.

Scenario I) the system works normally and the parameters of the filters are set as:

xˆ 0 = [ −0.9,0.4,0.4,0.4] ,



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

P0 = 10 −6 I 6 , Q k = 10−6 I 4 k = 1, 2,... , R k = 10−6 I 3 .

According to simulation results the tuning parameter κ has not significant effect on the performance of UKF and the κ=-3 has been selected in all simulations. The results showed the difference between RMSEs obtaining from UKF for different values κ=-1,-2,-3,-3.9 are negligible. Table 3 lists the RMSE that is generated by the CKF, the UKF and the EKF when the system is subjected to different measurement noise levels. It can be perceived, the EKF has the worse RMSE than other filters. Since the EKF works based on linearization and the more difference between linearized model and real one causes more error in state estimation. Note that, both CKF and UKF avoid linearization and works based on exploiting sample points naming cubature points and sigma points. Table 3 demonstrates in some cases the CKF outperforms and in some cases the converse happens. Since in this case study the parameters of UKF have been adjusted well this filter has acceptable state estimation. Although the difference between the RMSE of CKF and the RMSE of UKF is not large, by choosing improper tuning parameter for UKF it may leads to worse estimation or instability while the CKF has fixed parameters in generating cubature points and stay on a fix acceptable performance. The RMSEs results of the second state to fourth state are written in Tables 4 to 6. By conducting the obtained results the same discussion for the first state can be extended to other estimation of states. To have a visual comparison, one of the simulation results is illustrated in Figure 3. For brevity only the state estimation of first state is shown. Figure 3 is obtained when Q=10-6 and R=10-2. As shown in this Figure the CKF and UKF outperform the EKF and have almost


Page 16 of 38

Page 17 of 38

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


same accuracy in estimation. The same condition for setting the parameters of filters is considered.

Scenario II) in this scenario the input vector at a specific time suddenly has been set to zero and after a period the input vector would get the last value before changing. The input is set as u k = [ 0 0 0 0 0] when 100 < k < 150. Hence, in this interval of samples the states are T

differ from their steady states and the performance of state estimation regard to sudden changes can be conducted. For instance the estimation of the first state for Q=10-6 and R= 102

is illustrated in Figure 4. Although all filters can follow the sudden change in the system,

the different performances cause different RMSEs showing in Tables 7 to 10. As shown in Tables 7 to 10, among the mentioned filters the EKF has the worst performance. According to the discussion in the Scenario I the performance of the CKF and UKF are rather same. Note that, because of the sudden change, the RMSEs of the second Scenario are a bit bigger than first Scenario. However, the differences of RMSEs between first and second Scenarios are negligible and it shows the mentioned filters can follow the sudden changes with the same performances.

Case study II) This case study investigates the performances of UKF and CKF regard to a pH process with nonlinear measurement function taking place within a continuous stirred tank reactor (CSTR) in which a neutralization reaction takes place between a strong acid (HA) and a strong base (BX). In this case study experimental data are considered to evaluate the mentioned nonlinear filters. 39 The state vector is defined as:

x1 =  A−  , x2 =  B +  , x3 =  X −  x = [ x1 x2 x3 ]

T

where  A−  ,  B +  and  X −  are the concentrations of the acid, base and the buffer agent, respectively. The process equations are given as follows: 17 ACS Paragon Plus Environment


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

x&1 =

1

(x θ

1,i

x&2 = − x&3 = −

1

θ 1

θ

− x1 ) −

Page 18 of 38

1 x1q B V

(47)

x2 +

1 ( x2,i − x2 ) qB V

(48)

x3 +

1 ( x3,i − x3 ) qB V

(49)

F ( x, ξ ) ≡ ξ + x2 + x3 − x1 −

Kw

ξ

−

x3

1 + ( K xξ / K w )

=0

(50)

where ξ = 10− pH and θ =V / q A . The term qA and qB represent the inlet flow rates. x1,i , x2,i and x3,i are the inlet concentrations of species, and K x and Kw are the dissociation

constant of the buffer and water, respectively. The system parameters are summarized in Table 11. By solving equation (50) pH is given by: 1/3 1/3   3 3  pH = − log   q + q 2 + ( r − p 2 )  +  q − q 2 + ( r − p 2 )  +     

where

q = p3 +

a = 1,

b=

Kw + x3 + x2 − x1 , Kx

c = ( x2 − x1 − K x )

 p 

Kw , Kx

(51)

d =−

K w2 , Kx

p=−

b , 3a

bc − 3ad c , and r = . 6a 2 3a

The output equation (51) illustrates a severe nonlinearity. As discussed in ref

22

the

strongest nonlinear behavior happens when pH ∈ [ 4,6.5] . The interval of pH in our experimental data is pH ∈ [3.69,11.79] . The sample rate of experimental data is ∆t = 1s . Hence, the process equation described by (47)-(49) has been discretized by Rung-Kutta method with ∆t = 1s . The equation (51) denotes the nonlinear measurement equation. The objective of presenting this simulation is to compare the CKF and the UKF performances in a nonlinear environment under experimental


Page 19 of 38

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


data. In addition, in the case of UKF five different value of tuning parameter κ is considered to achieve better comparison. The parameters of the filter are set as:

xˆ 0 = [ 0.6 0.2 0.4] ×10−3 , P0 = 10−3 I3 , Q k = 10−6 I3 , k = 1, 2,... ,

R k = 10−1 k = 1, 2,... . To evaluate the performance of nonlinear filters the RMSE for output is written in Table 12. In this case study choosing a proper value for tuning parameter κ leads to high performance of the UKF. For a fair comparison, we try to adjust the tuning parameter κ leading to the high performance of UKF. The experimental data include acid flow in liters, base flow in liters as input vector and measurement of pH as output. According to the RMSEs which are written in Table 12, the UKF results in high performance when the tuning parameter is κ=0. Figure 5 shows the output of pH process and a comparison between the CKF and the UKF state estimations when κ=0. To compare the results more precisely, the zoomed areas 1 and 2 in Figure 5 are shown in Figure 6 and Figure 7. In addition, the state estimation for different tuning parameters of κ for the UKF algorithm is shown in Figures 6 and 7. Figure 6 demonstrates that the CKF has better state estimation than UKF. Furthermore, when pH belongs to the high interval of nonlinearity in measurement equation, the state estimation by CKF has proper outcome. However, in the highest performance of UKF, i.e. κ=0, the estimation by UKF is not acceptable as much as CKF one. In addition, the performance of UKF depends on choosing a proper tuning parameter κ, where κ=1 and κ=2 the UKF output estimation is far from true data. Nevertheless, the CKF is free of adjusting parameters and has the fix cubature



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 38

points, which is an advantage for the CKF. Figure 7 illustrates the zooming area 2, where the pH is in the soft nonlinearity interval. This Figure demonstrates the better estimation of CKF in this area than the UKF. Consequently, Table 12, Figures 5, 6, and 7 verifies the advantages of CKF and show the higher performance of CKF in practical application with experimental data.

Case study III) A model of CSTR under constraints is considered. In this system, concentrations are needed to be estimated while the estimation of concentrations should be positive because the negative concentration has not physical meaning. Assume the reversible gas phase takes place in a batch reactor.40

A ⇌ B+C

2B ⇌ C

(52)

The first principle model for a constant volume, well mixed, isothermal batch reactor is given as follows:

dC A = − k1C A + k2CB Cc dt dCB = k1C A − k 2CB Cc − 2k3CB2 + 2k4CC dt dCC = k1C A − k 2CB Cc + k3C B2 − k4 CC dt where

[ k1

(53)

k2 k3 k4 ] = [ 0.5 0.05 0.2 0.01] and C j indicates the concentration of

species. The system state vector is denoted as x = [C A C B Cc ] and the measurement model T

is

the

pressure

which

is

assumed

to

be

y = [ RT RT RT ] x .

For

this

case

RT = 3.84 mol atm/L . The state space model is written corresponding to concentrations and

the nonnegative estimation of states is imposed at each time step. Hence the constraint for all step times can be written as:

xi > 0 i = 1, 2,3.

(54) 20 ACS Paragon Plus Environment

Page 21 of 38

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


The model described by equation (53) has been discretized using fourth order RungeKutta method with ∆ t = t k +1 − t k = 0.25 . Moreover, it is considered that the system is corrupted by Gaussian noises both in the state and measurements. The parameters used to generate true data are assumed as: x 0 = [0.5 0.05 0]T , Q = 0.012 I 3 , R = 0.252.

The initial parameters used in both CCKF and CKF are:

xˆ 1,0 = [0 0 0]T , xˆ 2,0 = [0 0 4]T Q k = 0.012 I 3

k = 1, 2,... ,

Rk = 0.252 k = 1, 2,... ,

Pk = 0.52 I 3 .

In this case study, first we try to compare CCKF and CKF and evaluate the performances of each filter with regard to constraint. Next, a comparison between CCKF and MHE algorithm, which is a conventional algorithm for constrained state estimation, has been drawn. Two initial guesses for the initial state xˆ 0 are assumed, and the simulation results show that the convergence of CKF to true states depends on choosing initial states. Nevertheless, with both initial guesses the CCKF converges to true states. Further, the effect of process and measurement noise covariance is investigated and the results for each state are written in Tables 13 to 15. First, the initial state xˆ 1,0 is chosen and the results are shown in Figure 8. The true state is illustrated by continuous black line, the blue thin dashed line denotes the estimation results of CKF and the red thick dashed line shows the estimated states by CCKF. As shown in this figure the CKF has not been converged to true state and the estimation of the second and third 21 ACS Paragon Plus Environment


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

states are negative that is unacceptable. However, the CCKF estimates the state truly and able to tolerates constraint well. Figure 9 shows the results when the initial state xˆ 2,0 is assumed. Although both CKF and CCKF converge to true states finally, the CKF estimates the concentrations negative in an interval of time. For instance, the third state is estimated negative before 20-th samples of time. In addition, as shown in Figure 9 the convergence rate of CCKF is faster than CKF. At each step time for each state six cubature points are generated that some of them are out of admissible area. The behavior of the generated cubature points and projected cubature points in the CCKF algorithm for the first state when the initial state xˆ 2,0 is chosen has been studied. The result for the first five steps time is illustrated in Figure 10. In this interval, the second, third, and sixth cubature points are the same with their projections and they have not stood out of the admissible area. The first, fourth, and fifth cubature points stood out of admissible area and the projected ones that used to estimate the predicted state are shown in Figure 10. Note that, the boundary line is x1=0. This figure shows the information of constraint is imposed to estimate the predicted state and subsequently to calculate the predicted covariance and other steps of the CCKF algorithm. This case study focuses on state estimation under positive constraints, hence to compare the CCKF and CKF performance the RMSE of each state has been written in Tables 13 to 15. To make a fair comparison, process and measurement noise are assumed to be varying, and the results have been obtained over 100 independent Monte Carlo runs. As shown in Tables 13 to 15, it is clear that the estimation of CCKF is more accurate than CKF because the constrained information has not been taken into account in the CKF algorithm. This case study is investigated in ref

41

and moving horizon estimation algorithm has

been proposed to this state estimation problem. To evaluate the performance of the proposed CCKF, it is compared with MHE algorithm. To draw a comprehensive comparison different 22 ACS Paragon Plus Environment

Page 22 of 38

Page 23 of 38

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


initial states are considered. At first, Initial parameters for CCKF and MHE algorithms are selected as:

xˆ = [0.5 0 0]T Q k = 0.012 I 3

k = 1, 2,... ,

Rk = 0.252 k = 1, 2,... ,

Pk = 0.52 I 3 .

Note that in this case initial states are chosen close to the true states. Figure 11 illustrates the estimated state by CCKF and MHE algorithms. As shown in this figure, both CCKF and MHE are able to estimate the state subject to constraint truly. Table 16 shows the RMSE and computation time of CCKF compared to MHE algorithms. According to this table the MHE algorithm has the better estimation accuracy than CCKF. Although the MHE algorithm converges to the true state faster, the computation burden of MHE is much more than CCKF. In the next case initial guess is selected as xˆ = [2 1 1]T , which is far from the true state vector compared to previous one. The results for both CCKF and MHE algorithms are shown in Figure 12. Moreover, the RMSE and computation time of CCKF compared to MHE algorithms, for this case, are summarized in Table 17. The obtained results show that for this case the CCKF not only has the better estimation but also has the less computational effort. For another initial guess when xˆ = [0 0 4]T the MHE algorithm diverges while the CCKF, as shown in Figure 9, converges to the true state. Therefore, the MHE algorithm has the high performance when the initial guess of state vector was close to the true state. Moreover, note that, the difference of RMSEs between CCKF and MHE, when initial states are close to the true states, is not major and both algorithms have acceptable performance. In addition, the MHE needs complex and heavy numerical computation in each step time, while the CCKF has less computational efforts than MHE and has acceptable performance for state estimation.



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

Consequently, to have a tradeoff between computational efforts and accuracy of estimation, for linear constraints and nonlinear constraints, that can be linearized, the CCKF can be used as an alternative to the MHE algorithm. 6. Conclusion In this paper, the application of the CKF in three chemical processes which one of them has been investigated by experimental data is assessed. It has been shown that using cubature points leads to more accurate estimation than using sigma points, which are utilized in the UKF. Moreover, the EKF has larger estimation error than the CKF and UKF. It is obvious that the CKF algorithm needs less computational effort than the UKF because the CKF exploits 2n sample points and the UKF exploits 2n+1 sample points to state estimation. In addition, the CKF and the UKF have been evaluated by experimental data and the obtained results demonstrate the higher performance of CKF. Consequently, the CKF can be successfully applied to nonlinear systems for state estimation purposes. Next, constrained version of CKF (CCKF) is proposed, and its performance has been investigated by simulation. Consequently, due to the acceptable performance of CCKF and lower computational effort, it can be used as an alternative to MHE to handle linear constraints and soft nonlinear constraints in nonlinear state estimation problems. References

(1)

Anderson, B.; Moore, J., Optimal Filtering. Prentice-Hall, Englewood Cliffs, NJ:

1979. (2)

Davis, M. H. A., Linear estimation and stochastic control. London: 1977.

(3)

Jacobs, O. L. R.; Jacobs, O., Introduction to control theory. Oxford University Press

Oxford, UK: 1993; Vol. 2. (4)

Ghil, M.; Malanotte-Rizzoli, P., Data assimilation in meteorology and oceanography.

Adv. Geophys 1991, 33, 141-266.


Page 24 of 38

Page 25 of 38

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


(5)

Bavdekar, V. A.; Deshpande, A. P.; Patwardhan, S. C., Identification of process and

measurement noise covariance for state and parameter estimation using extended Kalman filter. Journal of Process Control 2011, 21, (4), 585-601. (6)

Gopalakrishnan, A.; Kaisare, N. S.; Narasimhan, S., Incorporating delayed and

infrequent measurements in Extended Kalman Filter based nonlinear state estimation.

Journal of Process Control 2011, 21, (1), 119-129. (7)

Dorsey, A. W.; Lee, J. H.; Russell, S. A., An extended Kalman filter formulation for

systematic transfer of information from batch to batch. Industrial & engineering chemistry

research 2003, 42, (8), 1753-1760. (8)

Llinas, J.; Hall, D. L.; Liggins, M. E., Handbook of Multisensor Data Fusion: Theory

and Practice. CRC Press: 2009. (9)

Wilson, D.; Agarwal, M.; Rippin, D., Experiences implementing the extended

Kalman filter on an industrial batch reactor. Computers & chemical engineering 1998, 22, (11), 1653-1672. (10)

Athans, M.; Wishner, R.; Bertolini, A., Suboptimal state estimation for continuous-

time nonlinear systems from discrete noisy measurements. IEEE Transactions on Automatic

Control 1968, 13, (5), 504-514. (11)

Gelb, A., Applied optimal estimation. The MIT press: 1974.

(12)

Lu, M.; Qiao, X.; Chen, G., Parallel computation of the modified extended kalman

filter. International journal of computer mathematics 1992, 45, (1-2), 69-87. (13)

Ho, Y.; Lee, R., A Bayesian approach to problems in stochastic estimation and

control. IEEE Transactions on Automatic Control 1964, 9, (4), 333-339. (14)

Alspach, D.; Sorenson, H., Nonlinear Bayesian estimation using Gaussian sum

approximations. IEEE Transactions on Automatic Control 1972, 17, (4), 439-448. (15)

Arulampalam, M. S.; Maskell, S.; Gordon, N.; Clapp, T., A tutorial on particle filters

for online nonlinear/non-Gaussian Bayesian tracking. IEEE Transactions on Signal

Processing 2002, 50, (2), 174-188. (16)

Gordon, N. J.; Salmond, D. J.; Smith, A. F. Novel approach to nonlinear/non-

Gaussian Bayesian state estimation, In IEE Proceedings on Radar and Signal Processin, 1993, 107-113. (17)

Julier, S.; Uhlmann, J.; Durrant-Whyte, H. F., A new method for the nonlinear

transformation of means and covariances in filters and estimators. IEEE Transactions on

Automatic Control 2000, 45, (3), 477-482.



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

(18)

Julier, S. J.; Uhlmann, J. K., Unscented filtering and nonlinear estimation.

Proceedings of the IEEE 2004, 92, (3), 401-422. (19)

Arasaratnam, I.; Haykin, S., Cubature Kalman Filters. IEEE Transactions on

Automatic Control 2009, 54, (6), 1254-1269. (20)

Wan, E. A.; Van Der Merwe, R. The unscented Kalman filter for nonlinear

estimation, In Adaptive Systems for Signal Processing, Communications, and Control Symposium (AS-SPCC), 2000; 153-158. (21)

Akin, B.; Orguner, U.; Ersak, A. State estimation of induction motor using unscented

Kalman filter, In Proceedings of IEEE Conference on Control Applications, 23-25 June, 2003, 915-919. (22)

Romanenko, A.; Santos, L. O.; Afonso, P. A., Unscented Kalman filtering of a

simulated pH system. Industrial & engineering chemistry research 2004, 43, (23), 75317538. (23)

A. Romanenko, J. A. A. M. C., The unscented ﬁlter as an alternative to the EKF for

nonlinear state estimation: a simulation case study. Computer and Chemical Engineering 2004, 28, 347–355. (24)

Mandela, R.; Kuppuraj, V.; Rengaswamy, R.; Narasimhan, S., Constrained unscented

recursive estimator for nonlinear dynamic systems. Journal of Process Control 2012. (25)

Kandepu, R.; Foss, B.; Imsland, L., Applying the unscented Kalman filter for

nonlinear state estimation. Journal of Process Control 2008, 18, (7–8), 753-768. (26)

Pakki, K.; Chandra, B.; Gu, D. W.; Postlethwaite, I. Cubature information filter and

its applications, In American Control Conference, 2011, 3609-3614. (27)

Chandra, K. P. B.; Da-Wei, G.; Postlethwaite, I., Square Root Cubature Information

Filter. IEEE Sensors Journal 2013, 13, (2), 750-758. (28)

Li, Q.; Song, Y. Cubature MCL: Mobile robot Monte Carlo Localization based on

Cubature Particle Filter, In 31st Chinese Control Conference (CCC), 2012, 5141-5145. (29)

Xiaojun, T.; Jianli, W.; Kai, C. Square-root adaptive cubature Kalman filter with

application to spacecraft attitude estimation, In 15th International Conference on Information Fusion, 9-12 July 2012, 1406-1412. (30)

Pesonen, H.; Piche; x; R. Cubature-based Kalman filters for positioning, In 7th

Workshop on Positioning Navigation and Communication (WPNC) 11-12 March 2010, 4549. (31)

Zarei, J.; Shokri, E., Robust Sensor Fault Detection Based on Nonlinear Unknown

Input Observer. Measurement 2014, 48, 355-367. 26 ACS Paragon Plus Environment

Page 26 of 38

Page 27 of 38

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


(32)

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, (2), 246-258. (33)

Zavala, V. M.; Laird, C. D.; Biegler, L. T., A fast moving horizon estimation

algorithm based on nonlinear programming sensitivity. Journal of Process Control 2008, 18, (9), 876-884. (34)

Zavala, V. M., Stability analysis of an approximate scheme for moving horizon

estimation. Computers & Chemical Engineering 2010, 34, (10), 1662-1670. (35)

Simon, D., Kalman ﬁltering with state constraints: a survey of linear and nonlinear

algorithms. IET Control Theory & Applications 2009, 4, (8), 1303-1318. (36)

Alouani, A.; Blair, W., Use of a kinematic constraint in tracking constant speed,

maneuvering targets. IEEE Transactions on Automatic Control 1993, 38, (7), 1107-1111. (37)

Simon, D.; Chia, T. L., Kalman filtering with state equality constraints. IEEE

Transactions on Aerospace and Electronic Systems, 2002, 38, (1), 128-136. (38)

Romanenko, A.; Castro, J. A. A. M., The unscented filter as an alternative to the EKF

for nonlinear state estimation: a simulation case study. Computers & chemical engineering 2004, 28, (3), 347-355. (39)

McAvoy, T. J.; Hsu, E.; Lowenthal, S., Dynamics of pH in controlled stirred tank

reactor. Industrial & Engineering Chemistry Process Design and Development 1972, 11, (1), 68-70. (40)

Haseltine, E. L.; Rawlings, J. B., Critical evaluation of extended Kalman filtering and

moving-horizon estimation. Industrial & engineering chemistry research 2005, 44, (8), 24512460. (41)

Qu, C. C.; Hahn, J., Process monitoring and parameter estimation via unscented

Kalman filtering. Journal of Loss Prevention in the Process Industries 2009, 22, (6), 703709.



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

bounds

projected points

admissible area

points out of admissible area Figure 1. Admissible area and point projecting for second order systems x k ∈ 2.

Figure 2. Schematic of a non-adiabatic CSTR.


Page 28 of 38

Page 29 of 38

True CKF UKF EKF

CA concentraion

-0.985 -0.99 -0.995 -1 -1.005 -1.01 -1.015 0

10

20

30 40 Time (sample)

50

60

70

Figure 3. Scenario I) state estimation of x1.

True CKF UKF EKF

-0.975 -0.98 -0.985 -0.99

A

C 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


-0.995 -1 0

20

40

60

80

100 120 Time (sample)

Figure 4. Scenario II) state estimation of x1.


140

160

180

200


12 11

pH measurement

10

True CKF UKF κ =0

Zooming area 2

9 8 7 6 5 4 3 0

Zooming area 1 200 400 600

800 1000 1200 TIme (sample)

1400

1600

1800

2000

Figure 5. The state estimation of pH measurement correspond to experimental data. Zooming area 1 True CKF UKF κ =-2 UKF κ =-1 UKF κ =0 UKF κ =1 UKF κ =2

12

pH Measurement

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 38

10 8

6 4

2 0

10

20

30

40

50 60 Time (sample)

70

80

90

100

Figure 6. Zooming area 1 of Figure 5, CKF and UKF estimations correspond to different κ.


Page 31 of 38

Zooming area 2 True CKF UKF κ =-2 UKF κ =-1 UKF κ =0 UKF κ =1 UKF κ =2

pH Measurement

11.4

11.2

11

10.8

10.6 1200

1210

1220

1230

1240 1250 1260 Time (sample)

1270

1280

1290

1300

Figure 7. Zooming area 2 of Figure 5, CKF and UKF estimations correspond to different

True state CKF CCKF

0.5 0 -0.5 0

CB concentration

CA concentration

κ.

CC 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


50

100

150 True state CKF CCKF

1 0 -1 0

50

100


2 0 -2 0

50

100

150

Time (sample)

Figure 8. The states estimation that are obtained by CCKF and CKF algorithms correspond to first initial guess.


CC 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

CB concentration CA concentration


1.5 1 0.5 0

Page 32 of 38

True state CKF CCKF 20

40

60

80

100

1.5

120

140

True state CKF CCKF

1 0.5 0 20

40

60

80

100

120

140


0 -1

20

40

60 80 100 Time (sample)

120

140

Figure 9. The states estimation that are obtained by CCKF and CKF algorithms correspond to second initial guess.


First cubature point First projected cubateure point Fouth cubature point Fourth projected cubature point Fifth cubature point Fifth projected cubature point

1

0.5

0

-0.5

-1 0.5

1

1.5

2

2.5 3 Time (step)

3.5

4

4.5

5

C

A

Figure 10. Evaluating cubature points before and after using projection function.

True state MHE CCKF

0.4 0.2 0

10

20

30

40

50

60

C

B

0.4 0.2 0

C

True state MHE CCKF 10

20

30

40

0.5

0

50

60

True state MHE CCKF

C

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


Cubature points and project cubature points for first state

Page 33 of 38

10

20

30 Time (step)

40

50

60

Figure 11. Estimated states by CCKF and MHE algorithms when initial guess of state vector is close to the true states.


C

1.5 1 0.5 0

Page 34 of 38


20

30

40

50

60

70

C

B

1

C

True state MHE CCKF

0.5 0

C

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

A


10

20

30

40

50

0.8 0.6 0.4 0.2 0

60

70


20

30 40 Time (sample)

50

60

70

Figure 12. Estimated states by CCKF and MHE algorithms when initial guess of state vector is far from the true states.


Page 35 of 38

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


Table 1: Model parameters Parameter

Unit

Expression

value

FR ,ref

3.333 × 10 −2

p1

VR

s −1

p2

k0

4.08×107

s −1

p3

Ea RTR ,ref

25.347

-

p4

( hAT )i VR ( ρC p ) f

6 . 63 × 10 − 1

s −1

1.45

-

p5

−( ∆H R )C A ,ref ( ρC p ) f TR ,ref

p6

(hAT ) i VW ( ρC p )W

5.97

s −1

p7

(hAT ) e VW ( ρC p )W

5.97

s −1

1 . 67 × 10 − 1

FJ ,ref

p8

s −1

VJ

(hAT )e VJ (ρCp ) J

p9

s −1

1.33

Table 2: Steady states Low-temperature Unstable

High-temperature stable

stable

x1

-0.0140582

-0.37748

-0.97640

x2

0.0068168

0.18304

0.47345

x3

0.0061321

0.16465

0.42590

x4

0.0054473

0.14627

0.37834



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

Table 3: Scenario I) RMSE computation of x1 over 100 independent Monte Carlo run for varying process noise and measurement noise levels. Q 10-6I4 10-4I4 -3 -2 -3 R 10 I3 10 I3 10 I3 10-2I3 -4 -3 -4 6.33×10 6.3×10 6.32×10 6.4×10-3 CKF -4 -3 -4 6.89×10 5.8×10 6.82×10 5.9×10-3 UKF -4 -3 -4 9.83×10 9.8×10 9.85×10 9.9×10-3 EKF

Table 4: Scenario I) RMSE computation of x2 over 100 independent Monte Carlo run for varying process noise and measurement noise levels. 10-4I4 Q 10-6I4 -3 -2 -3 R 10 I3 10 I3 10 I3 10-2I3 -4 -3 -4 6.7×10 6.64×10 6.6×10-3 6.61×10 CKF -4 -3 -4 6.5×10 6.94×10 6.5×10-3 6.91×10 UKF -4 -3 -4 9.8×10 9.9×10 9.89×10 9.9×10-3 EKF

Table 5: Scenario I) RMSE computation of x3 over 100 independent Monte Carlo run for varying process noise and measurement noise levels. Q 10-6I4 10-4I4 -3 -2 -3 R 10 I3 10 I3 10 I3 10-2I3 -4 -3 -4 7.88×10 3.1×10 7.89×10 3.0×10-3 CKF -4 -3 -4 6.67×10 3.2×10 6.66×10 3.2×10-3 UKF -4 -3 -4 8.29×10 4.0×10 8.29×10 3.9×10-3 EKF

Table 6: Scenario I) RMSE computation of x4 over 100 independent Monte Carlo run for varying process noise and measurement noise levels. Q 10-6I4 10-4I4 -3 -2 -3 R 10 I3 10 I3 10 I3 10-2I3 -4 -3 -4 6.47×10 6.5×10 6.49×10 6.4×10-3 CKF 9.18×10-4 6.2×10-3 9.17×10-4 6.1×10-3 UKF -4 -3 -4 9.88×10 9.9×10 9.91×10 9.8×10-3 EKF

Table 7: Scenario II) RMSE computation of x1 over 100 independent Monte Carlo run for varying process noise and measurement noise levels. Q 10-6I4 10-4I4 -3 -2 -3 R 10 I3 10 I3 10 I3 10-2I3 -4 -3 -4 6.5×10 6.55×10 5.63×10 6.3×10-3 CKF -4 -3 -4 8.87×10 7.8×10 7.2×10 6.1×10-3 UKF -4 -3 -4 18×10 9.9×10 18×10 10.8×10-3 EKF

Table 8: Scenario II) RMSE computation of x2 over 100 independent Monte Carlo run for varying process noise and measurement noise levels. Q 10-6I4 10-4I4 -3 -2 -3 R 10 I3 10 I3 10 I3 10-2I3 -4 -3 -4 6.73×10 6.5×10 6.72×10 6.6×10-3 CKF -4 -3 -4 9.2×10 9.7×10 7.38×10 6.5×10-3 UKF -4 -3 -4 15×10 11×10 15×10 10.3×10-3 EKF


Page 36 of 38

Page 37 of 38

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


Table 9: Scenario II) RMSE computation of x3 over 100 independent Monte Carlo run for varying process noise and measurement noise levels. Q 10-6I4 10-4I4 -3 -2 -3 R 10 I3 10 I3 10 I3 10-2I3 -4 -3 -4 12×10 3.0×10 6.80×10 4.2×10-3 CKF -4 -3 -4 9.8×10 3.6×10 6.74×10 5.3×10-3 UKF -4 -3 -4 18×10 5.3×10 10×10 6.5×10-3 EKF

Table 10: Scenario II) RMSE computation of x4 over 100 independent Monte Carlo run for varying process noise and measurement noise levels. Q 10-6I4 10-4I4 -3 -2 -3 R 10 I3 10 I3 10 I3 10-2I3 -4 -3 -4 8.23×10 7.0×10 6.49×10 6.4×10-3 CKF -4 -3 -4 9.0×10 9.17×10 6.1×10-3 9.12×10 UKF 16×10-4 12×10-3 9.91×10-4 9.8×10-3 EKF

Table 11: The values of process parameters. Process Parameter description Value Parameters x1,i acid inlet concentration 3.2×10-3 mol/L x2,i base inlet concentration 5.0×10-2 mol/L x3,i buffer inlet concentration 2.5×10-3 mol/L Kx buffer dissociation constant 10-7 mol/L Kw water dissociation constant 10-14 mol2/L2 V reactor volume 1100 L

Table 12: RMSE computation of output. 0.3113 CKF κ= -2 0.7131 κ=-1 0.6205 κ=0 0.5159 UKF κ=1 0.6697 κ=2 1.0742

Table 13: RMSE computation of x1 over 100 independent Monte Carlo run for varying process noise and measurement noise levels. Q 0.001I3 0.01I3 R 0.01 0.1 0.01 0.1 0.0630 0.0631 0.0654 0.0659 CCKF 0.3079 0.3084 0.3095 0.307 CKF




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


Table 16: RMSE and the computation time of states over 100 independent Monte Carlo run for CCKF and MHE algorithms when initial guess of state vector is close to the true states. RMSE Computation time (second) x1 x2 x3 0.0752 0.0250 0.0236 0.321 CCKF 0.0417 0.0263 0.006 5.759 MHE

Table 17: RMSE and the computation time of states over 100 independent Monte Carlo run for CCKF and MHE algorithms when initial guess of state vector is far from the true states. RMSE Computation time (second) x1 x2 x3 0.1838 0.1170 0.1264 0.318 CCKF 0.2103 0.1296 0.1626 5.882 MHE


Page 38 of 38

Nonlinear and Constrained State Estimation Based on the Cubature

Recommend Documents