Kalman Filtering

575

Chapter 41

Kalman Filtering 41.1 Introduction Linear regression and non-linear regression methods to fit a model through a number of data points have been discussed in Chapters 8, 10 and 11. In these regression methods all data points are collected first followed by the estimation of the parameters of a postulated model. The validity of the model is checked by a statistical evaluation of the residuals. Generally, the same weight is attributed to all data points, unless a weighing procedure is applied on the residuals, e.g., according to the inverse of the variance of the experimental error (see Section 8.2.3.2). In this chapter we discuss an alternative method of estimating the parameters of a model. There are two main differences from the regression methods discussed so far. First, the parameters are estimated during the collection of the data. Each time a new data point is measured, the parameters of the model are updated. This procedure is called recursive regression. Because at the beginning of the measurement-estimation process, the model is based on a few observations only, the estimated model parameters are imprecise. However, they are improved as the process of data collection and updating proceeds. During the progress of the measurement-estimation sequence, more data points are measured leading to more precise estimates. The second difference is that the values of the model parameters are allowed to vary during the measurement-estimation process. An example is the change of the concentrations during a kinetics experiment, which is monitored by UV-Vis spectrometry. Multicomponent analysis is traditionally carried out by measuring the absorbance of a sample at a number of wavelengths (at least equal to the number of components which are analyzed) and calculating the unknown concentrations (see Chapter 10). These concentrations are the parameters of the model. The basic assumption is that during the measurement of the absorbances of the sample at the selected wavelengths, the concentrations of the compounds in the mixture do not vary. However, if the measurements are carried out during a kinetics experiment, the concentrations may vary in time, and as a result the composition of the sample varies during the measurements. In this case, we cannot simply estimate the unknown concentrations by multiple linear regression as explained in Chapter 10. In order to estimate the concentrations of the

576

compounds as a function of time during the data acquisition, two models are needed, the Lambert-Beer model and the kinetics model. The Lambert-Beer model relates the measurements (absorbances) to the concentrations. This model is called the measurement model {A =J{a,c)). The kinetics model describes the way the concentrations vary as a function of time and is called the system model (c = J{k,t)). In this particular instance, the system model is an exponential function in which the reaction rate k is the regression parameter. The terms * systems' and 'states' are associated with the Kalman fdter. The system here is the chemical reaction and its state at a given time is the set of concentrations of the compounds at that time. The output of the system model is the state of the system. The measurement and system models fully describe the behaviour of the system and are referred to as the state-space model. Thus the system and measurement models are connected. The parameters of the measurement model (in our example, the concentrations of the reactants and the reaction product) are the dependent variables of the system model, in which the reaction rate is the regression parameter and time is the independent variable. Later on we explain that system models may also contain a stochastic part. It should be noted that this dual definition implies that two sets of parameters are estimated simultaneously, the parameters of the measurement model and the parameters of the systems model. In this chapter we introduce the Kalman filter, with which it is possible to estimate the actual values of the parameters of a state-space model, e.g., the rate of a chemical reaction from the evolution of the concentrations of a reaction product and the reactants as a function of time which in turn are estimated from the measured absorbances. Let us consider an experiment where a flow-through cell is connected to a reaction vessel, in which the reaction takes place. During the reaction, which may be fast with respect to the scan speed of the spectrometer, a spectrum is measured. At r = 0 the measurement of the spectrum is started, e.g., at 320 nm in steps of 2 nm per second. Every second (2 nm) estimates of the concentrations of the compounds are updated by the Kalman filter. When arrived at the end of the spectral range (e.g. 600 nm) we may continue the measurement process as long as the reaction takes place, by reversing the scan of the spectrometer. As mentioned before, the confidence in the estimates of the model parameters improves during the measurement-estimation process. Therefore, we want to prevent a new but deviating measurement to influence the estimated parameters too much. In the Kalman filter this is implemented by attributing a lower weight to new measurements. An important property of a Kalman filter is that during the measurement and estimation process, regions of the measurement range can be identified where the model is invalid. This allows us to take steps to avoid these measurements affecting the accuracy of the estimated parameters. Such a filter is called the adaptive Kalman fdter. An increasing number of applications of the Kalman filter

577

has been published, taking advantage of the formulation of a systems model which describes the dynamic behaviour of the parameters in the measurement model and exploiting the adaptive properties of the filter. In this chapter we discuss the principles of the Kalman filter with reference to a few examples from analytical chemistry. The discussion is divided into three parts. First, recursive regression is applied to estimate the parameters of a measurement equation without considering a systems equation. In the second part a systems equation is introduced making it necessary to extend the recursive regression to a Kalman filter, and finally the adaptive Kalman filter is discussed. In the concluding section, the features of the Kalman filter are demonstrated on a few applications.

41.2. Recursive regression of a straight line Before we introduce the Kalman filter, we reformulate the least-squares algorithm discussed in Chapter 8 in a recursive way. By way of illustration, we consider a simple straight line model which is estimated by recursive regression. Firstly, the measurement model has to be specified, which describes the relationship between the independent variable x, e.g., the concentrations of a series of standard solutions, and the dependent variable, 3;, the measured response. If we assume a straight line model, any response y^ is described by: yi = ^0 + ^1 ^i + ^i

bg and b^ are the estimates of the true regression parameters % and Pj, calculated by linear regression and e^ is the contribution of measurement noise to the response. In matrix notation the model becomes: y. = xj h + e^ where x- is a [2x1] column vector 1

and b is [2x1] column vector

h A. A recursive algorithm which estimates PQ and pj has the following general form: New estimate = Previous estimate + Correction

578

After each new observation, the estimates of the model parameters are updated (= new estimate of the parameters). In all equations below we treat the general case of a measurement model with p parameters. For the straight line model /? = 2. An estimate of the parameters b based ony - 1 measurements is indicated by b(/ - 1). Let us assume that the parameters are recursively estimated and that an estimate b(/ - 1) of the model parameters is available fromj - 1 measurements. The next measurement y(j) is then performed at x(/), followed by the updating of the model parameters to b(/). The first step of the algorithm is to calculate the innovation (I), which is the difference between measured y(j) and predicted response y(j) at x(/). Therefore, the last estimate b(/ - 1) of the parameters is substituted in the measurement model in order to forecast the response y(j), which is measured at x(j):

y(j) = bo(j-l) + b,(j-l)x(j)or >;(/•) = xT(/*)b(/'-l)

(41.1)

The innovation /(/) (not to be confused with the identity matrix I) is the difference between measured and predicted response at x(j). Thus /(/) = y(j) - yij)The value of the innovation is used to update the estimates of the parameters of the model as follows: b(7) = b ( 7 - l ) + k ( ; ) / ( j ) px\

px\

pxl

(41.2)

1x1

Equation (41.2) is the first equation of the recursive algorithm. k(/) is a/7xl vector, called the gain vector. Looking in more detail at this gain vector, we see that it weights the innovation. For k equal to the null vector, b(/) = b(/ - 1), leaving the model parameters unadapted. The larger the k, the more weight is attributed to the innovation and as a consequence to the last observation. Therefore, one may intuitively understand that the gain vector depends on the confidence we have in the estimated parameters. As usual this confidence is expressed in the variancecovariance matrix of the parameters. Because this confidence varies during the estimation process, it is indicated by P(/) which is given by:

nj)-One expects that during the measurement-prediction cycle the confidence in the parameters improves. Thus, the variance-covariance matrix needs also to be updated in each measurement-prediction cycle. This is done as follows [1]: P(7)=P(;-l)-k(;)xT(7)P(7-l) pxp

pxp

pxl

Ixp

pxp

(41.3)

579

Equation (41.3) is the second equation of the recursive filter, which expresses the fact that the propagation of the measurement error depends on the design (X) of the observations. Once the fitting process is complete the square root of the diagonal elements of P give the standard deviations of the parameter estimates. Key factor in the recursive algorithm is the gain vector k, which controls the updating of the parameters as well as the updating of the variance-covariance matrix P. The gain vector k is the most affected by the variance of the experimental error r(j) of the new observation y(j) and the uncertainty PO - 1) in the parameters b(/ - 1). When the elements of PO* - 1) are much larger than r(/), the gain factor is large, otherwise it is small. After each new measurement j the gain vector is updated according to eq. (41.4) k(7) = P ( y - l ) x ( y ) ( x T ( y ) P ( y - l ) x ( ; ) + r ( ; ) ) - i pxl

pxp

pxl

\xp

pxp

px\

(41.4)

1x1

The expression x^(/)P(/ - l)x(/) in eq. (41.4) represents the variance of the predictions, y(j), at the value x(j) of the independent variable, given the uncertainty in the regression parameters P(/). This expression is equivalent to eq. (10.9) for ordinary least squares regression. The term r(j) is the variance of the experimental error in the response y(j). How to select the value of r(j) and its influence on the final result are discussed later. The expression between parentheses is a scalar. Therefore, the recursive least squares method does not require the inversion of a matrix. When inspecting eqs. (41.3) and (41.4), we can see that the variancecovariance matrix only depends on the design of the experiments given by x and on the variance of the experimental error given by r, which is in accordance with the ordinary least-squares procedure. Typically, a recursive algorithm needs initial estimates for b(0) and P(0) to start the iteration process and an estimate of r(j) for all j during the measurementestimation process. When no prior information on the regression parameters is available b(0) is usually set equal to zero. In many textbooks [1], it is recommended to choose for P(0) a diagonal matrix with very large diagonal elements, expressing the large uncertainty in the chosen starting values b(0). As explained before, r(j) represents the variance of the experimental error in observation y(j). Although the choice of P(0) and r(j) introduces some arbitrariness in the calculation, we show in an example later on that the gain vector (k) which fully controls the updating of the estimates of the parameters is fairly insensitive for the chosen values r(/) and P(0). However, the obtained final value of P depends on a correct estimation of r(j). Only when the value of r(j) is equal to the variance of the experimental error, does P converge to the variance-covariance matrix of the estimated parameters. Otherwise, no meaning can be attributed to P. In summary, after each new measurement a cycle of the algorithm starts with the calculation of the new gain vector (eq. (41.4)). With this gain vector the variance-

580

Iteration step

Fig. 41.1. Evolution of the innovation during the recursive estimation process (iteration steps) (see Table 41.1).

covariance matrix (eq. (41.3)) is updated and the estimates of the parameters are updated (eq. (41.2)). By monitoring the innovation sequence, the stability of the iteration process can be followed. Initially, the innovation shows large fluctuations, which gradually fade out to the level expected from the measurement noise (Fig. 41.1). An innovation which fails to converge to the level of the experimental error indicates that the estimation process is not completed and more observations are required. However, P(0) also influences the convergence rate of the estimation process as we show with the calibration example discussed below. By way of illustration, the regression parameters of a straight line with slope = 1 and intercept = 0 are recursively estimated. The results are presented in Table 41.1. For each step of the estimation cycle, we included the values of the innovation, variance-covariance matrix, gain vector and estimated parameters. The variance of the experimental error of all observations y is 25 10"^ absorbance units, which corresponds to r = 25 10"^ au for a l l / The recursive estimation is started with a high value (10^) on the diagonal elements of P and a low value (1) on its off-diagonal elements. The sequence of the innovation, gain vector, variance-covariance matrix and estimated parameters of the calibration lines is shown in Figs. 41.1^1.4. We can clearly see that after four measurements the innovation is stabilized at the measurement error, which is 0.005 absorbance units. The gain vector decreases monotonously and the estimates of the two parameters stabilize after four measurements. It should be remarked that the design of the measurements fully defines the variance-covariance matrix and the gain vector in eqs. (41.3) and (41.4), as is the case in ordinary regression. Thus, once the design of the experiments is chosen

581 TABLE 41.1 Recursive estimation of the parameters BQ and Z?^ of a straight line (see text for symbols; y are the measurements). InitiaUsation: bT(0) = [0 0]; r = 25 10-^; P^^(0) = P2,2(^)= 1000; Fi 2CO) = PiA^) = 1 J

x(/')

y(j)

y(j)

I(j)

h(j)

1

m

1 0.1

0.101

0

0.101

0.09999 0.01010

0.99000 0.09998

2

1 0.2

0.196

0.102

0.094

0.0060 0.9500

3

1 0.3

0.302

0.291

0.011

4

1 0.4

0.403

0.401

5

1 0.5

0.499

6

1 0.6

0.608

7

1 0.7

0.703

0.706

8

1 0.8

0.801

9

1 0.9

10 1 1.0

P(/") 9.8990 -98.99

-98.99 989.9

-1.0000 9.99995

1.37 10"" -8.09 lO""

-8.16 10-" 5.249 10"'

-0.0024 1.0072

-0.7307 5.2029

-6.0 10"' -2.61 10-"

-2.62 10-" 1.303 10"'

-0.002

-0.0032 1.0138

-0.5254 3.0754

3.73 10-' -1.26 10-"

-1.26 10-" 5.06 10^

0.504

-0.005

-0.0012 1.0042

-0.4085 2.0236

2.68 10-' -7.4 10"'

-7.41 10"' 2.49 10""

0.601

0.007

-0.0035 1.0138

-0.334 1.433

2.10 10-' -4.88 10'

-4.89 10"'

-0.003

-0.0025 1.0104

-0.2838 1.0694

1.7 10"' -3.5 10'

-3.5 10"' 8.78 10"'

0.806

-0.005

-0.0014 1.0064

-0.2468 0.8290

1.5 10' -2.6 10'

-2.6 10"' 5.84 10'

0.897

0.904

-0.007

-0.0002 1.0015

-0.2185 0.6618

1.27 10' -2.02 10'

-2.02 10"' 4.08 10"'

1.010

1.002

0.0082 -0.0014 1.0060

-0.1962 0.5408

1.13 10' -1.62 10'

-1.62 10' 2.97 10"'

1.41 10^

(x(j),j = 1, ..., n), one can predict how the variance-covariance matrix behaves during the iterative process and decide on the number of experiments or decide on the design itself. This is further explained in Section 41.3. The relative insensitivity of the estimation process for the initial value of P is illustrated by repeating the calculation with the diagonal elements P(l,l) and P(2,2) of P set equal to 100 instead of equal to 1000 (see Table 41.2). As can be seen from Table 41.2 the gain vector, the variance-covariance matrix and the estimated regression parameters rapidly converge to the same values. Also using unrealistically high values for the experimental error (e.g. r(j) = 1) does not affect the convergence of the gain factors too much as long as the diagonal elements of P remain high. However, we also

582

Iteration number

Fig. 41.2. Gain factor during the recursive estimation process (see Table 41.1). lt+4 1E+2

b 1E+0 \ \

P(2,2)

1E-2

\pair'"---~.^___^ 1E-4

1

8

9

10

Iteration number

Fig. 41.3. Evolution of the diagonal elements of the variance-covariance matrix (P) during the estimation process (see Table 41.1).

observe that P no longer represents the variance-covariance matrix of the regression parameters. If we start with a high r(j) value with respect to the diagonal elements of P(0) (e.g. 1:100), assuming a large experimental error compared to the confidence in the model parameters, the convergence is slow. This is indicated by comparing the innovation sequence for the ratio r(j) to P(0) equal to 1:1000 and 1:100 in Table 41.2. In recursive regression, new observations steadily receive a lower weight, even when the variance of the experimental error is constant (homoscedastic). Consequently, the estimated regression parameters are generally not exactly equal to the values obtained by ordinary least squares (OLS).

583 Q.
i

1

^—"

o 1 Q)

C

"(0 o 0) (A

0.8 0.6

:

0.4

-

•

slope

1 / /

0.2 u

intercept

^ / ^ ,,. 1

1

,1

1

,

1

, 1.

,

1

1

1

10 Iteration number

Fig. 41.4. Evolution of the model parameters during the iteration process (see Table 41.1).

This straight line example exemplifies the fact that by exploiting the recursive approach alternative calibration practices are possible. The standard practice in calibration is to measure calibration standards prior to the analysis of the unknown samples and to construct a calibration line. Thereafter one starts the analysis of a sample batch. To check whether the analysis system performs well, each batch of samples is preceded by a check sample with a known concentration. The results of the check sample are plotted in a control chart (see Chapter 7). When the value exceeds the warning or action limits, one may decide to re-establish the calibration line by remeasuring a new set of calibration standards. This procedure disregards the knowledge on the calibration factors already available from the previous calibration run. By applying the recursive algorithm, it is in principle sufficient to remeasure a limited set of calibration standards and to update the estimates of the calibration factors already available from the previous calibration step. However, as can be seen from the example given in Table 41.1, this would mean that the weight attributed to the new measurements steadily decreases (small gain vector). This leads to the undesirable situation where the calibration factors are hardly updated. However, it should be bome in mind that the reason to decide to recalibrate was that during the measurement of the unknown samples, calibration factors may be changed. As a result, the uncertainty in the calibration factors were judged to be unacceptably large, requiring recalibration. Expressed in terms of a system, the system state (here slope and intercept) has not been observed for some time. Because this system state varies in time (in an unknown deterministic or stochastic way), the uncertainty about the system state increases and the state has to be observed again. This uncertainty is taken into account in the Kalman filter by the system equation, which is discussed in the next section. One of the methods to

584 TABLE 41.2 Influence of initial values of the diagonal values of the variance-covariance matrix P(0) and the variance of the experimental error on the gain vector and the Innovation sequence (see Table 41.1 for the experimental values, y) Pu(0)

r=\

r=25 10-6 kl 1

Pi,i(0) == P2 2(0) = 100

= P2,2<[0)=1000

k2

kl

I

r=l

r = 25 10"^ k2

kl

I

0.99 0.10

0.101

0.98 0.10

0.101

0.99

k2

k2

kl

I

0.11

0.101

0.98 0.11

I 0.101

2

-1.0

10.0

0.094

-0.75 8.31

0.094

-1.0

9.99

0.094

0.00 3.33

0.095

3

-0.73 5.20

0.011

-0.62 4.76

0.035

-0.66

4.99

0.011

-0.33 3.31

0.105

4

-0.53 3.08

0.002

-0.48 2.94

0.011

-0.50

2.99

0.002

-0.37 2.48

0.069

0.0005

-0.40

1.99 -0.045

-0.34 1.80

0.034

-0.33 1.42

-0.010

-0.33

1.42

0.007

-0.30 1.34

0.034

-0.003

-0.28

1.07

-0.001

-0.28

1.07 -0.003

-0.27 1.03

0.016

b

-0.41 2.02

-0.005

6

-0.33

1.43

0.007

7

-0.28

1.07

-0.39

1.98

8

-0.24 0.83

-0.005

-0.25 0.83

-0.003

-0.25

0.83 -0.004

-0.24 0.81

0.009

9

-0.22 0.66

-0.007

-0.22 0.66

-0.006

-0.22

0.66 -0.007

-0.21 0.65

0.003

10

-0.19 0.54

0.008

-0.20 0.54

0.009

-0.20

0.54

-0.19 0.54

0.016

Pi.i

P2.2

Pu

PKI

^2.2

0.008

P2,2

Pu

P2,2

1

9.9

989.9

9.9

989.9

0.988

98.8

1.96

98.8

2

1.4 E-4

5.2 E-3

4.2

166

1.2 E-4

4.9 E-3

1.96

65.5

3

6.0 E-5

1.3 E-3

2.23

47.5

5.7 E-5

1.2 E-3

1.64

32.8

4

3.7 E-5

5.0 E-4

1.47

19.6

3.7 E-5

4.9 E-4

1.27

16.5

5

2.7 E-5

2.5 E-4

1.09

9.89

2.7 E-5

2.5 E-4

1.01

9.01

6

2.1 E-5

1.4 E-4

0.86

5.67

2.1 E-5

1.4 E-4

0.82

5.37

7

1.7 E-5

8.8 E-5

0.71

3.56

1.8 E-5

8.8 E-5

0.69

3.43

8

1.5 E-5

5.8 E-5

0.60

2.37

1.5 E-5

5.9 E-5

0.59

2.31

9

1.3 E-5

4.1 E-5

0.52

1.66

1.3 E-5

4.1 E-5

0.52

1.63

10

1.1 E-5

3.0 E-5

0.47

1.21

1.1 E-5

3.0 E-5

0.46

1.19

assess this uncertainty as a function of time is by monitoring (e.g., in a Shewhart control chart) the value of the calibration factors obtained by a traditional calibration procedure, and modelling them by, e.g., an autoregressive model (see Section 20.4). However, because analytical chemists are unfamiliar with this type of filter, and the recursive approach is not simple, Kalman filters are hardly applied in calibration.

585

41.3 Recursive multicomponent analysis At this point we introduce the formal notation, which is commonly used in literature, and which is further used throughout this chapter. In the new notation we replace the parameter vector b in the calibration example by a vector x, which is called the state vector. In the multicomponent kinetic system the state vector x contains the concentrations of the compounds in the reaction mixture at a given time. Thus x is the vector which is estimated by the filter. The response of the measurement device, e.g., the absorbance at a given wavelength, is denoted by z. The absorbtivities at a given wavelength which relate the measured absorbance to the concentrations of the compounds in the mixture, or the design matrix in the calibration experiment (x in eq. (41.3)) are denoted by h^. Because zij) = zij) + v(/) and zij) = h^(j)x(j - 1 ) , the measurement equation given by eq. (41.1) is transformed into: z{j) = h^j)x(j^ \xp

1x1

l)+v(7) for; = 0,1, 2,... pxl

(41.5)

1x1

where v(j) includes the random measurement error (not to be confounded with r(j) which is the variance of the random measurement error) and the systematic error caused by the bias in model parameters, x(/ - 1). The equations for the variancecovariance update and Kalman gain update are adapted accordingly. The overall sequence of the equations is given in Table 41.3. It should be noted that here the state vector x is time-invariant. TABLE 41.3 Kalman filter algorithm equations for time-invariant system states Initialisation: x(0), P(0) Kalman gain update k ij) = P ( ; - l)h (;) (h'^ (;•) P(y - l)h (;) + K;))"^ px\

pxp pxl

Ixp

(41.6)

pxp pxl

Variance-covariance update PO; = PO"-l)-kWh'^0-)PO'-l) pxp

pxp pxl

(41.7)

Ixp pxp

Measurement zij) State update x(j)=xU-l)-^k(j)(z{j)-h^ (j)xU-l)) pxl pxl pxl 1x1 pxl 1x1

(41.8)

586 TABLE 41.4 Absorptivities of CI2 and Br2 in chloroform [2] Wavenumber cm-'/103

Absorptivities (>a 03) CI2

Absorbance of sample Br2

4.5

168

0.0341

24

8.4

211

0.0429

26

20.0

158

0.0335

28

56

30

0.0117

30

100

32

71

22

4.7 5.3

0.0110 0.0080

By way of illustration we work out the two-component analysis of chlorinebromine by UV-Vis spectrometry. In Table 41.4 the absorptivities of CI2 and Br2 in chloroform at six wavenumbers are given [2], with the absorbances measured for a supposedly unknown sample with concentrations (CQ^ = 0.1 and c^^. = 0.2). Concentrations of CI2 and Br2 are estimated by recursive regression. With this example we want to demonstrate that the evolution of the estimations of the unknown concentration depends on the design of the experiments, here the set of wavelengths and the order in which these wavelengths are selected. A measure for the precision of the estimation process is the variance-covariance matrix P. As said before, this precision can be evaluated before any measurement has been carried out. The two following wavelength sequences are evaluated — sequence A: 22-32 10^ cm"^ in ascending order, and sequence B: 22, 32, 24, 30, 26, 28 10^ cm~^ The measurement equation for the multicomponent system is given by:

za;=h^(/>o'-i) where z(j) is the predicted absorbance at the 7th measurement, using the latest estimated concentrations x(j - 1) obtained after the (/* - l)th measurement. h^O) contains the absorptivities of CI2 and Br2 at the wavelength chosen for the jth measurement. Step 1. Initialisation (/ = 0) P(0) =

"1000

1

1

1000

r(j) = 1 for ally

x(0) =

587

Step 2. Update of the Kalman gain vector k(l) and variance-covariance matrix P(l) k(l) = P(0)h(l)(hT(l)P(0)h(l) + i r ' P(l) = P(0)-k(l)h'^(l)P(0) This gives for design A:

k(l) =

fiooo L 1

0.0045 1000 1 1 [0.0045" +1 |[([0.0045 0.168] 1 1000 0.168 1000 [ 0.168

r

[0.1596" ' L 5-744 _

P(l) =

"1000 1

1 ~

1000 1 "0.1596" [0.0045-0.168] 1 1000 _ 5.74^t_

K)00^

999.2 -25.8 -25.8 34.88

Step 3. Predict the value of the first measurement z(l) z(l) = [0.0045 0.168]

To' 0

=0

Step 4. Obtain the first measurement: z(l) = 0.0341. Step 5. Update the predicted concentrations: x(l) = x(0) + k(l) (z(l) - 2(1)) x(l) =

"0" 0

+

"0.16'

[0.0341-0] =

_5.7_

0.0054 0.196

Step 6. Return to step 2. These steps are summarized in Tables 41.5 and 41.6. The concentration estimates should be compared with the true values 0.1 and 0.2 respectively. For design B the results listed in Table 41.7 are obtained. From both designs a number of interesting conclusions follow. (1) The set of selected wavelengths (i.e. the experimental design) affects the variance-covariance matrix, and thus the precision of the results. For example, the set 22, 24 and 26 (Table 41.5) gives a less precise result than the set 22, 32 and 24 (Table 41.7). The best set of wavelengths can be derived in the same way as for multiple linear regression, i.e. the determinant of the dispersion matrix (h^h) which contains the absorptivities, should be maximized.

588 TABLE 41.5 Calculated gain vector k and variance covariance matrix P\ f 11(0) = P2.2(^) = 1000; P^j^O) = ^2.1(0) = 1 > ''= 1 Step

Wavelength

k(«)

P(n)

A:l 1

22

2 3

kl

^2,2

^1.2-^2,1

34.88

-25.88

0.16

5.74

999.2

24

1.16

2.82

995.8

14.73

-34.12

26

9.37

1.061

859.7

12.98

-49.53

4

28

13.17

245.0

11.4

-18.11

5

30

7.11

-0.52

71.4

10.5

-5.6

6

32

3.71

-0.25

52.6

10.4

-4.3

-0.67

TABLE 41.6 Estimated concentrations (see Table 41.5 for the starting conditions) Step

Wavelength

x2 CL

Br.

1

22

0.0054

0.196

2

24

0.0072

0.2004

3

26

0.023

0.2022

4

28

0.0802

0.1993

5

30

0.0947

0.1982

6

32

0.0959

0.198

TABLE 41.7 Calculated gain vector, variance covariance matrix and estimated concentrations (see Table 41.5 for the starting conditions) Step

Wavelength

k{n) k\

x(«)

P(A2)

kl

Pu

^2.2

^1,2-^2,1

JCl

x2

999.2

34.88

-25.88

0.0054

0.196

166.2

34.43

-6.42

0.0825

0.194

166.2

13.81

-6.54

0.0825

0.198

62.6

13.68

-2.86

0.0938

0.197

1

22

0.16

2

32

11.76

3

24

0.016

4

30

6.24

5

26

0.59

1.560

62.1

10.4

-4.1

0.0941

0.198

6

28

2.82

0.069

52.2

10.4

-4.3

0.0955

0.198

5.744 -0.27 2.859 -0.22

589 TABLE 41.8 Concentration estimation with the optimal set of wavelengths (see Table 41.5 for the starting conditions) Step

Wavelength

k(n) k\

x(n)

P(^^) k2

Pli

P22

^1,2-^2,1

x\

x2

999.2

34.88

-25.88

0.0054

0.196

1

22

2

32

11.76

-0.27

166.2

34.43

-6.42

0.0825

0.1942

3

30

6.244

-0.18

62.6

34.34

-3.42

0.0940

0.1939

0.16

5.744

(2) From the evolution of P in design B (Table 41.7), one can conclude that the measurement at wavenumber 24 10^ cm~^ does not really improve the estimates already available after the sequence 22, 32. Equally the measurement at 26 10^ cm~^ does not improve the estimate already available after the sequence 22,32,24 and 30 10^ cm~^. This means that these wavelengths do not contain new information. Therefore, a possibly optimal set of wavenumbers is 22,32 and 30. Inclusion of a fourth wavelength namely at 28 10^ cm~^ probably does not improve the estimates of the concentrations already available, since the value of P converged to a stable value. To confirm this conclusion, the recursive regression was repeated for the set of wavelengths 22, 32 and 30 10^ cm"^ (see Table 41.8). Thyssen et al. [3] developed an algorithm for the selection of the best set of m wavelengths out of n. Instead of having to calculate 10^^ determinants to find the best set of six wavelengths out of 300, the recursive approach only needs to evaluate a rather straightforward equation 420 times. The influence of the measurement sequence on the speed of convergence is well demonstrated for the four-component analysis (Fig. 41.5) of a mixture of aniline, azobenzene, nitrobenzene and azoxybenzene [3]. In the forward scan mode a quick convergence is attained, whereas in the backward scan mode, convergence is slower. Using an optimized sequence, convergence is complete after less than seven measurements (Fig. 41.6). Other methods for wavelength selection are discussed in Chapters 10 and 27.

41.4 System equations When discussing the calibration and multicomponent analysis examples in previous sections, we mentioned that the parameters to be estimated are not necessarily constant but may vary in time. This variation is taken into account by

590

a

1 "

K

M

1 .6

•

"

K

X K M

M

K

X

«""

1.2 •".

X

K *

X

X X

r>

0.8

Z

*-

rzi

xa^ie"" -4

_

AA 1

1

14

1

1

28

1

1

42

1

1

56

1-

>

<

70

—> k

78

—> k

0.8 +

0.4 \

Fig. 41.5. Multicomponent analysis (aniline (jci), azobenzene fe), nitrobenzene fe) and azoxybenzene (JC4)) by recursive estimation (a) forward run of the monochromator (b) backward run (k indicates the sequence number of the estimates; solid lines are the concentration estimates; dotted lines are the measurements z).

591

l.B f

1.2

e.4 I ps..

X2K10 XltUB

M

28

42

56

70

—> k

Fig. 41.6. Multicomponent analysis (see Fig. 41.5) with an optimized wavelength sequence.

the system equation. As explained before, the system equation describes how the system changes in time. In the kinetics example the system equation describes the change of the concentrations as a function of time. This is a deterministic system equation. The random fluctuation of the slope and intercept of a straight line can be described by a stochastic model, e.g., an autoregressive model (see Chapter 20). Any unmodelled system fluctuations left are included in the system noise, w(j). The system equation is usually expressed in the following way: x(/*) = F ( / j - l ) x ( / - l ) + w(/')

(41.9)

where F(/j - 1) is the system transition matrix, which describes how the system state changes from time ti_^ to time /^. The vector w(/) consists of the noise contributions to each of the system states. These are system fluctuations which are not modelled by the system transition matrix. The parameters of the measurement equation, the h-vector and system transition matrix for the kinetics and calibration model are defined in Table 41.9. In the next two sections we derive the system equations for a kinetics and a calibration experiment. System state equations are not easy to derive. Their form depends on the particular system under consideration and no general guidance can be given for their derivation.

592 TABLE 41.9 Definition of the state parameter (x), h-vector and transition matrix (F) for two systems State parameters (x)

h-vector

1. Calibration

Slope and intercept of [ 1 c] where c is the concentthe calibration line at ration of the calibration time t standard measured at time t

2. 1 st order kinetics monitored by Uv-Vis A -> B

Concentrations of A and B at time t

Absorbance coefficients of A and B at the wavelength of the reading at time t

Transition matrix (F) time constant of the the variations of slope and intercept 1st order reaction rate

41.4.1 System equation for a kinetics experiment Let us assume that a kinetics experiment is carried out and we want to follow the concentrations of component B which is formed from A by the reaction: A -> B. For a first-order reaction, the concentrations of A (= jCj) and B (= x^^) as a function of time are described by two differential equations: dxj Idt = -k^x^ djC2/dr = fcjjCj

which can be rewritten in the following recursive form: jc,(r + 1) = (1 - A:i)jci(0 + w,(0

(41.10)

x^{t + 1) = k^x^{t) + JC2(0 + ^2(0 w indicates the error due to the discretization of the differential equations. These two equations describe the concentrations of A and B as a function of time, or in other words, they describe the state of the system, which in vector notation becomes: x ( r + l ) = Fxft) + w(0 with x\t^\)^[x,{t^\)x^{t+\)] w^(r+l) = [vi;i(r+1)^2(^+1)] and the transition matrix equal to: "l-A:,

0"

F

Jr —

k^

1

(41.11)

593

41.4.2 System equation of a calibration line with drift In this section we derive a system equation which describes a drifting calibration line. Let us suppose that the intercept x,(/ + 1) at a time; + 1 is equal to x,(/) at a time j augmented by a value a(/'), which is the drift. By adding a non-zero system noise w„ to the drift, we express the fact that the drift itself is also time dependent. This leads to the following equations [5,6]: x,(7+l) = x,0') + aO') a(/ + l) = a(/) + w„(/ + l) which is transformed into the following matrix notation:

.a(7 + l)

0 1 1 •^i(y) + 0 1 a(7) WaCi+l)

(41.12)

A similar equation can be derived for the slope, where P is the drift parameter: X2(;+l)

.P0"+1).

0

1 ip2(;) .0 iJIpCy).

(41.13)

H ' p ( ; + l)

Equations (41.12) and (41.13) can now be combined in a single system model which describes the drift in both parameters: \x,{i+\)

x^ii+i) a(/+l)

.PO'+l).

10

1 0 " \x\{j)^

0

0 1 U2(;) + 0 0 1 0 «(;•) H'aCy + i) >vp(; + l) 0 0 0 1 ^PO).

0

10

or x(/ + 1) = Fx(/) + w(7' + 1) with

F=

10 10" 0 10 1 0 0 10 0 0 0 1

Xi

andx =

a

.P.

F describes how the system state changes from time tj to tj^^.

(41.14)

594

For the time invariant calibration model discussed in Section 41.2, eq. (41.14) reduces to: "1

OTA:, •^iO')1

0

lJ[jC2lU)}

where x^ = intercept and X2 = slope.

41.5 The Kalman filter 41.5.1 Theory In Sections 41.2 and 41.3 we applied a recursive procedure to estimate the model parameters of time-invariant systems. After each new measurement, the model parameters were updated. The updating procedure for time-variant systems consists of two steps. In the first step the system state x(/ - 1) at time ti_^ is extrapolated to the state x(j) at time tj by applying the system equation (eq. (41.15)) in Table 41.10). At time tj a new measurement is carried out and the result is used to TABLE41.10 Kalman filter algorithm equations Initialisation: x(0), P(0) State extrapolation: system equation x(/!/- 1) = F(/)xO- 11/'- 1) + wO)

(41.15)

Covariance extrapolation P(/V- 1) = F(/-)P(/- II7- 1)F''(/) + Q 0 - 1)

(41.16)

New measurement: z(j) Measurement equation z(J) = h'(j)x(j\j - 1) + v(/') forj = 0„ 1 „ 2,... Covariance update P(/'iy) = P(/"iy - 1) - mh'ijWiiy

-1)

(41.17)

Kalman gain update k(/) = POV- l)h(/)(h''(/-)POV- l)hO) + rO))-' State update x(J\J) = ^(J\J-^) + mKz(j)-h'(j)x(j\j-

D)

(41.18)

(41.19)

595

update the state x(j) (eq. (41.19) in Table 41.10). In order to make a distinction between state extrapolations by the system equation and state updates when making a new observation, a double index (j\j - 1) is used. The index (j\j - 1) indicates the best estimates at time tj, based on measurements obtained up to and including those obtained at point tj^^. Equations (41.15) and (41.19) for the extrapolation and update of system states form the so-called state-space model. The solution of the state-space model has been derived by Kalman and is known as the Kalman filter. Assumptions are that the measurement noise v(j) and the system noise w(;) are random and independent, normally distributed, white and uncorrected. This leads to the general formulation of a Kalman filter given in Table 41.10. Equations (41.15) and (41.19) account for the time dependence of the system. Eq. (41.15) is the system equation which tells us how the system behaves in time (here inj units). Equation (41.16) expresses how the uncertainty in the system state grows as a function of time (here inj units) if no observations would be made. Q(/ - 1) is the variance-covariance matrix of the system noise which contains the variance of w. The algorithm is initialized in the same way as for a time-invariant system. The sequence of the estimations is as follows: Cycle 1 1. Obtain initial values for x(OIO), k(0) and P(OIO) 2. Extrapolate the system state (eq. (41.15)) to x(llO) 3. Calculate the associated uncertainty (eq. (41.16)) P( 110) 4. Perform measurement no. 1 5. Update the gain vector k(l) (eq. (41.18)) using P(IIO) 6. Update the estimate x(llO) to x(lll) (eq. (41.19)) and the associated uncertainty P(lll) (eq. (41.17)) Cycle 2 1. Extrapolate the estimate of the system state to x(2ll) and the associated uncertainty P(2I1) 2. Peform measurement no. 2 3. Update the gain vector k(2) using P(2I1) 4. Update (filter) the estimate of the system state to x(2l2) and the associated uncertainty P(2I2) and so on. In the next section, this cycle is demonstrated on the kinetics example introduced in Sections 41.1 and 41.4. Time-invariant systems can also be solved by the equations given in Table 41.10. In that case, F in eq. (41.15) is substituted by the identity matrix. The system state, x(/), of time-invariant systems converges to a constant value after a few cycles of the filter, as was observed in the calibration example. The system state,

596

x(/), of time-variant systems is obtained as a function of y, for example the concentrations of the reactants and reaction products in a kinetic experiment monitored by a spectrometric multicomponent analysis. 41.5.2 Kalman filter of a kinetics model Equation (41.11) represents the (deterministic) system equation which describes how the concentrations vary in time. In order to estimate the concentrations of the two compounds as a function of time during the reaction, the absorbance of the mixture is measured as a function of wavelength and time. Let us suppose that the pure spectra (absorptivities) of the compounds A and B are known and that at a time t the spectrometer is set at a wavelength giving the absorptivities h^(0- The system and measurement equations can now be solved by the Kalman filter given in Table 41.10. By way of illustration we work out a simplified example of a reaction with a true reaction rate constant equal to k^ =0.1 min"^ and an initial concentration jCi(O) = 1. The concentrations are spectrophotometrically measured every 5 minutes and at the start of the reaction after 1 minute. Each time a new measurement is performed, the last estimate of the concentration A is updated. By substituting that concentration in the system equation x^(t) = x^(0)txp(-k^t) we obtain an update of the reaction rate k. With this new value the concentration of A is extrapolated to the point in time that a new measurement is made. The results for three cycles of the Kalman filter are given in Table 41.11 and in Fig. 41.7. The "c ik (0 I i •G

1?

(Q

g>

1

O C

o (Q

08

C

8

0.6

oC o 0.4 0.2 o

0

_l

I

I

l_

_J

I

\

o

o

o

L_

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 • time

Fig. 41.7. Concentrations of the reactant A (reaction A ^ B) as a function of time (dotted line) (CA = 1, CB = 0); • state updates (after a new measurement), O state extrapolations to the next measurement (see Table 41.11 for Kalman filter settings).

597 TABLE41.il The prediction of concentrations and reaction rate constant by a Kalman filter; (JCJCO) =l,k^=0.l Time

Concentrations

Wavelengthy

Discrete

Continuous

Absorptivities

Absorbance

A

z

B

z

min-^

Estimate of the concentration of A State update State extra(41.19)

1 ^-

polation (41.15)

1)

A

B

A

B

0

1

0

1

0

1

0.90

0.10

0.90

0.10

2

0.82

0.18

0.81

0.19

0.82

3

0.74

0.28

0.73

0.27

0.74

4

0.67

0.33

0.66

0.34

5

0.61

0.39

0.59

0.41

6

0.55

0.45

0.53

0.47

0.57

7

0.50

0.50

0.48

0.52

0.52

8

0.45

0.55

0.43

0.57

0.48

9

0.41

0.59

0.39

0.61

0.43

x(OIO) 1.00 ->1

10

1

9.1

10

jc(lll)0.91

;c(llO)

0.90

0.67 ->2

5

2

3.8

3.97

3.87

3.81

jc(2l2) 0.63

4211)

jc(3l2)

0.62

0.40

10

0.37

0.63

0.35

0.65

11

0.33

0.67

0.31

0.69

0.34

12

0.29

0.71

0.28

0.72

0.31

13

0.27

0.73

0.25

0.75

0.28

14

0.25

0.75

0.23

0.77

0.26

15

0.22

0.78

0.21

0.79

^ 3

->4

2

7

5

2

3.10

3.15

A:(3I3) 0.38

;c(4l4)0.16

x(4l3)

0.23

16

0.20

0.80

0.18

0.82

0.15

17

0.18

0.82

0.17

0.83

0.13

18

0.16

0.84

0.15

0.85

0.12

19

0.15

0.85

0.14

0.86

0.11 0.10

')State extrapolations in italic are obtained by substituting the last estimate of k^ equal to -\n{x(j\j))/t in the system equation (41.15).

dashed line in Fig. 41.7 shows the true evolution of the concentration of the reactant A. For a matter of simplicity we have not included the covariance extrapolation (eq. (41.16)) in the calculations. Although more cycles are needed for a convincing demonstration that a Kalman filter can follow time varying states, the example clearly shows its principle.

598

41.5.3 Kalman filtering of a calibration line with drift The measurement model of the time-invariant calibration system (eq. (41.5)) should now be expanded in the following way: z(j) = h'iDxij - 1) + v(j) for; = 0, 1, 2,...

(41.20)

where h^(/) = [lc(/)0 0] c(j) is the concentration of the analyte in theyth calibration sample

x'{j-l)

= [b, b, a p ]

The model contains four parameters, the slope and intercept of the calibration line and two drift parameters a and p. All four parameters are estimated by applying the algorithm given in Table 41.10. Details of this procedure are given in Ref. [5].

41.6 Adaptive Kalman filtering In previous sections we demonstrated that a major advantage of Kalman filtering over ordinary least squares (OLS) procedures is that it can handle timevarying systems, e.g., a drifting calibration parameter and drifting background. In this section another feature of Kalman filters is demonstrated, namely its capability of handling incomplete measurement models or model errors. An example of an incomplete measurement model is multicomponent analysis in the presence of an unknown interfering compound. If the identity of this interference is unknown, one cannot select wavenumbers where only the analytes of interest absorb. Therefore, solution by OLS may lead to large errors in the estimated concentrations. The occurrence of such errors may be detected by inspecting the difference between the measured absorbance for the sample and the absorbance estimated from the predicted concentrations (see Chapter 10). However, inspection of PRESS does not provide information on which wavelengths are not selective. One finds that the result is wrong without an indication on how to correct the error. Another type of model error is a curving calibration line which is modelled by a straight line model. The size and pattern of the residuals usually indicate that there is a model error (see Chapter 8), after which the calculations may be repeated with another model, e.g., a quadratic curve. The recursive property of the Kalman filter allows the detection of such model deviations, and offers the possibility of disregarding the measurements in the region where the model is invalid. This filter is the so-called adaptive Kalman filter.

599

41.6.1 Evaluation of the innovation Before we can apply an adaptive filter, we should define a criterion to judge the validity of the model to describe the measurements. Such a criterion can be based on the innovation defined in Section 41.2. The concept of innovation, /, has been introduced as a measure of how well the filter predicts new observations:

/(/•)=zij) - h^ox/ -1)=zij) - i(j) where zij) is the jth measurement, x(/ - 1) is the estimate of the model parameters after 7 - 1 observations, and h^(/) is the design vector. Thus /(/) is a measure of the predictive ability of the model. For the calibration example discussed in Section 41.2, x(/ - 1) contains the slope and intercept of the straight line, and h^(/) is equal to [1 c(j)] with c(j) the concentration of the calibration standard for the jth calibration measurement. For the multicomponent analysis (MCA), x(/ - 1 ) contains the estimated concentrations of the analytes after j - I observations, and h^(/) contains the absorptivities of the analytes at wavelength;. It can be shown [4] that the innovations of a correct filter model applied on data with Gaussian noise follows a Gaussian distribution with a mean value equal to zero and a standard deviation equal to the experimental error. A model error means that the design vector h in the measurement equation is not adequate. If, for instance, in the calibration example the model was quadratic, h^ should be [1 c(j) c(/)^] instead of [1 c(j)]. In the MCA example h^(/) is wrong if the absorptivities of some absorbing species are not included. Any error in the design vector h^ appears by a non-zero mean for the innovation [4]. One also expects the sequence of the innovation to be random and uncorrelated. This can be checked by an investigation of the autocorrelation function (see Section 20.3) of the innovation. 41.6.2 The adaptive Kalman filter model The principle of adaptive filtering is based on evaluating the observed innovation calculated at each new observation and comparing this value to the theoretically expected value. If the (absolute) value of the innovation is larger than a given criterion, the observation is disregarded and not used to update the estimates of the parameters. Otherwise, one could eliminate the influence of that observation by artificially increasing its measurement variance r(/), which effectively attributes a low weight to the observation. For a time-invariant system, the expected standard deviation of the innovation consists of two parts: the measurement variance (r(/)), and the variance due to the uncertainty in the parameters (P(/)), given by [4]:

600

OJU) = [r(j) + h^(/*)P(/' - l)mr'

(41.21)

As explained before, the second term in the above equation is the variance of the response, £(/), predicted by the model at the value h(/) of the independent variable, given the uncertainty in the regression parameters P(/ - 1) obtained so far. This equation reflects the fact that the fluctuations of the innovation are larger at the beginning of the filtering procedure, when the states are not well known, and converge to standard deviation of the experimental error when the confidence in the estimated parameters becomes high (small P). Therefore, it is more difficult to detect model errors at the beginning of the estimation sequence than later on. Rejection and acceptance of a new measurement is then based on the following criterion: if |/0)|>3a,(/): reject otherwise: accept Adaptation of the Kalman filter may then simply consist of ignoring the rejected observations without affecting the parameter estimates and covariances. When using eq. (41.21) a complication arises at the beginning of the recursive estimation because the value of P depends on the initially chosen values P(0) and thus is not a good measure for the uncertainty in the parameters. When large values are chosen for the diagonal terms of P in order to speed up the convergence (high P means large gain factor k and thus large influence of the last observation), eq. (41.21) overestimates the variance of the innovation, until P becomes independent of P(0). For short runs one can evaluate the sequence of the innovation and look for regions with significantly larger values, or compare the innovation with r(j). By way of illustration we apply the latter procedure to solve the multicomponent system discussed in Section 41.3 after adding an interfering compound which augments the absorbance at 26 10^ cm"^ with 0.01 au and at 28 10^ cm"^ with 0.015 au. First we apply the non-adaptive Kalman filter to all measurements. The estimation then proceeds as shown in Table 41.12. The above example illustrates the self adaptive capacity of the Kalman filter. The large interferences introduced at the wavelengths 26 and 28 10^ cm~^ have not really influenced the end result. At wavelengths 26 and 28 10^ cm"^ the innovation is large due to the interferent. At 30 10^ cm"^ the innovation is high because the concentration estimates obtained in the foregoing step are poor. However, the observation at 30 10^ cm~^ is unaffected by which the concentration estimates are restored within the true value. In contrast, the OLS estimates obtained for the above example are inaccurate (jCj = 0.148 and X2 = 0.217) demonstrating the sensitivity of OLS for model errors.

601 TABLE 41.12 Non-adaptive Kalman filter with interference at 26 and 28 10-3 cm-^ (see Table 41.5 for the starting conditions) Step 7

Wavelength

0

Innovation

^1

^2

Absorbance

CI2

Br2

Measured^)

Estimated^)

0

0

1

22

0.054

0.196

0.0341

0

0.0341

2

24

0.072

0.2004

0.0429

0.0414

0.0015

3

26

0.169

0.211

0.0435

0.0331

0.0104

4

28

0.313

0.204

0.0267

0.0158

0.0110

5

30

0.099

0.215

0.0110

0.0409

-0.030

6

32

0.098

0.215

0.0080

0.0082

-0.0002

'^Values taken from Table 41.4. 2)Calculated with the absorbtivity coefficients from Table 41.4.

The estimation sequence when the Kalman filter is adapted is given in Table 41.13. This illustrates well the adaptive procedure, which has been followed. At 26 10-^ cm~^ the new measurement is 0.0104 absorbance units higher than expected from the last available concentration estimates (x^ = 0.072 and X2 = 0.2004). This deviation is clearly larger than the value 0.005 expected from the measurement noise. Therefore, the observation is disregarded and the next measurement at 28 10^ cm~^ is predicted with the last accepted concentration estimates. Again, the difference between predicted and measured absorbance (= innovation) cannot be explained from the noise and the observation is disregarded as well. At 30 10^ cm"^ the predicted absorbance using the concentration estimates from the second step is back within the expectations, P and k can be updated leading to new concentration estimates x^ = 0.098 and X2 = 0.1995. Thereafter, the estimation process is continued in the normal way. The effect of this whole procedure is that the two measurements corrupted by the presence of an interferent have been eliminated after which the measurement-filtering process is continued.

41.7 Applications One of the earliest applications of the Kalman filter in analytical chemistry was multicomponent analysis by UV-Vis spectrometry of time and wavelength independent concentrations, which was discussed by several authors [7-10]. Initially, the spectral range was scanned in the upward and downward mode, but later on

602 TABLE 41.13 Adaptive Kalman filter with interference at 26 and 28 cm-' (see Table 41.5 for the starting conditions) Step 7

Wavelength

0 1

22

Innovation

^1

X2

Absorbance

CI2

Br2

Measured')

Estimated^)

0

0

0.054

0.196

0.0341

0

0.0341 0.0015 0.0104

2

24

0.072

0.2004

0.0429

0.0414

3

26

0.0331

28

-

0.0435

4

-

0.0267

0.0100

0.0166 -0.002 -0.00001

5

30

0.0980

0.1995

0.0110

0.00814

6

32

0.0979

0.1995

0.0080

0.00801

')Values taken from Table 41.4. 2)Calculated with absorbtivity coefficients from Table 41.4.

optimal sequences were derived for faster convergence to the result [3]. The measurement model can be adapted to include contributions from linear drift or background [6,11]. This requires an accurate model for the background or drift. If the background response is not precisely described by the model the Kalman filter fails to estimate accurate concentrations. Rutan [12] applied an adaptive Kalman filter in these instances with success. In HPLC with diode array detection, threedimensional data are available. The processing of such data by multivariate statistics has been the subject of many chemometric studies, which are discussed in Chapter 34. Under the restriction that the spectra of the analytes should be available, accurate concentrations can be obtained by Kalman filtering in the presence of unknown interferences [13]. One of the earliest reports of a Kalman filter which includes a system equation is due to Seelig and Blaunt [14] in the area of anodic stripping voltametry. Five system states — potential, concentration, potential sweep rate and the slope and intercept of the non-Faradaic current — were predicted from a measurement model based on the measurement of potential current. Later on the same approach was applied in polarography. Similar to spectroscopic and chromatographic applications, overlapping voltamograms can be resolved by a Kalman filter [15]. A vast amount of applications of Kalman filters in kinetic analysis has been reported [16,17] and the performance has been compared with conventional non-linear regression. In most cases the accuracy and precision of the results obtained from the two methods were comparable. The Kalman filter is specifically superior for detecting and correcting model errors.

603

The Kalman filter is particularly well-suited to monitor the dynamic behaviour of processes. The measurement procedure in itself can be considered to be a system which is observed through the measurement of check samples. One can set up system equations, e.g., a system equation which describes the fluctuations of the calibration factors. Only a few applications exploiting this capability of a Kalman filter have been reported. One of the difficulties is a lack of system models, which describe the dynamic behaviour of an analytical system. Thyssen and coworkers [17] demonstrated the potential of this approach by designing a Kalman filter for the monitoring of the calibration factors. They assembled a so-called self-calibrating Flow Injection Analyzer for the determination of chloride in water. The software of the instrument included a system model by which the uncertainty of the calibration factors was evaluated during the measurement of the unknown samples. When this uncertainty exceeded a certain threshold the instrument decided to update the calibration factors by remeasuring one of the calibration standards. Thyssen [18] also designed an automatic titrator which controlled the addition of the titrant by a Kalman filter. After each addition the equivalence point (the state of the system) was estimated during the titration.

References 1. 2. 3.

4. 5. 6.

7. 8. 9. 10.

11.

D. Graupe, Identification of Systems. Krieger, New York, NY, 1976. Landolt-Bornstein, Zalen Werte und Funktionen. Teil 3, Atom und Molekular Physik. Springer, Berlin 1951. P.C. Thijssen, L.J.P. Vogels, H.C. Smit and G. Kateman, Optimal selection of wavelengths in spectrophotometric multicomponent analysis using recursive least squares. Z. Anal. Chem., 320(1985)531-540. A. Gelb (Ed.), Applied Optimal Estimation. MIT Press, Cambridge, MA, 1974. G. Kateman and L. Buydens, Quality Control in Analytical Chemistry, 2nd Edn. Wiley, New York, 1993. P.C. Thijssen, S.M. Wolfrum, G. Kateman and H.C. Smit, A Kalman filter for calibration,evaluation of unknown samples and quality control in drifting systems: Part 1. Theory and simulations. Anal. Chim. Acta, 156 (1984) 87-101. H.N.J. Poulisse, Multicomponent analysis computations based on Kalman Filtering. Anal. Chim. Acta, 112 (1979) 361-374. C.B.M. Didden and H.N.J. Poulisse. On the determination of the number of components from simulated spectra using Kalman filtering. Anal. Lett., 13 (1980) 921-935. T.F. Brown and S.D. Brown, Resolution of overlapped electrochemical peaks with the use of the Kalman filter. Anal. Chem., 53 (1981) 1410-1417. S.C. Rutan and S.D. Brown, Pulsed photoacoustic spectroscopy and spectral deconvolution with the Kalman filter for determination of metal complexation parameters. Anal. Chem., 55 (1983) 1707-1710. P.C. Thijssen, A Kalman filter for calibration, evaluation of unknown samples and quality control in drifting systems: Part 2. Optimal designs. Anal. Chim. Acta, 162 (1984) 253-262.

604 12.

13.

14. 15.

16. 17.

18.

19.

S.C. Rutan, E. Bouveresse, K.N. Andrew, P.J. Worsfold and D.L. Massart, Correction for drift in multivariate systems using the Kalman filter. Chemom. Intell. Lab. Syst., 35 (1996) 199-211. J. Chen and S.C. Rutan, Identification and quantification of overlapped peaks in liquid chromatography with UV diode array detection using an adaptive Kalman filter. Anal. Chim. Acta, 335(1996) 1-10. P. Seelig and H. Blount, Kalman Filter applied to Anodic Stripping Voltametry: Theory. Anal. Chem., 48 (1976) 252-258. C.A. Scolari and S.D. Brown, Multicomponent determination in flow-injection systems with square-wave voltammetric detection using the Kalman filter. Anal. Chim. Acta, 178 (1985) 239-246. B.M. Quencer, Multicomponent kinetic determinations with the extended Kalman filter. Diss. Abstr. Int. B 54 (1994) 5121-5122. M. Gui and S.C. Rutan, Determination of initial concentration of analyte by kinetic detection of the intermediate product in consecutive first-order reactions using an extended Kalman filter. Anal. Chim. Acta, 66 (1994) 1513-1519. P.C. Thijssen, L.T.M. Prop, G. Kateman and H.C. Smit, A Kalman filter for caHbration, evaluation of unknown samples and quality control in drifting systems. Part 4. Flow Injection Analysis. Anal. Chim. Acta, 174 (1985) 27-40. P.C. Thijssen, N.J.M.L. Janssen, G. Kateman and H.C. Smit, Kalman filter applied to setpoint control in continuous titrations. Anal. Chim. Acta, 177 (1985) 57-69.

Recommended additional reading S.C. Rutan, Recursive parameter estimation. J. Chemom., 4 (1990) 103-121. S.C. Rutan, Adaptive Kalman filtering. Anal. Chem., 63 (1991) 1103A-1109A. S.C. Rutan, Fast on-line digital filtering. Chemom. Intell. Lab. Syst., 6 (1989) 191-201. D. Wienke, T. Vijn and L. Buydens, Quality self-monitoring of intelligent analyzers and sensor based on an extended Kalman filter: an application to graphite furnace atomic absorption spectroscopy. Anal. Chem., 66 (1994) 841-849. S.D. Brown, Rapid parameter estimation with incomplete chemical calibration models. Chemom. Intell. Lab. Syst., 10 (1991) 87-105.