Nonlinear Parameter Estimation - American Chemical Society

ent methods for the actual estimation process. ... w|xk,/3] = M[xkl/?. (1). In (1) xk is the kth vector of indepen- dent and dependent variables and /...
1 downloads 0 Views 4MB Size
instrumentation

Nonlinear Parameter estimation can also consider vector measure­ ments where each measurement vector consists of two or more individual measurements taken at the same time. For the nonlinear case, the model becomes a function of the parameters and this is denoted as the scaler quan­ tity

Thomas A. Brubaker Department of Electrical Engineering Colorado State University Fort Collins, Colo. 80523

Kelly R. O'Keefe Department of Chemistry Colorado State University Fort Collins, Colo. 80523

ω[χφ] = ΐ[χκβ].

This is a second paper in a two-part sequence on the application of param­ eter estimation in analytical chemis­ try. Linear parameter estimation was treated in the first paper (1), and this will serve as background for the present discussion. In this paper non­ linear parameter estimation is de­ scribed with emphasis given to gradi­ ent methods for the actual estimation process. Two examples, one from spec­ troscopy and one from kinetics, are used to demonstrate how nonlinear parameter estimation can be applied. As mentioned previously, the para­ mount problem in successfully using parameter estimation in chemistry is the development of a parametric ana­ lytical model or a valid experimental data model for the chemistry being studied. For the linear case, the model is linear in the parameters and for scaler measurements the model out­ put as a function of the independent variables, possibly past dependent variables, and the parameters is given as

«1**01 = M[*k]0.

(1)

In (1) Xk is the kth vector of indepen­ dent and dependent variables and β is the parameter vector. The term ω[χΐ(ΡΊ is the model output for a spe­ cific vector Xk which implies known or measured values for the independent and dependent variables. The term M[xk] is a one row matrix whose ele­ ments are the independent and depen­ dent variables or are functions of these variables. Carefully note that we 0003-2700/79/A351-1385S01.00/0 © 1979 American Chemical Society

(2)

In both (1) and (2) the vector notation is simply a shorthand way of saying the model is a function of a set of in­ dependent and dependent variables and a set of parameters. An example is an output that is the sum of two expo­ nential decays, such as might be ob­ served for the phosphorescence of a binary mixture, where the model is given by o>[kAt,/3] = /3 3 e-ft kAt + /34β-Α*Δ« (3) where At is the increment of observa­ tion time. For both the linear and nonlinear models, the most common error crite­ rion describing how well the model fits the data is least squares where the summed square error is given as m

S - Σ \y[xk}-u\xUW.

(4)

k= l

In (4), y[xk] is the measurement for the known or measured values of the kth vector of independent and depen­ dent variables. In vector form the summed square error is given by S - (y - ω)'(.y - ω)

(5)

where y is a vector of all m measure­ ments and ω is a vector of all model outputs. As shown in the previous paper, for the linear case ω = Ηβ

(6)

where Η is a matrix of all of the M's given in (1). The resulting optimal es­ timate found by minimizing (5) is ρ· = (Η

Since Qi is quadratic in β about the point ft, Q; is bowl-shaped and has a single minimum. Since S(ft) is a con­ stant, taking the derivative of Qi with respect to β and setting the result to zero will give the bowl minimum. Doing this yields

H = qi + H|(j8 - ft) = 0.

(14)

Assuming ft is the ith estimate and letting ft+i = β = ith plus one esti­ mate (14) now gives

ft+1 =

ft-(Hi)-'qi.

(15)

Equation (15) is known as the Newton or Newton-Raphson method. This method of gradient selection is effec­ tive as long as Hi is positive definite so that we know the algorithm is pointing toward a minimum. This is generally true if one is close to the minimum. However, over the total parameter

(17)

where Jj is a matrix of first derivatives with elements

Qi = S(ft) + (qdHfi - ft) + i (j8 - ft)*Hi09 - ft). (13)

(16)

_

afjU.rf)

(18)

oft,

This selection of Aj is useful for chem­ ical parameter estimation because it is relatively easy to compute on an itera­ tive basis. In addition, if the model chosen is correct, Ji'J; is proportional to the Hessian as the error approaches the minimum, thus, in this region the algorithm behaves in a theoretic manner. Examples. In the first example, we will demonstrate nonlinear parameter estimation using the simple sampled process x(kAT) = e-

-kAT

(19)

where τ is the time constant of the ob­ served signal and Δ Τ is the sampling interval. Such signals commonly occur in time-resolved fluorescence and phosphorescence spectroscopy and the task is to estimate r when the data are corrupted by noise. Using this model with the collected data y, at intervals Δ Τ , the summed square error is

:

A

2

4

6

6 9 Estimate Number

Time (Seconds) Figure 1. Ideal single exponential decay

A

12

15

Figure 2. Convergence of estimates for time constants of single exponential decay

1386 A · ANALYTICAL CHEMISTRY, VOL. 51, NO. 13, NOVEMBER 1979

0.7

α 0 35 ο .Û

Actual

Simulated

0

2

3 Error Surface

Milliliters of Picrate Figure 3. Actual and simulated photometric titration curves with kinetic complications

S = Σ (y(iAT) -

β-

ίΔΤ/ΐ 2

) ·

(20)

i=l

For the case of τ = 2 seconds, the ideal data is given in Figure 1. The parame­ ter estimate is 2.015 seconds with ad­ ditive noise with zero mean and vari­ ance 0.1. This corresponds to a signal to noise ratio of about 3. The conver­ gence of the Marquardt Algorithm for an initial guess of 1.0 seconds is shown in Figure 2. Note that the convergence is dependent on the initial guess and the amount of noise power, but that in this case, convergence to a final esti­ mate occurs in less than 9 iterations. The second example involves the es­ timation of kinetic parameters from a photometric titration curve. If a titrand, A, is titrated at a constant rate with titrant Τ to produce a product, P, it is the sole absorbing species at the measurement wavelength, and in ad­ dition the titration reaction is A + T:

P.

Figure 4. Global error surface for estimate of kinetic parameters

task is to estimate kj and k2 so as to obtain good agreement between simu­ lated titration curves and those exper­ imentally measured. When 3.70 Χ 10" 4 mol L" 1 creati­ nine was titrated with 1.19 X 10 - ' 2 mol L " 1 picrate, both in 0.200 mol L" 1 NaOH, the data labeled actual in Fig­ ure 3 were obtained. Pairs of initial guesses for ki in the range of 0 to 30 L m o l - 1 sec - 1 and for k2 in the range 0 to 3 s e c - ' were supplied to the simplex parameter estimation algorithm; the differential equation was simulated for the guesses; and summed square errors were calculated using equation 4. Upon convergence, values of kj and k> were obtained that yield the titra­ tion curve labeled simulated in Figure 3. Repeated titration/estimation gave estimates of 4.67 ± 0.11 L m o l - 1 s e c - 1 for k! and 0.00172 ± 0.00012 s e c - 1 for k2.

Figure 4 shows the summed square error surface for pairs of ki and k 2 in the parameter space. The complex na­ ture of this surface results from the nonlinearity of the model and estima­ tion structure.

Thomas A. Brubaker received his Ph.D. in Electrical Engineering in 1963 from the University of Arizona. Since that, time he has taught at the University of Missouri, Columbia, and is presently a professor of electri­ cal engineering at Colorado State University. His research interests are digital systems, digital signal process­ ing, and their applications to instru­ mentation areas such as chemistry.

Kelly R. O'Keefe received his Ph.D. in Analytical Chemistry from the Uni­ versity of Illinois in 1975. He joined the faculty of Colorado State Univer­ sity in the same year. His research in­ terests are in the areas of multiele­ ment atomic absorption spectrome­ try, kinetic methods of analysis, and applications of parameter estimation methods in analytical chemistry.

References

(1) T. Brubaker, R. Tracy, and C. Pomernacki, "Linear Parameter Estimation," Anal. Chem., 50, 1017A~1024A (1978). (2) S. N. Deming and S. L. Morgan, "Sim­ plex Optimization of Variables in Ana­ lytical Chemistry," Anal. Chem., 45, 278A-283A(1973). (3) Y. Bard, "Nonlinear Parameter Esti­ mation," Academic Press, New York, N.Y., 1974. (4) D. Marquardt, "An Algorithm for Least-Squares Estimation of Nonlinear Parameters," J. Soc. Ind. Appl. Math, 11, No. 2(1963). This work was supported by Contract 8403203 between Lawrence Livermore Laboratory and Colorado State University.

(21)

The differential equation that de­ scribes the rate of reaction at any time in the titration is d[A] dt

-kalAjjrCtt/V - C a Vi/V + [A]) + k2|CaVi/V-[A]|-r[A]/V

(22)

where r is the rate of titrant addition in Lsec - 1 , V; is the initial volume of the titrand in L, C t and C a are the ti­ trant and titrand concentrations, re­ spectively, in mol L - 1 , V is the volume at time t, and d[A]/dt is the rate of ap­ pearance of species A. If all of these parameters in a titration are known, digital or analog techniques may be used to simulate the titration curve, i.e., produce an absorbance versus time plot. If, on the other hand, the ki­ netic parameters are unknown, our

1388 A · ANALYTICAL CHEMISTRY, VOL. 51, NO. 13, NOVEMBER 1979