1001
Anal. Chem. 1993, 65, 1061-1068
Kalman Filter-Optimized Simulation for Step Voltammetry Robert S. Bear, Jr. and Steven D. Brown* Department of Chemistry and Biochemistry, University of Delaware, Newark, Delaware 19716
The technlque-spectflcnature of electrochemlcalslmulatlons llmlts thelr appllcatlon for modellng the klnetlcs of heterogeneous charge-transfer reactlons. llme-consumlng curvematching algorlthmswed for flttlng dmulatedvoltammograms toexoerimentaldata are an addltbnal Ilmltatlon. I n thls paper, a gemrallzedmethod for slmulatlngthe v o t t a m t r k response of a heterogeneouselectron-transfer reactlon to a potentlalstep sequence Is presented. By coupllng the slmulatlon wlth an extended Kahnan fllter, a recurdve optlmlzatlon of the dmulatlon Is performed to flt experlmental data. Estlmates of the heterogenous charge-transfer parametersare obtalned In a dngle paso through the simulation. The reductlons of Cr(H20)sa+In 0.5 M NaC104-0.01 M HC104and Zn( 11) In 1.0 M KN03 are presented as example appllcatlons.
INTRODUCTION Electrochemical are widely used to model systems for which no closed-form solution exists. The use of these methods to obtain kinetic parameter estimates from experimental voltammograms requires a means to optimize the simulationmodel so that an optimal fit to the experimental data is generated. Manual optimization, done by visual matching of a simulation to the experimental voltammogram, is routinely used. While manual optimization is easy to implement, it is extremely time consuming because many trial simulations must be run to produce the candidate simulation models. Another disadvantage is that the parameter estimates obtained from such visual matches lack statistical validity. Simulations can also be coupled to nonlinear regressions to obtain parameter estimates through repetitive simulation runs.3-5 Because of the large number of iterations required for the estimates to converge, these methods also use an extensive amount of computer time to fit a single voltammogram. A more computationally efficient approach is to perform the optimization of the simulation model to best fit the experimentaldata recursively. Brown et al.6 used an extended Kalman filter to recursively optimize an explicit finite difference (EFD) simulation for linear sweep voltammetry. Little computational overhead is used by the filter, so that the parameter estimation is performed in essentiallythe same time required for the simulation. Estimates of the standard rate constant and transfer coefficient for a simple heterogeneous electron-transfer reaction are readily obtained by this method. Kalman filter-based optimizations of cyclic voltammetric simulations have also been reported and compared with nonlinear optimizations of simulation models
by simplex and Marquardt methods? While applications of this optimization and fitting method have thus far been limited to fitting data obtained by linear sweep and cyclic voltammetry, in principle a K h a n filter-optimized simulation can be used with any current-potential-time response. This study was focused on developing an extension to work previously reported on the Kalman filter methods that could be applied to a wide range of voltammetric techniques. Specifically,techniquesthat are defined as a series of potential steps or pulses were considered. An EFD simulation that models the voltammetric current-potential-time response produced by a potential-step sequence was developed. The simulation has been generalized, so that any potential-step sequencecan be used. Additionally,a correctionfor the effects of the double layer on the reaction kinetics was incorporated into the model. Optimization of the simulation model is performed in a single pass through the data to be fitted by the use of an extended K h a n filter. This generalized approach allows parameter estimates to be made independent of the potential-step protocol used to collectthe voltammetric data. To demonstrate the performance of the generalized approach, four different potential-step sequences were used in both synthetic and experimental studies. Using this method, double-layer-correctedkinetic parameter estimates were obtained for two typical heterogeneous electron-transfer reactions. These parameter estimates were in good agreement with published values, regardless of the specific potentialstep sequence used.
THEORY Simulation Model. In this investigation, the example of a simple heterogeneous electron-transfer reaction is considered. The reaction takes place at the surface of a hanging mercury drop electrode (HMDE), which is under potentiostatic control in an unstirred solution of the reactant species. The simple heterogeneous reduction reaction is represented as O+ne-*R (1) where 0 and R represent the oxidized and reduced species and n is the number of electrons. This reaction and the diffusional transport processes in spherical coordinates are described by the following set of equations and their initial and boundary conditions:
* Author to whom correspondence should be addressed.
(1)Feldberg, S.W.Electroanal. Chem. 1969,3, 199-296. (2)Britz, D. Digital Simulation in Electrochemistry; SpringerVerlag: Berlin, 1988. (3) Arena, J. V.; Rusling, J. F. Anal. Chem. 1986,58,1481-1488. (4)Gampp, H.Anal. Chem. 1987,59,2456-2460. ( 5 )Evans, D.H.; Jimenez, P. J.; Kelly, M. J. J . Electroanal. Chem. 1984,163,145-157. (6)Brown, T.F.;Caster, D. M.; Brown, S. D.Anal. Chem. 1984,56, 1214-1221. 0003-2700/93/0365-1061$04.00/0
(4)
(7) Lavangnini, I.; Pastore, P.; Magno, F. Anal. Chin. Acta 1989,223, 193-204.
0 1993 American Chemical Society
1062
ANALYTICAL CHEMISTRY, VOL. 05, NO. 8, APRIL 15, 1993
[O] = [o]b,&
at all r, when t = 0
[Rl = 0; at all r; when t = 0
-
(6)
(7)
[O]= [O]bulk; at r QJ,when t > 0 (8) [R] =O; atr-,-, whent>O (9) where DOand D Rare the diffusion constants for the oxidized and reduced species, r and t are the radial and time coordinates, and r, is the radius at the electrode surface. Butler-Yolmer kinetics6 are assumed, where the potentialdependent forward and reverse rate constants in eq 2 are related to the heterogeneous standard rate constant (k,)and transfer coefficient (a)by
?:
k, = k, exp(- -(E
Kalman Filter. In this work the simulation described above is used as the nonlinear measurement model for an extended Kalman filter. Descriptions of the Kalman filter have been given elsewhere in the chemical literature,6JOJl and a brief summary of the extended Kalman filter algorithm is given in the Appendix. The state vector (x),a vector comprised of the set of parameters to be estimated, consists here of a concentration adjustment coefficient (y), the heterogeneous standard rate constant, the transfer coefficient, and the instrumental offset
($7. x=k]
-EO/))
(12)
(10) (11)
where F, R, T,E, and Eo' have their usual electrochemical definitions. An EFD algorithm is adequate to model the simple kinetic scheme considered here. The EFD simulation is also convenient because it is easily coupled to the discrete, extended Kalman filter. The electrochemicalsystem described above was simulated using an EFD method with even time and spacegrids,ls2 modified for use with potential-step techniques. To correctly model the response produced in a potential-step experiment, it is necessary to calculate the flux present at the end of each time increment. An estimate of the flux is first calculated at the beginning of each time increment. From this estimate, the concentrationprofiles are updated, changing the concentration gradients present at the interface. Using the new concentration gradients, a calculation of the flux is then run to provide a more accurate estimate of the current at the end of the time increment. A five-point approximation was used in this study to obtain the desired concentration gradienta.9 This scheme provides an increase in the accuracy of the simulation for modeling potential-step experiments, while retaining the simplicity of the explicit method. To make the simulation and kinetic parameter estimation method described here independent of the specificpotentialstep protocol, the applied potential must be defined externally as a time series of potential values. Techniques like staircase voltammetry, where the spacing of data points is determined by the width of the potential step, require each time step in the data set to be subdivided into a larger number of simulation time increments. These subdivisionsare required for the EFD simulation to properly model the changing voltammetricresponse across a potential step. Voltammetric protocols that require the potential to the sampled at intervals across a single potential step can be handeled by defining the step as a series of time increments with the same potential value. Each potential step appears as a separate chronoamperometric experiment, with the current being sampled at intervals throughout its duration. As the number of sampling intervals for a single potential step is increased, the number of time-increment subdivisions required in the simulation is decreased. When the number of data points across the step is large, the time increment of the simulation can be set equal to the spacing of experimental data points. Using this approach, staircase voltammetry, square-wave voltammetry, normal-pulse voltammetry, or any other voltammetric protocol that can be describedby a series of discrete potential steps can be modeled with a single simulation. ( 8 ) Bard, A. J.; Faulkner, L. R. Electrochemicul Methods; John Wiley and Sons, Inc.: New York, 1980. (9) Britz, D. Anal. Chim. Acta 1987, 193, 277-285.
While the concentrations of each of the test solutions are known, the incorporation of a concentration adjustment parameter allows for compensation of small uncertainties in both the concentration and electrode surface area. Additionally, the incorporation of a offset term in the state vector prevents minor offsets from interfering with the kinetic parameter estimation. The Butler-Volmer kinetic model used is based on the assumption that the kinetic parameters, contained in the state vector, are time and potential invariant. Because of this property, the system dynamics (eq A l ) can be simplified. When the additional assumption that the model noise is zero is made, the system dynamics matrix is reduced to I, the identity matrix. The measured response for the system-the current passed through the working electrode-is expressed in terms of the dimensionless flux calculated by the simulation. This relationship defines the nonlinear measurement model for the extended Kalman filter (eq A3). Coupling of the simulation to the filter is accomplished through6
h(k) = nFA'YIOIbulkJ(k,,a) &/At + C in terms of the filter states, this relation becomes
(13)
where A is the electrode surface area and J is the flux calculated in the EFD simulation. In eqs 13 and 14, the index k refers to the kth data point in the experimental voltammogram contained the vector z in eq A3. Filtering an experimental voltammogram with the EFD simulation as the filter measurement model amounts to the estimation of the parameters in the state vector at each data point k, so that the measurement model best fits the experimental data. The state parameters control the simulation, and so the process of generating an optimal fit to the experimental data simultaneously optimizes the simulation. This process is done recursivelyand the simulationparameters are adjusted at each data point; convergence of these parameters to stable values is usually within a single pass through the experimental voltammogram. All that is needed to begin the recursive fitting and optimization is an experimental voltammogram and an initial guess for the state parameters, as well as a guess of their error covariance. As is demonstrated below, an accurate initial guess is helpful in speeding convergence but is not essential to obtain satisfactory estimates in most cases. The output of the filter consists of the state vector, containing the optimal parameters for the simulation model, and the estimated error covarianceof the states. Because, in essence, the operation of filtering treats (10)Brown, S. D. Anal. Chim. Acta 1986, 181, 1-26. (11)Rutan, S. C. J. Chemom. 1990,4, 103-121.
1063
ANALYTICAL CHEMISTRY, VOL. 65, NO. 8, APRIL 15, 1993 101
I
I
I
0
0.5
I
0
1
Time (S)
0.5 Time (S)
I
Time (S)
Time (S)
:.o.:k 30
30
0
I
a"
-0.2 I
I
-0.1
a"
-0.2
1
1
0
3 8
0.5 Time (S)
1
0
I
0.5 Time (S)
1
Flgurr 1. Expanded view of the applled potential functlons. Only the first second Is dlsplayed for clarity. The complete applied potential functions consisted of 2500 time steps: (a) conventlonal staircase voltammetry; (b) iarge-amplltude staircase voltammetry;(c) randomly perturbed Staircase voltammetry;(d) square-wavevoltammetry. The time Interval between sampling points was 0.001 s. the electrochemical data as a time series, with each point k of the experimental voltammogram z processed one by one, it is possible, given sufficient computing power, to implement the method discussed here in real time.
EXPERIMENTAL SECTION Preparation of Solutions. All solutions were prepared from reagent-grade chemicals,and distilled-deionized water was used throughout. The stock solutions of Cr(Hz0)c3+and Zn(I1) were made slightly acidic by the addition of a small amount of dilute "03. The supporting electrolyte used for the Cr(Hz0)e3+ reduction was 0.5 M NaC104-0.01M HC104,and that used for the Zn(I1) studies was 1.0 M KN03. Instrumentation. All data were collectedusing a homemade potentiostat, which was interfaced to a 80286-based microcomputer. Interfacing software was written in the TURBO PASCAL programming language. Potential-stepsequences,corresponding to the different voltammetric techniques, were supplied to the potentiostat as an ASCII listing of the potential values. A Metrohm E410 hanging mercury drop electrode (HMDE) was used as the working electrode. Enough Hg was dispensed to give a drop radius of 0.49 mm. A platinum wire was used as the counter, and a Ag/AgCl reference electrode was used. All simulations and Kalman filter routines were written in the MATLAB language (The Mathworks, Inc., 21 Eliot St., South Natick, MA 01760) and were run on an Apple Macintosh IIcx and a 80486-based microcomputer. Potential Functions. Experiments were performed using four types of potential-step sequences,all of which were generated in the MATLAB environment and saved as ASCII files. The first sequence corresponded to a conventional staircase voltammetry experiment, which was designed to provide a close approximation to a linear sweep. The potential-step size was restricted to be less than 1 mV, and the response current was sampled at the end of each step. A modified form of the staircase experiment was also used. In this staircase technique, the potential steps were of much larger amplitude (severalmillivolts) and the response current was sampledat evenlyspaced increments across the entire width of the step. The third experimental technique performed was a modified form of square-wave voltammetry, in which the response current was sampled across each pulse. In the final technique, a small-amplitude, random pulse sequencewas superimposedonto a large-amplitudestaircase experiment. Figure 1 illustrates each of the applied potential functions considered in this paper. The respective simulated response currents for a system exhibiting quasi-reversible kinetics are displayed in Figure 2. Each potential step in the conventional
10-
0
1
Time (S)
2
-101 0
I 1
2
Time (S)
Flgurr 2. Simulatedresponse currentsfor each of the applied potential functions shown in Flgure 1. The parameters used to generate these responses were P' = -0.200 V, k, = 1 X lo-* cm s-l,a = 0.5, At = 0.001 s, n = 1, A = 0.025 cm2,analyte & = 0.001 M, and & = 4=5 X cm2s-l: (a) conventional Staircase voltammetry; (b) large-amplltudestalrcasevoltammetry;(c)randomly perturbed staircase voltammetry; (d) squarawave voltammetry. staircase voltammogram has a single data point, which corresponds to the current present at the end of the step. In the two large-amplitude staircase techniques each potential step corresponds to 50 time increments, with the current being measured at the end of each. The remaining technique, square-wave voltammetry, incorporated 25 time increments to represent each potential step. Procedure. The cell assembly was contained in a tempered bath maintained at 25.0 A 0.1 "C. Each of the solutions was purged with NBfor a period of 15 min prior to data collection. Avoltammogram of the supporting electrolyteblank was recorded for each technique. Three replicates were recorded for each voltammetric run; the mean was used for the filtering routines. The size of the time steps was chosen with the constraint that an adequate suppression of the charging current was obtained for all data points.
RESULTS AND DISCUSSION The performance of the Kalman filter-optimized simulation method was characterized through a series of evaluations. These evaluations, described in detail below, were conducted in three distinct areas. First, the adequacy of the simulation model was validated through comparison of the optimized simulation model with theoretical results. Next, synthetic voltammograms were used to study the performance of the Kalman filter at fitting an optimized EFD simulation model under different conditions. Convergence properties, accuracy of estimates, and precision of estimates were all evaluated in these studies. In the final series of evaluations, the Kalman filter was used to obtain charge-transfer parameter estimates for two typical heterogeneous electron-transfer reactions studied by several pulse voltammetric methods. Validation of the Simulation. The accuracy of the simulation model was validated through comparisons with theoretical results for two simple systems. The first performance benchmark used was a simulation of the Cottrell experiment, a system for which an analytical solution exists. As the response current for a system exhibiting Cottrell-like behavior is readily calculated, this comparison provided a simple and accurate method to evaluate the performance of the simulation. The EFD pulse simulation method converged to a relative error of less than 0 . 1 % after 250 time steps, and an error of less than 0 . 4 5 % relative to the Cottrell current was obtained after only 50 time steps.
1064
ANALYTICAL CHEMISTRY, VOL. 65, NO. 8, APRIL 15, 1993
The second comparison evaluated the accuracy of simulation when a staircasevoltammetric experiment for a chargetransfer system displaying reversible behavior was modeled. A numerical calculation of the staircasevoltammetry response done with high accuracy was used as a reference for these comparisons.12 The accuracy of the simulationswas examined as a function of the size of the potential steps used and the number of simulation time increments per step required to adequatelyreproducethe staircase response. With a staircase voltammogram generated using 10-mV potential steps, 11 simulation time increments were required to match the reference to less than 1% relative error and an error of less than 0.55% was obtained using 20 increments. As the size of the potential steps was decreased,the number of increments requiredto achieve a set error level also decreased. Reducing the step size to 2.5 mV yielded an error of less than 0.75% using 10 simulation time incrementa. In the limit that the potential-step size was reduced to 0.5 mV, only five time increments were required to reduce the simulation error to less than 1% relative to the reference response. Performance of the Filter. The quality of the parameter estimates from the Kalman filter was examined over a wide range of values of k, and a. A series of 500 staircase voltammogramswas generated, using different combinations of k, and a. Values of k, were varied over several orders of magnitude, while those of a ranged from 0.05 to 1.0. All other simulation parameters were held constant throughout the series. The resulting set of simulations covered the range of system kineticsfrom Nernstian to totally irreversiblebehavior. These simulations were characterized by the Matauda and Ayabe A parameter's (15) which is often used as a measure of the reaction kinetics. In eq 15, v is the scan rate. White noise was added to each of the simulated voltammograms to give a peak signal-to-noise ratio of 150. The filter was used to obtain estimates of the kinetic parameters for each of these data seta. In this series of fits,multiple passes through the data were used. The initial guesses supplied to the filter for the first pass were set to contain no a priori information; a value of 0.5 was input for a,and k, was set essentially to zero. State vector estimates obtained from the first pass through the filter were then used as initial guesses for the next pass. Convergence of the parameter estimates was usually obtained within two or three passes, depending on the noise level present on the voltammogram to be fitted. When the reaction kinetics displayed quasi-reversible or irreversible behavior, the estimated parameter values were in excellent agreement with those used to generate the simulations. Typical errors were less than 1%relative to the true value for system where A was less than 1. In fitting voltammograms with reversible kinetic behavior, where A took on values greater than 10, erratic results were obtained because the voltammogram shape became independent of the state parameters k, and a. Because the estimation of parameters requires that the states of the model be observable from the data, kinetic parameter estimates for reversible systems cannot be obtained. The two or three passes through the filter required for convergence of the parameter estimates were a function of the noise level present on the data sets, as noted above. When no noise was present, only a single pass through the filter was required. A single-passKalman filter fit required only a slight increase in computation time over that used by the simulation itself; the exact magnitude of the increase was observed to be (12) Kalapathy, U.; Tallman, D. E. Anal. Chem. 1992,64,2693-2700. (13) Matsuda, H.; Ayabe, Y. 2.Elektrochem. 1955,59,494-503.
8
.
i
* .
-10-
hdA)
Flgws 3. Percent relative error In ksestimates vs the log of Matsuds and Ayabe's A parameter.I3 All results were obtalned from a single pass through the Kalman filter, using randomtyselected Initial guesses for k, and a.
a function of the simulation used. For example, 76.2 s was required (i486based microcomputer! to generate a single EFD simulation of a typical staircase voltammogram, while the Kalman filter required only 78.9 s to optimally simulate the same voltammogram. Other simulation methods, not discussed in this paper, gave similar computational overhead times when coupled to the Kalman filter (e.g., a CrankNicolson simulation with expanding space grid2 required 41.3 s to simulate a typical staircase voltammogram while 47.9 s was required to optimize this simulation with the Kalman filter). Even with a substantial noise level and no a priori information about the values of k, and a,the two or three passes required by the Kalman filter optimization represent a substantial savings in computer time over other methods based on nonlinear regression. To evaluate the estimation accuracy and precision for the extended Kalman filter, Monte Carlo studies were employed. When a linear Kalman filter is used, the error covariance matrix (P) provides a direct means to evaluate estimation errors. In the extended Kalman filter, the error covariance matrix propagated is only an approximation to the true error covariance. As a result the accuracy of the filter is dependent on the trajectory of the state ~ector.1~ The only way to properly assess both the accuracy and precision of the parameter estimates is through replicate fittinge with different initial guesses for the state vector. This approach was employed in the second set of synthetic data studies. A series of voltammograms was generated, using a value of CY = 0.5 and varying k,over severalorders of magnitude. Band-limited white noise was added to each of the voltammograms to give a peak signal-to-noiae ratio of 1OOO. The voltammograms were fit using randomly selected values of k, and a as initial guesses. These initial guesses were centered about the values used to generate the simulations with a relative standard deviation of 25%. Parameter estimates were obtained using a single pass through the K h a n filter. Errors in the k, estimates obtained from this series of fits are displayed in Figure 3. In systems with quasi-reversible and irreversible kinetics, the estimates for k, and a obtained by the filter converged to the values used to generate the voltammograms. From these experiments it was concluded that, even in the presence of noise, accurate estimates can be obtained in a single pass when the initial guesses supplied to the filter are close to the true values. No real convergence from the initial (14) Gelb, A., Ed. Applied Optimal Estimation; MIT Press: Cambridge, MA, 1974.
ANALYTICAL CHEMISTRY, VOL. 05, NO. 8, APRIL 15, 1993
Table I. Relative Parameter Estimation Error and MAD Values for Different Voltammetric Methods method staircase large-amplitude staircase randomly perturbed staircase staircase large-amplitude staircase randomly perturbed staircase
ks
k,
a
a
170
0.77
3.2
0.090
0.69
1700 1700
0.30 0.12
0.80 0.71
0.074 0.012
0.22 0.19
1700
0.53
0.63
0.12
0.17
SINa error,* % MAD, % error,* % MAD, 76 170 2.2 4.8 0.41 1.4 170 0.78 3.6 0.092 1.0
*
Signal-to-noiseratio. Calculated as ((median(estimates)- true)/ true) X 100%.
distribution of parameter values was observed for reversible systems. The precision of the parameter estimates obtained from Kalman filter fitting of the optimized simulation model to data was also examined with a Monte Carlo method. A simulated staircase voltammogram for a charge-transfer system exhibiting irreversible kinetic behavior ( k , = 1X 10-5, a = 0.5; DO= D R= 5 X 10-6cm2s-l; scan rate, 1V s-l), consisting of 500 potential steps, was used as the input data set for the fitting routine. Random noise was added to the synthetic data set prior to each fit to give a peak signal-to-noise ratio of 1700. Randomly selected initial guesses for the rate constant and transfer coefficient were input to the filter for a series of 250 trials. The distribution of initial guesses was centered about the values of the kinetic parameters used to generate the simulation, with a median absolute deviation (MAD) of 25 7% relative to the true values. Fitting was done by using a single pass of the filter through the data to produce an optimized simulation model along with estimates of k , and a. Median values for the rate constants and transfer coefficients estimated by the filter were within 0.03% and 0.004 % , respectively, of the true values. Median absolute deviation values for the estimates of k , and a were 1.1% and 0.19% of the estimated values, respectively. Modeling Generalized Response Surfaces. The voltametric protocol independent performance of the Kalman filter was demonstrated through fitting voltammograms produced by three distinctly different potential-step sequences: conventionalstaircase, large-amplitudestaircase, and the randomly perturbed staircase. The following parameter values were used to generate the simulations: k , = 1.0 X 1W cm s-l; a = 0.5; n = 1; DO = DR = 5.0 X 10-6 cm2 5-1; sampling interval, 0.001 s. Each of the simulations consisted of a total of 500 data points, which corresponded to a scan rate of 2.0 V 5-1 for the conventional staircasemethod. The two large-amplitude pulse methods consisted of a total of 50 potential steps, with 10 data points each. Using the synthetic voltammograms, two series of Monte Carlo studies were performed to investigate the quality of parameter estimates. Different peak signal-to-noiseratios, 170and 1700, were used in the two series. At both signal-to-noise ratios, 50 fits were performed for voltammograms generated from each of the three experimental techniques, with random selection of the initial guesses for k , and a. The results are displayed in Table I. Very accurate parameter estimates were obtained for each of the three voltammetric methods in both series. In the lower signalto-noise series both the relative estimation errors and the relative spread of the estimates decreased as the random characterof the voltammogram increased. Thissame decrease in the relative spread of the estimates was observed in the higher signal-to-noise series. The decrease in the range of values for the parameter estimates is indicative of better
1065
convergence in the fits. This trend is due to an increase in the accessible qualitative informationpresent in the responses generated by random potential-step sequences.lSl8 The increase in the observablekinetic informationwas the primary motivation for applying such an unconventional potentialstep protocol. Even with the restricted amplitude of the random component used in this study, the voltammograms with more random character yielded better convergence in the parameter estimates. The increase in the observed error for the median k , and a estimates in the high signal-to-noise, randomly perturbed staircase experiment may be a function of the specific potential-step sequence used and should not pose a problem for the general implementation of these methods. Corrections for Double-Layer Effects. One important component present in the study of real electrochemical kinetics that was not previously incorporated into the Kalman filter methods was the effect of the electrochemical double layer on the observed heterogeneous charge-transferkinetics. In the present study, a Frumkin correction was applied. The observed rate constant is related to the corrected value by19 (16) where kobs and kcor, are the observed and corrected rate constants, a is the transfer coefficient,2 is the charge present on the oxidized species, 02 is the average potential present at the outer Helmholtz plane, and the remaining symbols have their usual electrochemical meanings. The value of a2 is obtained for a nonspecificallyadsorbed 1:lelectrolyte from Guoy-Chapman-Stern theory as20 0, =
(g)
sinh-’(
rn
’”)
(8ecJITc)
(17)
where qrnis the charge present per unit area of the electrode surface at a given potential, c is the electrolyte concentration, e is the dielectric constant for the solvent, E, is the permitivity of a vacuum, Z is the charge present on the electrolyte, and the remaining symbols have their usual meaning. In the present study this correction was built into the simulation model within the filtering routine. Hence, supplying the Kalman fiiter with the oxidized species charge and appropriate values for 0 2 as a function of applied potential resulted in double-layer-corrected parameter estimates of the heterogeneousrate constant. Thiscorrectionwas usedin the analysis of the experimental data discussed below. Performance of the Method with Experimental Data. The examination of experimental voltammetric data was required both to demonstrate the adequacy of the generalized step voltammetry simulationand to evaluate the performance of the Kalman filter optimization in the presence of possible small, but unmodeled, interferences and nonidealities. Two systems that have been previously characterized were used to study the performance of the method with real data. The first reaction studied experimentally was the reduction of Cr(H20)63+in 0.5 M NaC104-0.01 M HC104. This reaction has been used previously as a model system for heterogeneous charge-transferkinetics by many authors.21-24 Values of qrn,25 diffusion constants for the oxidized and reduced species (5.6 (15)Ichise, M.;Nagayanagi, Y.; Kojima,T.J.Electroanal. Chem. 1974, 49,187-198. (16)Yamagishi, H.J.Electroanal. Chem. 1987,235, 117-129. (17)Yamagishi, H.J.Electroanal. Chem. 1990,283, 115-123. (18)Yamagishi, H.J. Electroanal. Chem. 1991,297, 125-131. (19)Delahay, P. Double Layer and Electrode Kinetics; Wiley: New York, 1965. (20)Mohilner, D. M.Electroanal. Chem. 1966, 1, 241-409. (21)Parsons, R.;Passeron, E. J. Electroanal. Chem. 1966, 12, 524529.
ANALYTICAL CHEMISTRY, VOL. 65, NO. 8, APRIL 15, 1993
1066
h
-I
-0.7
-1.1 I
0
0.5
1
1.5
2 TIME (SI
2.5
3
3.5
1
4
Table 11. Kinetic Parameter Estimates for Cr(HzO)es+ Reduction in 0.5 M NaC104-0.01 M HClO4 method standard staircase large-amplitude staircase randomly perturbed staircase square wave mean of all methods polarography median of several methods
k,, cm 8-1
ref
a
(1.91 f 0.23) X 10-7 0.504 f 0.012" c (2.02 f 0.06) X le7a 0.493 f 0.003" c (1.72 f 0.15) X le7 0.496 f 0.01@ c (I
(2.01 f 0.17) X le7 (1.9 f 0.1) X le7 1.29 x le7 5.2 X
0.481 f 0.011" 0.49 f O.0lo 0.47u 0.61b
e c 23 6,23, 24
a Values corrected for double-layer effects. Values not corrected for double-layer effects. This work.
-5
1
0
0.5
1
1.5
2
2.5
3
3.5
4
TIME (5)
Flgurr 4. Experlmental data for the reduction of a 1.8 mM solution of Cr(H20)a3+In 0.5 M NaC104-0.0 1 M HCIO,, using a randomly perturb& stalrcase technlque. The potentiel-step sequence (vs SCE) Is shown In the top plot. A sampling Interval of 0.001 s was used. The lower plot dlsplays the voltammetric response (solid line),the Kalman flltergenerated fit (dashed line),and the Innovations for the flt (dotted line). Parameter estimates for the fit are glven In Table 11.
10-6 and 7.9 X 10-6 cm2 s-l, respectively29, and the value of E"' (-0.650V vs SCE6) were obtained from the literature. Data sets were recorded for each of these systems using each of the four potential-step sequences discussed above. For each sequence, at least two different step durations were used, and for each, three replicate runs were obtained as noted above. Based on the above values and the experimental parameters pertaining to each of the voltammograms, estimates of the standard rate constant and transfer coefficient for the C r ( H ~ 0 ) 6 reduction ~+ were made using a single pass through the fiiter. Filtering of replicate data sets, using initial guesses based on previously published values, yielded a narrow range of estimates. This observation indicated that an adequate level of convergence was achieved within a single pass. In runs when the initial guesses were substantially inaccurate (i-e.,in excess of 1order of magnitude for the case of the rate constant), multiple passes through the filter were required to achieve acceptable convergence. For these situations, the parameter estimates obtained from the previous pass were supplied to the filter as initial guesses for the present pass. An example fit obtained using the randomly perturbed staircase experiment is presented in Figure 4. A summary of results for fitting the experimental voltammograms produced with each of the four potential-step sequences is given in Table 11. All of the results presented in the table were obtained using a single pass through the filter. Values for the double-layer-corrected charge-transfer parameters reported here are in good agreement with those previously published, also listed in Table 11. The reduction of Zn(I1) at the HMDE in 1.0 M KN032630 was the second reaction used to test the performance of the X
(22) Weaver, M. J.; Liu, H. Y.; Kim, Y. Can. J.Chem. 1981,59,19441953. (23) Anson, F. C.; Rathjen, N.; Frisbee, R. D. J. Electrochem. SOC. 1970,177,477-482. (24) Zielinska-Ignaciuk, M.; Galus, Z. J.Electroanal. Chem. 1974,50, 41-53. (25) Wroblowa, H.; Kovac, Z.; Bockris, J. OM. Trans. Faraday SOC. 1965,61,1523-1548. (26) Go, W. S.;ODea, J. J.;Osteryoung, J. J.ElectroanaL Chem. 1988, 255, 21-44. (27) Galus, Z. CRC Crit. Reu. Anal. Chem. 1975, 4, 359. (28) ODea, J.;Osteryoung, J.; Osteryoung, R. A. J.Phys. Chem. 1983, 87,3911-3918. (29) ODea, J.; Osteryoung, J.;Lane, T. J.Phys. Chem. 1986,90,27612764. (30) Go, W. S.; Osteryoung, J. J. ElectroanaL Chem. 1987,233, 275281.
Table 111. Kinetic Parameter Estimates for Zn(I1) Reduction in 1.0 M KNO, voltammetric method k,, cm s-1 a standard staircase large-amplitude staircase randomly perturbed staircase square wave mean of all methods square wave square wave
ref
(3.19 f 0.03) X 10-4 a 0.1678 f 0.0082" c (2.72 f 0.13) X 10-4 0.1360 f 0.0057u c (2.93 f 0.18) X lo4 (3.12 f 0.21) X 10-4 (3.0 f 0.2) X lo4 a (2.64 f 0.16) X lo4 (4.6 f 0.3) X
*
a
(I
0.1630 f 0.0071* c 0.1846 f 0.0079" c 0.16 f 0.02" C (0.20 f 0.02)(I 24 0.23 f O.Olb 28
Values corrected for double-layer effects. Values not corrected for double-layer effects. This work.
Kalman filter method. Experimental voltammograms were collected using each of the four potential-step sequences. For each voltammetric technique, three replicate runs were recorded, using two different step durations. Values of a2as a function of applied potential for the 1.0M KN03 solution,26 diffusion constants (Dz,,(II) = 6.7 X 1o-S cm2s-1 26 and Dz,.,(H~) = 1.9 X cm2 s-* 27), and an estimated value for the reversible half-wave potential (-1.OOO f 0.001 V vs S C E 9 were obtained from the literature. Table I11 lists the results obtained from fitting the experimental voltammograms for the Zn(I1) reduction using a single pass through the fiitm. The values obtained in this study were also in excellent agreement with the previously published results.26 An example fit for the Zn(I1) system is shown in Figure 5,which displays the results of fitting a large-amplitude staircase voltammogram. The fits obtained for the Zn(I1) data were very good but are of slightly lower quality than those obtained in the Cr(&0)e3+ studies. Estimating ED' for Quasi-Reversible Kinetics. Fitting data for systems exhibiting quasi-reversiblebehavior requires a very good estimate of E O ' , due to the sensitivity of the shape of the response curves to estimates of E"' in this kinetic regime. If the value of Eo' input to the filter is not in good agreement with the true Eo' for the data, the simulation model optimized by the filter will not match the shape of the experimental voltammogram and the fit (and parameter estimates) will be poor. In the synthetic studies described above, the values of Eo' used to generate the data sets were known exactly, resulting in filtering runs that yielded both excellent parameter estimates and fit quality. It was necessary to confirm that the literature value of the reversible half-wave potential for Zn(I1) in 1.0 M KN03,26when used in the filter optimization, yielded a simulation model that best fit the experimental voltammetric data presented here. A series of fits was performed using one of the randomly perturbed staircase voltammetric data sets for the Zn(1I) reduction described above. Because the version of the filter discussed here does not estimate the value of EO', any
ANALYTICAL CHEMISTRY, VOL. 85, NO. 8, APRIL 15, lQg3 -0.8 I
-1.2j 0
1
0.1
0.2
0.4
0.3
0.5
0.6
0.7
0.8
0.5
0.6
0.7
0.8
TIME (s) 15'
I
-5 I 0
0.1
0.2
0.3
0.4
i
TIME (s)
Figure 5. Experimentaldata for the reduction of a 0.62 mM solution of Zn(I1) in 1.0 M KN03, using a largeamplitude staircase technique. The potentiaCstep sequence (vs SCE) Is shown In the top plot. A sampling interval of 0.0005 5 was used. The lower plot displays the voltammetricresponse(solid line),the Kalmanfllter~nerated fit (dashed line), and the innovatlons for the fit (dotted line). Parameter estimates for the fit are listed in Table 111.
simplex optimization was also constrained, so that the search was confined to the region -1.030 < E"' < -1.OOO V vs SCE, where it was reasonable to expect the true value of the formal potential. Only a few iterations were required by the simplex to estimate the single parameter. Several runs of the optimization were performed,each with different initial values for E"'. The minimum was found at an E"' value of -1.0107 V vs. SCE. In summary, a Kalman filter-optimized simulation and fitting method for estimating kinetic parameters from step voltammetry data, independent of any specificpotential-step sequence, was demonstrated. The simulation used as the filter measurement model was an EFD aigorithm modified to model the voltammetricresponse of a heterogeneouschargetransfer reaction in a potential-step experiment. Corrections for double-layer effects were built into the simulation model for studies involvingexperimental data. Becausethe Kalman filter optimization is computationally efficient, a single pass through the data with the filter uses only slightly more computer time than a single, stand-alone EFD simulation, but also yields estimates of parameters used in the simulation model and estimates of the error covariance of these parameters. Using this K h a n fiiter-basedfitting method, accurate charge-transferparameter estimates were obtained for a wide range of reaction kinetics, with little or no a priori information. While the EFD simulationused here was adequate to model a simple heterogeneous charge-transfer reaction, more sophisticated algorithms should be employed for accurate modeling of systems with complex kinetic schemes. Additionally, the use of an expanding space grid will increase the computational efficiency of the simulation. These extensions to the present method are currently being investigated in this laboratory and will be the topic of future papers.
i: 0'
-1.02
-1.018
-1.016
-1.014
-1.012
-1.01
-1.008
-1.006
-1.004
1
-1.002
1067
ACKNOWLEDGMENT
This work was supported by Grant DE-FG02-86ER13542 from the Division of Chemical Sciences,Office of Basic Energy Sciences, US. Department of Energy.
Value of EO'Input to Filter
Figure 6. Sum of squared innovationsVI the value of P' input to the Kalmanfilter for the Zn(1I ) reduction In 1.O M KNO,, using the randomly perturbedstalrcase experiment. All other parameter values were held constant for the series of fits.
adjustments to this parameter must be made external to the fitting process. The value of E"' used in the filter model was varied over a severalmillivolt range about the reference value, while all other input parameters were held constant. Parameter estimates and the filter innovations for each fitlo were recorded. Figure 6 displays a plot of the sum of squared innovations obtained as a function of the value of E"' used. A minimum was observed for an E"' value of -1.011 V vs SCE, which was in good agreement with that calculated from the literature value for the reversible half-wave potential26 and the known diffusion constanta.26-27 While the above series of fits represents a large amount of computer time, the estimation of the E"' value can be made more computationally efficient. The extended Kalman filter optimization of the simulation model described above was run inside of an external simplex optimization to estimate the best value for E"' by minimizing the sum of squared innovations for the fit. In this way, an optimal formal potential was determined, and the optimal simulation model for that formal potential was generated and fit to the experimental voltammogram. This dual optimization was tested on the Zn(I1) voltammetric data discussed above. As before, all ather input parameters were held constant. The
APPENDIX The system model for the extended Kalman filter is given by
i ( k ) = f(x(k),k) + w(k)
(Al) where f(x(k),k) is a function describing the nonlinear system dynamics, w is the model noise vector, and the index k is that of the present time step. The model noise is assumed to be normally distributed with zero mean and covariance Q ( k ) .
-
w(k) N(O,Q(k)) The nonlinear measurement model is
(A2)
z ( k ) = h(x(k),k)+ u(k) (A3) where z ( k ) is an observation from the observation vector z, h(x(k),k) is a nonlinear function of the state vector, and u(k) is the measurement noise. A normal distribution with zero mean and covariance R ( k ) is assumed for the measurement noise.
-
u(k) N ( O 3 ( k ) ) (A4) An additional assumption is made with respect to the model and measurement noise. E[wG)u(k)]= 0 for all j and k (A51 The initial guess for the state vector is specified as normally
1068
ANALYTICAL CHEMISTRY, VOL. 65. NO. 8. APRIL 15, 1993
*
distributed about with a covariance given by the initial error covariance matrix (Po) x(0) = N*,P,)
(A61
P(0) = E[(x(O)- X,)(X(O) - XdT1
(A7)
From the nonlinear system model above, a linearized system dynamics matrix is obtained from a series expansion about the present trajectory of the state vector.14
The measurement function is treated in an analogousmanner as
P(klk-1) =
F(x(k-l),k-l)P(k-llk-l)FT(X(k-l),k-l) + Q(k-1) ( A l l ) where Q is the covariance matrix for the model noise. From the current estimate of the state vector, the value of the nonlinear measurement model is calculated. With the linearized measurement model, the Kalman gain matrix is calculated. K(k) = P(klk-l)HT(x(klk-l))[H(x(klk-l)) X P(k(k-l)HT(x(klk-l)) + R(k)]-' (A12) whereR(k) is the measurement noise covariancematrix. Based on the discrepancy between h(x(klk-l),k) and the present measurement t(k), a correction, weighted by the Kalman gain matrix, is applied to x(k(k-1) yieldingthe updated state vector x(klk). x(klk) = x(k)k-1)
The vector x(klk-1) is the present estimate of x based on all information up to but not including the current time step. There are five steps in the extended K h a n algorithm. At the beginning of each iteration, an estimate of the state vector x(klk-1) is projected, baaed on the previous value x(k-llk-1). x(k)k-1) = f(x(k-l),k-1)
+ W(k-1)
The error covariance matrix is propagated as
W0)
+ K(k)[z(k) - h(x(k)k-l),k)]
x(klk-1) + K(k)v(k) (A13) where v(k) is the filter innovations. The error covariance matrix is then updated P(k(k) = [I-K(k)H(x(klk-l))]P(klk-l) (A14) where I is the identity matrix. The sequence of eqs AlO-Al4 is repeated for each point of the data set being filtered.
RECEIVEDfor review September 8, 1992. Accepted January 11, 1993.
