A Moving Window Formulation for Recursive Bayesian State

Article pubs.acs.org/IECR

A Moving Window Formulation for Recursive Bayesian State Estimation of Systems with Irregularly Sampled and Variable Delays in Measurements† Vinay A. Bavdekar Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Alberta, Canada T6G 2G6

Jagadeesan Prakash Department of Instrumentation Engineering, Madras Institute of Technology Campus, Anna University, Chennai, 600004, India

Sachin C. Patwardhan Department of Chemical Engineering, Indian Institute of Technology Bombay, Powai, Mumbai, 400076, India

Sirish L. Shah* Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Alberta, Canada T6G 2V4 S Supporting Information *

ABSTRACT: The time delay involved between sampling and obtaining measurements of certain quality variables is a common scenario in various process applications. Further, this delay is not fixed and can vary for various reasons. Moreover, certain measurements may be sampled at irregular time intervals. The state estimation algorithms available in the literature have been developed for the scenario where the measurements are sampled regularly or are available after a fixed time delay. In this work, a recursive moving window Bayesian state estimator formulation is proposed to utilize such measurements with variable time delays to compute the state estimates. The length of the moving window ensures that the algorithm utilizes all the available measurements (delayed or otherwise) for computing the state estimates. In practice, it may also become necessary to account for the physical bounds on the states. A constrained version of the moving window recursive state estimator is also developed to yield state estimates that are consistent with their respective bounds and constraints. The efficacy of the unconstrained moving window state estimator is demonstrated by application on the benchmark Tennessee Eastman simulation case study and an experimental two-tank heater−mixer setup, while the efficacy of the constrained moving window state estimator is demonstrated by simulation of a benchmark gas-phase batch reactor system.

1. INTRODUCTION In several chemical and biochemical processes, measurements of various quality variables (such as product concentration, melt viscosity, polymer average molecular weight, dry weight, etc.) are often available after a certain time due to inherent delays in the measurement. Moreover, there are certain primary variables (such as biomass concentration) or efficiency parameters that are difficult to measure online or costly to measure at regular intervals. Measurements of such variables may be available through offline laboratory analysis and/or calculations with irregular measurement delays. Thus, only a subset of the measurements, such as pressure, temperature, level, etc., are available online and at regular intervals. The conventional nonlinear Bayesian state estimation algorithms such as the extended Kalman filter (EKF),1 unscented Kalman filter (UKF),2 ensemble Kalman filter (EnKF),3 and particle filters4 have been developed under the assumption that the measurements are available at regular intervals. There is a rich literature for state estimation of nonlinear systems5−7 for the scenario © 2014 American Chemical Society

where the measurements are sampled at regular intervals, without any time delay. State estimation using delayed and irregularly sampled measurements has received relatively less attention in the literature. The accuracy of the state estimates obtained can be improved if the delayed measurements can be incorporated into the update step of the state estimation algorithms. Modified versions of the EKF have been proposed in the literature to incorporate the delayed measurements for computing the state estimates. Gopalakrishnan et al.9 provide a comprehensive review of different modifications in the EKF algorithm that incorporate the delayed and irregularly sampled measurements. These modifications have been broadly classified into two different categories: (a) methods that are Received: Revised: Accepted: Published: 13750

March 5, 2014 July 29, 2014 August 7, 2014 August 7, 2014 dx.doi.org/10.1021/ie5009585 | Ind. Eng. Chem. Res. 2014, 53, 13750−13763

Industrial & Engineering Chemistry Research

Article

order hold) until the measurement becomes available at a subsequent sampling interval. Thus, in this approach the size of the measurement vector and measurement noise covariance matrix are constant in the MHE objective function. The second approach proposed is a “variable structure” approach,14,15 in which only those measurements that are available at a particular sampling instant are used to compute the optimal estimate of the states. In this approach, the size of the measurement vector and noise covariance matrix, in the MHE objective function, will change depending on the availability of the measurements. The advantage of the fixed structure approach is that, in certain applications, it does not cause large fluctuations in the estimation error covariance, when the slow or delayed measurements become available. This is, however, a possibility with the variable structure approach, which may result in a nonsmooth decay of the estimation error.14 However, computing the “arrival cost” in the MHE algorithm is an important issue and an open research problem, especially for nonlinear systems. Various approaches proposed13,15−17 involve the use of a Bayesian state estimator to approximate the arrival cost, which directly affects the accuracy of the state estimates obtained. The state estimation algorithms for use of delayed measurements are developed under the assumption that the time delay in obtaining the measurement is constant. There is very little work reported for the situation where the time delay in obtaining the measurement is variable and not constant. In practice, however, different measurements might be delayed by different time lengths and the delay in each measurement can be variable. In this work, a generalized framework, based on the moving window formulation, is proposed for state estimation of nonlinear processes by incorporating measurements that are irregularly sampled and/or available with a time delay that may vary. The proposed framework can be used with any nonlinear Bayesian state estimator. The main contributions of this work are as follows. A moving window formulation of the Bayesian state estimators is proposed to utilize the delayed and/or irregularly sampled measurements for reconstructing the estimates of the process states. The proposed approach combines the advantages of the MHE formulation and those of the Bayesian state estimation approach. The moving window formulation allows incorporating measurements with variable time delays and/or irregular sampling intervals, while the Bayesian approach yields the minimum mean squared state estimates, along with their statistical moments. For processes operating close to constraints on states or parameters, the moving window formulation can use constrained Bayesian state estimators18−22 to account for the state and/or parameter constraints. The proposed formulation falls into the second category of state estimation methods based on the classification by Gopalakrishnan et al.9 The rest of the paper is organized as follows. Section 2 describes the assumptions involved in constructing the process model used for simulations and state estimation. The proposed moving window state estimator formulation for utilizing the delayed and/or irregularly sampled measurements to reconstruct the state estimates is described in section 3. In section 4, the efficacy of the proposed moving window state estimator is demonstrated on the Tennessee Eastman challenge problem, in which the none of the states are directly measured and a large number of measurements are sampled irregularly and are available after a certain delay. The efficacy of the constrained moving window state estimator is demonstrated on the two-

based on state augmentation and (b) methods that provide a procedure to fuse the delayed measurements based on their availability. Gudi et al.8 propose to augment the model equations with extra integrating states equal to the number of sampling instants required to obtain the delayed measurement. One drawback of this approach, which has been pointed out by Gudi et al.8 themselves, is the possibility of loss of observability of the augmented system. Moreover, it is difficult to handle variable time delays, since the number of integrated states are set equal to the number of sampling instants by which the measurements are delayed. Mutha et al.10,11 propose a fixed-lag smoothing based EKF to recompute the state trajectories for the sampling intervals equal to the delay in obtaining the measurements. In this algorithm, the states between the current sampling instant and previous sampling instants equal to the delay are smoothed once the delayed measurement becomes available. This is achieved by running two EKFs in parallel: one which computes the filtered state based on measurements available without delay and the second which preempts the Kalman gain for the sampling instants between when the delayed measurement is sampled and when it becomes available. Once the delayed measurement becomes available, the states are filtered using only the delayed measurement and smoothed using the latest nondelayed measurements. The principal disadvantage in this approach is the storage and computation cost associated with smoothing. Since the covariance matrices and states have to be stored, the costs increase with the state dimensions and with the delay involved in obtaining the slow measurement. Moreover, in the proposed approach, the linearization of the state and measurement equations is carried out using the available state estimates and are not corrected based on the smoothed estimates. This may lead to incorrect computation of the covariance matrices and Kalman gain, thereby leading to a possibility of obtaining inaccurate state estimates. Prasad et al.12 propose to use two EKFs in parallel. However, this algorithm does not involve smoothing of the states, which reduces the storage and computation costs. The first EKF computes the state estimates based on the fast-rate measurements. When a delayed measurement becomes available, the second EKF updates the states corresponding to the sampling instant at which the delayed measurement was sampled. These corrected states are then passed on to the first EKF, which then recomputes the trajectory of the filtered state estimates based on the available fast rate measurements. An important limitation of this approach is that assumes that all the delayed measurements have an identical time delay. Thus, the number of EKFs running in parallel increases with the different time delays involved with different measurements. The moving horizon estimation (MHE) algorithm13 was developed to account for the constraints on the states or parameters. A constrained maximum likelihood estimates (MLE) based optimization problem is formulated to obtain the estimates of the states and/or parameters, using a moving window of past measurements. The MHE algorithm can also be modified to incorporate the delayed measurements for state estimation. This can be achieved by setting the length of the moving window at least equal to the maximum anticipated delay in availability of the measurements.14 Krämer and Gesthuisen14 propose two approaches to formulate an MHE for multirate and/or delayed measurements. The first approach is a “fixed structure” MHE, in which the last prediction error of the slow and/or delayed measurement is held fixed (a zero13751

dx.doi.org/10.1021/ie5009585 | Ind. Eng. Chem. Res. 2014, 53, 13750−13763


Article i k = ki,l the measurement ys,k is available, due to the measurement delay, di, it actually represents a measurement taken at sampling instant k − di. Here, di denotes the number of sampling instants by which the measurement is delayed. Thus, defining ti = ki,l − di, the measurement model for the ith irregularly sampled variable is represented as follows:

state gas-phase reactor. The results obtained by applying the moving window state estimator for state estimation on an experimental case study are presented in section 5, and the conclusions derived from this study are presented in section 6.

2. PROCESS MODEL FOR STATE ESTIMATION The first-principles model for a continuous nonlinear process can be described using the following nonlinear equations dx = f(x , u , t ) dt

ys,i t = h s, ti(x s, ti) + vtii i

where i = 1, 2, ..., ns. It is further assumed that wk, vf,k, and vs,k are mutually uncorrelated, independent, and identically distributed random variables. It may be noted that the value of di may vary with time and each element in ys may have a different di.

(1)

yT(t ) = h(x , t )

where x ∈ n denotes the states of the process, u ∈ m denotes the manipulated inputs, and yT ∈ r denotes the process measurements obtained. The measurement vector yT is further partitioned as yT = [ y Tf y sT ]T , where yf is the subsequence of measurements that are not delayed and ys is the subsequence of measurements that are available after a certain time delay, which may be variable. It may be noted that both types of measurementsdelayed or otherwisemay be sampled at irregular intervals. The process simulations and modeling for state estimation are carried out under the following simplifying assumptions Assumption 1. The measurement noise is modeled as a zeromean white noise sequence with a known distribution. Assumption 2. The manipulated inputs are piecewise constant over the smallest sampling interval, T: u (t ) = u k

∀ tk ≤ t < tk + T

3. STATE ESTIMATION FOR SYSTEMS WITH IRREGULAR AND DELAYED MEASUREMENTS 3.1. Brief Review of Conventional State Estimation Algorithms. This section includes a brief review of nonlinear Bayesian state estimation algorithms for the case where measurements are sampled regularly and are available without any delay. These recursive Bayesian state estimators can be viewed as minimum mean square error (MMSE) estimators for nonlinear systems. The state estimation consists of two main steps: (a) prediction and (b) measurement update or correction of the first two moments of the distribution of the states. The prediction step propagates the mean of the states as follows:

(2)

(5)

yk̂ | k − 1 = h(x̂k | k − 1) + vk

(6)

K k = Pεe , k Pee , k −1

(7)

where Pεe,k is the cross covariance between the prediction error of the states and the innovations and Pee,k is the covariance of the innovations. These matrices can be obtained either through linearization of the state and measurement equations (EKF) or as sample statistics (UKF or EnKF). The corrected states are then obtained as follows:

kT

∫(k−1)T f(x(τ), uk−1) dτ + wk−1

= F(x k − 1, uk − 1, wk − 1) yf, k = h f, k(x k) + vf, k ys, k = h s, k(x k − di) + vs, k − di

x̂k | k − 1 = F(x̂k − 1 | k − 1, uk − 1, wk − 1)

This can be achieved by either directly propagating the updated mean of states obtained at the previous sampling instant (as in the EKF), or as a sample mean of the predictions obtained by drawing samples of the updated mean of the states and process noise (as in the UKF or EnKF) from the previous sampling instant. The key step in MMSE state estimation is computation of the observer gain matrix in the measurement update step and correcting the state estimates using this Kalman gain, the measurements and the predicted states. The observer gain matrix is computed as

Assumption 3. The choice of the sampling interval, T, is small enough that the variations in the unmeasured disturbances can be approximated as piecewise constant functions. The unmeasured disturbances are assumed to evolve as a white noise sequence with a known distribution, whose source is unknown. Their effect on the states is modeled as additive in the process dynamics. Assumption 4. The irregular measurements are available at sampling intervals, which are integral multiples of the smallest sampling time, T. Further, the time delay involved in obtaining the delayed measurements is assumed to be an integral multiple of the smallest sampling time, T. Under assumptions 1−4, a discrete time representation of the true process can be obtained as follows: xk = xk−1 +

(4)

(3)

n

where wk − 1 ∼  is a zero-mean white noise sequence, with a known distribution. yf,k represents the subset of those measurements which are available instantaneously but may be sampled at irregular intervals. The subset of irregularly sampled and delayed measurements is represented in a different manner. Consider the ith irregularly sampled and delayed measurement element, yis, and let {ki,l: l = 1, 2, ...} represent the subsequence of sampling instants {k = 0, 1, 2, ...} for which the measurements of ys,k become available. Now, even though at

ek = yk − yk̂ | k − 1

(8)

x̂k | k = x̂k | k − 1 + K kek

(9)

3.2. Moving Window Formulation for State Estimation Using Delayed and/or Irregularly Sampled Measurements. To utilize multirate sampled measurements with variable time delays in updating the states, the concept of a moving window is used, which is defined as follows: Let dmax = maxi di,max, where di,max represents the maximum anticipated delay in obtaining the ith (i = 1, 2, ..., r) measurement. Then, the minimum length of the moving window is defined as N = dmax + 1. At the kth sampling instant (k > N), consider the set 13752



Article

of N past measurements, defined as YN = {yk−N, ..., yk}. It may be noted that the size of the measurement vectors can change within the window depending upon the delayed measurement becoming available at the kth sampling instant. Suppose a delayed measurement of the output y becomes available at the kth sampling instant, after a delay of di samples, then the vector yk−di is appended to include this new measurement. The concept of the moving window is demonstrated using the schematic shown in Figure 1. The schematic shows

over the window [ k − N k ] starting from the initial state x̂k−N−1|k−N−1. The steps involved in implementing the EnKF are as follows: (l) Np samples of the process and measurement noise {w(l) j−1, vj : l = 1, 2, ..., Np}; j = {k − N, ..., k} are drawn randomly from their respective known distributions. In the prediction step, the ensemble of the states x̂j−1|j−1: j = {k − N, ..., k} are propagated as follows: x̂(jl|)j − 1 = F(x̂(jl−) 1 | j − 1, u j − 1, w(jl−) 1)

(10)

ŷ (jl|)j − 1 = hj(x̂(jl|)j − 1) + v(jl)

(11)

ŷ(l) j|j−1

v(l) j

It should be noted that the dimensions of and will change depending on the number of measurements available at the jth sampling instant. The ensemble of the state and measurement predictions are then used to estimate the sample mean and covariance as follows: x̅ j | j − 1 =

yj̅ j − 1 = Figure 1. Schematic representation of the moving window of measurements.

Pεe , j =

representations of measurement vectors for two consecutive time windows, [ k − N − 1 k − 1] and [ k − N k ], which contain three different types of measurements: (i) y1 represents the irregularly sampled measurements, which are available after a certain time delay, (ii) y2 represents the irregularly sampled measurements available without any delay, and (iii) y3 represents the regularly sampled measurements, available without any time delay. As can be seen from the schematic, when moving from window [ k − N − 1 k − 1] to window [ k − N k ], the measurement vectors yj belonging to the common subset of measurements in YN−1 and YN can have different dimensions depending on the new and/or delayed measurements become available at the kth sampling instant, i.e., yj ∈ q, where q = {1, 2, ..., r}. The choice of the window length ensures that the length of the measurement vector does not change at the jth sampling instant in the past such that j < k − N. In this work the ensemble Kalman filter (EnKF) is used for developing the moving window state estimation framework. The EnKF algorithm is a combination of Monte Carlo based filters, which draw random samples of the states and disturbances from their respective distributions, and minimum mean square error (MMSE) estimators, in which the update step is based on minimizing the estimation error covariance. Thus, the EnKF provides an advantage that it can be used for state estimation in the presence of disturbances with nonGaussian distributions along with the computational efficiency resulting from the requirement to compute only the first two moments of the distribution. The details of the moving window state estimation algorithm are explained using the EnKF framework. However, it should be noted that the moving window framework can be used with any recursive Bayesian state estimator, with appropriate modifications. At each sampling instant k, it is proposed to implement the EnKF

Pee , j =

Np

1 Np

∑ x̂(jl|)j − 1

1 Np

∑ ŷ(jl|)j − 1

l=1

(12)

Np l=1

1 Np − 1 1 Np − 1

(13)

Np

∑ [εj(l)][e(jl)]T l=1

(14)

Np

∑ [e(jl)][e(jl)]T l=1

(15)

where εj(l) = x̂(jl|)j − 1 − x̅ j | j − 1

(16)

e(jl) = ŷ (jl|)j − 1 − yj̅ | j − 1

(17)

The Kalman gain matrix and the ensemble of updated states are calculated as follows: K j = Pεe , jPee , j−1

(18)

Ψ(jl|)j − 1 = yj − ŷ (jl|)j − 1

(19)

x̂(jl|)j = x̂(jl|)j − 1 + K jΨ(jl|)j − 1

(20)

The dimensions of Pεe , j ∈ n × q , Pee , j ∈ q × q may change with j, depending on q, the number of measurements available. The resultant dimensions of the Kalman gain matrix are K j ∈ n × q and the dimension of the innovations vector is

Ψ(jl|)j − 1 ∈ q × 1. The updated states are computed as the mean of the updated ensemble:

1 x̅ j | j = Np

Np

∑ x̂(jl|)j l=1

(21)

To summarize, the steps involved in the implementation of the state estimator at the kth sampling instant are as follows: (i) Obtain the measurements that become available at the kth sampling instant and, depending on the measurement delay in 13753



Article

the individual channel, update the measurement vectors in the window [ k − N k ]. (ii) Implement the EnKF described above (eqs 10−21) over the window [ k − N k ]. In order to implement the moving window state estimator, the following variables need to be stored: statistical properties of the posterior of the states at the beginning of the window, the measurements, and the manipulated inputs for the sampling instants in the moving window. For the moving window EnKF, the dimensions of the variables that need to be stored are n × Np for the ensemble of the states, maximum of r × N for the measurements, and m × N for the manipulated inputs. When going from the kth to the (k + 1)th sampling instant, the inputs and the measurements from the (k − N + 1)th instant to the kth instant and the ensemble of the states x̂(l) k−N+1|k−N+1 (l = 1, 2, ..., Np) are retained. It can be seen that the computation time of the algorithm is also dependent on the length of the moving window, N, which denotes the largest delay in obtaining a measurement after sampling. Thus, as N increases, the computation time increases. It may also be noted that if none of the measurements are delayed (dmax = 0), then the length of the moving window reduces to 1, and the algorithm reduces to the conventional state estimation algorithms. 3.3. Constrained Moving Window State Estimator for Delayed and Irregular Measurements. In many practical problems of interest, it becomes necessary to account for the physical bounds on the states and parameters being estimated. Various constrained state estimation formulations, such as constrained EKF,18 constrained UKF,19−21 and constrained EnKF,22 have been proposed for the scenario in which measurements are regularly available. These approaches can be incorporated into the moving window formulation to accommodate delayed and/or irregular measurements. The moving window EnKF algorithm described in section 3.2 is modified as follows, to accommodate for the constraints on the states/parameters. If the states or process noise is constrained, it becomes necessary to generate the random particles of the states and noise which are consistent with these bounds. The concept of “truncated distributions” can be used to generate such random particles. A truncated distribution is a conditional distribution, conditioned on the bounds of its corresponding random variable. For example, given a probability density function f(ζ), and its cumulative distribution function ϕ(ζ) defined over (−∞,∞), the truncated probability density function for (a ≤ ζ ≤ b) is defined as follows:22 f (ζ ) f (ζ |a ≤ ζ ≤ b) = ϕ(b) − ϕ(a)

η=

b

f (ζ |a ≤ ζ ≤ b) =

1 ϕ(b) − ϕ(a)

P1/2

(25)

⎡ s11 ⎢ ⎢s = ⎢ 21 ⎢⋮ ⎢⎣ sn1

0 s22 ⋮ sn2

... 0 ⎤ ⎥ ... 0 ⎥ ⎥ ⋱ ⋮ ⎥ ... snn ⎥⎦

(26)

truncated univariate normal distributions 5T(0, 1|t L, j ≤ t j ≤ t H, j) can be defined where j−1

t L, j =

xL, j − x̅ j − ∑r = 1 sjr tr sjj j−1

t H, j =

xH, j − x̅ j − ∑r = 1 sjr tr sjj

(28)

The above transformation requires that the random samples be drawn recursively. Thus, t 1 is drawn first from 5T(0, 1|t L,1 ≤ t1 ≤ t H,1) by setting the limits [tL,1, tH,1]. Then, the limits [tL,2, tH,2] are set and t2 is drawn from 5T(0, 1|t L,2 ≤ t 2 ≤ t H,2). This procedure is repeated until tn is drawn. Thus, a sample from the truncated normal distribution 5T(x̅ , P), can be expressed as follows: x = x̅ + P1/2T T = [ t1 t 2 ··· tn ]T (29) At the kth sampling instant, the constrained moving window EnKF is implemented over the window [ k − N k ]. The prediction step in the constrained state estimator is identical to the one in the unconstrained moving window EnKF (eqs 10 and 11). At the jth sampling instant (j ∈ [ k − N k ]), any particle of the predicted states, which violates the constraints, is projected back onto the constrained space20

f (ζ ) d ζ = 1

(l) x̂(jl|),c j − 1 = Pr[x̂ j | j − 1]

In particular, the expression for a truncated univariate normal distribution is as follows:

(30)

∀ l = 1, 2, ..., Np. The projection operator is applied as follows Pr[x̂(jl|)j − 1] = x̂(jl|)j − 1

η2 −2

( )

1 σ(2π )1/2 ϕ(ηb) − ϕ(ηa)

b − ζ̅ σ

(27)

(23)

5T(ζ ̅ , σ 2|a ≤ ζ ≤ b) =

ηb =

⎡ xL,1 − x1̅ ⎤ ⎤ ⎡ xH,1 − x1̅ ⎢ ⎥ ⎥ ⎢ s s 11 11 ⎢ ⎥ ⎥ ⎢ ⎡ ⎤ t 1 ⎢x − x − s t ⎥ ⎥ ⎢x − x − s t H,2 2 21 1 L,2 2 21 1 ⎢ ⎥ ̅ ̅ ⎢ ⎥ ⎥ t2 ⎥ ⎢ ⎢ ⎢ ⎥ ⎥≤ s22 s22 ≤⎢ ⎢ ⎥ ⎢ ⎥ ⎥ ⋮⎥ ⎢ ⎢ ⋮ ⋮ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ t ⎢ ⎥ ⎥ ⎣ n⎦ ⎢ n−1 n−1 ∑ ∑ x x s t x x s t − − − − ⎢ L, n ⎢ H, n ̅n ̅n r = 1 nr r ⎥ r = 1 nr r ⎥ ⎢⎣ ⎥⎦ ⎥⎦ ⎢⎣ snn snn

b

exp

a − ζ̅ , σ

Using a transformed random vector T, such that

(22)

∫a

ηa =

Consider a random variable vector x, which has a truncated normal distribution 5T(x̅ , P), bounded in [xL, xH]. In this work, the following approach is used to generate random samples from a given truncated normal distribution.22 Since P is a symmetric, positive definite matrix, its Cholesky decomposition is obtained as

such that

∫a

ζ − ζ̅ , σ

(24)

where 13754

∀ xL ≤ x̂(jl|)j − 1 ≤ xH

(31)

= xL

∀ x̂(jl|)j − 1 ≤ xL

(32)

= xH

∀ x̂(jl|)j − 1 ≥ xH

(33)



Article

Figure 2. Schematic of the Tennessee Eastman process.

Koläs et al.23 have shown that the above approach is the optimal least-squares solution, if the constraints are on individual states. The measurement update step is modified as follows. The sample covariance of the predicted state particles at the jth sampling instant, in the moving window, is obtained as Pj | j − 1 =

1 Np − 1

The solution of the above optimization problem yields the updated state particles, x̂(l) j|j (l = 1, 2, ..., Np), such that they are consistent with their respective constraints. The updated mean of the state estimates is computed using eq 21. Remark 1. The dimensions of yj and hj(·) change according to the number of measurements available at the jth. Consequently, the size of the weighting matrix (measurement noise covariance), Rj, also changes to reflect the number of measurements available. Remark 2. When the measurement noise model is linear, i.e., yk = Ckxk + vk, then the constrained optimization problem defined by eqs 36−39 can be reduced to the following quadratic optimization problem:22

Np T c c (l),c ∑ [x̂(jl|),c j − 1 − x̅ j | j − 1][x̂ j | j − 1 − x̅ j | j − 1] l=1

(34)

where x̅ cj | j − 1

1 = Np

Np

∑ x̂(jl|),c j−1 l=1

x̂ j | j = arg min(x j)T Hj(x j) + fTj (x j) xj

(35)

where

For the lth particle, at the jth sampling instant, measurement update is now obtained as a solution of the following constrained optimization problem: x̂(jl|)j = min[(εj(l))T Pj | j − 1−1(εj(l)) + (e(jl))T R j−1(e(jl))] xj

(36)

εj(l) = x j − x̂(jl|),c j−1

(37)

e(jl)

= yj − [hj(x j) +

Hj = [Pj | j − 1−1 + CTj R j−1Cj]

(41)

−1 T fTj = −2[(x̂(jl|),c + (yj − v(jl))R j−1Cj] j − 1) Pj | j − 1

(42)

subject to

v(jl)]

xL ≤ x j ≤ xH

(38)

(43)

4. SIMULATION CASE STUDIES 4.1. Tennessee Eastman Problem. The proposed moving window EnKF is applied for state estimation on the benchmark Tennessee Eastman (TE) process challenge problem.24 The TE

Subject to

xL ≤ x j ≤ xH

(40)

(39) 13755



Article

sampled at every 3 s, as they are used to compute the output of their respective PI controllers. The values of these measurements sampled at every 0.05 h are used for the purpose of state estimation. The measurements of the gas-phase and product concentrations are sampled at every 0.1 h (every two sampling instants) and 0.25 h (every five sampling instants) respectively,24 thereby introducing an irregularity in the sampling of the measurements. Further, there is a variable delay between 0.15 and 0.25 h (three to five sampling instants) in the availability of the product stream measurements. The gas-phase concentration measurements are delayed by 0.1 h (two sampling instants). Thus, there are six regularly sampled measurements, which are not delayed, and 10 irregularly sampled and delayed measurements. The random unmeasured process disturbances are simulated by adding a zero-mean Gaussian white noise to the state dynamics. Similarly, the measurements are corrupted with a zero-mean Gaussian white noise at the times for which the measurements are sampled. The values of the noise covariance matrices are taken from the work by Kadu et al.21 and are given in Table 2.

process is highly nonlinear and involves measurements sampled at different sampling rates. Also, there is a variable time delay in obtaining some of the measurements. The plant produces two products (G and H) using four reactants (A, C, D, and E). The process produces two more side products, of which one is an inert (B) and the other is a byproduct (F). A schematic of the TE process is shown in Figure 2. The plant consists of five major unit operations: reactor, condenser, vapor/liquid separator, recycle compressor, and stripper. The four reactants enter the reactor, which produces gaseous products. These gaseous products, along with the unconverted reactants, pass through a condenser followed by a vapor/liquid separator, where the condensed products get separated from the unconverted reactants. The unconverted reactants are passed through a recycle compressor and sent back into the reactor. The products (G and H) are obtained from the bottoms of the stripper. A purge stream, from the top of the vapor/liquid separator, removes the inert (B) and byproduct (F) from the process streams. The nonlinear mechanistic model of the TE process developed by Ricker and Lee25 is used for the purpose of simulations and state estimation. The steady state of the process is taken at the base case described by Ricker and Lee.25 The model consists of 26 states and 10 manipulated inputs. The states and their values for the base case are given in Table 1 of the Supporting Information. The process is open loop unstable and is first stabilized by closing the unstable loops using proportional−integral (PI) controllers. The reactor level, stripper level, separator level, and reactor pressure are locally regulated using PI controllers, which now makes the entire process open loop stable. The tuning parameters of the PI controllers, taken from the work by Karra et al.,26 are given in Table 2 in the Supporting Information. While carrying out state estimation, the PI loops on these variables are closed and their set points can be manipulated. The smallest sampling time for the state estimation is 0.05 h (3 min). However, the PI controllers take action at every 3 s. The fast rate of action by the PI controllers is required to maintain the stability of the process. The 16 measurements that are chosen for the state estimation problem are given in Table 1. Of these, the reactor pressure (Pr), reactor liquid level (VLr), separator liquid level (VLs), and stripper bottom level (VLp) measurements are

Table 2. TE Process: Process States and Their Steady State Values

variable

description

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Pr VLr Ps VLs VLp Pm 100yA9 100yB9 100yC9 100yD9 100yE9 100yF9 100yG9 100yH9 100yG11 100yH11

reactor pressure reactor liquid level separator pressure separator liquid level stripper bottom level stripper pressure A in purge B in purge C in purge D in purge E in purge F in purge G in purge H in purge G in product H in product

unit kPa % kPa % % kPa mol mol mol mol mol mol mol mol mol mol

% % % % % % % % % %

value

Q R x̂0

10−6 × I26 square of std dev given in Table 1 1.05 × xss

4.1.1. Results. The maximum delay in obtaining the measurements is for the product stream concentrations, which is five sampling instants. Hence the length of the moving window is N = 5 + 1 = 6. State estimation is carried out using the proposed moving window EnKF algorithm, described in section 3.2, with Np = 500. The simulation is carried out for a period of 50 h of operating time for the TE process. The initial conditions of the process states and the manipulated inputs are the steady state of the base case.25 A 5% error is introduced in the initial states for the moving window EnKF as given in Table 2. A +10% step change was introduced in F1, the feed flow rate of reactant A, and maintained at 10 h < t ≤ 40 h. Further, a −10% step change was introduced in F2, the feed flow rate of reactant D, at 25 h < t ≤ 40 h. The values of F1 and F2 are brought back to their steady state values and maintained at the same values at 40 h < t ≤ 50 h. The simulations were carried out in Matlab 2013a. The process equations were solved using the explicit Euler method, with a time step of 3 s as suggested by Ricker and Lee.25 On a computer equipped with an Intel Core i7-4770K 3.50 GHz processor and 16 GB of RAM, the simulation time for completing a 50 h run of state estimation for the TE problem using the proposed moving window EnKF algorithm was approximately 8 h 10 min. Out of the 26 states, the figures for the trajectories of the estimation errors of four states are shown here. Two of the states, namely NAr (kilomoles of A in reactor) and NEr (kilomoles of E in reactor), are unmeasured, and the estimation errors for these are shown in parts a and b, respectively, of Figure 3. The other two, namely NGp (kilomoles of G in the product stream) and NHp (kilomoles of H in the product stream), have indirect measurements, which are irregularly sampled and delayed, and the estimation errors for these are shown in parts a and b, respectively, of Figure 4. From Figures

Table 1. TE Process: List of Measurements no.

parameter

std dev 0.30 0.50 0.30 1.00 1.00 0.30 0.25 0.10 0.25 0.10 0.25 0.025 0.05 0.05 0.5 0.5 13756



Article

Figure 3. (a) Error in estimates of NAr and (b) error in estimates of NEr in the reactor.

Figure 4. (a) Error in estimates of NGp and (b) error in estimates of NHp in the stripper.

3 and 4, it can be seen that the estimation errors of the unmeasured states and the states with delayed and irregular measurements come close to and remain near zero, despite the mismatch in the initial conditions and noise. The reduction of the estimation error to values close to zero indicates that the moving window EnKF is able to track the states with reasonably good accuracy. In order to benchmark the results of the simulations under multirate sampling and variable time delays, a comparison of the state estimates is presented for three cases: case A, where only those measurements that are available without delay are used (six measurements); case B, in which the delayed and irregularly sampled measurements are used along with the measurements available without delay and state estimation is carried out using the proposed moving window EnKF algorithm (16 measurements); and case C, where all 16

measurements are available at regular time intervals and without any delay (16 measurements). State estimation for cases A and C is carried out using the conventional EnKF algorithm,3 as the length of the moving window is reduced to N = 1. It should be noted that cases A and C represent two extreme scenarios of measurement availability. In one case, none of the delayed measurements are available, and the other represents an ideal scenario where all the measurements are regularly sampled and measurement delays are absent. The process and simulation conditions are identical to those mentioned previously. Similarly, the conditions for the state estimator are also identical to those used in the case where measurements are delayed. To make conclusions independent of specific noise realizations, 10 Monte Carlo runs of the TE problem were performed to obtain the root mean square error (RMSE) of the 13757



Article

the three cases. The RMSEs of the state estimates reported are the mean of the RMSE values obtained in the 10 runs. When the RMSEs of the state estimates are compared for case A and case B, it can be seen that 16 out of the 26 states have a lower or comparable RMSE when the proposed moving window state estimation algorithm is used along with the delayed and/or irregularly sampled measurements (see Table 3), while for the remaining 10 states there is a deterioration in performance of the proposed state estimator (see Table 4). Further improvement can be obtained over case A and the performance can be moved closer to case C if the sampling frequency of the gas and liquid phase concentration measurements is increased. However, the sampling frequency of these concentration variables used in case B is chosen consistent with the sampling frequency specifications given in the TE problem.24 4.2. Gas-Phase Reactor. Consider the following irreversible gas-phase reaction in a well-mixed, constant volume isothermal batch reactor:6

state estimates for all three cases. Each run was for 50 h of the TE process. The initial conditions of the process and the state estimator were maintained at those given in Table 1 (in the Supporting Information) and Table 2 (in this text), respectively. The input changes, too, were identical in each run. However, the realizations of the state and measurement noise were different in each run. For every state, the RMSE is computed as follows: εx , j = xj − xĵ | j

RMSE =

1 Nd

(44) Nd

∑ (εx ,j)2 (45)

j=1

where Nd indicates the number of data points, xj is the true value of the state, and x̂j|j is the estimated value. The RMSE values of the state estimates are presented in Tables 3 and 4 for Table 3. TE Process: RMSEs of State Estimates for Case A, Only Measurements That Are Not Delayed; Case B, Delayed and Irregular Measurements; and Case C, All Measurements Regular and Not Delayeda

The first-principles model for this reaction is given by the following differential equations:

RMSE

a

no.

state

case A

case B

case C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

NFr NFs NFm NAs NAm NGp NDs NGs NAr NCm NCs NCr NHs NDr NGm NHm

0.0715 0.0421 0.1422 0.1171 0.1279 0.0386 0.0009 0.0456 0.0192 0.3844 0.3229 0.0523 0.0320 0.0011 0.0086 0.0045

0.0539 0.0318 0.1079 0.0967 0.1070 0.0332 0.0008 0.0408 0.0177 0.3609 0.3034 0.0496 0.0318 0.0012 0.0087 0.0052

0.0051 0.0035 0.0074 0.0189 0.0218 0.0331 0.0009 0.0402 0.0066 0.0318 0.0254 0.0051 0.0293 0.0012 0.0046 0.0048

dpA

a

case B

case C

1 2 3 4 5 6 7 8 9 10

NBs NBm NBr NHp NGr NEs NHr NEm NEr NDm

0.1040 0.1197 0.0174 0.0233 0.1354 0.0197 0.1266 0.0453 0.0245 0.0085

0.1106 0.1279 0.0186 0.0256 0.1491 0.0223 0.1504 0.0549 0.0333 0.0128

0.0135 0.0155 0.0046 0.0255 0.1151 0.0109 0.1036 0.0238 0.0161 0.0131

= k1pA 2

(48)

⎡ pA ⎤ PA + B = [1 1]⎢ ⎥ ⎣ pB ⎦

(49)

dt

where k1 = 0.16 is the reaction rate constant and x = [ pA pB ]T are the partial pressures of A and B, respectively. The total pressure PA+B, sum of the two partial pressures, is available as a measurement. For the process simulations and state estimation, the sampling time is chosen as T = 0.1. The measurement of the total pressure is delayed by two sampling instants. Further, the measurements are also sampled at irregular sampling intervals, which are integral multiples of the sampling time T. For the simulations, this scenario is created by randomly dropping samples of the pressure measurement. The measurements may be unavailable for a maximum of two consecutive sampling instants. The initial state of the process is x 0 = [3 1]T. The process is affected by zero-mean random unmeasured disturbances, modeled as additive in the state dynamics, with covariance Q = diag[10−6 10−6 ]. The measurements are corrupted with a zero-mean Gaussian random noise, with covariance R = 0.01. 4.2.1. Results. It has been demonstrated in earlier works6,22 that use of unconstrained state estimation algorithms yields state estimates that are not consistent with their physical limits, even when the measurements are available at regular intervals and without any delay. Hence, in this case, the moving window constrained EnKF (C-EnKF) is used for state estimation. The number of particles chosen for the C-EnKF was Np = 50. Since the measurements are delayed by two sampling instants, the length of the moving window is N = 3. The state estimator was initialized with the following conditions:

RMSE case A

(47)

dpB

Table 4. TE Process: RMSEs of State Estimates for Case A, Only Measurements That Are Not Delayed; Case B, Delayed and Irregular Measurements; and Case C, All Measurements Regular and Not Delayeda state

= −2k1pA 2

dt

The table shows states showing improvements in case B over case A.

no.

(46)

2A → B

The table shows states showing deterioration in case B over case A.

13758


Industrial & Engineering Chemistry Research ⎡ 0.1 ⎤ x̂ 0 = ⎢ ⎥ ⎣ 4.5⎦

⎡ 36 0 ⎤ P0 = ⎢ ⎥ ⎣ 0 36 ⎦

Article

Table 5. Average RMSE Values of pA and pB (50) RMSE (delayed measurements) RMSE (regular measurements)

The process dynamics equations were solved using the standard adaptive Runge−Kutta methods in Matlab 2013a. On a computer equipped with an Intel Core i7-4770K 3.50 GHz processor and 16 GB of RAM, one simulation run using the proposed constrained moving window EnKF algorithm required approximately 17 s. Figures 5 and 6 compare the

pA

pB

0.1415 0.1353

0.2217 0.1925

estimates for 50 runs of this scenario are also reported in Table 5. Similar to the previous example, the RMSE of the state estimates obtained for the scenario when measurements are delayed and irregularly sampled is higher than those obtained for the scenario when the measurements are regularly sampled and available without delay.

5. EXPERIMENTAL CASE STUDY The proposed moving window state estimator was implemented on the data collected from the benchmark heater− mixer experimental setup27 available at the Automation Laboratory, Department of Chemical Engineering, Indian Institute of Technology Bombay. 5.1. Heater−Mixer System. The heater−mixer setup consists of two tanks in series, the schematic of which is shown in Figure 7. Cold water is introduced in tank 1, in which Figure 5. Estimates of pA.

estimates of pA and pB with their respective true values. From Figures 5 and 6, it can be seen that the moving window CEnKF is able to track the states with reasonable accuracy, and also maintain them within their physical bounds, despite the mismatch in the initial value of the state estimates and the irregularly sampled and delayed measurements. The average RMSE values of the state estimates, for 50 runs of the moving window C-EnKF on the same set of data, are given in Table 5. The simulation is also run for the scenario where the measurements are available without any delay and are sampled regularly. Hence, the length of the moving window is reduced to N = 1, which is the conventional C-EnKF.22 It should be noted that the simulations cannot be carried out for the scenario corresponding to case A in section 4.1.1 since there is only one measurement in this case and discarding it would result in a scenario with no measurements. The initial process conditions and state estimator conditions are identical to those in the previous scenario. The RMSE values of the state

Figure 7. Schematic of the experimental two-tank heater−mixer setup.

Figure 6. Estimates of pB. 13759



Article

level loop is closed using a PI controller between h2 and I3. The flow dynamics is perturbed by introducing a 2 mA signal in the set point to h2. The controller output (I3) was available and recorded at every sampling instant. 5.2. Obtaining the Delayed and Irregularly Sampled Measurements. The scenario where certain measurements are delayed and/or irregularly sampled is created as follows: The measurement of water temperature in tank 1 (T1) is delayed by three sampling instants (15 s). Further, the measurement of T1 is sampled at irregular intervals. The entries of tank 2 temperature (T2), too, are randomly dropped to reflect a scenario of irregularly sampled measurements. However, T2 measurements are assumed to be available without any delay. The level in tank 2 (h2) is measured at every sampling instant and is also available without any delay. The unavailability of T1 and T2 measurements is restricted to a maximum of two consecutive sampling instants (10 s). However, for the purpose of evaluating the estimator performance, all the measurements are separately recorded at every sampling instant. It may be noted that the input−output data used in this work are identical to the data used by Bavdekar et al.28 5.3. Results. The moving window EnKF, with Np = 100, was implemented for state and parameter estimation of the heater−mixer setup. The values of the process and measurement noise covariance matrices are taken from the work by Bavdekar et al.28 The values of Q and R and the estimate of the initial states and error covariance are given in Table 6. Since the delay in obtaining the measurement T1 is three sampling instants, the length of the moving window is N = 4.

it is heated using a 4 kW heating coil. The hot water from tank 1 overflows into tank 2, in which it is mixed with another cold water stream that is introduced separately. The water in tank 2 is heated using another 3 kW heating coil. The contents in both tanks are well agitated by using two variable speed agitators. The water exits from the bottom of tank 2 under the action of gravity. The heat inputs, Q1 and Q2, to both tanks can be manipulated independently by changing the current to the thyristor power control units, which can be varied between 4 and 20 mA. The cold water flow rates, F1 and F2, to tank 1 and tank 2, respectively, can be independently manipulated through pneumatic control valves. The water temperatures in each tank (T1 and T2) and the water level in tank 2 (h2) are the measured states. The cold water flow rates, F1 and F2, and the heating inputs, Q1 and Q2, are available as manipulated inputs. The temperature of the inlet water stream (Tc) is a measured disturbance. The setup is interfaced with a computer using a data acquisition module (Advantech 5000 series hardware) and LabView and Matlab software. The first-principles model developed for this system28 is described by the following differential equations: V1

α1Q 1(I1) dT1 = F1(D)(Tc − T1) + dt ρCp

A 2h 2

(51)

dT2 = F1(D)(T1 − T2) + F2(I3)(Tc − T2) dt 1 + [α2Q 2(I2) − 2πr2h2U (T2 − Ta)] ρCp (52)

A2

dh 2 = F1(D) + F2(I3) − 10−3k h2 dt

Table 6. Heater−Mixer Setup: Covariance Matrices and Initial State Estimates

(53)

The parameters α1 and α2 are introduced to account for the changes in heat transfer correlations due to unknown variations in heat loss to the atmosphere from tank 1 and tank 2, respectively. The dynamics of these two variables and the inlet cold water temperature (Tc) are modeled as a random walk model α1, k = α1, k − 1 + wα1, k − 1 (54) α2, k = α2, k − 1 + wα2 , k − 1

(55)

Tc, k = Tc, k − 1 + wTc , k − 1

(56) 2

0.200 1.052 ⎤ ⎥ − 0.133 0.181 ⎥ − 0.002 − 0.006 ⎥ ⎥ 0.782 0.005 3.564 − 0.100 − 0.098 ⎥ − 0.133 − 0.002 − 0.100 0.022 0.043 ⎥ ⎥ 0.181 − 0.006 − 0.098 0.043 0.208 ⎦

R

10−3

⎡177.598 ⎢ − 32.410 ×⎢ ⎢ 0.212 ⎢ ⎣ 70.117

− 32.410 0.212 70.117 ⎤ ⎥ 69.587 0.024 − 5.949 ⎥ 0.024 0.004 0.040 ⎥ ⎥ − 5.949 0.040 147.081⎦

x̂0

[65.0 52.0 0.290 36.0 0.99 0.99]T

P0

diag[0.01 0.01 0.001 0.01 0.002 0.002]

Q

where wα1 ∼ 5(0, qα ), wα2 ∼ 5(0, qα ), and wTc ∼ 5(0, qT ). 1

10−3

⎡137.494 ⎢ ⎢ 0.156 ⎢− 0.690 ×⎢ ⎢ 0.775 ⎢ 0.200 ⎢ ⎣1.052

c

The process noise is treated to arise from unknown sources and is, hence, modeled as additive in the state dynamics. The values of the process and measurement noise covariance matrices (Q and R) are taken from the work by Bavdekar et al.28 The flow rates (F1 and F2) are modeled as a function of the percent of the input current range (D and I3) to the respective control valves. The heat terms (Q1 and Q2) are modeled as a function of the percent of the input current range (I1 and I2) to the respective thyristor power control modules. The algebraic equations that related the milliAmpere signals to the model inputs are given in the Appendix. The sampling time (T) for this experiment was chosen as 5 s. The measurement data are generated by perturbing the system over the entire range of the heat inputs (Q1 and Q2). A pseudo random binary signal (PRBS) of 2 mA was introduced in the current signal D (to manipulate F1). To stabilize the level dynamics in tank 2, the

0.156

− 0.690 0.775

44.388 − 0.018 0.782 − 0.018 0.009 0.005

The state estimates obtained are compared with their corresponding measurements. For the first 1000 s of the experimental data, the estimates of T1 (with delayed and irregularly sampled measurements) are shown in Figure 8. Similarly, estimates of the irregularly sampled state, T2, are shown in Figure 9. The variations in the parameters α1 and α2 are shown in Figure 10. The root mean square error (RMSE) values between the filtered states and their obtained measurements are given in Table 7. For the experimental case, the RMSE for each state estimate is computed as RMSE =

13760

1 Nd

Nd

∑ (yj − xĵ | j)2 j=1

(57)



Article

Figure 8. Heater−mixer setup. Delayed and irregularly sampled tank 1 temperature: comparison of measured and estimated values.

Table 7. Heater−Mixer Setup: RMSE of Error between State Estimates and Measured Values state RMSE

T1 0.7301

T2 0.4425

h2 0.0061

accuracy, despite the delayed and irregularly sampled measurements.

6. CONCLUSIONS In this work, a moving window Bayesian state estimation approach has been developed, which is capable of utilizing irregularly sampled measurements having time varying measurement delays to carry out state estimation. The proposed approach can be used with any unconstrained or constrained recursive Bayesian state estimator. The ability to handle time varying measurement delays in multiple measurements is a distinguishing feature of the proposed approach with reference to the available literature. The difficulties arising from time varying measurement delays are handled by selecting the length of the moving window greater than the maximum possible delay in any of the measurements. In contrast to some of the approaches available in the literature, which require storing the trajectories of the state estimates and associated covariance matrices until a delayed measurement becomes available, the only storage required in the proposed approach is that of the statistical properties of the state estimates at the beginning of the moving window and the inputs and measurements corresponding to the sampling intervals covered by the window. It may be noted that the computation time for this proposed method is a function of the window size, which is dictated by the maximum delay in the measurements. The efficacy of the proposed moving window state estimation scheme has been demonstrated by conducting simulation and experimental studies on three benchmark problems in the literature, viz., the Tennessee Eastman (TE) problem, the gas phase reactor system, and the heater−mixer experimental setup. The ensemble Kalman filter (EnKF) and its constrained version (C-EnKF) have been used for formulating the moving window state estimation schemes. For the TE problem, the performance of the moving window EnKF was

Figure 9. Heater−mixer setup. Irregularly sampled tank 2 temperature: comparison of measured and estimated values.

Figure 10. Heater−mixer setup. Variation in α1 and α2.

Figures 8−10 and the low RMSE values (when compared to the process values) indicate that the proposed moving window EnKF is able to track the process states with reasonable 13761



Article

like to acknowledge WestGrid and Compute/Calcul Canada for providing computational resources to carry out simulations on the Tennessee Eastman problem.

compared with the following two extreme scenarios: (a) when none of the delayed measurements are used to compute the state estimates and (b) when all the measurements are sampled regularly and delays are absent. The simulation results indicate that the use of the delayed measurements in estimation improves the state estimates for the majority of the states when compared with the scenario that ignores delayed measurements. The performance of the state estimator can be further improved if the frequency of measurement sampling is increased, bringing it closer to the performance of the ideal case, where all measurements are available regularly. Qualitatively similar results were obtained for the gas-phase reactor system when the constrained moving window EnKF was used for estimation. In the case of the experimental study, the proposed moving window EnKF was found to track the process states with a reasonable accuracy, despite the delayed and irregularly sampled measurements.

■

(1) Muske, K. R.; Edgar, T. F. Nonlinear State Estimation. In Nonlinear Process Control; Henson, M. A., Seborg, D. E., Eds.; Prentice Hall: Upper Saddle River, NJ, 1997. (2) Julier, S. J.; Uhlmann, J. K. Unscented filtering for nonlinear state estimation. Proc. IEEE 2004, 92 (3), 401−422. (3) Evensen, G. Data AssimilationThe Ensemble Kalman Filter; Springer-Verlag: Heidelberg, 2007. (4) Gordon, N. J.; Salmond, D. J.; Smith, A. F. M. Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEE Proc. Radar Signal Process. 1993, 140 (2), 107−113. (5) Daum, F. Nonlinear filters: Beyond the Kalman filter. IEEE Trans. Signal Process. 2005, 50 (2), 174−188. (6) Rawlings, J. B.; Bakshi, B. R. Particle filtering and moving horizon estimation. Comput. Chem. Eng. 2006, 30 (10−12), 1529−1541. (7) Patwardhan, S. C.; Narasimhan, S.; Prakash, J.; Gopaluni, R. B.; Shah, S. L. Nonlinear Bayesian state estimation: A review of recent developments. Control Eng. Pract. 2012, 20 (10), 933−953. (8) Gudi, R. D.; Shah, S. L.; Gray, M. R. Adaptive multi-rate state and parameter estimation strategies with application to a bioreactor. AIChE J. 1995, 41 (11), 2451−2464. (9) Gopalakrishnan, A.; Kaisare, N. S.; Narasimhan, S. Incorporating delayed and infrequent measurements in extended Kalman filter based nonlinear state estimation. J. Process Control 2011, 21 (1), 119−129. (10) Mutha, R. K.; Cluett, W. R.; Penlidis, A. A new multiratemeasurement-based estimator: Emulsion copolymerization batch reactor case study. Ind. Eng. Chem. Res. 1997, 36 (4), 1036−1047. (11) Mutha, R. K.; Cluett, W. R.; Penlidis, A. On-line nonlinear model-based estimation and control of a polymer reactor. AIChE J. 1997, 43 (11), 3042−3058. (12) Prasad, V.; Schley, M.; Russo, L. P.; Bequette, B. W. Product property and production rate control of styrene polymerization. J. Process Control 2002, 12 (3), 353−372. (13) Rao, C. V.; Rawlings, J. B. Constrained process monitoring: moving-horizon approach. AIChE J. 2002, 48 (1), 97−109. (14) Krämer, S.; Gesthuisen, R. Multirate state estimation using moving horizon estimation. Proceedings of the 16th IFAC World Congress, 2005; International Federation of Automatic Control: Laxenburg, Austria, 2005. (15) López-Negrete, R.; Biegler, L. T. A moving horizon estimator for processes with multi-rate measurements: A nonlinear programming sensitivity approach. J. Process Control 2012, 22 (4), 677−688. (16) Qu, C. C.; Hahn, J. Computation of arrival cost for moving horizon estimation via unscented Kalman filtering. J. Process Control 2009, 19 (2), 358−363. (17) Ungarala, S. Computing arrival cost parameters in moving horizon estimation using sampling based filters. J. Process Control 2009, 19 (9), 1576−1588. (18) Vachhani, P.; Rengaswamy, R.; Gangwal, V.; Narasimhan, S. Recursive estimation in constrained nonlinear dynamical systems. AIChE J. 2005, 51 (3), 946−959. (19) Vachhani, P.; Narasimhan, S.; Rengaswamy, R. Robust and reliable estimation via unscented recursive nonlinear dynamic data reconciliation. J. Process Control 2006, 16 (10), 1075−1086. (20) Kandepu, R.; Foss, B.; Imsland, L. Applying the unscented Kalman filter for nonlinear state estimation. J. Process Control 2008, 18 (7−8), 753−768. (21) Kadu, S. C.; Bhushan, M.; Gudi, R. D.; Roy, K. Modified unscented recursive nonlinear dynamic data reconciliation for constrained state estimation. J. Process Control 2010, 20 (4), 525−537. (22) Prakash, J.; Patwardhan, S. C.; Shah, S. L. Constrained nonlinear state estimation using ensemble Kalman filter. Ind. Eng. Chem. Res. 2010, 49, 2242−2253.

■

ALGEBRAIC RELATIONS FOR THE HEATER−MIXER SETUP The algebraic equations that relate the model inputs to the milliamp signals are as follows. 1. Heat to tank 1: Q 1(I1) = −0.0026I13 + 0.0584I12 + 53.579I1 W

(58)

2. Heat to tank 2: Q 2(I2) = −0.0035I2 3 − 0.782I2 2 + 78.82I2 W

(59)

3. Flow to tank 2: F2(I3) = ( −0.1115I3 4 + 17.2404I33 − 817.6I32 + 18858I3) × 10−10 m 3/s

(60)

4. Flow to tank 1: F1(D) = (1.2676D3 − 108.02D2 + 8166D) × 10−10 m 3/s (61)

The variables I1, I2, I3, and D are expressed as percentage values of the current range (0−100%) to the final control elements.

■

ASSOCIATED CONTENT

S Supporting Information *

Table of steady state values of the process states of the Tennessee Eastman problem; table of PI controller parameters used for stabilizing the unstable loops in the TE problem. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. † A condensed version of this work was presented at the 18th IFAC World Congress, Milan, Italy, in August 2011.

■

ACKNOWLEDGMENTS The authors would like to acknowledge the financial support from the NSERC-Matrikon-Suncor-iCORE Industrial Research Chair program in Computer Process Control at the University of Alberta for carrying out this work. The authors would also 13762



Article

(23) Koläs, S.; Foss, B. A.; Schei, T. S. Constrained nonlinear state estimation based on the UKF approach. Comput. Chem. Eng. 2009, 33 (8), 1386−1401. (24) Downs, J. J.; Vogel, E. F. A plant-wide industrial process control problem. Comput. Chem. Eng. 1993, 17 (3), 245−255. (25) Ricker, N. L.; Lee, J. H. Nonlinear modeling and state estimation for the Tennessee Eastman challenge process. Comput. Chem. Eng. 1995, 19 (9), 983−1005. (26) Karra, S.; Shaw, R.; Patwardhan, S. C.; Noronha, S. Adaptive model predictive control of multi-variable time varying systems. Ind. Eng. Chem. Res. 2008, 47 (8), 2708−2720. (27) Srinivasrao, M.; Patwardhan, S. C.; Gudi, R. D. Nonlinear predictive control of irregularly sampled multirate systems using nonlinear blackbox observers. J. Process Control 2007, 17 (1), 17−35. (28) 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. J. Process Control 2011, 21 (4), 585−601.

13763


A Moving Window Formulation for Recursive Bayesian State

Recommend Documents