Probability-Grid Filtering for Chemical Batch Reactors

Copyright © IFAC 12th Triennial World Congress, Sydney, Australia, 1993

PROBABILITY-GRID FILTERING FOR CHEMICAL BATCH REACTORS P. Terwlesch and M. Agarwal TCL, EidgenOssische Technische lIochschule , CH-8092 Zurich, Switzerland

distributions are a priori known to be not normal. For example, modeling the purity of a substance through a normal distribution, with a nominal value of 99% and a standard deviation of 5%, would allow " purity" to exceed 100%, and does not incorporate the physically obvious constraint information. Or, the possibility of competing reaction mechanisms corresponding to two different values of a rate constant , depending on presence or absence of a non-detectable impurity, dictates an initial distribution with two maxima. Also, due to the ever-changing operating point, batch-reactor models typically involve significant nonlinearities coupled with poor initial estimates of states and parameters . A nonlinear filter is advocated for such situations, which overcomes the limitations of the EKF at the expense of higher computation requirement. For weakly nonlinear systems with initial distributions close to the normal, the EKF is advantageous on account of its computational efficiency. The nonlinear filter uses a grid discretization of the state probability distribution, which is updated every sampling instant using Bayes' rule. Consequently, it is referred to here as the Probability- Grid Filter (PGF). A vector representation of the grid, for reasons of computational efficiency, is responsible for the method being termed the p-vector approach in the literature [2] . The advantage of PGF is that the grid allows good representation of non-normal distributions, and the linearization errors are kept small as linearization is performed at each grid point rather than only at the expected value as in the EKF. Consequently, better estimates and quicker convergence are achieved, albeit at the cost of increased computational effort.

Abstract: Use of a state-space grid-form of the Bayes' rule is advocated over the popular Extended Kalman Filter (EKF) for state estimation given nonnormal probability distributions or strongly nonlinear processes, as is conunon in batch reactors. The EKF cannot account for a non-normal initial distribution or the distortion of a normal distribution through nonlinear model, and is subject to relatively large error due to linearization around only the expected state value at the current time . A grid discretization of the state space allows for arbitrary initial distributions, faithful processing of the distribution through non linear model, and smaller linearization error as linearization is performed around each grid point at the current time. The Probability-Grid Filter (PGF) used here thus exhibits faster convergence to better estimates and narrower distributions, increased robustness to model-plant mismatch, reduced sensitivity to tuning parameters, and extenuated susceptibility to instability. A simulation example of a typical batch reactor is used to demonstrate the benefits of the PGF over the EKF when the process is significantly nonlinear or the initial distribution is significantly non-normal. The benefits of the PGF come at the cost of increased computational effort. Keywords: batch reactor, state estimation, nonlinear filtering

1

Introduction

For safe and optimal operation of chemical batch reactors, estimates of on-line unmeasured quantities such as concentration are useful, as well as their probability distributions. When the dynamic state model is not perfect or parameters or initial states are inaccurately known, then a filter is required that also uses the measurements and a measurement model.

2

Probability-Grid Filtering

The following description of probability-grid filtering is based mainly on a review by Sorenson [2], who applied this technique to target-tracking problems.

For this purpose, the Extended Kalman Filter (EKF) has received most attention and many successful applications. A key assumption in the EKF is that the initial and the current state uncertainties have to be described by normal distributions, characterized by means and variances. This limitation precludes the EKF from incorporating possibly contrary a priori information about the initial distribution, as well as from processing even a normal initial distribution correctly through the system nonlinearities. 1 Moreover, for nonlinear models and poor initial estimates, the error due to linearization around the estimated rather than the true value can lead the EKF to instability.

Consider a discrete-time model of the plant,l with the following structure: XHl Ylo

f(XIo)

=

h(XIo)

+ Wlo + Vlo

(1) (2)

where k is discrete time; Xlo state; Ylo measurement ; additive noises on state and measurement , respectively; and f, h known nonlinear functions. In absence of better information,' Wlo and Vlo are assumed to be zero-mean white noises described by Wlo, Vlo

lUse of a discrete-time model is not a major restriction since measurements are almost always sampled and continuous models can be converted using explicit discretization techniques. 'The algorithm allows use of any noise distribution.

Batch reactors involve situations where the initial 1 A normal distribution processed through a nonlinear system will yield a non-normal distribution.

57

normal probability-density functions N(w, 0, Q) and N(v, 0, R) with known (or user-specified) Q and R, respectively.4 Here, N(Ot,/3, P) stands for the normal probability-density function of a variable Ot with mean /3 and covariance P.

The estimated value i. of the state Xr. is calculated for the PGF as: m

i. = E{x.}:::::

.=1

The measurement update based on a new measurement YIo, yields the conditional posterior distribution for the state Xlo using Bayes' rule:

where E { ... } denotes the expectation operator.

3

(3) The calculation of p(XIo 1 YIo-t} from p(XIo-1 1 YIo-t} available in the previous instant is treated below. The conditional density for the measurement, p(YIo 1 Xlo), in eq. 3 can be expressed as N(YIo,h(xlo),R).

"10.,

Cl

''V

3.1

(4)

Simulated Process

'

"10.)

Consider an isothermal batch reactor in which a second-order reaction produces a substance B from A. B reacts further giving C in a consecutive firstorder reaction.

p(YIo 1

= 1, ... , m. The normalizing constant used for numerical scaling.

Vi

Simulation Example

This section illustrates the main features of the PGF by comparing its performance to that of the EKF for a nonlinear process. Two case. are presented. First, the case of a normal initial distribution is examined in order to isolate the effect on the two filters of the system nonlinearity. The second case demonstrates the ability of the PGF to benefit from incorporation of a priori knowledge such as physically constrained irUtial estimates and multi-modal initial distributions.

Since computer representation of a probability distribution that cannot be expressed analytically requires discretization, the range of possible values for the state vector Xlo is represented by grid points i = 1, ... , m. In this grid representation, the measurement update of eq. 3 becomes

P(Il •• 1Y.) =p(" •• 1y.-t}N(Yr.,h("r..),R)

I>(" •. )Ilr..

Cl

can be

The dynamic update of P("r.-1,. 1 Y.-1) to p("r.. 1 Y.-1) uses the model in eq. 1 as described in lazwinski [1], p. 34 and 86:

" •• =f("._l,.) p

(Illo. 1Yr.-1, Q = 0)

Vi=l, ... ,m

(5)

= P("1o-1,.I Y.-d IIJ(" •• )11

(6)

The reactions are described by the following differential equations:

d[A]

where this update step is conditional upon Q = 0, and IIJ("r.;)11 is the absolute value of the determinant of the lacobian for the inverse mapping f-1. The lacobian of the inverse can be obtained either as the derivative of an explicit inversion of the state equations or as the inverse of the forward lacobian using the Inverse-Function Theorem.

d[B]

m

(9)

dt

It is assumed that the only available measurement is (10)

If the user specifies Q t:- 0, the following convolution is to be performed additionally:5

P("Io. 1Yk-1,Q t:- 0)

(8)

dt

For the simulation, the following values have been used unless otherwise stated: 1:1 = 0.02, kl = 0.01, Cl = 5, Cl = 10, [A](t = 0) = 1 and [B](t = 0) = O. Zero-mean white noise with co variance 0.04 is added to the measurement tI in eq. 10.

= (7)

The goal is to estimate the concentrations [A] and [B], with prior knowledge of equations 8-10 as well as

;=1

At k = 0, the conditional distribution P("Io. 1y.-t} is to be replaced using the a priori known distribution P(llo.)' At each sampling instant, the entire update comprises use of equations 4, 5,6 and 7. Compared to the other equations, eq. 7 involves relatively expensive computation. Each of the m grid points " •• requires m calls to the normal-density function, leading to a total of m 2 function evaluations. This means that the PGF is computationally expensive for Q t:- o. This expense could be reduced by a choice of more efficient grid-point loactions, that yield similar accuracy using a smaller number of points. 6

of k l , Cl, Cl, and the measurement-noise covariance, but with 1:1 and the irUtial concentration not being accurately known.

3.2

Model used for Estimation

To allow estimation of both concentrations and the unknown rate constant k 1 , the state equations (8, 9) are extended with zero dynamics for k 1 • Writing Zl = [A], Zl [B] and Z3 = k1 a first-order Euler discretisation with ilt = 1 yields:

=

However, noise models providing an accurate statistical description of the true process are rarely known. 4 Time-varying f., h., Q. and Rio are admissible. 5The user may choose to specify Q "I 0 when, for example, the state noise or the model-mismatch is expected to be non-negligible. 6 For example, a grid contraction f(!Jr._1,.) --+ !J •• could be incorporated instead of eq. 5, although choice and implementation of a suitable grid transformation is non-trivial.

Zl,Io+1

Zl,. -

Zl,.H

Z2,.

Z3,.+1

tI.

=

ilt[Z3,.zt.]

+ ilt[Z3,.Ztr. -

k 2 z 2 ,.]

(11)

Z3,Io

C1Z~,.

+ ClZ~,r.

k 2 , Cl, C2 and the measurement noise covariance are known at true values. This system is fully observable for almost any state.

58

ii[?;: ; ; o

~O

100

1.$0

200

~~p-J-~" h_'V_/ H :--1

1 ].50

0

SO

100

1.$0

ZOO

~~ :~b;'" -Jj' H: o

so

100

1$0

200

2.SO

0

0.015

~t>-'-

Figure 1: Concentrations [A] and [B] and measurement y. The inverse mapping r- l is obtained in a straightforward manner by solving eq. 11 for x •. The absolute value of the lacobian's determinant, evaluated at the new state, is thus computed to be:

-G.01

0

$0

100

_

"

j

:

150

200

i---

.

__

100

Figure 2: Estimation error - PGF, - - EKF.

150

z. -

200

Xi,

i

=

1 250

1 2$0

I, '" 3.

extra degrees of freedom in its distribution representation_ Then, after k ::::: 30 until about k ::::: 160, the PGF performs better because at each instant it linearizes around not only the expected value but rather each individual value of the state grid. Consequently, the linearization error is smaller. After k ::::: ISO, EKF gives good estimates, having converged towards the true values so that the linearization error is small. Note that although EKF estimates are close to the true values also at k ::::: 30, no convergence is achieved due to a high filter gain caused by the still large uncertainty of the estimates. That the large uncertainty is not due to excessively high Po was confirmed by use of other Po tunings, which could not give any improvement.

As this direct computation may sometimes be hard to perform, one can compute IIJ(x.+dll also through the Inverse-Function Theorem using the identity 1 det(AA- l ) = det(A)det(A- l ). For the state noise, Q = 0 is the default setting but is increased when mismatch is treated.

=

3.3

.50

210

Case of Normal Initial Distribution

In the first case to be examined, a normal initial distribution is assumed for both the filters. For PGF: p(xo) = N(xo,xo,P o) where Po is the initial covariance matrix in EKF: Po = diag{10- 1 , 2.5 . 1O- 3 ,1O- 4 }. This is to investigate the different processing by the two filters of the initial distribution through the nonlinear process model. Both filters start with incorrect initial estimates Zl,O = 0.90, Z2,O = 0.05, Z3,O = 0.03 and use Q = 0 and R = 0.04. The PGF collocates each state with 20 points T covering Zi ,O ± 2.5J[P o]... Concentrations [A], [B] and measurement y for the simulated run are shown in fig. 1. The corresponding estimation errors are given in fig . 2. It can be seen that overall the PGF performs better than the EKF, with the exception of the tail-end period starting at k ::::: 160 where the PGF gives marginally inferior estimates, being limited by the coarseness of the used grid . However, this is not a major drawback, as its bias is already small and could be improved at will at the expense of more computation.'

It is important to note that the PGF performs significantly better during the time interval (namely k ::::: 30 to k ::::: 160) that constitutes the critical phase for batch optimization. Prior to this phase, not enough measurement information is available to switch from the off-line-optirnized input trajectory to an on-line updated one, although the end-result (e .g., final yield) is maximally sensitive to input corrections. After the critical phase, although more accurate estimates of the states are available, the end-result often becomes relatively insensitive to input corrections. Fig. 3 illustrates the evolution with time of the marginal probability distributions p(Zi) for both the PGF and the EKF. The nine plots are organized such that within each row time k is constant, and each column represents a state Zi. On each plot a dashed verticalline indicates the location of the true value. Note that due to the employed point-mass approximation,i not the area underneath each curve but rather the sum of the discrete points representing each curve is unity. While both filters start out with the same initial distribution, the EKF comes close to the true value at k = 25 but overshoots it afterwards, whereas the PGF does not fully approach the true value at k = 25 but rides it at k = 50. As expected, the PGF distribu-

In the beginning, until k ::::: 30, the PGF does not yet exhibit its superiority as not enough measurements are available for it to take advantage of the T A collocation with (12,12,10) points for Zl, Zl and Z3, respectively, has yielded similar estimation quality for the exrunple presented here, whereas in other examples 15 to 20 collocation points were found to be a good compromise between computational effort and accuracy. sin the (20,20,20)-collocation-point case, 0.2 seconds per sampling time were needed on a VAX 6000 computer, whereas the (12,12,10) collocation required only 0.05 seconds.

iTo allow for a comparison of the EKF density to the PGF point-mass approximation, the EKF distribution also is represented by its point-mass approximation. In each plot the EKF curve is depicted for only the limited Zirange used by the grid of the PGF.

59

"0

k=O

10

. . ..

u

I

:

.~

U

u~"L ........

t.I

I

:u ~~. ~ ~

...

u

...

uu

:[J.,. =[IJ.k=~ I ...

:

&AUUU

x,

u

I

......

JL

UIUUO'.U

x,

L
U

I

: U

:[J" :J.L

I I

1..

"'DJI [JI WI ~I ~.. L. .J\A .~

.: ...

x,

k-20

::

.

u

k-20

U

j"~

..

,"

:[J =40

UD~&AlUI

[[l I

... )V\... ...

.. r\

....

k=O

I

u

k-20

...

u

='r . 1ll u[J :UU . ......,

...

U

lA...

...

..

=40

Ut

k=40

..lL

. . , '. . . . . . .

x.

Figure 3: Evolution of marginal probability p(Zi). PGF (solid), EKF (dotted), true value (dashed).

UU

......

x,

UD

GAl

IL06

x,

Figure 4: Influence of non-normal initial distribution. Evolution of marginal probability p(Zi). PGF (solid), EKF (dotted), true value (dashed) .

tion is narrower and asymmetrical due to its better processing of the system nonlinearity.

tributions are almost centered around the true values. The EKF's estimates are poorer and have a larger variance. Although starting with the same mean value, the EKF performs significantly worse in estimating the rate constant Z3. Note that what appears to be final bias on the estimates actually vanishes for k > 150 with both filters.

In the above run, since the initial estimates are relatively close to the true values, the nonlinearities in the state equations do not severely disadvantage the EKF. The crucial nonlinearity responsible for the linearization error resides in the measurement equation (10). This observation is confirmed by replacing, in both the process and the estimation model, the nonlinear measurement equation by the linear form y = c(1 - Zl - Zl) for which the EKF and the PGF yield practically identical estimates (not shown) . This shows that when linearization errors are small and the probability distribution is close to normal, then the performance of the PGF reduces to that of the EKF, which is computationally advantageous. On the other hand, even with a linear measurement equation, when initial estimates are poorer so that the state nonlinearities become significant, then the PGF retains its superior performance (not shown).

3.4

I I I

a. . . .

UW I k-2~ UDJIk-2~ lA

••

.. .. ~

I

&. . . . . .

k=O

U[] =O

:

u

:

k=O

I

Model-plant mismatch in the form of a side-reaction unknown to the estimator has also been investigated. The results can not be presented here due to space constraints.

4

Conclusion

A nonlinear filtering technique has been applied to the state estimation problem in chemical batch reactors. Linearization is performed at each point on a discrete grid used to represent the probability distributions, thereby considerably reducing the linearization error compared to the extended Kalman filter, although at the expense of more computational effort. In cases with poor initial estimates, strong model nonlinearities, or non-normal initial distributions, the probability-grid filter is shown to outperform the popular extended Kalman filter by giving a narrower state distribution closer to the true value.

Case of Non-normal Initial Distribution

To illustrate the incorporation of a priori knowledge into the PGF, Zl and Z2 are now interpreted as purities, so that they are constrained as Zl,Z2 E [0,1). Also, the rate-constant distribution p( Z3) is assumed to have two maxima, which can arise, e. g., when the reaction is sensitive to presence or absence of a nondetectable impurity. While the PGF can use this type of information to improve the estimates, the EKF is unable to incorporate it. The first row of fig. 4 shows the initial distributions and the starting values (dashed). Note that the true values for Zl and Zl are intentionally simulated on the constraint boundaries, which is the worst case for the PGF. Even then, however, the PGF performs visibly better than the EKF . While at k = 20 the impact of the bi-modal initial distribution for Z3 pervades all states,10 at k 40 the second peak has vanished already and the PGF dis-

Acknowledgment: First author acknowledges financial support from the Schweizerischer N ationalfonds zur Forderung der wissenschaftlichen Forschung.

References [1) A.H. lazwinski. Stocha6tic Proce66e6 and Filtering Theory. Academic Press, 1970. (2) H.W. Sorenson. Recursive estimation for nonlinear dynamic systems . In J.C . Spall, editor, Baye6ian analY6i6 of time 6erie6 and dynamic model" M. Dekker, 1988.

=

lOThis can be interpreted as hypothetical system dynamics corresponding to the incorrect mode of the "'3distribution, which has not yet been entirely rejected.

60