ARTICLE pubs.acs.org/IECR
Numerical Analysis of Alternating Current Voltammetry: Identifiability, Parameter Selection and Experimental Design A. Vikhansky,*,† S. M. Matthews,‡ and A. C. Fisher‡ † ‡
School of Engineering and Materials Science, Queen Mary, University of London, London E1 4NS, U.K. Department of Chemical Engineering and Biotechnology, University of Cambridge, Cambridge CB2 3RA, U.K. ABSTRACT: This paper describes a numerical analysis of the parametric identifiability of electrochemical systems. First, we analyze global identifiability of the entire set of parameters in a single alternating current voltammetry experiment and examine the effect of different waveforms (square, sawtooth) on the accuracy of the identification procedure. The analysis of global identifiability is equivalent to finding a global optimum of a specially designed function. The optimization problem is solved by a random search method, and a statistical analysis of the obtained solution allows for selection of a subset of the parameters (or their linear combinations), which can be identified. Finally, we discuss optimization of the waveform for better identifiability.
kf
A þ ne- h B kb
∂t CA - D 1=2 ∂xx CA ¼ 0,
∂t CB - D -1=2 ∂xx CB ¼ 0
ð2Þ
where D is the ratio of the diffusion coefficient of A to the diffusion coefficient of B. The forward and backward rate constants are given by kf ¼ k0 expðR½EðtÞ - E0 - Ru IÞ kb ¼ k0 expð- ð1 - RÞ½EðtÞ - E0 - Ru IÞ
ð3Þ
where E(t) is the applied potential, E0 is the reversible potential, k0 is the rate constant (i.e., Damk€ohler number), Ru is the uncompensated resistance, and I = Ic þ If, the total current consisting of capacitive current Ic and faradaic current If. The capacitive current is given by dE ð4Þ Ic ¼ Cdl dt where Cdl is the capacitance of the double layer. The system has the following boundary conditions: x ¼ 0: If ¼ D1=2 ∂x CA ¼ -D-1=2 ∂x CB ¼ -kf CA þ kb CB x ¼ ¥: CA ¼ C¥ ; CB ¼ 0 ð5Þ In a typical voltammetric experiment, eqs 2-5 are solved with different sets of parameters, namely, E0, R, k0, Ru, Cdl, and D, until the predicted total current coincides with the experimental voltammogram. In other words the task is to find a set of the parameters providing solution to the following problem: Ipredicted - Imeasured
)
min
ðE0 ,R,k0 ,Ru ,Cdl ,D Þ
ð6Þ
ð1Þ
Here A and B are the solvent and the solute. For the sake of simplicity hereafter all equations are given in dimensionless form (for the dimensional form of the equations and their nondimensionalization, r 2011 American Chemical Society
see, e.g., refs 1 and 2). The transport equations for the concentrations CA and CB read
)
1. INTRODUCTION In a typical voltammetric measurement, the current response to an applied potential is analyzed to provide quantitative information on parameters such as the following: analyte concentration, diffusion coefficients, charge-transfer rate, and the structure of the electrical double layer, etc. Aside from fundamental information about the nature of interfacial electron-transfer processes, voltammetry can also quantify the kinetics of chemical reactions occurring in solution. Early measurements focused on the use of relatively simple voltammetric waveforms such as the potential step, linear sweep, and cyclic voltammetry. However, to gain insight on more complex processes, a number of more convoluted waveforms have been adopted. These include the imposition of a higher frequency perturbation with a variety of forms such as a sine wave, square wave, and sawtooth. Recent developments in voltammetric methodologies have enabled the numerical and experimental studies of some of these waveforms; these have shown that the imposition of a variety of periodic perturbations to a direct current (dc) sweep are all part of a single category of techniques. This is due to the fact that each waveform can be described as a combination of a series of sine waves (equation). As a result, there is a series of similar response characteristics which enables the development of a unified analytical methodology to be applied to all waveforms. This paper outlines the development of a numerical method, which allows for the analysis of identifiability of the parameters of the system and provides a quantitative measure for the efficiency of different waveforms. The system under consideration is semiinfinite mass transport to a planar electrode, where a reduction or oxidation reaction occurs:
Received: May 7, 2010 Accepted: December 8, 2010 Revised: November 24, 2010 Published: January 28, 2011 2831
dx.doi.org/10.1021/ie101048b | Ind. Eng. Chem. Res. 2011, 50, 2831–2838
Industrial & Engineering Chemistry Research where the difference between the measured and the predicted currents is taken in any norm. Although this approach is intuitively obvious, its implementation does not necessarily bring the desired results. First, the solution of (6) might be not unique; second, a small component of experimental noise can cause a large deviation of the estimated parameters. Therefore, an identifiability test has to be performed before starting the experiments in order to ensure that a given experimental design allows for the identification of the parameters with sufficient accuracy. If the model is unidentifiable under given experimental conditions, the number of the estimated parameters should be reduced or, alternatively, the experimental design should be changed in order to guarantee the identifiability. The uniqueness of the inverse problem for partial differential equations has been the subject of intensive research during the past 7 decades (see, e.g., refs 3-5). Unfortunately, the vast majority of the results are dedicated to the systems described by a single parabolic equation and cannot be directly applied to electrode dynamics. Although there are a number of other analytical methods based on Taylor series expansion,6,7 local state isomorphism, and differential algebra,8 which can be used for identifiability analysis of the models described by ordinary differential and differential-algebraic equations, the applicability of these methods to partial differential equations is questionable. As an alternative to the above-mentioned methods, the authors of refs 9 and 10 used a direct approach; namely, they explicitly built the M-dimensional output surface (for details see eq 7 in section 2) in the N-dimensional space of the measured parameters. The system is identifiable if the output surface does not intersect itself. Certainly, this approach works only if M is not too large (in ref 9, e.g., M = 2 and N = 3). A statistical approach to the problem of parameters estimation has been adopted in refs 11-14, where the parameters of the model are treated as random. A bimodal distribution of the identified parameters or their broad distribution is signaling that the system is poorly identifiable. A heuristic approach based on recognition of patterns of behavior has been proposed in refs 1 and 2. It is shown that different patterns are sensitive to variation of a single parameter and the method is based on a consecutive estimation of the parameters. The uniqueness of the result is verified by comparison of the predicted and measured higher harmonics. Mathematically, this approach is equivalent to solution of (6) by a point-relaxation method. Although this method converges for sufficiently smooth functions, due to numerical imperfectness it can stop far from the optimum.15 Therefore, in spite of the significant progress in this field, there is a need for further numerical investigation based on more rigorous principles. In the present study we exploit an identifiability test,16 which requires the solution of an unconstrained optimization problem. The paper is organized as follows. In section 2 we describe the numerical methods of identifiability, formulate our identifiability test function, and discuss numerical implementation of the method. Section 3 presents the results of the numerical experimentation. We analyze the identifiability of the parameters in a single alternating current (ac) voltammetry experiment and discuss the methods for parameter selection and experimental design. Finally, we consider identifiability of the entire set of parameters in a series of experiments with varying scan rate. The paper closes with a brief summary.
ARTICLE
2. IDENTIFIABILITY TEST AND NUMERICAL IMPLEMENTATION An identifiability test aims to answer the question of whether the parameters of a model can be uniquely determined from experimental data.17 We consider a set of outputs yi (i = 1, ..., M) measured in an experiment (e.g., time series of the current). It is assumed that yi is a function of an N-dimensional set of parameters θj ∈ Θ (M G N) and an R-dimensional set of operating conditions uk (e.g., applied potential), which can be controlled by the experimentalist: y ¼ Yðθ;uÞ ð7Þ Therefore, for a given u the outputs form an N-dimensional output surface9,10 with parametrization θj. The mathematical model Y( 3 ; 3 ) is identifiable under the operating conditions u if for every y there exists only one θ, which satisfies eq 7. In the other words ð8Þ Yðθ1 ;uÞ ¼ Yðθ2 ;uÞ S θ1 ¼ θ2 Hereafter if we assume that the operating conditions u are fixed, we remove them from the formulas if it is unambiguous. To proceed further, we renormalize the parameter set as follows. We assume that the vector of parameters θ̂ = (E0,R,k0, E θ̂i E θ̂max . Then Ru,Cdl,D) belongs to a hypercube, i.e., θ̂min i i we introduce a new parameter’s set 0 E θi E 1 such that ^ min ^ max þ ð1 - θi Þθ ^ ¼ θi θ ð9Þ θ i
i
i
In the present work we slightly modify the identifiability test16 by defining the following optimization problem: ( ) Fðθ1 - θ2 Þ 2 η ¼ max Hðθ1 ,θ2 Þ ð10Þ θ1 , θ2 jYðθ1 Þ - Yðθ2 Þj2 where at the beginning we assume for simplicity that F(θ1-θ2) = |θ1 - θ2|2, although later we will discuss another form of the function F. Equation 10 implies that jθ1 - θ2 j E η j Yðθ1 Þ - Yðθ2 Þj
ð11Þ
and if we associate ΔY = |Y(θ1) - Y(θ2)| with an experimental error, the identifiability criteria η provides an upper bound for the error of the identification procedure Δθ = |θ1 - θ2|; i.e., in the worst case the relative error of the parameters identification is η times higher than the relative error of the measurements. Since a significant amount of research literature is devoted to sensitivity analysis,9,10,18-20 it is instructive to demonstrate how the identifiability test (10) relates to the linearization-based sensitivity analysis. If the maximum of (10) is attained at θ*1 ≈ θ*, 2 Taylor series expansion of numerator and denominator of (10) yields ( T ) ∂Y ∂Y 2 2 ðθ1 - θ2 Þ ðθ1 - θ2 Þ jθ1 - θ2 j = η max θ1 ,θ2 ∈Θ ∂θ ∂θ ¼
1 λmin ð12Þ
where λmin is the smallest eigenvalue of the Gramian matrix T ∂Y ∂Y G¼ ð13Þ ∂θ ∂θ 9,10
The solution of the identifiability problem crucially depends on our ability to find the global optimum of (10). The computation 2832
dx.doi.org/10.1021/ie101048b |Ind. Eng. Chem. Res. 2011, 50, 2831–2838
Industrial & Engineering Chemistry Research
ARTICLE
Figure 1. Typical waveforms used in alternating current (ac) voltammetry. The insertions show the ac components of the signal: (a) eq 18; (b) square wave; (c) forward sawtooth; (d) backward sawtooth.
procedure is as follows. Equations 2-5 are solved for two sets of parameters x = (θ1,θ2) by a finite difference backward differentiation method on an expanding grid and then the identifiability criteria H(x) H(θ1,θ2) is evaluated. In the present work we use the simultaneous perturbation gradient approximation (SPSA) method21 in order to solve the optimization problem. It is a random search method, which uses a finite-difference estimate of the gradient of the objective function based only on two measurements of the objective function as ∂H Hðx þ cΔÞ - Hðx - cΔÞ g ∂xi 2cΔi
ð14Þ
where the random value Δi = (1 with probability 1/2 for each outcome. Then the new pair of parameter sets are generated as x ¼ x þ ag R1
ð15Þ γ
Parameters a = a0/(k þ A) and c = c0/(k þ A) are decreasing functions of the number of the iteration k, and the common values for R1 and γ are 0.602 and 0.101, respectively, and A, a0, and c0 are user-specified constants.21 The intrinsic noise of the method prevents it from being stuck at a local optimum, and it is shown that under quite general conditions the method can be used for global optimization.22 To ensure convergence to the global maximum, we have used a search with multiple starting points. Each run of the SPSA algorithm consists of 2 103 steps, and we performed between 20 and 40 independent runs. Since each iteration of the algorithm requires four solutions of the direct problem, eqs 2-5 have been solved (1.6-3.2) 105 times. The entire procedure takes about 2 h of a Pentium 4 processor running at 2.4 GHz.
3. RESULTS AND DISCUSSION 3.1. Identifiability of a Low-Resistance System. In the
present work we consider the triangular waveform used in direct
current (dc) cyclic voltammetry: 8 t > > < 2A0 T , Edc ðtÞ ¼ t > > : 2A0 1 - T ,
t 1 < T 2 t 1 G T 2
ð16Þ
where A0 and T are the amplitude and the period, respectively. A sum of sine waves of higher frequency is superimposed onto it: N X 20πnt ð17Þ An sin EðtÞ ¼ Edc ðtÞ þ T n¼1 as it is shown in Figure 1. We start with an “easy’’ case, namely, when the effect of the resistance and, therefore, double layer capacitance can be neglected. Then we can assume that Ru = 0, Cdl = 0, and the identifiability has to be checked only with respect to four parameters. In the present study we experiment with the ac waveforms which have been proposed recently. A waveform consisting of sinusoidal waves with amplitude inversely proportional to the frequency has been used in ref 2: An ¼ n-1=2
ð18Þ
Early, square, forward, and backward sawtooth waveforms have been examined in refs 23 and 24. These signals can be expanded as Fourier series involving the sine functions with amplitudes: 4 An ¼ , n ¼ 1,3,::: πn ð19Þ 2ð-1Þn An ¼ , n ¼ 1,2,3,::: πn where the first formula corresponds to the square waveform and “-’’ and “þ’’ are for forward and backward sawtooth waveforms, respectively. In our calculations we set T = 1 and A0 = 14. 2833
dx.doi.org/10.1021/ie101048b |Ind. Eng. Chem. Res. 2011, 50, 2831–2838
Industrial & Engineering Chemistry Research
ARTICLE
Figure 2. Total current for θ̂1 (see Table 1) (left) and the difference between the outputs produced by θ̂1 and θ̂2 (right) in a low-resistance system with waveform (18). C¥ = 1.
Table 1. Solution of the Optimization Problem (10) under the Conditions of Figure 2 ̂ θ 1
̂ θ 2
θ̂min
θ̂max
E0
6.49
6.62
6.0
10.0
R
0.681
0.587
0.3
0.7
k0
8.09
10.0
1.0
10.0
D
3.22
2.88
0.25
4.0
The results of the calculations with waveform (18) and the parameter’s sets are presented in Figure 2 and Table 1. The results for the second set of parameters θ2 are not shown because they are undistinguishable from those for θ1. Note that the presented pair θ1,θ2 is not unique; i.e, there are many other combinations of the parameters, which produce similar values of H. It means that the identifiability test does not look for one “bad’’ combination of the parameters but provides a reasonable estimate for the error under the conditions, which are likely to be met in a real experiment. Since the applied potential consists of a slowly varying component and a high-frequency sinusoidal part, it is natural to compare the high- and low-frequency components of the measured current separately. A Fourier transform has been used in refs 1 and 2 in order to separate different patterns of the signal behavior. In the present work we apply a filter in order to separate the effect of high harmonics; i.e., we use a uniform slidingwindow filter as follows: Z _ _ 80 tþT=160 IðτÞ dτ, I 0 ðtÞ ¼ IðtÞ - I ðtÞ I ðtÞ ¼ ð20Þ T t-T=160 and the expression for |Y(θ1) - Y(θ2)|2 reads 2R T 3 _ RT _ ½I 0 ðθ1 Þ - I 0 ðθ1 Þ2 dt ½ I ðθ1 Þ - I ðθ1 Þ2 dt 1 5 _ max 4R0T , R0 _ 4 ½I 0 ðθ1 Þ þ I 0 ðθ1 Þ2 dt T ½ I ðθ1 Þ þ I ðθ1 Þ2 dt 0
0
rate (1/T)) values of k0, and variation of the scan rate might improve the identifiability further. To assess the effect of the applied potential on the identifiability, we performed calculations with square, forward, and backward sawtooth waveforms. The results of the calculations are presented in Figure 3. The computed values of the identifiability criteria for these three waveforms are η = 7.50, 7.66, and 6.73. The computed currents are presented in Figure 3, and the obtained results demonstrated that although the backward sawtooth waveform corresponds to the smallest value of η, all four examined signals provide very close levels of accuracy. Note that the identifiability criteria characterize the amplification of the experimental error under the worst case scenario and the actual error of the identification procedure might be lower. 3.2. Optimization of the Applied Potential. It has been noted earlier that the output is a function of the set of parameters θ and the set of operating conditions u (the applied potential). Therefore, it is natural to apply an optimization procedure in order to find an applied potential, which guarantees better identifiability. Taking (10) into account one can formulate the following minimax problem: min max Hðθ1 ,θ2 ;uÞ u θ1 ,θ2
i.e., our task is to minimize the worst-case value of H. The solution of (22) is a saddle point of H(θ1,θ2;u). Due to the structure of the SPSA algorithm it is equally suitable for finding of maxima (minima) and saddle points. However, in the present work we are interested only in that saddle points, which correspond to the global maxima of H(θ1,θ2;u) for the given u. In order to exclude the suboptimal saddle points from consideration we modify the SPSA algorithm as follows. We generate a number (20-40) of pairs of the parameter sets xi = (θ1i,θ2i). Then we simultaneously perturb xi and the vector of operating conditions u = (A1,A2,...,AN). The component of the gradient with respect to xi is calculated according to eq 14, while the gradient in the operating condition’s space is estimated as
ð21Þ The computed value of the identifiability criteria is η = 8.68; i.e., in the worst case the relative error of the parameters identification is 8.68 times higher than the relative error of the measurements. Note that the error in (21) is an integral of the measured signal, and one can expect that the experimental noise is partially smoothed out. Therefore, the obtained result indicates that when the effect of resistance is negligible, the parameters of the model can be estimated with a reasonable accuracy in a single experiment. The less identifiable set of parameters correspond to high (in comparison with the scan
ð22Þ
gu ¼
Hðx i þ cΔi ;u þ cu Δu Þ-Hðxi - cΔi ;u - cu Δu Þ ð23Þ 2cu Δu
where i* = arg max(H(xi;u)) and cu, Δu have the same meaning as in (14). Finally, u is updated as follows u ¼ u-au gu
ð24Þ
The results of the calculations are presented in Figures 4 and 5 and Table 2. The computed value of η is 5.72; i.e., there is a modest improvement of around 30% in comparison to signal (18) and 15% in comparison to the backward sawtooth waveform. It is 2834
dx.doi.org/10.1021/ie101048b |Ind. Eng. Chem. Res. 2011, 50, 2831–2838
Industrial & Engineering Chemistry Research
ARTICLE
Figure 3. Total current for θ̂1 (solid line) and θ̂2 (dashed line) (left) and the difference between the outputs produced by θ̂1 and θ̂2 (right) for (top to bottom) square, forward, and backward sawtooth waveforms. C¥ = 1.
Figure 4. Amplitudes (left) and the total applied voltage (right) of the optimized experimental design.
worth noting that the optimal applied potential presented in Figure 4 is not unique: there are many other signals of similar intensity, which provide similar values of η. Our experience with different waveforms and the gradients ∂H/∂u estimated by the SPSA method indicate that the identifiability is not very sensitive to the particular form of the applied potential; i.e., every voltage with a frequency comparable to the characteristic reaction rate k0 yields nearly the same level of identifiability.
3.3. Parameters Selection. If the number of parameters is too large, the model becomes unidentifiable in a single experiment. The example is a system where the effect of resistance and the double-layer capacitance cannot be neglected. The value of identifiability criteria obtained in our calculations is η = 50, and the system is unidentifiable. An analysis of the obtained results indicates that effects of different parameters are mutually compensated, e.g., an 2835
dx.doi.org/10.1021/ie101048b |Ind. Eng. Chem. Res. 2011, 50, 2831–2838
Industrial & Engineering Chemistry Research
ARTICLE
Figure 5. Total current for θ̂1 (solid line) and θ̂2 (dashed line) (see Table 2) (left) and the difference between the outputs produced by θ̂1 and θ̂ 2 (right) for the optimized potential presented in Figure 4. C¥ = 1.
Table 2. Solution of the Optimization Problem (22) under the Conditions of Figure 5 θ̂1
θ̂2
θ̂min
θ̂max
E0
7.41
7.43
6.0
10.0
R
0.691
0.345
0.3
0.7
k0
10.0
8.20
1.0
10.0
D
1.07
1.82
0.25
4.0
Table 3. Columns of Matrix S S1
S2
S3
S4
S5
S6
θE
-0.0502
-0.0396
0.0099
0.9555
-0.2874
0.0183
θR
-0.6827
-0.7254
-0.0547
-0.0679
-0.0081
0.0040
θk0 θRu
-0.6197 -0.0996
0.5569 0.0621
0.4902 0.1759
-0.0860 0.1671
-0.2394 0.6177
-0.0276 0.7388
θCdl
0.0920
-0.0428
-0.1670
-0.2128
-0.6808
0.6731
θD
0.3590
-0.3955
0.8353
-0.0430
-0.1228
-0.0048
i.e., θ~ is the new reduced set of parameters. If we redefine the function F as ~1 - θ ~ 2 j2 ¼ Δ θ ~ T Δθ ~ ¼ ΔθT ðPT PÞΔθ Fðθ1 - θ2 Þ ¼ jθ ¼ tr½PðΔθΔθT ÞPT
ð26Þ
the projection matrix P can be determined from the solution of the following minimax problem: min max Hðθ1 ,θ2 ;PÞ P θ1 , θ2
Figure 6. Values of Δθk = θ1k - θ2k and ΔθR = θR1 u - θR2u.
increase of k0 can be canceled by a corresponding increase of the resistance R u. It has be noted earlier that the poorly identifiable set of parameters is not unique. Iterations of the SPSA method generate a significant number of pairs (θ1, θ2) with nearly maximum values of H as it is shown in Figure 6. This scatter plot does not provide full information, because different points correspond to different values of H; however, one can see that there is a correlation between two parameters. Therefore, although the entire set of the parameters is unidentifiable, one can expect that there is a (small) number of combinations of the parameters, which can be extracted from the experimental data with reasonable accuracy. In the present work we consider linear combinations of the ~ N (N ~ < N) parameters; namely, we postulate that there is an N matrix with orthonormal rows P such that ~ ¼ Pθ, θ
~ θ ¼ PT θ,
PPT ¼ I
ð25Þ
ð27Þ
i.e., the parameter selection can be considered as a part of experimental design.25 The SPSA algorithm described above can be used in order to solve (27); however, in the present work we use another approach, which does not require much of additional computations and can be implemented as a postprocessing of the identifiability analysis. The method is as follows. Consider the sequence of parameters (θ1i,θ2i) generated by the SPSA algorithm and keep only n points, which correspond to the 10% of the points with the highest values of H (the 10% threshold is chosen arbitrarily; use of the “top 5%’’ or “top 3%’’ does not change the results significantly). Then we define the following correlation matrix: C¼
1 X Δθi ΔθTi n i jΔY j2
ð28Þ
P Since tr[C] = (1/n) Hi and taking (26) into account, one can conclude that 1X Hðθ1i ,θ2i ;PÞ ð29Þ P ¼ arg minftr½PCPT g ¼ arg min n i i.e., the P solves the problem (27) on average. To find this matrix, we consider the matrix S = {s1,s2,...,sN} consisting of the normalized eigenvectors of the C, where i > i* implies that the corresponding eigenvalues λi E λi*. An example of this matrix is presented in Table 3. The first N0 columns of this matrix 2836
dx.doi.org/10.1021/ie101048b |Ind. Eng. Chem. Res. 2011, 50, 2831–2838
Industrial & Engineering Chemistry Research
ARTICLE
Figure 7. Total current for θ̂1 (see Table 4) for scan rates (top to bottom) T = 10, 1.0, and 0.1 (left) and the difference between the outputs produced by θ̂1 and θ̂2 (right). C¥ = 1.
correspond to the directions of the highest spread of the parameters; i.e., there are these combinations of the variation of the parameters, which have the smallest effect on the output and cannot be identified. Therefore, we have to project our parameters' space onto the subspace orthogonal to s1,...,sN0 ; i.e., the formula for matrix P reads P ¼ fsN 0 þ1 ,sN 0 þ2 ,:::,sN g T
Table 4. Solution of the Optimization Problem (10) under the Conditions of Figure 7 θ̂1
ð30Þ
To check our procedure, we performed an identifiability analysis according to eq 10, where the function F(θ1,θ2) has been calculated according to eq 26. The obtained values of η are η = 18.8 for N0 = 2 and η = 5.78 for N0 = 3. Therefore, there is a significant improvement in comparison to η = 50 obtained at N0 = 0. 3.4. Identifiability of the Entire Set of Parameters. Finally, we try to identify the experimental conditions, which allow for reliable estimation of the entire set of parameters. It has been noted earlier that a single experiment is not sufficient for this purpose. As a first attempt we performed calculations with varying scan rate; namely, we simultaneously solved eqs 2-5 for T = 10, 1.0, and 0.1, calculated the function HT(θ1,θ2) for each of the scan rates, and then H(θ1,θ2) in eq 10 was calculated
θ̂2
θ̂min
θ̂max
E0
7.99
7.83
6.0
10.0
R
0.436
0.698
0.3
0.7
k0
3.23
9.49
1.0
10.0
Ru Cdl
1.55 6.08 10-4
1.98 4.67 10-4
0.0 0.0
3.0 10-3
D
2.63
3.07
0.25
4.0
as H = min(HTi). The results are presented in Figure 7 and Table 4. The results for θ2 are not shown because they are undistinguishable from those for θ1. The computed value of the identifiability criteria is η = 23, and the identifiability of the system is poor. As one can see, there are sets of parameters with rather different values of k0, which produce very close outputs at three scan rates. It has been suggested in refs 1 and 2 that changing the concentration of the soluble phase might improve the identifiability. To verify this suggestion, we performed the calculations 2837
dx.doi.org/10.1021/ie101048b |Ind. Eng. Chem. Res. 2011, 50, 2831–2838
Industrial & Engineering Chemistry Research with three different scan rates and with C¥ = 1.0 and 0.1, which makes the total number of the virtual experiments equal to 6. There is only a modest improvement in comparison to the previous case: η = 16.4. The main problem is associated with the high uncertainty in the preexponential factor k0. It is in high contrast with the reversible case. In our other simulations we fixed k0 = 106 and R = 0.5 and examined the identifiability of E0, Ru, Cdl, and D with T = 10, 1.0, and 0.1 and C¥ = 1.0. The computed value of the identifiability criteria is η = 4.17, and the system is well-identified.
4. CONCLUSIONS In the present work we implemented a numerical test in order to analyze the identifiability of an electrochemical system. The method allows for the global identifiability analysis of a nonlinear system prior to undertaking any actual experimentation, The identifiability test is based on an unconstrained minimization problem, which has been solved by a random search method. It is shown that if the effect of resistance is negligible, the entire parameter’s set can be identified from a single ac voltammetry experiment, and an experimental design procedure is proposed in order to improve the identifiability. Comparison of the efficiency of different waveforms indicates that the identifiability is not very sensitive to the particular form of the applied potential, and optimized potential offers only a small improvement in comparison to the square wave and sawtooth waveforms proposed in the literature. If the effect of resistance cannot be neglected, several experiments with varying scan rates and varying concentrations of the soluble phase should be performed, or, alternatively, one can choose a subset of parameters, which can be identified with sufficient accuracy. In the present work we propose a statistical postprocessing procedure, which allows for selection of an identifiable linear combination of the parameters. ’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ REFERENCES (1) Sher, A. A.; Bond, A. M.; Gavaghan, D. J.; Harriman, K.; Feldberg, S. W.; Duffy, N. W.; Guo, S.-X.; Zhang., J. Resistance, capacitance, and electrode kinetic effects in fourier-transformed largeamplitude sinusoidal voltammetry: Emergence of powerful and intuitively obvious tools for recognition of patterns of behavior. Anal. Chem. 2004, 76, 6214–6228. (2) Tan, Y.; Stevenson, G. P.; Baker, R. E.; Elton, D.; Gillow, K.; Zhang, J.; Bond, A. M.; Gavaghan, D. J. Designer based fourier transforme dvoltammetry: A multi-frequency, variable amplitude, sinusoidal wave-form. J. Electroanal. Chem. 2009, 634, 11–21. (3) Alifanov, O. M.; Artyukhin, E. A.; Rumyantsev, S. V. Extreme Methods for Solving Ill-Posed Problems with Applications to Inverse Problems; Begell House: Redding, CT, USA, 1995. € (4) Ozisik, M. N.; Orlande, H. R. B. Inverse Heat Transfer; Taylor and Francis: Oxford, U.K., 2000. (5) Isakov, V. Inverse Problems for Partial Differential Equations; Springer: Berlin, 2006. (6) Yates, J. W.T.; Evans, N. D.; Chappell, M. J. Structural identifiability analysis via symmetries of differential equations. Automatica 2009, 45, 2585–2591. (7) Evans, N. D.; Chappell, M. J. Extensions to a procedure for generating locally identifiable reparameterisations of unidentifiable systems. Math. Biosci. 2000, 168, 137–159.
ARTICLE
(8) Walter, E. Identifiability of State Space Models; Springer: Berlin, 1982. (9) Zimmerman, W. B.; Rees, J. M.; Craven, T. J. Rheometry of nonNewtonian electrokinetic flow in a microchannel T-junction. Microfluid Nanofluid 2006, 2, 481–492. (10) Zimmerman, W. B. Electrochemical Microfluidics. Chem. Eng. Sci. 2010, in press. (11) Braumann, A.; Kraft, M. Incorporating experimental uncertaintiesinto multivariate granulation modelling. Chem. Eng. Sci. 2010, 65, 1088–1100. (12) Braumann, A.; Kraft, M.; Mort, P. Parameter estimation in a multi-dimensional granulation model. Powder Technol. 2010, 197, 196– 210. (13) Man, P.; Braumann, A.; Kraft, M. Resolving conflcting parameterestimates in multivariate population balance models. Chem. Eng. Sci. 2010, 65, 4038–4045. (14) Man, P.; Braumann, A.; Kraft, M. Statistical approximation of the inverse problem in multivariate population balance modeling. Ind. Eng. Chem. Res. 2010, 49, 428–438. (15) Glowinski, R.; Lions, J. L.; Tremolieres, R. Numerical Analysis of Variational Inequalities; North-Holland: Amsterdam, 1981. (16) Vikhansky, A. Numerical analysis of the global identifiability of reaction-diffusion systems. Comput. Chem. Eng. 2010, submitted for publication. _ om, K. J. On structural identifiability. Math. (17) Bellman, R.; Astr€ Biosci. 1970, 7, 329–339. (18) Vikhansky, A.; Kraft, M. A Monte Carlo methods for identification and sensitivity analysis of coagulation processes. J. Comput. Phys. 2004, 200 (1), 50–59. (19) Vikhansky, A.; Kraft, M.; Simon, M.; Schmidt, S.; Bart, H.-J. Population balance modelling of a droplet size distribution in a RDC: Aninverse problem approach. AIChE J. 2006, No. 52, 1441–1450. (20) Vikhansky, A.; Kraft, M. Two methods for sensitivity analysis of coagulation processes in population balances by a monte carlo method. Chem. Eng. Sci. 2006, 4966–4972. (21) Spall, J. C. Multivariate stochastic approximation using a simultaneous perturbation gradient approximation. IEEE Trans. Autom. Control 1992, 37, 332–341. (22) Maryak, J. L.; Chin, D. C. Global random optimization by simultaneous perturbation stochastic approximation. IEEE Trans. Autom. Control 2001, 53, 780–783. (23) Gavaghan, D. J.; Elton, D.; Oldham, K. B.; Bond, A. M. Analysis of ramped square-wave voltammetry in the frequency domain. J. Electroanal. Chem. 2001, 512, 1–15. (24) Gavaghan, D. J.; Elton, D.; Bond., A. M. A comparison of sinusoidal, square wave, sawtooth, and staircase forms of transient ramped voltam-metry when a reversible process is analysed in the frequency domain. J. Electroanal. Chem. 2001, 513, 73–86. (25) Chu, Y.; Hahn, J. Integrating parameter selection with experimental design using uncertainty for nonlinear dynamic systems. AIChE J. 2008, 54, 2310–2320.
2838
dx.doi.org/10.1021/ie101048b |Ind. Eng. Chem. Res. 2011, 50, 2831–2838
