Optimization-Based Sigma Points Selection for ... - ACS Publications

Jan 7, 2013 - set of suitable optimization problems. Kadu et al.9 proposed a modified URNDDR algorithm that can handle arbitrary nonconvex constraints...
0 downloads 0 Views 909KB Size
Article pubs.acs.org/IECR

Optimization-Based Sigma Points Selection for Constrained State Estimation Sachin C. Kadu,†,‡ Mani Bhushan,*,† and Kallol Roy§ †

Department of Chemical Engineering, Indian Institute of Technology Bombay, Mumbai India 400076 Research Reactor Maintenance Division, BARC, Mumbai, India 400085

§

ABSTRACT: Various deterministic sampling-based constrained state estimation techniques for nonlinear dynamical systems have been proposed in the literature. These algorithms mainly focus on correcting the estimated (or filtered) states such that they satisfy all of the given constraints. However, the aspect of correcting predicted states (after model propagation) such that they satisfy all of the constraints has not been given much attention in the literature. These states are obtained by propagating deterministically chosen points (called sigma points) through the process model. Most approaches either do not incorporate constraints during the generation of these sigma points or use arbitrary strategies for constraint incorporation. Further, the weights associated with these points are not systematically updated after incorporation of any state constraint. In this article, we propose an optimization-based sigma points selection approach that overcomes these limitations. Additionally, we incorporate the effect of constraints by directly working with truncated probability density functions. The efficacy of the proposed approach is demonstrated by comparing its performance against those of various strategies available in the literature on several simulation case studies.

1. INTRODUCTION In many engineering applications encountered in fault diagnosis, chemical process simulation and control, and so on, the underlying systems have nonlinear dynamics. Additionally, the process states are often subject to some constraints that arise from various physical considerations. For example, in a chemical process reactor, the species concentrations, system temperature, and system pressure are non-negative. Hence, the problem of constrained state estimation for nonlinear dynamical stochastic systems has been the subject of considerable research in the literature. The most popular estimator used for state estimation of nonlinear dynamical systems is the extended Kalman filter (EKF). It is an extension of the traditional Kalman filter (KF) and involves linearization of the nonlinear terms. Recently, a new class of derivative-free Kalman filters, namely, the unscented Kalman filter (UKF), was proposed by Julier et al.1 as an alternative to direct model linearization. The UKF relies on model propagation of deterministically sampled points [called sigma points (SPs)] to estimate mean and covariance, instead of using a linearized model. These formulations were originally presented for unconstrained state estimation. Since then, the applicability of the basic UKF algorithm has been demonstrated by several researchers for a variety of applications.2−5 Over the past few years, several researchers have proposed variations of the basic UKF algorithm to incorporate constraints during the state estimation procedure. Li and Leung6 proposed a constrained UKF algorithm to localize vehicles. Further, Ma et al.7 applied dual perceptually constrained UKF for enhancing speech degraded by colored noise. In their approach, the Kalman gain is constrained to handle the state constraints. In all of these approaches, sigma points generated during the prediction step cannot handle any type of constraints. Vachhani et al.8 proposed unscented recursive nonlinear dynamic data reconciliation © 2013 American Chemical Society

(URNDDR), which replaces the update step in the UKF by a set of suitable optimization problems. Kadu et al.9 proposed a modified URNDDR algorithm that can handle arbitrary nonconvex constraints. Even though these algorithms can handle any kind of constraints during the state estimation step, the method used for sigma points generation [labeled as the intervalconstrained unscented transform (ICUT) by Teixeira et al.10] incorporates only upper and lower bounds on the states by simply clipping them at their respective bounds. The mean (which is one of the sigma points) and covariance of the underlying probability distribution are not explicitly updated after this clipping. In addition, the weight updating procedure is arbitrary. Kandepu et al.11 proposed a simple procedure of projecting transformed infeasible sigma points onto the boundary of the feasible region to include state constraints in the UKF algorithm. However, the weights of the projected sigma points are not updated to reflect their asymmetric location after projection. The application of the manifold-constrained UKF for aerial vehicles was demonstrated by Sipos,12 who dealt with state constraints by using modified UT and a modified time update model. However, similarly to ref 11, in this case, the weights associated with sigma points are also not updated at each time instance to represent their asymmetric nature. Kolas et al.13 investigated several variants of the UKF algorithm for constraint handling for cases where the EKF leads to poor performance. In their approaches also, the weights are not updated to account for the asymmetric location of sigma points. Recently, Teixeira et al.10 investigated various forms of UKF algorithm for handling state interval constraints based on combinations of various Received: Revised: Accepted: Published: 1916

November 25, 2011 July 31, 2012 January 7, 2013 January 7, 2013 dx.doi.org/10.1021/ie2027487 | Ind. Eng. Chem. Res. 2013, 52, 1916−1926

Industrial & Engineering Chemistry Research

Article

methods available for prediction and filtering (or estimation) steps. The various forms proposed by Teixeira et al.,10 such as projected interval UKF (PIUKF) and truncated interval UKF (TIUKF), used ICUT for handling the constraints during the prediction step. The ICUT has certain limitations as discussed earlier. In the literature, several other constrained estimation approaches based on stochastic sampling were proposed such as that by Prakash et al.14 In this article, we do not review such approaches because we focus our attention on constrainthandling algorithms developed under a deterministic sampling framework. Similarly, articles related to constraint handling within the traditional EKF approach have also not been surveyed. The limitations of various UKF-based constraint-handling algorithms have already been discussed. To summarize, most of the algorithms have the following drawbacks: (i) correction of only filtered (or estimated) states but not the predicted (modelpropagated) states, (ii) arbitrary strategies for weight updating (after incorporating constraints) during the sigma points selection step, and (iii) no reflection of the effect of constraints on the underlying probability distribution (usually assumed to be Gaussian) of the estimated states. In the traditional Kalmanfilter-based constraint-handling area, Simon15 presented a probability density function (PDF) truncation method. The PDF truncation approach is used to truncate the mean and covariance of estimated states to account for state constraints. Further, Simon and Simon16 formulated constrained Kalman filtering by density function truncation for the turbofan engine health estimation problem. They showed that the truncated PDF approach provides a more accurate way of incorporating inequality constraints than other traditional approaches such as the projection approach to constrained filtering. They also pointed out that, in terms of root-mean-square (RMS) estimation error, the PDF truncation-based approach performs well compared to projection-based approach. In this article, we propose a constrained sigma points truncated UKF (CSPTUKF) algorithm that has the following features: It is based on deterministic sampling ideas as in UKF, but it borrows the idea of the PDF truncation approach from the KF literature to reflect the effect of constraints on the PDF of the estimated states. It also ensures that the sigma points and their corresponding weights can be systematically generated during the prediction step in the presence of any type of constraints. Further, the number of sigma points used in this algorithm is less than the number used in UKF. This article is structured as follows: In section 2, we first briefly summarize the basic UKF algorithm (for unconstrained systems) and also summarize various UKF-based constrained estimation approaches available in the literature. In section 3, we present our proposed CSPTUKF approach, which overcomes several of the drawbacks of the algorithms available in the literature. The utility of CSPTUKF is discussed in section 4 in terms of its application to an isothermal batch reactor,14 two interacting tanks in series,17 and a polymerization reactor18 simulation. In these case studies, the performance is compared against that of several techniques available in the literature (which were summarized in section 2). Finally, concluding remarks are presented in section 5.

2.1. UKF. The UKF was proposed by Julier et al.1 The fundamental component of the UKF is the unscented transformation (UT), which uses a set of deterministically chosen weighted sigma points to parametrize the mean and covariance of a probability distribution. Consider the following discrete-time nonlinear state-space model of the process and measurements with additive noises

xk + 1 = f (xk , uk) + γk

(1)

yk = h(xk , uk) + νk

(2)

where k represents the sampling instant, xk ∈ Rn represents the vector of states, uk ∈ Rm is the vector of known manipulated input variables, and yk ∈ Rq is the vector of measured variables. The process noise γk represents the unmodelled process dynamics in eq 1 and is assumed to be white and Gaussian with zero mean and covariance matrix Qk. The measurements are corrupted by noise vk, which is also assumed to be white and Gaussian with zero mean and covariance matrix Rk. Further, vk and γk are assumed to be uncorrelated. Functions f and h represent the process and the measurement model, respectively. For given estimated states x̂k(+) and associated covariance Pk(+), the UKF can be summarized as consisting of the following steps. 2.1.1. Prediction or Forecast Step. A set of 2n + 1 sigma points χi,k(+) is chosen symmetrically about x̂k(+) as follows ⎧ xk̂ ( +) i=0 ⎪ ⎪ x ̂ ( +) + (n + κ ) i = 1, ..., n ⎪ k χi , k ( +) = ⎨ [ Pk( +) ]i ⎪ ⎪ xk̂ ( +) − (n + κ ) i = n + 1, ..., 2n ⎪ ⎩ [ Pk( +) ]i − n

(3)

1/2

The term [(Pk(+)) ]i is the ith column of the matrix square root of Pk(+); that is, [(Pk(+))1/2] is a matrix such that [(Pk(+))1/2][(Pk(+))1/2]T = Pk(+). The associated weights for the given sigma points are specified as w0 =

κ , (n + κ )

wi =

1 2(n + κ )

i = 1, ..., 2n

(4)

where κ is a tuning parameter. The transformed set of sigma points is obtained by propagating each sigma point through the process model χi̅ , k + 1 ( −) = f (χi , k ( +), uk)

i = 0, 1, ..., 2n

(5)

The predicted mean and the state prediction error covariance matrix (SPECM) are then computed as 2n

xk̂ + 1( −) =

∑ wiχi̅ ,k+ 1 (−) i=0

(6)

2n

Pk + 1( −) =

∑ wi[χi̅ ,k+ 1 (−) − xk̂ + 1(−)] i=0

[χi̅ , k + 1 ( −) − xk̂ + 1( −)]T + Q k

2. UKF-BASED ALGORITHMS FOR CONSTRAINED STATE ESTIMATION In this section, we first briefly summarize the basic UKF approach for unconstrained state estimation. Then, we also summarize the various UKF-based constrained state estimation algorithms available in the literature.

(7)

To account for the effect of Qk, sigma points are regenerated from the computed x̂k+1(−) and Pk+1(−) values. The new predicted sigma points χi,k+1(−) are computed in a manner similar to that given by eq 3. These new predicted sigma points are used in the next step, which is the estimation step. The weights (wi) 1917

dx.doi.org/10.1021/ie2027487 | Ind. Eng. Chem. Res. 2013, 52, 1916−1926

Industrial & Engineering Chemistry Research

Article

respective bounds. The mean (which is the 0th sigma point) and covariance of the underlying probability distribution are not explicitly updated after this clipping. In addition, the weight updating procedure is also arbitrary. 2.3. Truncated UKF (TUKF). TUKF is a three-step algorithm consisting of prediction, estimation, and truncation steps. It is obtained by appending the probability density function (PDF) truncation procedure to the UKF algorithm.10 In the TUKF, the prediction step is given by eqs 3−7, and the estimation step is given by eqs 8−14 together with the PDF truncation step as described in ref 15. In the PDF truncation procedure, the probability density function (which is assumed to be Gaussian) computed during the estimation step is truncated to account for linear constraints on the states. In the TUKF, use of the PDF truncation approach enables direct incorporation of the effect of constraints on the parameters (mean, covariance) of the underlying probability distribution. Even though the TUKF algorithm can handle constraints during the state estimation step, the method used for sigma points generation (as given by eq 3) does not handle any type of the constraints on the state variables. Further, the weights associated with sigma points are not updated at any time step. The inability of the TUKF to handle any type of constraints during the sigma points generation step is overcome by the truncated interval UKF, which is an extension of the TUKF for handling intervaltype constraints during the state prediction step. It is highlighted in the next section. 2.4. Truncated Interval UKF. The truncated interval UKF (TIUKF) proposed by Teixeira et al.10 is also a three-step algorithm similar to the TUKF. The major difference between the TUKF and the TIUKF is in the sigma points generation step. In the TUKF, sigma points and associated weights are generated using the UT procedure as given by eqs 3 and 4. On the other hand, the TIUKF uses the ICUT procedure as in the URNDDR approach. The rest of the procedure is the same for both the TUKF and the TIUKF. Even though the TIUKF algorithm can handle constraints during the state estimation step, the method used for sigma points selection (ICUT) incorporates only upper and lower bounds on the states by simply clipping them at their respective bounds. In addition, the weight updating procedure is arbitrary. Next, we discuss another variant of the TUKF in which sigma points are generated using the minimal skew simplex points (MSSP) procedure as proposed by Julier and Uhlmann.21 This algorithm is not discussed in the literature, and we refer to it as the MSSPTUKF. 2.5. Minimal Skew Simplex Points TUKF (MSSPTUKF). The MSSPTUKF is a variation of the TUKF, where eqs 3 and 4 used for sigma points and their weight generation are replaced by the MSSP set as proposed by Julier and Uhlmann.21 The key idea in MSSP set is to capture the mean and covariance of given distribution with minimum number of sigma points. The MSSP are asymmetric in nature and have the property that they minimize the third-order moment. The MSSP set is generated in two steps as described next.21 2.5.1. Generation of MSSP Set from a Distribution Having Zero Mean and Identity Covariance Matrix. An MSSP set having zero mean and identity covariance matrix is generated using the algorithm given in ref 21. Let the corresponding sigma points be labeled χi, i = 0, 1, ..., n + 1. These minimal skew sigma points are asymmetrically located around the origin. They form a polytope in higher-dimension space in which all the sigma points are not equidistant from the origin.

associated with the sigma points remain the same at all time instances. 2.1.2. Estimation or Data Assimilation Step. The new predicted sigma points are propagated through the measurement model to predict the measurement sigma points Υi,k+1(−) and their mean ŷk+1(−) Υi , k + 1( −) = h(χi , k + 1 ( −), uk + 1)

i = 0, 1, ..., 2n

(8)

2n

yk̂ + 1 ( −) =

∑ wiΥi ,k+ 1(−) i=0

(9)

The covariance in innovations [Pee,k+1(−)] and the cross covariance between the state prediction error and the innovations [Pse,k+1(−)] are computed as 2n

Pee, k + 1( −) =

∑ wi[Υi ,k+ 1(−) − yk̂ + 1 (−)] i=0

[Υi , k + 1( −) − yk̂ + 1 ( −)]T + R k + 1

(10)

2n

Pse, k + 1( −) =

∑ wi[χi ,k+ 1 (−) − xk̂ + 1(−)] i=0

[Υi , k + 1( −) − yk̂ + 1 ( −)]T

(11)

The Kalman filter gain Kk+1, the updated mean state vector x̂k+1(+), and the state estimation error covariance matrix Pk+1(+) (SEECM), are then computed as Kk + 1 = Pse, k + 1( −)[Pee, k + 1( −)]−1

(12)

xk̂ + 1( +) = xk̂ + 1( −) + Kk + 1[yk + 1 − yk̂ + 1 ( −)]

(13)

Pk + 1( +) = Pk + 1( −) − Kk + 1Pee, k + 1( −)KkT+ 1

(14)

Here, + and − indicate the times immediately after and before the measurement at any time instance, respectively. The basic UKF algorithm discussed so far does not handle any type of constraints on the states. In the literature, various extensions of the UKF are available to handle state constraints.19 Next, we summarize some relevant constrained UKF formulations, which are later used for comparison with the performance of our proposed formulation. 2.2. Unscented Recursive Nonlinear Dynamic Data Reconciliation. In the unscented recursive nonlinear dynamic data reconciliation (URNDDR) formulation, sigma points and associated weights (generated by the UKF formulation) are updated to account for bounds on the state variables during the prediction step. Toward this end, sigma points and weights are first generated according to eqs 3 and 4. In case any of these sigma points violate a bound constraint on any state, a backstepping procedure is used to ensure that the bound is satisfied. The weights of the sigma points are then also modified to ensure that the sum of the weights adds up to 1. This procedure of generating sigma points and weights was labeled as the interval-constrained unscented transform (ICUT) by Teixeira et al.10 The estimation step of the UKF (eqs 10−14) is replaced by a set of optimization problems (one for each sigma point) that can incorporate convex constraints on the states. Details of the URNDDR approach are presented in refs 8, 20. Even though the URNDDR algorithm can handle convex constraints during the state estimation step, the method used for sigma points selection (ICUT) incorporates only upper and lower bounds on the states by simply clipping them at their 1918

dx.doi.org/10.1021/ie2027487 | Ind. Eng. Chem. Res. 2013, 52, 1916−1926

Industrial & Engineering Chemistry Research

Article

2.5.2. Transformation of MSSP Set to Capture Mean and Covariance of a Given Distribution. The following transformation is used to obtain the sigma points χi,k(+) having mean x̂k(+) and covariance Pk(+) from the sigma points χi belonging to the above generated MSSP set χi , k ( +) = xk̂ ( +) +

Pk( +) χi

For a given estimated state x̂k(+) and associated covariance Pk(+), the CSPTUKF can be summarized as consisting of the following steps. 3.1. Truncation Step. In this step, the mean x̃k(+) and the covariance P̃ k(+) of the truncated PDF are computed from the original (untruncated) Gaussian distribution N[x̂k(+),Pk(+)] after incorporation of the s constraints given in eq 16. The truncation algorithm for the constraints as given in eq 16 is described in ref 15. The truncated mean x̃k(+) and covariance P̃k(+) are used in prediction step for the next time instant. 3.2. Prediction or Forecast Step. In this step, the set of sigma points and their associated weights are computed by solving the optimization problem

i = 0, 1, ..., n + 1 (15)

Two variants of this algorithm are possible by tuning the weight associated with the mean (w0). w0 = 0 results in an MSSP set with n + 1 sigma points; otherwise, n + 2 sigma points are obtained. Apart from the sigma points generation step, the rest of the procedure for the MSSPTUKF is the same as that for the TUKF. The advantage of the MSSPTUKF lies in the use of fewer sigma points compared to the TUKF. However, it has drawbacks similar to those of the TUKF, namely, it cannot handle any type of constraints during the sigma points generation step.

min

χi , k (+), wi , k(+)

zkL, j ≤ ϕkT, jχi , k ( +) ≤ zkU, j 0 < wi , k( +) < 1

j = 1, ..., s ; i = 1, ..., Np

i = 1, ..., Np

Np

∑ wi ,k(+) = 1 i=1 Np

∑ wi ,k(+)χi ,k (+) = xk̃ (+) i=1 Np

∑ wi ,k(+)[χi ,k (+) − xk̃ (+)][χi ,k (+) − xk̃ (+)]T

= Pk̃ ( +)

i=1

(18)

In this proposed approach, at each kth time instant, a single optimization problem with Np × (n + 1) decision variables is solved. The decision variables in the given optimization problem are the Np updated sigma points vectors χi,k(+) and their associated weights wi,k(+). Here, Np = n + 1 is the minimum number of sigma points required to capture the mean and covariance of the given distribution.22 This formulation allows the selection of individual sigma points such that the sample mean of the sigma points is the truncated mean and the sample covariance is the truncated covariance. Also, unlike other approaches in the literature, the weights associated with each sigma point are not fixed a priori and are instead allowed to vary so as to satisfy the given constraints. As a result, weights associated with different sigma points can be different. This is especially relevant when the sigma points are located asymmetrically and, hence, the weights of different points need not be the same. The transformed set of sigma points is obtained by propagating each of the generated Np sigma points through process model as

Figure 1. Steps of CSPTUKF.

Consider the state description of the process as given by eqs 1 and 2 together with the following inequality constraints j = 1, 2, ..., s

(17)

subject to the constraints

3. PROPOSED ESTIMATION ALGORITHM: CSPTUKF In the previous section, we briefly discussed some existing approaches in the literature. As mentioned earlier, these approaches have the drawbacks of not incorporating constraint information during computation of the parameters of the probability density function or updating of the weights, as well as not ensuring that the predicted state estimates also satisfy the state constraints. In this section, we present our proposed optimization-based truncated UKF approach for constrained state estimation for nonlinear dynamical systems that overcomes these limitations. The approach is labeled as the constrained sigma points truncated unscented Kalman filter (CSPTUKF). It is a three-step algorithm consisting of truncation, prediction, and estimation steps and is summarized in Figure 1. These steps are next described in more detail.

zkL, j ≤ ϕkT, jxk ≤ zkU, j

F=0

χi̅ , k + 1 ( −) = f (χi , k ( +), uk)

(16)

where ϕTk,jxk represents the jth linear combination of the states on which inequality constraints (lower and upper bounds) have been placed and s is the total number of such constraints. zLk,j and zUk,j are the lower and upper bounds, respectively, for the jth inequality constraint at the kth time instant.

i = 1, ..., Np

(19)

The predicted mean and SPECM are then computed as Np

xk̂ + 1( −) =

∑ wi ,k(+)χi̅ ,k+ 1 (−) i=1

1919

(20)

dx.doi.org/10.1021/ie2027487 | Ind. Eng. Chem. Res. 2013, 52, 1916−1926

Industrial & Engineering Chemistry Research

Article

Np

Pk + 1( −) =

It is also to be noted that, in the absence of constraints, the CSPTUKF reduces to a two-step algorithm consisting of prediction and estimation steps. Further, for linear dynamical systems, it is equivalent to the KF. This is shown in Appendix A. It is also to be noted that, although we consider linear inequality state constraints in this work, our approach can be extended to handle nonlinear constraints as well. The optimization problems for sigma points generation can directly handle the nonlinear constraints, whereas linearized approximations of the nonlinear constraints can be used during PDF truncation. We now briefly discuss issues related to the choice of objective function in the optimization problem used to generate the sigma points and the computational effort required for the proposed approach. (i) In the implementation in our work, the objective function in formulation 17 is empty; that is, eq 17 represents a feasibility problem. However, different objectives can be used. These can be related to higher-order moments, such as minimizing the difference in some of the skew (third-order moment) components as computed by the sigma points and the skew components of the truncated probability density function. This would require computation of the skew of the truncated probability density function and a larger value of Np. The objective can also be based on any other heuristic or logic, for example, maximizing the minimum distance between the sigma points pairwise, minimizing the maximum difference between the weight of sigma points, and so on. In this work, we have not explored these aspects further. (ii) The proposed CSPTUKF formulation requires that a feasibility problem be solved at every time instant to generate the sigma points and associated weights. This feasibility problem involves bilinear and trilinear constraints even for linear constraints on states. However, it is expected that it will still be comparable to other single-step optimization-based constrained state estimation approaches such as the URNDDR, which involves solving quadratic programming (QP) problems for linear constraints and a linear measurement map, because the CSPTUKF involves solving only feasibility problems. At each time instance, in the CSPTUKF, two feasibility problems, each with (n + 1)2 decision variables, are solved. On the other hand, in the URNDDR, 2(n + q) + 1 optimization problems, each having n decision variables, are solved separately. Thus, for smaller n, the computational effort for the CSPTUKF is expected to be lower than that for the URNDDR, whereas for larger n, the reverse might hold. The exact performance, however, will be problemspecific. Further, in the literature, use of 2n + 1 sigma points is popular for the UKF because sets of sigma points that are symmetric around the mean can then be generated. In our work, because we are dealing with constrained state estimation problems and because constraints are unlikely to be symmetric around the mean, we chose to work with n + 1 sigma points. However, a user can also use more sigma points with our approach. Using more points should, in general, lead to better state estimation but at the cost of additional computational burden mainly involved with the optimization problems. In the next section, we compare our proposed formulation, namely, the CSPTUKF, with several algorithms from the literature presented earlier.

∑ wi ,k(+)[χi̅ ,k+ 1 (−) − xk̂ + 1(−)] i=1

[χi̅ , k + 1 ( −) − xk̂ + 1( −)]T + Q k

(21)

To account for the effect of the addition of process noise covariance Qk, the sigma points are regenerated from a Gaussian distribution with mean x̂k+1(−) and covariance Pk+1(−). Toward this end, the optimization formulation given by eqs 17 and 18 is solved with mean and covariance of x̂k+1(−) and Pk+1(−), respectively. Let these new predicted sigma points be labeled χi,k+1(−) and their weights be labeled wi,k+1(−). These sigma points and weights are used in the next (estimation) step. 3.3. Estimation or Data Assimilation Step. The new predicted sigma points are propagated through the measurement model to predict the measurement sigma points Υi , k + 1( −) = h(χi , k + 1 ( −), uk + 1)

i = 1, ..., Np

(22)

The predicted observation is calculated as Np

yk̂ + 1 ( −) =

∑ wi ,k+ 1(−)Υi ,k+ 1(−) i=1

(23)

The covariance in innovations and the cross covariance between the state prediction error and the innovations are calculated as Np

Pee, k + 1( −) =

∑ wi ,k+ 1(−)[Υi ,k+ 1(−) − yk̂ + 1 (−)] i=1

[Υi , k + 1( −) − yk̂ + 1 ( −)]T + R k + 1

(24)

Np

Pse, k + 1( −) =

∑ wi ,k+ 1(−)[χi ,k+ 1 (−) − xk̂ + 1(−)] i=1

[Υi , k + 1( −) − yk̂ + 1 ( −)]T

(25)

Then, the Kalman filter gain is calculated as Kk + 1 = Pse, k + 1( −)[Pee, k + 1( −)]−1

(26)

The updated mean of the state vector and SEECM are computed as xk̂ + 1( +) = xk̂ + 1( −) + Kk + 1[yk + 1 − yk̂ + 1 ( −)]

(27)

Pk + 1( +) = Pk + 1( −) − Kk + 1Pee, k + 1( −)Kk + 1T

(28)

This completes one cycle: going from time instant k to time instant k + 1. This procedure is then repeated for the next time step (k = k + 1). The proposed CSPTUKF ensures that the individual sigma points generated during the prediction step satisfy the given constraints. It overcomes the limitation of the ICUT algorithm used to handle interval constraints during the prediction step in various constrained state estimation algorithms such as the URNDDR, the MURNDDR, and various forms proposed in ref 10. In addition, the updating of weights during each time instance is not arbitrary. They are solutions of an optimization problem. The asymmetric location of the sigma points can result in the assignment of different weights to the points at each time instant. These features, namely, the ability to handle any type of constraints during sigma points generation and the systematic weight updating during the prediction step, are missing from all of the constraint-handling algorithms proposed in the literature.

4. SIMULATION RESULTS The performance of the proposed CSPTUKF is compared with that of other approaches discussed in this article, namely, the TUKF, TIUKF, MSSPTUKF, and URNDDR, on the following 1920

dx.doi.org/10.1021/ie2027487 | Ind. Eng. Chem. Res. 2013, 52, 1916−1926

Industrial & Engineering Chemistry Research

Article

problems: (i) an isothermal batch reactor,14 (ii) two interacting tanks in series,17 and (iii) a polymerization reactor.18 The sum of square root estimation error (SSREE) averaged over M different simulation trials is used for this comparison. The average SSREE (ASSREE) is defined as ASSREE =

1 M

M

starting from feasible state values were set to their respective bounds during model propagation. At any time instance, the proposed CSPTUKF approach uses four sigma points, whereas the TUKF and TIUKF use seven and the URNDDR uses nine. The number of sigma points used in the MSSPTUKF is either four or five depending on the value of w0. The states as estimated by the CSPTUKF, MSSPTUKF (w0 = 0), TUKF, TIUKF, and URNDDR approaches are presented in Figure 2. It can be seen from this figure that the states estimated

N

∑∑

epT, kep , k (29)

p=1 k=1

where ep,k is the vector of errors between the true states and the estimated states at time instance k for the pth simulation trial. N is the total number of time instances in each simulation trial, and M is the total number of simulation trials. All of the computations were performed in MATLAB version 7.4 running in a Linux environment on a P4 Intel Core 2 Quad processor on a 4 GB RAM machine. TOMLAB version 5.9 was used to solve optimization problems whenever required. 4.1. Case Study I: Isothermal Batch Reactor.14 Consider the following gas-phase reversible reactions taking place in an isothermal batch reactor kc1

A HooI B + C,

kc3

2B HooI C

kc 2

(30)

kc4

The governing equations for this system are x ̇ = ν Tr ,

y = [32.84 32.84 32.84]x

(31)

where ν and r are the stoichiometric and reaction-rate matrices, respectively, and are given by ⎡−1 1 1⎤ ν=⎢ , ⎣ 0 −2 1⎥⎦

⎡ kc1x1 − kc 2x 2x3 ⎤ ⎥ r=⎢ ⎢⎣ kc3x 2 2 − kc4x3 ⎥⎦

Figure 2. Concentration tracking by various state estimation algorithms for (a) reactant A, (b) product B, and (c) product C.

by all of the approaches satisfied the upper and lower bounds. This is because all of these approaches can handle bound-type constraints during the state estimation step. However, the same is not true for the state prediction step. The percentage constraint violations during sigma points generation for the various algorithms are presented in Table 2. This percentage was

(32)

Further, kci is the reaction rate constant for the ith reaction, and the set of rate constants is given by [0.5 0.05 0.2 0.01]T. The vector of states (x) consists of the concentrations of components A, B, and C denoted as x1, x2, and x3, respectively. The measurement is of the scaled sum of all concentrations. The various parameters used for the state estimation problem are summarized in Table 1 in consistent units. The lower bounds

Table 2. Comparison of the Performance of Various State Estimation Algorithms for an Isothermal Batch Reactor percentage constraint violation during sigma points generation

Table 1. Parameters for the State Estimation Problem of an Isothermal Batch Reactor symbol

description

value

x0

initial true states

Δt N M Q R x̂0(+)

integration time number of time instances number of simulation trials process noise covariance matrix measurement noise variance initial guess for state estimation

[0 0 1]T

P0(+)

initial SEECM

0.25I3×3

algorithm

[0.5 0.5 0]T

TUKF TIUKF MSSPTUKF (w0 = 0) MSSPTUKF (w0 = 0.00005) CSPTUKF URNDDR

0.1 200 10 10−6I3×3 0.0625

before model propagation

after model propagation

ASSREE

avg computation time (s) per time instance

0.2143 0 0.2750

0 0 0

0.6706 0.8228 0.6948

0.0292 0.0298 0.0198

0.22

0

0.6950

0.023

0 0.1667

0 0

0.6534 1.4258

0.0821 0.2148

calculated as the ratio of the number of sigma points violating the constraints to the total number of sigma points generated during simulation trials, namely, N × M × Np. Constraint violation for the TUKF and MSSPTUKF occurs because there is no mechanism in these methods to incorporate any constraints during sigma points generation. The ICUT procedure used in the TIUKF and URNDDR is capable of handling interval-type constraints as used in this case study. Howvere, it is to be noted that the ICUT procedure used for generating sigma points from a

(LBs) and upper bounds (UBs) for the concentrations of A−C are taken as [0.0129, 0.4758], [0.0743, 0.3575], [0.0243, 0.6601], respectively. The values of the LBs and UBs were chosen such that the true states satisfied these constraints. The MATLAB inbuilt function ode45 with a step size of 0.1 (which is the same as the sampling time) was used to integrate the process equations to obtain true state data. In addition, to keep the model realistic, any infeasible states resulting from integration of the state equations 1921

dx.doi.org/10.1021/ie2027487 | Ind. Eng. Chem. Res. 2013, 52, 1916−1926

Industrial & Engineering Chemistry Research

Article

given distribution {i.e., N[x̂k(+),Pk(+)]} assumes that the mean x̂k(+) of the distribution satisfies the given interval constraints.10 If this assumption does not hold, then the ICUT procedure can result in the generation of sigma points that violate the constraints. It can be seen from Table 2 that the TIUKF does not violate constraints during the sigma points generation step because x̃k(+) used during the ICUT satisfies the bound constraints at all times as a result of the PDF truncation step. However, the same is not true for the URNDDR. For the URNDDR, the percentage constraint violation can be nonzero before the model propagation step. This is due to the fact that, at the first time instance, the initial guess x̂0(+) used for simulation might not satisfy the constraints. Hence, constraint violation will be observed only during the first time instance if x̂0(+) does not satisfy the constraints. In the present case study, x̂0(+) does not satisfy the bound constraints. Hence, the percentage constraint violation is nonzero for the URNDDR during sigma points generation before model propagation. It should be noted that, ideally, infeasible sigma points should not be propagated through the model, as doing so can lead to problems during model propagation and the particular state estimation algorithm should be terminated. However, for the sake of comparison, even when a sigma point was infeasible, we propagated it through the process model. For the current case study, doing so did not create any problems during model propagation. Another aspect to note is that, even though the ICUT procedure used in the TIUKF and URNDDR is able to handle interval constraints during the prediction step, its weight updating procedure is arbitrary. On the other hand, the proposed CSPTUKF algorithm can handle any type of constraint during the state prediction step, and the weights are part of the solution of an optimization problem. The asymmetric nature of the sigma points results in the assignment of unequal weightage to individual sigma points. Further, these weights are updated at each time instance. The variation of the weights of individual sigma points (generated by the CSPTUKF) with time is shown in Figure 3. This variation of weights is shown for the sigma points generated before the model propagation step. A similar variation was observed (not shown here) for weights associated with sigma points generated after the model propagation step. The performances of the various approaches (ASSREE and computational effort) are also listed in Table 2. It can be seen that

the CSPTUKF performed better than the other approaches with acceptable computational efforts. 4.2. Case Study II: Two Interacting Tanks in Series.17 In the previous case study, only lower and upper bounds on states were imposed. We now present a case study involving two interacting tanks in series, with bounds on linear combinations of states as well. The process schematic is shown in Figure 4. In this

Figure 4. Two interacting tanks in series.

process, water enters through the first tank and leaves through the second tank. The water level in each tank depends on the water levels in both tanks; hence, the system is interacting in nature. The water level in each tank is considered as a state, and the water level in tank 1 is assumed to be measured. A nonlinear mathematical model is given as dh1 F c = i − 1 h1 − h2 dt A1 A1

(33)

dh 2 c c = 1 h1 − h2 − 2 h2 dt A2 A2

In the process, A1 = 28 cm2 and A2 = 32 cm2 are the crosssectional areas of tanks 1 and 2, respectively. c1 = 3.1449 cm2.5/s is the flow resistance coefficient between tanks 1 and 2, and c2 = 2.2590 cm2.5/s is the flow resistance coefficient at the exit of tank 2. The inlet flow rate Fi was chosen to be 5 cm3/s. The values of some of these parameters were taken from ref 23, in which a widely studied quadruple-tank process was presented. The various parameters used for the state estimation problem are summarized in Table 3. Lower and upper bounds on the two Table 3. Parameters for the State Estimation Problem of Two Interacting Tanks in Series symbol

description

value

x0

initial true states (cm)

[2.0 0.9]T

Δt N M Q R x̂0(+)

integration time (s) number of time instances number of simulation trials process noise covariance matrix (cm2) measurement noise variance (cm2) initial guess for state estimation (cm)

0.1 100 10 0.1I2×2 0.0625

P0(+)

2

initial SEECM (cm )

[1.5 0.5]T 35I2×2

states were taken to be [0, 4.5] (state 1) and [0, 3] (state 2). Additionally, lower and upper bounds were also imposed on h1 − h2 as 0 and 4.5, respectively. These constraints were chosen such that they were satisfied by the true states. The MATLAB inbuilt function ode45 with a step size of 0.1 s (which is the same as the sampling time) was used to integrate the process equations to obtain true state data. The states as estimated by the CSPTUKF and MSSPTUKF (w0 = 0) are shown in Figure 5. The corresponding results for the TUKF, TIUKF, and URNDDR are not shown because the sigma

Figure 3. Variation of weights [wi,k(+)] associated with individual sigma points. 1922

dx.doi.org/10.1021/ie2027487 | Ind. Eng. Chem. Res. 2013, 52, 1916−1926

Industrial & Engineering Chemistry Research

Article

of the proposed algorithm on a relatively larger problem. This case study involved styrene polymerization in a continuous stirred-tank reactor (CSTR). This is a highly nonlinear system having six states, namely, Ci, Cs, Cm, λ1, Mn, and Mw, and four measurements, namely, Cs, Mn, Mw, and Cm. The differential equations for the concentrations of initiator (Ci), solvent (Cs), and monomer (Cm) are FC ⎛F ⎞ dC i i ii = − ⎜ i + k i ⎟C i + ⎝V ⎠ dt V FC dCs FC i s i + FmCsm =− t s + dt V V FmCm m − FC dCm t m = −k pCmP + dt V

The first moment (λ1) of the polymer molecular weight distribution, the number-average molecular weight (Mn), and the weight-average molecular weight (Mw) of the dead polymer chains in the CSTR are given by

Figure 5. Water-level tracking by various state estimation algorithms for (a) tank 1 and (b) tank 2.

Fλ dλ1 PM m =− t 1 +β dt V (1 − α)

points generation step used in the TUKF was not able to handle any type of constraints, whereas the ICUT procedure used in the TIUKF and URNDDR is not able to handle the constraint 0 ≤ h1 − h2 ≤ 4.5. It was observed during the simulations that, even though individual sigma points satisfied the bound constraints on h1 and h2, the constraint h1 − h2 ≥ 0 was violated by the TUKF, TIUKF, and URNDDR. During the prediction step, when we propagated such infeasible sigma points through the model, we obtained heights that were complex numbers because of the (h1 −h2)1/2 term in the model equations. Hence, it was not possible to execute these algorithms without any modification during model propagation. Further, even though the MSSPTUKF also cannot handle any type of constraints (as demonstrated in the previous case study by showing nonzero constraint violation during the sigma points generation step), in the present case study, constraint violation by the MSSPTUKF was not observed. This occurred just by chance. In contrast to these approaches, our proposed formulation CSPTUKF can handle all of these constraints during sigma points generation step. The corresponding results for the case study are shown in Figure 5. The performances of the CSPTUKF and MSSPTUKF were also compared by calculating ASSREE and the average time required per time instance. The results are listed in Table 4. For this case study, their performances were comparable.

PM n dM n = (β M m − η) λ1(1 − α) dt dM w PM m = [ψ − βM w (1 − α)] dt λ1(1 − α)2

MSSPTUKF (w0 = 0) MSSPTUKF (w0 = 0.05) CSPTUKF

ASSREE (cm)

avg computation time (s) per time instance

2.5375 2.5375

0.01625 0.01942

2.5381

0.09536

(35)

where P = (2f Ciki/kt)1/2 is the concentration of live polymer and Ft = Fi + Fm. The expressions for α, β, η, and ψ are as follows α=

k pCm (k p + k fm)Cm + k fsCs + k tP

β = (k fmCm + k tdP + k fsCs)(2α − 2α 2) + k tcP

η = [(k fmCm + k tdP + k fsCs)α + 0.5k tcP]M n(1 − α) ψ = [(k fmCm + k tdP + k fsCs)(α 3 − 3α 2 + 4α) + k tcP(α + 2)]M m

The rate constants used in these equations are given by ki =

Table 4. Comparison of the Performance of Various State Estimation Algorithms for Two Interacting Tanks in Series algorithm

(34)

0.693 60 × 10[(A i / T ) + Bi ]

k p = A p exp( −Ep/RT )

k fm = A fm exp( −Efm /RT ) k fs = A fs exp( −Efs /RT )

k t = A t exp( −Et /RT )

k td = 0.15k t

This case study demonstrates the utility of our proposed formulation in the presence of constraints on linear combinations of state variables. It avoids the generation of infeasible sigma points, which is generally not true for several other approaches available in the literature. 4.3. Case Study III: Polymerization Reactor.18 The objective of the third case study was to demonstrate the efficacy

k tc = 0.85k t The steady-state operating conditions and the relevant parameters for simulation were taken from ref 18. The various parameters used for the state estimation problem are summarized in Table 5 and have consistent SI units as given in ref 18. The LBs and UBs for the states were taken as [0 0 0 0.01 1923

dx.doi.org/10.1021/ie2027487 | Ind. Eng. Chem. Res. 2013, 52, 1916−1926

Industrial & Engineering Chemistry Research

Article

Table 5. Parameters for the State Estimation Problem of a Polymerization Reactor symbol

description

value

x0

initial true states

[1.9854 × 10

Δt N M Q

integration time number of time instances number of simulation trials process noise covariance matrix

1 14400 5

diag[ 6 × 10−12 1.5 × 10−8 1.5 × 10−8 1.575 × 10−4 0.005 0.01]

R

measurement noise covariance matrix

diag[0.003 100 200 0.003]

x̂0(+)

initial guess for state estimation

[1.9854 × 10−3 5.475 1.9408 1.7 × 102 4.5 × 103 7 × 103 ]T

P0(+)

initial SEECM

diag[1.25 5 5 250 5000 5000]

−3

100 100]T and [1 15 10 500 25 000 40 000]T, respectively.18 The true states were generated using the Euler−Maruyama scheme,24 and the change in the value of the initiator flow rate (Fi) was assumed to be step change from 1.5414 × 10−4 to 5.4 × 10−5 kmol/s at the 3600th time instance. The states as estimated by the CSPTUKF are shown in Figure 6. It can be seen from the

sigma points through the model, we obtained states with complex numbers because of the square-root term in variable P (concentration of live polymer). Hence, it was not possible to execute these algorithms without any ad hoc modification during model propagation. Our proposed formulation, on the other hand, does not require any ad hoc modification in the presence of constraints because it avoids the generation of infeasible sigma points.

5. CONCLUSIONS Accurate state estimation is essential for several activities related to control, monitoring, and optimization in any process. Apart from ensuring low uncertainty, it is important to ensure that the states obtained at various steps of the estimation algorithm also satisfy the constraints on the states arising from the physics of the underlying process. In this work, an optimization-based sigma points generation approach is presented for constrained state estimation problems. In contrast to existing approaches in the literature, the proposed approach, labeled the CSPTUKF, can handle any type of constraints during the sigma points generation step. In addition, the weights associated with each sigma point are updated at each time instance to match the mean and covariance of the given distribution. Further, PDF truncation is used to obtain accurate estimates of the states and the corresponding uncertainties. The performance of the CSPTUKF is shown to be better than or comparable to that of other algorithms available in the literature in terms of the average sum of square root estimation errors (ASSREE) for three case studies. In terms of computational efforts, the proposed CSPTUKF involves solving feasibility problems and was found to result in acceptable computational times. The real utility of our approach, however, is in its ability to handle constraints on states. The applicability of other approaches, for several constrained state estimation problems, will involve an element of luck, in the sense that these approaches might not work at all depending on the constraints (such as linear constraints) and the process parameters, as was demonstrated in the case studies considered in this article.

Figure 6. Tracking of states by CSPTUKF for the polymerization reactor.

figure that the trends of true and estimated states change at about the 3600th time instant. This is because of the step change in Fi. A similar type of state variation was observed for the TIUKF and URNDDR. The results for these methods are not plotted, as it is difficult to visually make out the differences between the states estimated by various algorithms. However, the state estimation accuracy in terms of ASSREE is better for the CSPTUKF than for these algorithms, with an acceptable computational time. The corresponding results are reported in Table 6. It is to be noted that the results for the TUKF and MSSPTUKF could not be generated because the sigma points generation step used in these algorithms is not able to handle the constraints. It was observed during the simulations that, when we propagated such infeasible

APPENDIX A. EQUIVALENCE BETWEEN KALMAN FILTER AND CSPTUKF In this appendix, the equivalence between the Kalman filter and the CSPTUKF is demonstrated for a linear (in terms of both state and measurement equations) unconstrained system. It is to be noted that, for unconstrained problems, the CSPTUKF formulation reduces to two steps, namely, a prediction step and an estimation step. Next, we show that, for a linear unconstrained problem, these two steps are exactly the same for the KF and the CSPTUKF.

Table 6. Comparison of the Performance of Various State Estimation Algorithms for a Polymerization Reactor algorithm

ASSREE

avg computation time (s) per time instance

TIUKF CSPTUKF URNDDR

4.5083 × 103 3.0633 × 103 7.3225 × 103

0.0152 0.6903 0.6107

T

5.475 1.9408 1.7 × 102 4.5 × 103 7 × 103 ]

1924

dx.doi.org/10.1021/ie2027487 | Ind. Eng. Chem. Res. 2013, 52, 1916−1926

Industrial & Engineering Chemistry Research

Article

To account for the effect of Qk, the SPs are regenerated from the computed x̂k+1(−) and Pk+1(−) values using an optimization formulation similar to that in eqs 17 and 18 and satisfy equations similar to eqs A.3−A.5. These new predicted SPs χi,k+1(−) and their weights wi,k+1(−) are then used for the next step.

Consider the following process and measurement model xk + 1 = ϕk xk + Γkuk + γk

(A.1)

yk = Ckxk + νk

(A.2)

where ϕk ∈ Rn×n is the state transition matrix, Γk ∈ Rn×m is the matrix relating inputs to states, and Ck ∈ Rq×n is the measurement matrix. Other symbols are as before. At the kth time instance, the sigma points χi,k(+) and associated weights wi,k(+) generated by the optimization formulation in eqs 17 and 18 satisfy the equations

A.2. Estimation or Data Assimilation Step

The new predicted sigma points are propagated through the measurement model to predict the measurement sigma points Υi,k+1(−) and their mean ŷk+1(−) Υi , k + 1( −) = Ck + 1χi , k + 1 ( −)

Np

xk̂ ( +) =

Np

∑ wi ,k(+)χi ,k (+)

yk̂ + 1 ( −) =

(A.3)

i=1

∑ wi ,k+ 1(−)Υi ,k+ 1(−) = Ck+ 1xk̂ + 1(−) (A.9)

i=1

Np

Pk( +) =

i = 1, ..., Np

The covariance in innovations is computed as

∑ wi ,k(+)[χi ,k (+) − xk̂ (+)][χi ,k (+) − xk̂ (+)]T

Np

i=1

Pee, k + 1( −) =

(A.4)

∑ wi ,k+ 1(−)[Υi ,k+ 1(−) − yk̂ + 1 (−)] i=1

Np

∑ wi ,k(+) = 1.0,

0 < wi , k( +) < 1

[Υi , k + 1( −) − yk̂ + 1 ( −)]T + R k + 1

i = 1, ..., Np

Np

i=1

(A.5)

=

i=1

A.1. Prediction or Forecast Step

− Ck + 1xk̂ + 1( − )] [Ck + 1χi , k + 1 ( − ) − Ck + 1xk̂ + 1(− )]T + R k + 1

The transformed set of sigma points is obtained by propagating each sigma point through the process model as χi̅ , k + 1 ( −) = ϕk χi , k ( +) + Γkuk

∑ wi ,k+ 1(−)[Ck+ 1χi ,k+ 1 (−)

i = 1, ..., Np

⇒ Pee, k + 1( − ) = Ck + 1Pk + 1(− )CkT+ 1 + R k + 1

(A.6)

(A.10)

The predicted mean is computed as

The cross covariance between the state prediction error and innovations and the Kalman filter gain are calculated as

Np

xk̂ + 1( −) =

∑ wi ,k(+)χi̅ ,k+ 1 (−) i=1

=

Np

Np

∑ wi ,k(+)[ϕkχi ,k (+) + Γkuk]

Pse, k + 1( −) =

= ϕk xk̂ ( +) + Γkuk

[Υi , k + 1( −) − yk̂ + 1 ( −)]T

(A.7)

Np ⎧ ⎪ wi , k + 1( −)[χi , k + 1 ( −) − xk̂ + 1( −)] =⎨ ∑ ⎪ ⎩ i=1 ⎫ ⎪ [χi , k + 1 ( −) − xk̂ + 1( −)]T ⎬ CT ⎪ k+1 ⎭

The SPECM is computed as Np

Pk + 1( −) =

∑ wi ,k+ 1(−)[χi ,k+ 1 (−) − xk̂ + 1(−)] i=1

i=1

∑ wi ,k(+)[χi̅ ,k+ 1 (−) − xk̂ + 1(−)] i=1

[χi̅ , k + 1 ( −) − xk̂ + 1( −)]T + Q k

= Pk + 1( −)CkT+ 1

or

(A.11) Np

Pk + 1( −) =

∑ wi ,k(+)[ϕkχi ,k (+) + Γkuk − ϕkxk̂ (+)

Kk + 1 = Pse, k + 1( −)[Pee, k + 1( −)]−1

i=1

− Γkuk][ϕk χi , k ( +) + Γkuk − ϕk xk̂ ( +) − Γkuk]T

= Pk + 1( −)CkT+ 1[Ck + 1Pk + 1( −)CkT+ 1 + R k + 1]−1 (A.12)

+ Qk

The equations then used for the updated mean of the state vector (eq 27) and SEECM (eq 28) are same are those used in the KF and, hence, no separate derivations are required. This establishes the equivalence.

Np ⎧ ⎪ wi , k( +)[χi , k ( +) − xk̂ ( +)] = ϕk ⎨ ∑ ⎪ ⎩ i=1 ⎫ ⎪ [χi , k ( +) − xk̂ ( +)]T ⎬ ϕT + Qk ⎪ k ⎭



AUTHOR INFORMATION

Corresponding Author

= ϕk Pk( +)ϕkT + Q k

*E-mail: [email protected]. Tel.: +91-22-25767214. Fax: +91-22-25726895.

(A.8) 1925

dx.doi.org/10.1021/ie2027487 | Ind. Eng. Chem. Res. 2013, 52, 1916−1926

Industrial & Engineering Chemistry Research

Article

(20) Narasimhan, S.; Rengaswamy, R. Reply to comments on “Robust and reliable estimation via unscented recursive nonlinear dynamic data reconciliation” (URNDDR). J. Process Control 2009, 19, 719−721. (21) Julier, S. J.; Uhlmann, J. K. Reduced sigma point filters for the propagation of means and covariances through nonlinear transformations. Proc. Am. Control Conf. 2002, 887−892. (22) Julier, S. J.; Uhlmann, J. K. Unscented filtering and nonlinear estimation. Proc. IEEE 2004, 92 (3), 401−422. (23) Johansson, K. H. The quadruple-tank process: A multivariable laboratory process with an adjustable zero. IEEE Trans. Control Syst. Technol. 2000, 8 (3), 456−465. (24) Kloeden, P. E.; Platen, E. Numerical Solution to Stochastic Differential Equations; Springer: New York, 1999.

Present Address ‡

Also affiliated with Reactor Projects Division, BARC, Mumbai, India 400085. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS This work was partially funded by a research grant from CSIR, India. REFERENCES



(1) Julier, S. J.; Uhlmann, J. K.; Durrant-Whyte, H. F. A new approach for filtering non-linear systems. Proc. Am. Control Conf. 1995, 1628− 1632. (2) Qu, C. C.; Hahn, J. Computation of arrival cost for moving horizon estimation via unscented Kalman filtering. J. Process Control 2009, 19, 358−363. (3) Qu, C. C.; Hahn, J. Process monitoring and parameter estimation via unscented Kalman filtering. J. Loss Prev. Process Ind. 2009, 22, 703− 709. (4) Romanenko, A.; Castro, J. A. A. M. The unscented filter as an alternative to the EKF for nonlinear state estimation: A simulation case study. Comput. Chem. Eng. 2004, 28, 347−355. (5) Romanenko, A.; Santos, L. O.; Afonso, P. A. F. N. A. Unscented Kalman filtering of a simulated pH system. Ind. Eng. Chem. Res. 2004, 43, 7531−7538. (6) Li, W.; Leung, H. Constrained unscented Kalman filter based fusion of gps/ins/digital map for vehicle localization. In Proceedings of the IEEE Conference on Intelligent Transportation Systems; 2003; pp 1362−1367. (7) Ma, N.; Bouchard, M.; Goubran, R. A. Dual perceptually constrained unscented Kalman filter for enhancing speech degraded by colored noise. In Proceedings of the 7th International Conference on Signal Processing; 2004; pp 2522−2525. (8) Vachhani, P.; Narasimhan, S.; Rengaswamy, R. Robust and reliable estimation via unscented recursive nonlinear dynamic data reconciliation. J. Process Control 2006, 16, 1075−1086. (9) Kadu, S. C.; Bhushan, M.; Gudi, R.; Roy, K. Modified unscented recursive nonlinear dynamic data reconciliation for constrained state estimation. J. Process Control 2010, 20, 525−537. (10) Teixeira, B. O. S.; Torres, L. A. B.; Aguirre, L. A.; Bernstein, D. S. On unscented Kalman filtering with state interval constraints. J. Process Control 2010, 20, 45−57. (11) Kandepu, R.; Foss, B.; Imsland, L. Applying the unscented Kalman filter for nonlinear state estimation. J. Process Control 2008, 18, 753−768. (12) Sipos, B. J. Application of the manifold-constrained unscented Kalman filter. In Proceedings of the IEEE/ION Symposium on Position, Location and Navigation; 2008; pp 30−43. (13) Kolas, S.; Foss, B. A.; Schei, T. S. Constrained nonlinear state estimation based on the UKF approach. Comput. Chem. Eng. 2009, 33, 1386−1401. (14) Prakash, J.; Patwardhan, S. C.; Shah, S. L. Constrained nonlinear state estimation using ensemble Kalman filters. Ind. Eng. Chem. Res. 2010, 49, 2242−2253. (15) Simon, D. Optimal State Estimation: Kalman, H∞, and Nonlinear Approaches; John Wiley & Sons: New Jersey, 2006. (16) Simon, D.; Simon, D. L. Constrained Kalman filtering via density function truncation for turbofan engine health estimation. Int. J. Syst. Sci. 2010, 41 (2), 159−171. (17) Coughanowr, D. R. Process Systems Analysis and Control, 2nd ed.; McGraw-Hill: New York, 1991. (18) Vachhani, P.; Rengaswamy, R.; Gangwal, V.; Narasimhan, S. Recursive estimation in constrained nonlinear dynamical systems. AIChE J. 2005, 51 (3), 946−959. (19) Simon, D. Kalman filtering with state constraints: A survey of linear and nonlinear algorithms. IET Control Theory Appl. 2010, 4 (8), 1303−1318.

NOTE ADDED AFTER ASAP PUBLICATION This paper was published on the Web on January 23, 2013, with minor text errors. The corrected version was reposted on January 25, 2013.

1926

dx.doi.org/10.1021/ie2027487 | Ind. Eng. Chem. Res. 2013, 52, 1916−1926