- Email: [email protected]
The Kalman filter in quantitative multicomponent analysis
trends in analytical chemistry
203
vol. 17, no. 6, 1992
The Kalman filter in quantitative multicomponent analysis* G. Gauglitz, M. Mettler and S. Weiss Ttibingen, Germany
Quantitative multicomponent analysis in UV-VIS spectroscopy can allow to replace high-performance liquid chromatographic (HPLC) methods or support them by increasing detection certainty. However, commercial software offers the possibility to carry out problematic analyses, because the potential of the individual multicomponent procedures is seldom exploited to full advantage. Usually, only one method is offered, without any modification possibilities. Difficult analysis may be caused by high concentrations and interferences between the compounds of the mixture. Thus the electron systems of the molecules respectively the UV spectra are altered in comparison with the spectra of the pure components in high dilution. These non-linearities may be compensated by multicalibration techniques with respect to the interpretation of multicomponent spectra, a possibility also existing for ordinary least squares (the K matrix approach). Spectral similarities between the individual calibration spectra shift non-linearities into the range of characteristic differences between the components. It is therefore useful to make a selection of the examined wavelengths to exclude similar ranges from the calculation of the concentrations. When the more common least squares methods are used, impurities or unknown components lead to wrong parameters (concentrations), and the procedure does not provide for correcting the calibration model. On the assumption that the above-mentioned faults can becombined in a single component, the adaptive Kalman filtering makes a proposal for the unmodelled compound [l]. This proposition can be used to alter the model and run the procedure a second time. The non-adaptive form of the filter is
approximately equivalent to ordinary least squares; the ordinary Kalman filter is sometimes called “sequential least squares” [2]. In the presence of an unmodelled compound the adaptive form provides a better result, both with simulated and real UV spectra. However, the estimation of the missing compound is of limited usability for the refinement of the model when using real spectra. Usually only certain wavelength ranges are computed in good accordance to the spectrum of the missing compound. In this case, the residuals would provide a better estimation. Another problem is the often narrow optimal range of the starting values for the variance+zovariance matrix in the adaptive mode. The following implementation of the Kalman filter in its adaptive and non-adaptive form is part of a larger multicomponent analysis package written in C for the IBM-PC. The adaptive filter adjusts the variance vector at each iteration step as the filter progresses, the ordinary filter uses the measurement variance vector. The symbols used for the description of the algorithm are also used in the program listing. The number i indicates the actual iteration step, i-l the previous one. Matrices are represented by uppercase bold characters, vectors by lowercase ones.
Algorithm [3] Calculation of Kalman gain factors: g(i) = C(i-1) A(i) / [r(i) + A’(i) C(i-1) A(i)]
where r(i) = l/m [gv$_j)]-A’(i)
C(i-l)A(i)
j=l
in the adaptive form, respectively the measurement variance in the non-adaptive form. Error covariance update:
* Dedicated to Prof. E. Bayer on the occasion of his 65th birthday.
C(i) = [ 1 - g(i) A’(i)] C(i-1)
204
trends
X(i) = X(i-1) + g(i) [y(i) - A’(i) X(i-l)] :=
7992
where A = calibration matrix; A’ = transpose of calibration matrix; r = variance (measurement variance or adapted variance); y = measurement vector.
Calculation of actual parameters:
[y(i) - A’(i) X(i-l)]
inanalytical chemistry,vol. 77, no. 6,
innovation v
Definitions The matrices are defined as pointers to structures:
typedef
struct matrix { double *m-mat; m_sa; int m_za; int **m_pnt; double } Matrix;
/* /* /* /*
pointer to matrix number of colums number of rows pointer to rows
A matrix library was written for allocating memory and calculating with matrix structures. The ac-
memory
*/ */ */ */
cess to the data structure and its elements is done via macros:
/* access to matrix memory */ (a) ->m_mat #define Mmat(a) /* get number of rows (a)->m_sa #define Msa(a) */ /* get number of colums */ (a)->m_za #define Mza(a) #define Mat(a,b,c) *(Mmat(a) + Msa(a)*(b) + (c)j access to element i,j of matrix u */ /* e.g. Mat(u,i,j):
Dimensions of matrices
Other variables
m = number of mixtures, k = wavelengths
double cov_parm
rows,columns X m,k A k,m m,l g C m,m templ m,m temp2 l,m m_var k,l innovations k, 1
int direction int iter int window temporary matrix temporary matrix measurement noise vector
int i,_i,p,pp
starting values of variance covariance matrix of the digital filter( 1=forward, -l=backward) iteration counter smoothing window of the adaptive filter loop counters
Listing /***** initial guess of X(0) and g *****/ if(direction==l) j=O; else j=k-1; for(i=m;i--; ) /* initial state (cont.) { Mat(X,i,j) = 0.0; Mat(g,i,j) = 1.0; /* Kalman gain factors /* diag. elem. variace/cov. Mat(C,i,i) = cov_parm; 1 Mat(innovations,j,O) = Mat(y,j,O); /*****X********* iteration loop ******************/ for( p=O,pp=k-2 ; ++p
--pp
)
*/ */ */
205
trends in analytical chemistry, vol. 11, no. 6, 1992
/***** estimate for(i=m;i--; ) Mat(X,i,iter)
next
system
state
*****/
= Mat(X,i,iter-direction)
/***** error covariance matcopy(C,C old);
extrapolation
*****/
/***** update Kalman gain factors *****/ if(mode>0) { /* a d a p t e d v a r i a n c e */ adaptvariance(r,innovations,A,C_old, iter,window,m,k,direction); /* u s e r i n s t e a d of m _ v a r */ update_kgain(C old,A, t e m p 2 , g , m , M a t ( r , i t e r , 0),iter);
} else update_kgain(C_old,A, temp2,g,m,Mat(m_var, iter, 0 ) , i t e r ) ; /****** update state estimate ******/ for (i = m, v k = 0.0; i--; ) /* c o m p u t e i n n o v a t i o n */ v k += M a t ( A , iter, i) * M a t ( X , i , i t e r - d i r e c t i o n ) ; v k = M a t ( y , i t e r , 0) vk; Mat(innovations,iter, 0) = vk; /* r e t a i n i n n o v a t i o n */ for (i : m; i--; ) Mat(X,i,iter) += v k * M a t ( g , i , 0 ) ; /***** update error covariance *****/ for (i = m; i--; ) for (j = m; j--; ) Mat(templ,i,j) : (i == j) * 1 . 0 Mat(g,i,0) * M a t ( A , iter, j); matmul(templ,C_old,C); } /* i t e r a t i o n l o o p */ void
update
kgain(Matrix *C_old, M a t r i x *A, M a t r i x *temp, M a t r i x *g, int m, d o u b l e m _ v a r , int iter) { d o u b l e kdiv, s; int j, i ; for (i : m; i--; ) /* i' (i)C(i-l) { for (j --m, k d i v = 0.0; j--; ) k d i v += M a t ( A , iter, j) * M a t (C old, i,j) ; Mat(temp,0,i) = kdiv;
*/
} for (i = m, k d i v = m _ v a r ; i--; ) /* .... * i(i) k d i v += M a t ( t e m p , 0 , i ) * M a t ( A , iter, i); /* = k d i v for (i = m; i--; ) /* C ( i - l ) A ( i ) / kdiv { for (j = m, s = 0.0; j--; ) s += M a t ( C old, i,j) * M a t ( i , iter, j); Mat(g,i,0) = s / kdiv;
*/ */ */
} }
R eferences 1
S.C. Rutan and S.D. Brown, Anal. Chim. Acta, 160 (1985) 99-119. 2 S.D. Brown, Chemom. lnteL Lab. Syst., 10 (1991) 87-105. 3 S.D. Brown, AnaL Chim. Acta, 181 (1986) 1-26.
G. Gauglitz, M. Mettler and S. Weiss are at the Institut fEir Physikalische und Theoretische Chemie der Universit~t TE~bingen, Auf der Morgenstelle 8, T~bingen, Germany
203
vol. 17, no. 6, 1992
The Kalman filter in quantitative multicomponent analysis* G. Gauglitz, M. Mettler and S. Weiss Ttibingen, Germany
Quantitative multicomponent analysis in UV-VIS spectroscopy can allow to replace high-performance liquid chromatographic (HPLC) methods or support them by increasing detection certainty. However, commercial software offers the possibility to carry out problematic analyses, because the potential of the individual multicomponent procedures is seldom exploited to full advantage. Usually, only one method is offered, without any modification possibilities. Difficult analysis may be caused by high concentrations and interferences between the compounds of the mixture. Thus the electron systems of the molecules respectively the UV spectra are altered in comparison with the spectra of the pure components in high dilution. These non-linearities may be compensated by multicalibration techniques with respect to the interpretation of multicomponent spectra, a possibility also existing for ordinary least squares (the K matrix approach). Spectral similarities between the individual calibration spectra shift non-linearities into the range of characteristic differences between the components. It is therefore useful to make a selection of the examined wavelengths to exclude similar ranges from the calculation of the concentrations. When the more common least squares methods are used, impurities or unknown components lead to wrong parameters (concentrations), and the procedure does not provide for correcting the calibration model. On the assumption that the above-mentioned faults can becombined in a single component, the adaptive Kalman filtering makes a proposal for the unmodelled compound [l]. This proposition can be used to alter the model and run the procedure a second time. The non-adaptive form of the filter is
approximately equivalent to ordinary least squares; the ordinary Kalman filter is sometimes called “sequential least squares” [2]. In the presence of an unmodelled compound the adaptive form provides a better result, both with simulated and real UV spectra. However, the estimation of the missing compound is of limited usability for the refinement of the model when using real spectra. Usually only certain wavelength ranges are computed in good accordance to the spectrum of the missing compound. In this case, the residuals would provide a better estimation. Another problem is the often narrow optimal range of the starting values for the variance+zovariance matrix in the adaptive mode. The following implementation of the Kalman filter in its adaptive and non-adaptive form is part of a larger multicomponent analysis package written in C for the IBM-PC. The adaptive filter adjusts the variance vector at each iteration step as the filter progresses, the ordinary filter uses the measurement variance vector. The symbols used for the description of the algorithm are also used in the program listing. The number i indicates the actual iteration step, i-l the previous one. Matrices are represented by uppercase bold characters, vectors by lowercase ones.
Algorithm [3] Calculation of Kalman gain factors: g(i) = C(i-1) A(i) / [r(i) + A’(i) C(i-1) A(i)]
where r(i) = l/m [gv$_j)]-A’(i)
C(i-l)A(i)
j=l
in the adaptive form, respectively the measurement variance in the non-adaptive form. Error covariance update:
* Dedicated to Prof. E. Bayer on the occasion of his 65th birthday.
C(i) = [ 1 - g(i) A’(i)] C(i-1)
204
trends
X(i) = X(i-1) + g(i) [y(i) - A’(i) X(i-l)] :=
7992
where A = calibration matrix; A’ = transpose of calibration matrix; r = variance (measurement variance or adapted variance); y = measurement vector.
Calculation of actual parameters:
[y(i) - A’(i) X(i-l)]
inanalytical chemistry,vol. 77, no. 6,
innovation v
Definitions The matrices are defined as pointers to structures:
typedef
struct matrix { double *m-mat; m_sa; int m_za; int **m_pnt; double } Matrix;
/* /* /* /*
pointer to matrix number of colums number of rows pointer to rows
A matrix library was written for allocating memory and calculating with matrix structures. The ac-
memory
*/ */ */ */
cess to the data structure and its elements is done via macros:
/* access to matrix memory */ (a) ->m_mat #define Mmat(a) /* get number of rows (a)->m_sa #define Msa(a) */ /* get number of colums */ (a)->m_za #define Mza(a) #define Mat(a,b,c) *(Mmat(a) + Msa(a)*(b) + (c)j access to element i,j of matrix u */ /* e.g. Mat(u,i,j):
Dimensions of matrices
Other variables
m = number of mixtures, k = wavelengths
double cov_parm
rows,columns X m,k A k,m m,l g C m,m templ m,m temp2 l,m m_var k,l innovations k, 1
int direction int iter int window temporary matrix temporary matrix measurement noise vector
int i,_i,p,pp
starting values of variance covariance matrix of the digital filter( 1=forward, -l=backward) iteration counter smoothing window of the adaptive filter loop counters
Listing /***** initial guess of X(0) and g *****/ if(direction==l) j=O; else j=k-1; for(i=m;i--; ) /* initial state (cont.) { Mat(X,i,j) = 0.0; Mat(g,i,j) = 1.0; /* Kalman gain factors /* diag. elem. variace/cov. Mat(C,i,i) = cov_parm; 1 Mat(innovations,j,O) = Mat(y,j,O); /*****X********* iteration loop ******************/ for( p=O,pp=k-2 ; ++p
--pp
)
*/ */ */
205
trends in analytical chemistry, vol. 11, no. 6, 1992
/***** estimate for(i=m;i--; ) Mat(X,i,iter)
next
system
state
*****/
= Mat(X,i,iter-direction)
/***** error covariance matcopy(C,C old);
extrapolation
*****/
/***** update Kalman gain factors *****/ if(mode>0) { /* a d a p t e d v a r i a n c e */ adaptvariance(r,innovations,A,C_old, iter,window,m,k,direction); /* u s e r i n s t e a d of m _ v a r */ update_kgain(C old,A, t e m p 2 , g , m , M a t ( r , i t e r , 0),iter);
} else update_kgain(C_old,A, temp2,g,m,Mat(m_var, iter, 0 ) , i t e r ) ; /****** update state estimate ******/ for (i = m, v k = 0.0; i--; ) /* c o m p u t e i n n o v a t i o n */ v k += M a t ( A , iter, i) * M a t ( X , i , i t e r - d i r e c t i o n ) ; v k = M a t ( y , i t e r , 0) vk; Mat(innovations,iter, 0) = vk; /* r e t a i n i n n o v a t i o n */ for (i : m; i--; ) Mat(X,i,iter) += v k * M a t ( g , i , 0 ) ; /***** update error covariance *****/ for (i = m; i--; ) for (j = m; j--; ) Mat(templ,i,j) : (i == j) * 1 . 0 Mat(g,i,0) * M a t ( A , iter, j); matmul(templ,C_old,C); } /* i t e r a t i o n l o o p */ void
update
kgain(Matrix *C_old, M a t r i x *A, M a t r i x *temp, M a t r i x *g, int m, d o u b l e m _ v a r , int iter) { d o u b l e kdiv, s; int j, i ; for (i : m; i--; ) /* i' (i)C(i-l) { for (j --m, k d i v = 0.0; j--; ) k d i v += M a t ( A , iter, j) * M a t (C old, i,j) ; Mat(temp,0,i) = kdiv;
*/
} for (i = m, k d i v = m _ v a r ; i--; ) /* .... * i(i) k d i v += M a t ( t e m p , 0 , i ) * M a t ( A , iter, i); /* = k d i v for (i = m; i--; ) /* C ( i - l ) A ( i ) / kdiv { for (j = m, s = 0.0; j--; ) s += M a t ( C old, i,j) * M a t ( i , iter, j); Mat(g,i,0) = s / kdiv;
*/ */ */
} }
R eferences 1
S.C. Rutan and S.D. Brown, Anal. Chim. Acta, 160 (1985) 99-119. 2 S.D. Brown, Chemom. lnteL Lab. Syst., 10 (1991) 87-105. 3 S.D. Brown, AnaL Chim. Acta, 181 (1986) 1-26.
G. Gauglitz, M. Mettler and S. Weiss are at the Institut fEir Physikalische und Theoretische Chemie der Universit~t TE~bingen, Auf der Morgenstelle 8, T~bingen, Germany