1214
Anal. Chem. 1984, 56, 1214-1221
better or equally good results in half the recording time. Whenever the background cannot be identically recorded or when the capillary response is to be studied (6), the ADPP technique is preferable. Two conclusions can be drawn from the results independent of the technique. Firstly, when the background can be properly subtracted, the normal pulse polarography is presumed to be the pulse technique that gives the best precision. This technique will always produce a bigger or equal current response as compared to DPP. It will also be much less sensitive to changes in the degree of reversibility. Secondly, for pulse polarographic determinations the method of standard additions is expected to be superior to the use of calibration curves. This conclusion is induced from the better precision obtained when the solution is not removed from the DME. Naturally, straight calibration curves must be ascertained, but this can be done through standard additions to the proper supporting electrolyte.
ACKNOWLEDGMENT The authors thank I\ke Olin for helpful comments on the manuscript.
LITERATURE CITED (1) Shalgosky, H. 1.; Watling, J. Anal. Chlm. Acta 1962, 26, 66. (2) Rooney, R. C. Z . Anal. Chem. 1967, 224, 263. (3) Lingane, J. J. Anal. Chlm. Act8 1969, 4 4 , 411-424. (4) Bond, A. M.; Grabaric, B. S. Anal. Chem. 1979, 51, 337. (5) Stutts, K. J.; Dayton, M. A.; Wightman, R. M. Anal. Chem. 1982, 54, 998-1000. (6) Christie, J. H.; Jackson, L. L.; Osteryoung, R. A. Anal. Chem. 1976, 48, 242-247. (7) Bond, A. M. “Modern Polarographic Methods in Analytical Chemlstry”; Marcel Dekker: New York, 1980. (8) Chrlstie, J. H.; Osteryoung, R. A. J . Electroanal. Chem. 1974, 49, 30 1-31 1. (9) Parry, E. P.; Osteryoung, R. A. Anal. Chem. 1965, 3 7 , 1634-1637. (10) Baecklund, P.; Danielsson, R. Anal. Chlm. Acta, in press. (11) Baecklund, P.; Wlkmark, G. UUIC 82/A03, Uppsala Universlty, 1982. (12) Meites, L. “Polarographic Techniques”, 2nd ed.; Wlley: New York, 1965. (13) Brumleve, T. R.; Osteryoung, J. J . Phys. Chem. 1982, 86, 1794-1801.
RECEIVED for review September 28, 1983. Accepted March 15, 1984.
Estimation of Electrochemical Charge Transfer Parameters with the Kalman Filter Teri F. Brown
Department of Chemistry, University of Idaho, Moscow, Idaho 83843
Donna M. Caster and Steven D. Brown* Department of Chemistry, Washington State University, Pullman, Washington 99164-4630
A dlscusslon Is presented of the appllcatlon of the Kalman filter to systems where the model is neither linear nor In closed form. The system studied here Involves the modeling of linear scan voltammetric responses, where heterogeneous klnetlc parameters are estimated. Results are presented for synthetic data and for reduction of Cr(II1) and Pb(I1) in perchlorate medla. It is demonstrated, using state estimation, that the Cr( I I I ) reduction follows Volmerlan klnetlc law.
Recently, it has been recognized that modeling of responses obtained from analytical measurements has a number of important applications. Among these are smoothing, the ability to remove noise efficiently without distorting or biasing the signal, and deconvolution, the ability to accurately separate the effects of two processes whose responses couple with one another. Processes such as smoothing, filtering, and prediction come under the general heading of state and parameter estimation. In state and parameter estimation, computer algorithms are used to determine optimal estimates of either the dependent variables of a process (state estimation) or of the parameters of the mathematical model of the process (parameter estimation) (1-3). Both involve combining a mathematical model for a system with noisy measurements made on that system to generate estimates that are optimal in some statistical sense (least-squares, maximum likelihood, etc.) (4). One of the most commonly used estimators is the Kalman filter, a recursive linear estimator. Over the past few years, 0003-2700/84/0356-1214$01.50/0
this estimator has been applied to a number of chemical problems ranging from the removal of noise (5-7) and background (7)to the separation of overlapped responses in electrochemical measurements (7, 8) and in spectroscopy (9-11). These examples are direct applications of the Kalman filter to systems with linear models written in closed form equations. When the models describing the system are not linear, other filters must be used. This situation is important when performing state and parameter estimation in chemistry, since many systems are neither linear nor describable in closed form. An example of such a system is the response obtained in linear-scan voltammetry (LSV), when an analyte is reduced a t a planar electrode. This system is not generally solvable in closed form, although various approximate solutions have been given (12-15); such systems are usually described numerically by finite-differencesimulations (16-18).This paper demonstrates the use of the Kalman fiter in an extended form for state and parameter estimation in LSV systems obeying Butler-Volmer kinetics (19). The method presented here is applicable to any system with a well-defined but not necessarily linear nor solvable model and demonstrates the power of the Kalman filter for chemical state and parameter estimation.
THEORY The Kalman Filter. The fundamentals of the Kalman filter have been presented elsewhere (4, 7, 9), but certain aspects must be reviewed here for clarity. The filter combines a model of a system defined in terms of states of the system 0 1984 Amerlcan Chemical Society
ANALYTICAL CHEMISTRY, VOL. 56, NO. 8, JULY 1984
Table I. Equation for the Nonlinear Kalman Filter Defining Equations of the Nonlinear System System Dynamics X ( k ) = H ( X ( k ) ,X(k - 1)) + w ( k ) Measurement Function Z ( k ) = F ( X ( k ) ) =u ( k ) Noise Assumptions (IIIa) (IIIb) (IIIC)
v ( k )= N ( O , R ( k ) ) w ( k ) =N(O,Q(k)) E [ ~ ( k ) . w ( j=) ]0, for all k,j
Initial Conditions E [ X ( O ) I = x, E [ P ( O ) ]= E [ ( X ( O )- x,)zl = Po Equations for the Extended Kalman Filter Algorithm System Extrapolation X ( k I k - 1)~:H ( X ( k l k - l), X ( k - I l k - 1)) (VI P(klk - 1) = h ( X ( k - I l k - l ) ) . P ( k - I l k - 1). hT + Q ( k ) (VI) Estimate Update measurement Z ( k ) is taken, then G ( k ) =P ( k l k - l ) . f T ( X ( k l k- 1)). [ f P f T + R(k)]-’ (VII) X ( k Ik) = X ( k lk - 1) t G ( k ) . [ Z ( k-) (VIII) F ( X ( k Ik - 1 ) ) l P ( k Ik) = P(k lk - 1) - G ( k ) . f T ( X ( kIk - 1)). (1x1 P ( k l k - 1) where N ( A , B )= a white-noise process with mean A and autocovariance B E [ A ] = the expectation value of A X(k) = the vector of N state variables at time k Z ( k ) = the vector of M measurables at time k w ( k )= the vector of model noise at time k ~ ( k=) the vector of measurement noise at time k H ( k , k - 1)= the state transition matrix ( N x N ) between times k - 1 and k F ( k ) = the measurement function matrix ( M X N ) at time k where h ( X ( k - Ilk - 1))is the matrix whose ijth element is given bv X(t)=X(h Ih-1) X(-)=X(h-l Ik-1)
and where f ( X ( k 1 k - 1))is the matrix whose ijth element is given bv
(X) and their interactions with each other in time, coupled with noisy measurements on the system (Z) to yield statistically optimal estimates of the actual states of the model and their covariance matrix. In this discussion, the term “states” will include both states and parameters. The defining equations of the nonlinear filter are given in Table I; these equations are in state-space notation in matrix form. Equation I defines the model for system dynamics (how the states of the system change with time); eq I1 relates the states of the system to measurable quantities of the system; eq I11 lists necessary assumptions regarding the noise in the system; and eq IV states that an initial guess of the state variables and their covariance matrix is necessary to initiate the filter. If eq I and I1 were linear and the noise assumptions were satisifed, the Kalman filter algorithm would give estimates of the states and their covariance that would be statistically optimal in either the maximum likelihood or the least-squares
1215
sense. In fact, the results would be optimal, no matter what functional form of the cost function were used in the optimization (20-22). When the noise assumptions are not met as, for example, in signals containing low frequency noise or 60-Hz line noise, the estimates are not guaranteed to be optimal. However, in electrochemical systems, the filter has given results which are very nearly optimal in the presence of many nonwhite noise processes ( 7 , 2 3 ) . When eq I and I1 are nonlinear, it is necessary to use another form of the Kalman filter. The most common form, called the “linearized Kalman filter”, uses a Taylor’s series expansion of the nonlinear function about a trajectory and then truncates the expansion after the linear terms. When the trajectory used in the expansion is the present best estimate of X, the filter is called the “extended Kalman filter”. Use of this form of the filter results in a decrease in the optimality of the filter estimates. Here, the optimality is limited by two effects: the truncation error caused by expanding the function and then truncating the expansion after two terms, and the error due to the form of the nonlinearity. The truncation error can be estimated from the Taylor’s series expansion. For a function of a single variable x , the expansion about a is f(x) = f(a)
+ ( x - a)-dx df
I
]c=a
-1
(x - a ) 2 d2f ++ ... + 2! dXQx,a
When the expansion is truncated after the second term for the linearization, the truncation error is given as
where 0 5 0 5 1 (24). The truncation error will be large if the choice of trajectory is poor (i.e., ( x - a) is large) or if the second derivative of the function is large (i.e., the function changes its slope rapidly) (25). This dependence of the truncation error on the trajectory will make the extended Kalman filter sensitive to the accuracy of the initial guess. If the guess is very poor, the truncation error will cause the filter to converge very slowly, or even diverge. This behavior is not seen in the linear Kalman filter (4, 7). The second error is due to the difficulty of identifying (or distinguishing) two variables when they appear in certain groupings. Consider the function of two variables f(x,y) = axY (3) No matter how many measurements are made on f(x,y), it is not possible to estimate both x: and y simultaneously; only the product can be estimated. The same principle applies to parameter or state estimation using the Kalman filter. The most important factor governing the choice of the functional form of the model and the state variables is this need to have an identifiable system. The equations of the extended Kalman filter algorithm are presented in the bottom half of Table I. Equations V and VI, called the system extrapolation equations, calculate the best estimate of X and its covariance matrix, P,a t time k , based on only the system dynamics from time k - 1. After a measurement, Z, is taken, the next three equations, called the estimate update equations, are calculated. Equations VI11 and IX calculate the best estimates of X and P,a t time k , based on both system dynamics and the measurement. The gain matrix, G(k),calculates a variance-based weight to apply to the term in the estimate update equations due to the deviation of the system from the measurement. These equations look like the equations of the original Kalman filter, except for the matrices h ( x ) and f ( x ) . The definitions for these two
1216
ANALYTICAL CHEMISTRY, VOL. 56, NO. 8, JULY 1984
Table 11. The Potter-Schmidt Formulation of the Kalman Filter Kalman Algorithm Extrapolation X ( k I k - 1 ) = H ( k , k - l ) . X ( k - I l k - 1) P ( k l h - 1)= H(h,k - l ) . P ( k - Ilk - l ) * H T+ Q ( k )
Gain G ( k )= P(k Ik - l ) . F ( k ) . [ F . P . P T t R(h)]-’
(I1 (11)
with boundary conditions 1. a t t = O , x L O
(111)
Update X ( h Ik) = X ( k Ik - 1) t G ( k ) . Z ( k )- F ( k ) . X ( h Ik P ( k Ik) = P ( h Ik - 1)- G ( k ) . F T ( h ) . P ( kIk - 1)
-
1) (IV)
(VI
2. a t t > O , x - m
Potter-Schmidt Algorithm Extrapolation X(klk- 1)=H(k,k- 1).X(k- Ilk- 1 ) P”2(hlk - 1) = H ( k , k - 1).P1’2(k- 1 J k- l ) . H T
Gain
(VI) (VII)
3. a t t > O , x = O
(VIII) (VIIIa) (VIIIb) (VIIIC)
G ( h ) a.P”’(klh - l ) . S ( k ) where S(kJ = P T i z ( k 112 - l ) . F ( k ) l / a = S ( k ) . S ( k )+ R ( k ) b = (1 t Ja.R(12))-’
Update
1
(9) nFA where
X ( k Ik) = X ( h l h - 1)t G ( h ) . Z ( h )- F ( k ) . X J k llz - 1 ) (IX) P’”(k Ik) = P ” 2 ( h l h - 1)- a.b.P”’*S(k).S ( k ) (XI where P ’ l z is the upper triangular square root matrix of the covariance, P , given by P = P 1 i z . P T / 2and , PT/‘ is the lower triangular square root matrix of P and is the transpose of P’’’
matrices are given in Table I; they are simply the derivative terms in the truncated Taylor series linearizations. When the equations for the extended K h a n filter are used on systems which are ill-defined mathematically due to the form of the nonlinearity, the covariance matrix, P , can become negative indefinite. Although this is a theoretical impossibility, it often occurs when an ill-defined problem is calculated on a digital computer with finite word size (26,27). A series of alternate algorithms, called square root algorithms, were developed in the 1970s to overcome this problem. Of these the Potter-Schmidt algorithm is usually considered the best (4, 27,28); it has been used throughout the present work. This algorithm is presented in the bottom half of Table 11. The algorithm wes a Cholesky decomposition (29) prior to the f i s t time increment to calculate the square root of the covariance matrix in upper triangular form. The algorithm is then formulated to propagate the square root covariance matrix instead of the covariance matrix, and the equations calculating the gain and covariance update are reformulated to decrease the possibility of computer round-off causing the diagonal elements of the covariance matrix to become negative. These changes make the algorithm much more stable, without degrading its speed. Digital Simulation Model of Charge Transfer Kinetics. The system under study consists of a three-electrode electrochemical cell containing an electrolyte and one electroactive species. The electrolyte solution is treated as a semiinfinite fluid with the electrolyte-working electrode interface as the finite boundary. When the working electrode is a mercury drop, the drop can also be treated as a semiinfinite fluid with the interface as its finite boundary. Diffusion is treated as one-dimensional linear diffusion relative to the interface; neither migration nor convection is allowed. The electroactive species is presumed to be present initially in the oxidized form and its reduction at the working electrode is assumed to follow Butler-Volmer kinetics. The equations for this electrochemical system are (30)
0
+ ne-
kf kb
R
(4)
kb = k , exp
( E - EO’)
E = Ei - vt
]
(lob) (10c)
In these equations, Co and CR are the concentrations of oxidized and reduced species, in mol/L or mmol/cm3, Co* is the initial concentration of Cot CR* is the initial concentration of CR,DOand DR are diffusion constants in cm2/$,k , and k b are the forward and backward heterogeneous rate constants in cm/s, k8 is the standard rate constant in cm/s, A is the electrode area in cm2,i is the current in mA, v is the scan rate in V/s, a is the transfer coefficient, E is the potential a t time t in V, Eo’is the formal potential of the redox couple, Ei is the initial potential, and R, F , T, and n hold their usual meanings. These equations have not been solved analytically, but they have been solved numerically (13). The present work uses the approach introduced by Feldberg (16-18). In this approach, the equations are first converted to explicit finite difference equations and then made dimensionless by the following transformations:
Di At
Di* = 4x2 ki 4 t ki* = ax where At is the time increment of the simulation and Ax is the length of the simulation box. The solution phase is divided into NMAX boxes of length Ax. For a time increment, t , the equations calculate the amount of material that diffuses between each box and its two adjacent boxes. Boundary condition one gives the initial concentration in each box. Boundary condition two limits the number of boxes that need to be calculated during each iteration; since the perturbation starts in box one, the box nearest box one where CO(x,t) eo*,is taken as the present “infinite” boundary, and boxes beyond it need not be calculated, Boundary condition three, which only applies to box one, is the most complex because
-
ANALYTICAL CHEMISTRY, VOL. 56, NO. 8, JULY 1984
simulation model, and the Ci are constants
Table 111. Dimensionless Finite Difference Equations
Ax c1 = FAt
Initial Conditions f o ( x , O ) = 1.0 fR(X,O)
1217
= 0.0
Box 1
1 2Do*
c3
=-
c4
=2&*
1
Unfortunately, eq 15 is nonlinear in X and involves two variables, fo(l,t) and fR(l,t), which are not states in the digital filtering model but are states in the digital simulation model. If the right-hand side of eq 15 is denoted as B(X,t), the truncate9 Taylor's series expansion of B(X,t) about the trajectory, X, is given by
when
+ A t ) - 1.01 d 0.001 l f ~ ( X , t+ A t ) l d 0.001
lfo(x,t
N
and
then NBOX = x and proceed to next time step with boxes x > NBOX left unchanged
where N is the dimension of the vector X. Substitution of this expression into (15) gives a linearized form of eq 15
i(t) - B(X,t)
+ XT*b(X,t)= bT(X,t)*X
(18)
where it involves mass balance, charge transfer, and current flux. This equation states that a decrease in oxidized species at the interface results in a direct increase in the reduced species at the interface and that this change is equal to the charge transfer rate described by Butler-Volmer kinetics. Lastly, this equation shows that the current flux measured in the system is linearly related to the (net) rate of charge transfer. The equations for this system in dimensionless, finite difference form are given in Table 111. For heterogeneous charge transfer kinetics, the parameters of interest are the forward and backward rate constants, kf and kb. These are usually expressed in terms of the standard rate constant, k,, and the transfer coefficient, a. Consider a three-state model, where the states are given by ~ ( 1 =) &o*A
(12a)
Since each of these states is time independent, the problem is actually one of static parameter estimation, and the system dynamics matrix, H , reduces to the identity matrix
The measurement function equation relating the states of the system to the single measurable quantity, the current response, i, is given by
or in terms of the state variables as
i(t) = C,x(l)x(2) x
where i(t) is in milliamperes, fo(1,t) and fR(l,t) are the dimensionless concentrations in box one given by the digital
b,(X,t) = -
aB;:,t)lx=x
A comparison of this equation with the notation in Tables I and I1 shows that the measurement, Z(k), at time k = t , is reduced to the scalar Z(k) = i(t) - B(X,t)
+ XT*b(X,t)
(19)
where only a single measurable quantity (the current) is used, and that the matrix F(X(k)) corresponds to the vector b(X,t). The complete procedure is as follows: (1)Between measurements, the digital simulation model continually extrapolates the concentrations in each box using the equations in Table I11 and the present estimate of the parameters. (2) The parameters and their covariance matrix are extrapolated by the digital filter using eq VI and VI1 in Table
11. (3) When a measurement is taken, the digital filter uses eq VIII, IX, and X in Table I1 to generate an estimate update of X and P based on the measurement Z(k), as defined in eq 19, and the present simulated values of f&t) and fR(1,t). EXPERIMENTAL SECTION Reagents. Stock solutions of Pb2+and Cr(H20),9+were made approximately 1.0 X loF2M by dissolution of Pb(N03)2and Cr(C104)3.6H20,respectively, in water. Supporting electrolytes were made from reagent-grade NaOH (Mallinckrodt),NaC104(G. F. Smith), and 72% HCIOl (Alfa). None were purified further before use. Water used in making solutions was obtained from a Super-Q water purification unit; its resistivity exceeded 18 MQ. Cells were leached with reagent-grade HN03 between runs. Cells, Electrodes, and Equipment. Cells used were purchased from Bioanalytical Systems and had a capacity of 15 mL. The working electrode was a Metrohm hanging mercury drop electrode. The capillary was siliconized prior to use. A Model RE-6Ag/AgCl electrode, available from Bioanalytical Systems, was used as a reference. It had a potential of -41 mV relative to a saturated calomel electrode. The counterelectrode was a length of heavy gauge Pt wire. Solutions were deaereated by a stream of high-purity N2 gas which was scrubbed by a V(I1) solution and by distilled H20. Typically, 15-25 min of deaereation time were found to be adequate for the solutions containing lead, but 30 min of deareation was needed for the solutions containing chromium.
1218
ANALYTICAL CHEMISTRY, VOL. 56, NO. 8, JULY 1984
The computer used to control the voltammetric scans to collect current data and to extract electrochemical parameters consisted of an LSI-11/23 CPU with 128K words of semiconductor memory. In addition to the floppy disk data storage and graphics display capabilities, it was equipped with four parallel 1/0 ports, four 12-bit digital-to-analog converters,a real-time clock, and a 12-bit, 16-channelanalog-to-digitalconverter. One of the parallel ports was used to interface an Analog Devices 1136K 16-bit, digitalto-analog converter. The potentiostat used for these studies was an II3M EC225 unit, modified for external control with the computer. The 16-bit, digital-to-analogconverter was used to supply the initial potential, and a homebuilt interface was used to provide trigger signals for computer control of the EC225 linear ramp generator. The ramp was monitored by the computer, via an analog-to-digitalconverter, and the scan rates were determined by linear least-squaresfitting of the voltage ramp. Current responses were also monitored through another channel of the analog-to-digital converter. Computer Programs. Several programs were developed as part of this work. All ran on the LSI-11/23 under the RT-11 operating system. To generate synthetic voltammograms, a finite difference digital simulation program was written in FORTRAN. This program, called HETERO, allowed input of various electrochemical parameters and generated a simulated linear sweep voltammogram consisting of 512 sampled current responses obtained over a specified potential range. Simulated voltammograms obtained from HETERO were stored on floppy disk for later use; a typical simulation required about 2 min to run on the LSI-11/23. Noise was added to the data by using the FORTRAN program ADDNOZ,which is described elsewhere (23). Experimental voltammograms were obtained with a FORTRAN driver coupled to a MACRO-11 assembly language data acquisition and control routine. This program was similar to ones previously reported (7). Data were collected by the computer and stored on floppy disk. A background subtraction routine written in FORTRAN was used to process the experimental data prior to input to the filter. No other preprocessing was required. The Kalman filter was programmed in FORTRAN by use of the Potter-Schmidt algorithm. Program KFHET required input of the voltammogram, estimates of various electrochemical parameters, and initial guesses of the variance in the voltammogram and of the values and variances of the parameters to be fitted It output the corrected values of the parameters and the covariance matrix and plotted the original and fitted voltammograms. Further filtering could be performed if desired, but this option was not used here, since iteration did not improve the quality of the fit, and frequently led to the filter’s diverging. It was also possible to output the states’ variation with potential in addition to the other quantities discussed. The program required about 3-4 min to process a 512-point voltammogram.
RESULTS AND DISCUSSION Fitting of Synthetic Data. T o thoroughly examine the capabilities of the Kalman filter for state and parameter estimation, two studies were performed. In the first study, the filter was applied to digitally simulated voltammograms. Forty-five simulated linear sweep voltammograms were generated by the program HETERO. These voltammograms covered a wide range of values for the heterogeneous rate constant k,, the transfer coefficient, the simulated sweep rate, and diffusion constants. In this manner, values could be selected to represent charge transfer behavior ranging from strictly reversible to completely irreversible. The peak potentials and peak currents observed for these simulated systems were compared to the results of Nicholson and Shain (23)for the reversible and the totally irreversible systems, and to the results of Matsuda and Ayabe (12) for the quasi-reversible systems. In all cases, the peak parameters obtained by examination of the simulated voltammograms agreed with those in the literature. To these simulated voltammograms, noise was added to give a signal-to-noise ratio of 10. Three noise processes were used: random noise, band-limited noise, and 6 0 - H ~noise. These noise processes were used because they best represent the
Figure 1. The effect of the heterogeneous rate constant and transfer coefficient on the shape of a linear scan voltammogram. Curves 1 through 4 have k , = 1.0 cm/s and cy increasing from 0.2 to 0.8; curves 5 through 8 have k , = 1.0 X cm/s and a increasing from 0.2 to 0.8; curves 9 through 12 have k , = 1.0 X lo-’ cm/s and a increasing from 0.2 to 0.8. In all cases €’’ = -0.2500 V, D o = D, = 1.0 X cm2/s, and the scan rate = 1.0 V/s.
types of noise commonly encountered in electrochemical measurements (23). The random noise was generated in the time domain by use of a 12-point, Gaussian pseudorandom number generator; the noise was added to the simulated response point-by-point. Band-limited noise was added by using the same generator in the frequency domain. Here, noise was added to the data over a small range of frequencies, and the dc signal was used to establish the desired signal-to-noiseratio. Single frequency (60 Hz) noise was also added in the frequency domain in the same manner. In all cases, the structure of the noise was checked by examination of the autocovariance function and the Fourier power spectrum as described elsewhere (23). When no noise was present, the filter provided excellent fits to simulated irreversible voltammograms, regardless of the accuracy of the initial guess. Both a and k , are obtained with essentially no error. When noise was present, the fitting was found to depend somewhat on the initial guess. Results from studies of the fitting of noise-corrupted synthetic voltammograms are shown in Table IV. Judging from this study, fitting of voltammograms corrupted with band-limited noise is superior to fitting of those corrupted by either of the other two noise processes, and fits using initial guesses where k , is lower than the actual value are superior to those whose initial guess for k , is high. In most cases, the quality of the estimated cy and k, values is very good; the coefficients of determination for the fitting ranged a t or above 0.99. When k, is increased and the simulated voltammogram shows quasi-reversible or reversible behavior, accurate estimation of parameters by the Kalman filter is impossible, because the shape of the voltammogram becomes less dependent on cy and k,, as shown in Figure l. Thus, a wide range of values for cy and k , give equally acceptable fits, and the fitting errors become large. It should be noted, however, that this situation arises in all methods based on the modeling of the shapes of responses, whether parametric (31)or empirical (32)methods are used. In general, precise evaluation of parameters will not be possible unless changes in curve shape occur with relatively small changes in the parameters; for quasi-reversible systems, methods based on convolution voltammetry may be applied (31). Fitting of Model Charge Transfer Systems. The second study involved fitting of unfiltered experimental data. Two model systems were used to demonstrate the effectiveness of using the Kalman filter to extract experimental k , and cy
ANALYTICAL CHEMISTRY, VOL. 56, NO. 8, JULY 1984
rl-
I1
a k
1219
= -1.28l I I -0.68 I !-8,21 I I -t
L1,48
-1.00
-0.80 POTENTIAL
-0.40
Flgure 2. Experimental data for reduction of Cr(H20)t+ in 0.5 M NaC10,/0.01 M HCIO,, at a scan rate of 2.99 VIS. The Kaiman filer-generated fit to the data is also shown,and the residuals for the fit are plotted. Heterogeneous parameters obtained from this fit are reported in Table V.
.-e
-
0 .-e*
.zz.!. 0 0 0
""rn
0 0 0
.
parameters from voltammograms. One system involved the reduction of Cr(H20)6s+in 0.50 M NaC104/0.01 M HC104. This system has received detailed study as a model for double-layer effects in electron transfer reactions (33-38). Several groups have measured the transfer coefficient and heterogeneous rate constant; their results are reported in Table V. Because this system is totally irreversible at all scan rates convenient for study by us, we have selected it for detailed study with the filter. In general, data were taken, and fits were made, in duplicate. A typical fit obtained from the application of the Kalman filter is shown in Figure 2. As with the synthetic data, the quality of the fit is excellent. The lack of significant structure in the residuals supports the claim that the digital simulation model, as adjusted by the filter, provides an accurate description of the chromium(II1) reduction in this medium. In these fits, previously measured values of the electrode area, Cr concentration, and formal potential (EO'= -0.650 V vs. SCE) were used. Diffusion constants were taken from ref 34. Various initial guesses were employed, but all gave results within the error bounds indicated. Results from fits by the Kalman fiter to this system are summarized in Table V. From those results, variance-weighted means of a = 0.61 f 0.02 and k , = (5.2 f 0.8) X lo4 cm/s may be calculated for the fitted parameters; these values are not corrected for double-layer effects (39). The second system, Pb2+in 0.9 M NaC104/0.1 M NaOH, was of interest because its reversibility is strongly dependent on the scan rate, ranging from strictly reversible at scan rates below about 0.01 V s-l to quasi-reversible at scan rates above 1 V s-l. This system is a good test of the filter, since the changes in shape for the P b response with scan rate will be fairly small, in view of the system's degree of reversiblity. Here again, previously measured values for the electrode area and Pb concentration were employed. Diffusion constants from ref 40 were used. The formal potential was observed to be -0.706 V vs. the SCE. As before, various initial guesses were tried. As expected, the quality of the fitting depended on the relative reversibility of the Pb2+ system's response; fitting improved as the system became less reversible. Results from fitting are reported in Table VI, and a typical fit produced by the Kalman filter to the experimental data is shown in Figure 3. The accuracy of the parameters estimated by the filter remains quite good in general, but the error terms associated with those estimates increase markedly. The variance-weighted means are k, = (6.23 f 0.14) X cm/s and a = 0.37 i 0.03. State Estimation. Theory (36) and previous work (33) on the reduction of Cr(H20)sS+in perchlorate media have suggested that the transfer coefficient a may depend upon potential, contrary to the behavior expected from Butler-Volmer
1220
ANALYTICAL CHEMISTRY, VOL. 56, NO. 8, JULY 1984
Table V. Heterogeneous Parameters for Cr(II1) Reduction in 0.5 M NaC10,/0.01 M HClO, method LSVC LSV LSVC LSV LSV LSV LSVC LSVC NPP CVf Polg
scan rate, V/s 0.997 2.99 2.99 5.01 6.98 6.98 6.98 6.98
k , , a cm/s ffa CDb ref (5.06 i 0.13) x 0.60 i 0.01 0.9992 this work (4.81 i 0.19) x 0.62 f 0.01 0.9985 this work (4.83 i 0.19) x 0.62 t 0.01 0.9981 this work (5.67 t 0.19) x 0.611 f 0.004 0.9995 this work (5.20 i 0.22) x 0.618 i 0.004 0.9992 this work (4.58 i 0.09) x 0.627 i 0.004 0.9992 this work (6.01 i 0.19) x 0.610 i 0.004 0.9991 this work (5.67 i 0.06) x 0.613 i 0.004 0.9991 this work (3.5 i 0.2) x 0.63 0.01 36, 37 0.05 5.0 X 0.55 35 8.1 x 0.60 34 a Results not corrected for double-layer effects. Coefficient of determination. Linear sweep voltammetry; analysis of wave using Kalman filter and digital simulation. Normal pulse polarography; analysis of Tafel plots. e Measured for 0.5 M NaC10,/0.005 M HClO,. Calculated from E"' = -660 mV (vs. SCE), -log K at -900 mV = 2.40 i 0.02, T = 25.0 "C. Measured by cyclic voltammetry; analysis of Tafel plots. 1M NaC10, supporting electrolyte. g Manual polarography, used to generate Tafel plots.
--
Table VI. Heterogeneous Parameters for Pb(I1) Reduction in 0.9 M NaClOJO.1 M NaOH method scan rate, Vis k,, cmis 01 CDa ref LSVb 4.99 (5.91 i 0.15) x 10.' 0.43 I0.03 0.9953 this work LSVb 4.99 (6.05 i 0.18) x l o - * 0.45 i 0.03 0.9954 this work LSVb 6.95 (5.94 i 0.12) x 10.' 0.39 t 0.03 0.9980 this work LSVb 6.95 (6.04 c 0.13) x 0.40 i 0.03 0.9981 this work LSVb 6.95 (6.61 i 0.13) x lo-' 0.9959 this work 0.32 i 0.02 LSVb 6.95 (6.72 i 0.13) x 10.' 0.9963 this work 0.33 i 0.02 FI (6.0 t 0.3) x 10.' 0.51 i 0.01 40 LSVd 1.00 (6.5 i 0.5) x lo-' 0.38 i 0.03 32 Coefficient of determination. Linear scan voltammetry; anal sis of wave using Kalman filter and digital simulation. Faradaic impedence; analysis of impedence vs. frequency data. 2Linear scan voltammetry; analysis of wave using k-nearest neighbor pattern classifier.
nCo' a
k, I N N O V A T I O N SEOUENCE O A I G I N R L VOLIANNOGARM
B Lp
I ?l,48
F -1.28
-1.00
-8.00
; -8.60
-hi18
: -8.20
t,BB
, 8.28
POTENIIAL Figure 3. Experimental data for reduction of Pb(I1) in 0.9 M NaCiO,/O.l M NaOH, at a scan rate of 6.95 V/s. The Kalman filter-generated fit to the data is also shown, and the residuals for the fit are plotted.
-'1.48
-1.28
-1.88
-8.88 -0.68
-0.48
-0.28
0.88
8.28
POTENTIAL
Heterogeneous parameters obtained from this fit are reported in Table VI.
Figure 4. Evolution in time of the states included in the Kalman filter for a fit to an experimental voltammogram. Also shown is the original voltammogram obtained for the reduction of Cr(H,O)z+ in 0.5 M NaCIO,/O.Ol M HCIO., at 2.99 V/s. The maximum value for each of the plotted variables has been normalized to 1.0 for ease in comparison. I n this case, a poor initial guess was provided to the filter.
law but consistent with that expected from Marcus theory (38); others (34, 36) have disputed this claim, but the possibility remains open, since stringent tests are not easy to design and since the predicted dependencies with potential are small (36, 38). Observation of the dependence of a on potential here is a simple extension of the estimation of a by the filter, since a calculation is made for the parameter a at each measured point, (Ek, i k ) . In this case, the assumption that X is a constant is removed, and plots of the dependence of the state variable 4 3 ) (which is equal to a ) on potential may be obtained. Figure 4 shows a plot of the dependence of the state variables on potential for the reduction of Cr(H20),3+in perchlorate a t a scan rate of 2.99 V s-l. Also included in this figure are the innovations sequence and the experimental voltammogram.
All are scaled to one for comparison of their respective dependences on potential. All of the states change rapidly as the filter endeavors to track the measurement. Once the measurement and filter results coincide, changes in state variables appear small. To determine whether a changes with potential, the filter's best estimate of the states are entered as the initial guesses. Thus, the filter need not change any states to track the data, unless those states are changing with potential. Figure 5 shows the dependence of a on potential for a poor initial guess and with a guess coinciding with the filter's best estimate. When an accurate guess is given, little change in a occurs with potential. For this voltammogram, the maximum change in a is 0.002, well within the experimental error in data acquisition. It appears that any change in a is attributable directly to data collection errors; a similar
ANALYTICAL CHEMISTRY, VOL. 56, NO. 8, JULY 1984
: I + POOR 6UCSS
6DOD N E S S
a
y1
Pl,M - I , ? #
-1,N
-i.SS
-8.61
-0.41
-0.28
0.88
I
8.28
POTENTIAL Figure 5. Evolution in time of the state x ( 3 ) , corresponding to the transfer coefficient, for the voltammogram shown in Figure 4. Results are shown for a poor guess of x ( 3 ) and for the guess corresponding to the optimal x(3).
conclusion was reached by Anson et al. (34) by other means.
CONCLUSIONS It has been demonstrated here that the Kalman filter can be coupled to a nonlinear, nonclosed-form model to produce results of high quality. In this example, use of the Kalman filter allows extraction of kinetic parameters on a single linear sweep voltammogram in a few minutes, with results being of equal (or possibly superior) quality to those obtained by the more tedious generation and analysis of Tafel plots. In this manner, systems unsuited to analysis by conventional means, those which decompose quickly in solution, for example, can be studied directly, since the data are obtained quickly. Additionally, the opportunity to follow the change in electrochemical parameters with potential offers a valuable tool to studies of electrode kinetics, since few tests of adequate sensitivity exist at present. Of course, the model discussed here is not the only one suited to coupling with the extended version of the Kalman filter. Other finite-difference models can also be combined with the filter. Subsequent resports will discuss results from these studies. Registry No. Cr(Hz0)63+, 14873-01-9;Pb, 7439-92-1.
LITERATURE CITED (1) Beck, J. V.; Arnold, K. J. "Parameter Estimation in Engineerlng and Science"; Wiley: New York, 1977.
1221
Brubaker, T. A.; Tracy, R.; Pomernacki, C. L. Anal. Chem. 1978, 50, 1017A. Brubaker, T. A.; O'Keefe, K. R. Anal. Chem. 1979, 51, 1385A. Geib, A., Ed. "Applied Optimal Estimation"; MIT Press: Cambridge, MA, 1974. Seelig, P. F.; Blount, H. N. Anal. Chem. 1978, 48, 252. Seelig, P. F.; Blount, H. N. Anal. Chem. 1979, 57, 327. Brown, T. F.; Brown, S. D. Anal. Chem. 1981, 53, 1410. Brown, T. F.; Caster, D. M.; Brown, S. D. NBS Spec. Pub/. 1981, No. 678, 163. Poulisse, H. N. J. Anal. Chim. Acta 1979, 772, 361. Didden, C. B. M.; Pouiisse, H. N. J. Anal. Left. 1980, 73 (Al), 921. Poulisse, H. N. J.; Engelen, P. Anal. Lett. 1980, 73 (A14), 1211. Matsuda, H.; Ayabe, Y. 2.Elektrochem. 1955, 59, 494. Nicholson, R. S.; Shah, I. Anal. Chem. 1964, 36,706. Oldham, K. B.; Spanier, J. J . Nectroanal. Chem. 1970, 26, 331. Imbeaux, J. C.; Saveant, J. M. J . Nectroanal. Chem. 1973, 44, 169. Feldberg, S. W. Nectroanal. Chem. 1989, 3, 199. Bard, A. J.; Faulkner, L. R. "Electrochemical Methods"; Wiley-Interscience: New York, 1960; Appendix B. Feldberg, S. W. I n "Computers in Chemistry and Instrumentation"; Mattson, J. S., Mark, H. B., Jr., MacDonaid, H. C., Eds.; Marcel Dekker: New York, 1972; Vol. 2, Chapter 7. Bard, A. J.; Faulkner, L. R. "Electrochemical Methods"; Wiley-Interscience: New York, 1980; Chapter 3. Kalman, R. E. J . Basic Eng. 1980, 82. 34. Kalman, R. E.; Bucy, R. S. J . BasicEng. 1961, 83, 95. Kalman, R. E. I n "Proceedings of the First Symposium on Engineering Applications of Random Function Theory and Probabllity"; Bagdanoff, J. L., Kozin, F., Eds.; Wiley: New York, 1963; pp 270-388. Brown, T. F. Ph.D. Dlssertatlon; University of Washington, Seattle, WA, 1962. Franklin, P. "Dlfferentiai and Integral Calculus"; Dover: New York, 1969; Chapter 16. Bryson, A. E.; Ho, Y.-C. "Applied Optimal Control"; Halstead Press: New York, 1975. Kamlnski, P. G.; Bryson, A. E.; Schmidt, S. F. IEEE Autom. Control 1971, AC-76, 727. Potter, K. C. in "Astronomical Guidance"; Battin, R. H., Ed.; McGrawHili: New York, 1964. Bierman, 0. J. IEEE Aerosp. Electron. Sys. 1973, AES-9, 28. Fadeeva, V. N. "Computational Methods of Linear Algebra"; Dover: New York, 1959; pp 81-84. MacDonald, D. D. "Transient Techniques in Electrochemistry"; Plenum Press: New York, 1977; Chapter 6. Toman, J. J.; Brown, S. D., submitted to J . Electroanal. Chem. DePalma, R. A.; Perone, S. P. Anal. Chem. 1979, 57, 829. Parsons, R.; Passeron, E. J . Electroanal. Chem. 1966, 72, 524. Anson, F. C.; Rathjen, N.;Frisbee, R. D. J . Electrochem. SOC. 1970, 777, 477. Zieiinska-Ignacluk, M.; Galus, 2. J. €lectroanal. Chem. 1974, 50, 41. Tyma, P. D.; Weaver, M. J. J . Electroanal. Chem. 1980, 7 1 7 , 195. Weaver, M. J.; Lin, H. Y.; Kim, Y. Can. J . Chem. 1981, 59, 1944. Weaver, M. J. Isr. J . Chem. 1979, 18, 35. Frumkin, A. N. 2.Phys. Chem. 1933, 764, 121. Biegler, T.; Laitinen, H. A. Anal. Chem. 1965, 37, 572.
RECEIVED for review April 29,1983. Accepted February 27, 1984.
