.3;- V vi>v.vu ‘ao. rgt... . 3 .21.!) 3.4:; .5. : 1...: .1? a .3: .. . to. x? :l»..: ...01. THESIS ' Illlll'llllllllllllllllllllllllll’llll 3 1293 00885 2174 This is to certify that the dissertation entitled Multi-Component Kinetic Determination with the Extended Kalman Filter presented by Brett Michael Quencer has been accepted towards fulfillment of the requirements for Ph.D. degree in Analytical Chemistry W/Kévt Major professor Date July 13, 1993 MSU is an Affirmative Action/Equal Opportunity Institution 0- 12771 LIBRARY E Michigan State niversity fl ‘— -___k PLACE iN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE g fifi‘ MSU is An Affirmative ActiorVEquai Opportunity institution ckiMmS—p. 1 , _ . M 7 MULTICOMPONENT KINETIC DETERMINATIONS WITH THE EXTENDED KALMAN FILTER By Brett Michael Quencer A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Chemistry 1993 ABSTRACT MULTICOMPONENT KINETIC DETERMINATIONS WITH THE EXTENDED KALMAN FILTER By Brett Michael Quencer This dissertation describes a new method of multicomponent analysis based on kinetic methods and a chemometric data processing technique known as the extended Kalman filter. The resolving power of this technique is based on differences in both the kinetics and spectra of a series of parallel reactions. There have been many previous methods developed to allow kinetic determination of more than one component in a mixture simultaneously. Unfortunately, all of these previous techniques have suffered from some severe limitations which arise due to invalid assumptions made when developing the model for the system. The extended Kalman filter is shown to be free from these limitations, and to have improved applicability to multicomponent determinations, including the ability to determine a greater number of components than previously. A series of simulations have been carried out in order to test the validity of this method of data analysis and to find conditions under which the extended Kalman filter method is likely to fail. The results of the simulation studies are described in detail. It is shown that as long as there are some spectral differences, this method is applicable even in a series of three parallel reactions where all three rate constants are equal. It is shown that the method is applicable even under conditions of severe Spectral and kinetic overlap. A rigorous test of the method was performed on a chemical system which has severely overlapped spectra and kinetics. The system consists of a series of metal species which react simultaneously with 4-(2-pyridylazo)resorcinol, a complexing agent, to form a series of closely related complexes. It is shown by comparison the extended Kalman filter method greatly improves the accuracy, precision, and utility of multicomponent kinetic methods. The new method is also shown to have a greatly reduced dependence on a number of factors which severely affect previous techniques. These include the previous assumption of rate constant invariance between runs and use of invalid system models. To Lisa and Alyssa, for making the drive worthwhile iv ACKNOWLEDGMENTS It is nice to finally get a chance to thank the people who have made my years at Michigan State University more enjoyable and truly worthwhile. First and foremost, I would like to thank Dr. Stan Crouch for his guidance and assistance whenever a stumbling block arose in my research. There are several other people at Michigan State who have made the graduate school experience more pleasant. Stephen Medlin, who always was there with a kind word and encouragement; Jim Ridge, for always being willing to help with electronics problems, and always good for a laugh; Larry Bowman, for his knowledge and willingness to impart it; Dave McLane, for his many interesting philosophical discussions, some of them relating to chemistry, and for his stress relief treatments (the eighteen hole kind); Tony Hsieh, for his interest and involvement; Dana Spence, for being a pretty good guy for a first year; and of course, Eric Hemenway, for always being willing to go to lunch, and the occasional help with computer problems. On a personal note, I would like to thank my wife, Lisa, for her support and encouragement, and our daughter, Alyssa, who is a joy to be around. I certainly would not have ever even made it to graduate school were it not for my parents, Mike and Mary Quencer, and their instilling in me a love of learning. Finally, I would like to thank my in-laws, Ed and June Goldsmith, for always being supportive and understanding of the demands of graduate school. vi TABLE OF CONTENTS CHAPTER 11. LIST OF TABLES _________________________ LIST OF FIGURES ________________________ INTRODUCTION __________________________ Multicomponent Kinetic Determinations-An Historical Perspective ___________________ A. Introduction B. Principles of Multicomponent Kinetic Methods_ 1. Methods for Large Rate Differences 2. Traditional Methods for Small Rate Differences 3. Modern Computer-Based Methods C. Applications of Multicomponent Methods 1. Graphical Extrapolation Methods Proportional Equations Based Methods__ Single-Point Method Computer Methods Multiwavelength Methods Miscellaneous Methods 9999’!" D. Conclusions LIST OF REFERENCES vii PAGE xi xiii 31 32 34 37 37 41 43 43 47 CHAPTER PAGE 111. The Kalman Filter in Analytical Chemistry ______________________________ 5 2 A. Introduction 52 B. The Extended Kalman Filter Algorithm 55 C. Applications 60 1. Linear Kalman Filter 60 2. Extended Kalman Filter 62 3. Adaptive Kalman Filter 63 4. Kalman Filter Networks 64 D. Conclusions 64 LIST OFREFERENCES 66 IV. Simulations Studies of Two-Component Systems ________________________________ 7 1 A. Kinetic Model 71 l. Two-Component Mixture Following Irreversible Kinetics 7 1 2. Use of the Model in the Extended Kalman Filter 7 4 B. Experimental 7 6 1. Simulated Data 7 6 2. Generation of Random Noise 7 7 3. Data Processing 7 8 C. Results and Discussion 80 1. Error vs. Percent Completion 8 2 2. Spectral Overlap 8 9 3. Effect of Random Noise 9 6 4. Effect of Varying Estimated Rate Constant 101 5. Other Considerations 108 viii CHAPTER VI. LIST OF REFERENCES Simulations of Complex Chemical Reactions A. Kinetic Models 1. General Irreversible Reactions 2. Reversible Models 3. Use of Models with the Kalman Filter_ B. Experimental C. Results and Discussion 1. Three-Component Irreversible System___ 2. Four-Component Irreversible System___ 3. Eight-Component Irreversible System __ 4. Reversible Systems D. Conclusions Simultaneous Kinetic Determination of Lanthanides with the Extended Kalman Filter __________________________________ A. Introduction B. Experimental Section 1. Apparatus 2. Reagents 3. Procedure 4. Computational Aspects C. Results and Discussion 1. Single Component Studies 2. Two- and Three-Component Studies D. Conclusions ix PAGE 110 111 111 113 115 117 118 118 123 128 131 135 137 137 140 140 140 140 141 141 141 154 165 CHAPTER LIST OF REFERENCES APPENDIX 1. Program NRNDAT LIST OF REFERENCES APPENDIX 11. Program EKFNRN PAGE 168 169 176 177 TABLE 2.1 2.2 2.3 2.4 2.5 2.6 4.1 4.2 5.1 LIST OF TABLES Applications utilizing graphical extrapolation methods Applications utilizing the method of proportional equations Modern applications utilizing the single- point method Applications of computer methods to multi- component kinetic determinations Applications utilizing multiwavelength detection Miscellaneous applications Values of molar absorptivity times path- length (L mol‘l) used for simulations unless otherwise noted. [C1]0=0.OOOI M, [C2]o=0.0001 M, [R]o=0.01 M, k1=50 5'1, k2=lO-50 8'1 Molar absorptivities for spectral overlap study. Concentrations and rate constants appear in Table 4.1 Values of molar absorptivity for the three component system studied. [C1]o =[C2]o =[C3]o= 1 x 10‘4 M, [R]o =0.0l M, k1=100 S'l,k2=10- 100 S'l,k3=10-100 8'1. xi PAGE 33 35 38 39 42 44 81 90 119 TABLE 5.2 5.3 5.4 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 Final estimated concentration, and error in the estimated concentration for the four-comp- onent case with moderate spectral overlap.__ Final estimated concentration, and error in the estimated concentration for the four-comp- onent case with greater spectral overlap. Estimated concentration of each component and error in components for eight-comp- onent case. Conditions for spectral acquisition of single component data Conditions for spectral acquisition for two and three component mixtures Final extended Kalman filter estimates of cone- entration for mixture of lanthanum and praseodymium. Final extended Kalman filter estimates of cone- entration for mixture of lanthanum and neodymium. Final extended Kalman filter estimates of conc- entration for mixture of praseodymium and neodymium. Final extended Kalman filter estimates of cone- entration for mixture of lanthanum, praseo- dymium, and neodymium. Final extended Kalman filter estimates of conc- entration for mixtures of lanthanum, praseo- dymium, and neodymium. Comparison of the method of proportional equations and the extended Kalman filter_ xii PAGE 126 127 130 143 155 156 157 158 159 161 163 FIGURE 3.1 4.1 4.2 4.3 4.4 4.5 4.6 4.7 LIST OF FIGURES Pictorial representation of true absorbance (0) versus Kalman filter estimate (-1-) Histogram for added noise. 20000 total points, theoretical standard deviation=0.005_ Error in estimated concentration of [C1]o vs. % completion of the reaction. [C1]o=1x10'4 M, k=100 SCC‘l Error in estimated C1 as a function of concen- tration ratio and rate constant ratio. C2 constant for all experiments Error in estimated C2 as a function of concen- tration ratio and rate constant ratio. C2 constant for all experiments Error in estimated C] as a function of concen- tration ratio and rate constant ratio. C1 constant for all experiments Error in estimated C2 as a function of concen- tration ratio and rate constant ratio. C1 constant for all experiments Estimated concentration of C2 as a function of k1/k2 with no spectral overlap. (Solid line)- True concentration, (---)-5% error bars. xiii PAGE 53 79 83 85 86 87 88 92 FIGURE 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 4.16 Estimated concentration of C2 as a function of k1lk2 with medium spectral overlap. (Solid line)-True concentration, (---)-5% error bars. Estimated concentration of C2 as a function of k1/k2with one wavelength. (Solid line)-True concentration, (---)-5% error bars. Estimated concentration of C2 as a function of k1/k2with total spectral overlap. (Solid line)-True concentration, (---)-5% error bars. Error in the average estimated concentration of C2 as a function of rate constant ratio and the standard deviation of the noise Standard deviation of the mean of C2 as a function of rate constant ratio and the standard deviation of the noise Absorbance and estimated absorbance for k'1/k'2=l.20 and standard deviation of the noise of 0.10. (Solid line)-True absorbance, (+)-estimated absorbance. Error in estimated concentration of C1 as a function of the initial estimates for the rate constants when k1/ k2 =5.00. Error in estimated concentration of C2 as a function of the initial estimates for the rate constants when k1/ k2 =5.00. Error in estimated concentration of C1 as a function of the initial estimates for the rate constants when k1/ k2 =2.00. xiv PAGE 93 94 95 98 99 100 102 103 104 FIGURE 4.17 4.18 4.19 5.1 5.2 5.3 5.4 5.5 5.6 6.1 6.2 Error in estimated concentration of C2 as a function of the initial estimates for the rate constants when k1/k2 =2.00. Error in estimated concentration of C] as a function of the initial estimates for the rate constants when k1/k2=1.00. Error in estimated concentration of C2 as a function of the initial estimates for the rate constants when k1/k2=1.00. Error in estimated value of C1 versus k2 and k3 Error in estimated value of C2 versus k2 and k3 Error in estimated value of C3 versus R2 and k3 Synthetic spectra of eight components. (solid line)-products, (dotted line)-common reagent. Error in estimated C1 versus ratio of forward to reverse rate constants. Error in estimated concentration of C1 (0) and C2 (+) versus k1/k2. Structure of 4-(2-pyridylazo)resorcinol (PAR) Molar absorptivities of lanthanides and 4- (2-pyridylazo)resorcinol (PAR). (o)-PAR, (+)-La(PAR) complex, (_)-Pr(PAR) complex, (---)-Nd(PAR) complex. XV PAGE 105 106 107 120 121 122 124 132 134 137 144 FIGURE 6.3 6.4 6.5 6.6 6.7 6.8 6.9 Estimated initial concentration of 2.88 x 10'5 M lanthanum with respect to wavelength. True concentration shown by solid line. 8.5 ms scans. Estimated initial concentration of 1.42 x 10'5 M praseodymium with respect to wavelength. True concentration shown by solid line. 8.5 ms scans. Estimated initial concentration of 1.39 x 10'5 M neodymium with respect to wavelength. True concentration shown by solid line. 8.5 ms scans. Estimated concentration using 1 pixel at 495 nm of 1.44 x 10'5 M lanthanum with respect to time. Error bars 21:] standard deviation. (_)- True concentration, (---)-10% error in conc- entration. 10.5 ms scans. Estimated concentration using 1 pixel at 495 nm of 1.42 x 10‘5 M praseodymium with respect to time. Error bars :tl standard deviation. (_)- True concentration, (---)-10% error in conc- entration. 10.5 ms scans. Estimated concentration using 1 pixel at 495 nm of 1.39 x 10'5 M neodymium with respect to time. Error bars 1:] standard deviation. (_)- True concentration, (---)-10% error in conc- entration. 10.5 ms scans. Residual absorbance values using 1 pixel at 495 nm for a reaction of 1.39 x 10'5 M neodymium with respect to time. xvi PAGE 145 146 147 149 150 151 153 FIGURE PAGE 6.10 Estimated concentrations of praseodymium (o) and neodymium (+) versus time using two wavelengths. True concentrations: [Pr]=1.06 x 10'5 M, [Nd]=l.04 x 10'5 M, [PAR]=2.00 x 10’4 M. 164 xvii CHAPTER I Introduction Curiosity is the very basis of education and if you tell me that curiosity killed the cat, I say only that the cat died , nobly. Arnold Edinborough Reaction-rate methods of analysis are becoming increasingly popular techniques in analytical chemistry. This can be demonstrated by an increasing rate of appearance of articles in the literature relating to the subject. Some examples of this are given in Chapter 11. These methods offer several distinct advantages over equilibrium-based methods, including speed and selectivity. Perhaps one of the more interesting applications of reaction-rate methods is that of multicomponent analysis. The goal of multicomponent techniques is to obtain accurate estimators of the true concentrations of two or more components in a sample. Of course, in a perfect world, this task would not be very difficult because, ideally, each component in the mixture would exhibit a unique response from the detector. It would also be ideal if there were no overlap of the kinetic responses of any of the components in the sample. The world in which we live is not ideal, and most samples of analytical interest exhibit overlapped detection and kinetic behavior. How then, does one go about resolving these overlapped responses? There have been many attempts at providing methods to do just this. Chapter II discusses previous techniques which have been utilized to solve this problem. The limitations to these previous techniques are discussed, and justification is given for the development of a new technique, namely, the extended Kalman filter, which is, of course, the subject of this dissertation. Chapter III introduces the Kalman filter, and its nonlinear form, the extended Kalman filter. The uses of the various forms of the Kalman filter to solve chemical problems since its inception are described. At this point, the extended Kalman filter algorithm is presented, and its applicability to the current problem is discussed. The advantages of using this approach for multicomponent kinetic determinations are put forth. In order to test the advantages and limitations of using the extended Kalman filter for multicomponent kinetic determinations, an extensive series of simulation studies were performed. In Chapters IV and V, these simulation studies are described, and the results are discussed. In Chapter IV, a model for a two-component irreversible system to be studied is developed, while Chapter V develops models for more complex chemical systems. The ways in which these models are applied for use in the extended Kalman filter are described. The general kinetic model used for the irreversible simulation studies is that of a pseudo-first-order reaction, described as: C+R—L-)P (1.1) where the reaction is first-order with respect to species C, species R is in excess, k is the full second-order rate constant, and P is the product. Under pseudo-first-order conditions, the pseudo-first-order rate constant, which can be denoted by k', is equal to the product of the true rate constant k, and the concentration of species R. In Chapters IV and V, the models used to calculate the simulated data are described, as is the model used in the extended Kalman filter algorithm to process the data. Of course, in real situations, noise is present, so in order to present a more realistic model, noise is added to the simulated data. The method for generation of this noise is described in Chapter IV, and its characteristics are discussed. The results of the simulation studies are described, and the influence of various parameters on the accuracy of the filter estimates are analyzed. Finally, some conclusions are drawn about what conditions it is possible to apply the filter to a chemical system, and under what conditions application of the filter would not be advisable. Simulation studies are useful because they allow one to vary parameters with a degree of accuracy and precision not possible in a real experimental system. However, simulated data behave in a predictable manner, which may lead one to let biases creep into the conclusions made about a technique. Also, because real systems have unknown effects, such as synergism, impurities which interfere, ionic strength variations, and pH effects which are not modeled, it is necessary to test the technique on a real chemical system. Chapter VI describes a chemical system used to test the version of the extended Kalman filter developed for this research. The chemical system used was that of the general complexing agent, 4-(2- pyridylazo)resorcinol, which complexes a large number of metal species. The utility of this system as a good test case for the extended Kalman filter is described, as are some of its uses in process chemistry. The experimental procedures used for this study are described, and the results are analyzed. The results of the estimates of concentration by the extended Kalman filter technique are compared to the results obtained by other analysis methods under the same conditions. The advantages of using this technique are discussed, as are the disadvantages. The applicability of the Kalman filter method to even more complicated reaction schemes is described. In Chapter VII, conclusions on the utility of the technique to a wide variety of reactions are made, and the future direction of this work, as the author sees it, is discussed. One area where this technique should find great applicability is for principal component analysis using a network of extended Kalman filters. A possible extension of the current work to principal components analysis is described in some detail. Since simulation studies consist of such a large part of the current work, the program, NRNDAT, used to calculate the artificial data for the case of a general number of components following non- reversible kinetics is presented in Appendix I. This program is referred to in Chapters IV through VI. The program first calculates -the true values for the absorbance at each wavelength for each point in time in the reactions, and then adds a random noise value to the true absorbance in order to simulate the noise characteristics of an experimental system. Finally, Appendix 11 contains the program listing for EKFNRN, the extended Kalman filter algorithm used to process the data for a general number of components following non-reversible reactions. In other cases of different types of reacting systems described in this work, this program was modified in order to work with the parameters desired from those systems. However, the general algorithm remained the same for all systems studied, with the only changes being made to the system model, the linearized measurement function, and the number of components analyzed. The areas which were changed when using different models are described. The new technique described in this dissertation should become a powerful method in analytical chemistry. It has overcome most of the limitations inherent in previous approaches to problems of this nature. CHAPTER II Multicomponent Kinetic Determinations-An Historical Perspective Nothing in life is to be feared. It is only to be understood. Marie Curie A. Introduction Kinetic methods of analysis have become increasingly popular in recent years. This increase in popularity can be seen by a number of international conferences devoted to kinetics in analytical chemistryl, and by two recent books on the subject2»3. There have also been several recent review articles which deal with kinetic methods in analytical cheniistry4t5»6t7:3. Kinetic-based determinations offer several advantages over equilibrium methods. Among these advantages are those of simplicity, speed, and precision9. Crouch“), asserts that one of the major reasons for the increased interest in kinetic methods is the increased use of automation in kinetics-based procedures. The use of automation in kinetic methods can be traced to 1960 with the introduction of an automated, variable time method for glucose“. 7 The history of developments from this point until 1988 is summarized in reference 2. Of course, kinetic methods have some disadvantages when compared to equilibrium-based methods as well. Since the analytical results are dependent on the measurement of a time-dependent quantity which is influenced by a rate constant, any factor which affects the value of the rate constant can potentially affect the accuracy and precision of a kinetic technique. Rate constants, and therefore most kinetic methods, are dependent on factors such as pH, ionic strength, and temperature. Therefore, in order to obtain good results with kinetic based methods, one has to control experimental conditions more so than one usually would in an equilibrium-based method. Kinetic methods take the time dependence of a chemical reaction or an instrumental response into account in order to obtain analytical information. Although historically analytical chemists have given little or no recognition to kinetic methods, many popular techniques are, in fact, kinetic in nature. Pardue4 defines a kinetic method as any analytical procedure in which the measurement step is influenced by a transient process and proposes that a majority, rather than a minority of modern analytical methods are kinetic in nature. By this definition, any chromatographic, atomic, or mass spectrometric technique, as well as a number of other analytical techniques are kinetic in nature, and can be properly classified as kinetic methods. Multicomponent kinetic methods, sometimes referred to as differential kinetics, are one of the more interesting facets of reaction-rate methods of analysis. In the historical context, the term 8 "multicomponent" is somewhat of a misnomer, since most of the techniques developed to date consist of the resolution of two, or at most, three components. These techniques might better be called dual-component kinetic methods. The general goal of multicomponent kinetic methods is to determine the kinetic parameters or the analytical concentrations of the reactants in a system. A general mixture amenable to multicomponent kinetic methods can be described as follows: Cl +R-—:1—>P, ( 2.1) C2 + R—3—->P2 where C1 and C2 are two different species which react with a common reagent, denoted by R, with two different rate constants, R1 and k2 to form similar, but not identical products, P1 and P2. Reagent concentrations are usually maintained such that pseudo-first-order kinetics apply, in which case equation 2.1 reduces to: C‘ 4513‘ (2.2) c2 —-£z—>P, where k'l =k1[R]t and k'z =k2 [R]t. By assuring that pseudo-first-order kinetics are followed, previous researchers were able to obtain analytical information using the model given in equation 2.2 above by a variety of methods. The biggest problem that occurs in all of the previous techniques, however, is that as the ratio of the two rate constants k'l/k'z becomes closer to unity, the errors in the estimated concentrations of the components in the mixture get larger and 9 larger. The most common methods for the kinetic resolution of mixtures are described below, followed by a discussion of several recent techniques that offer the possibility of improved determinations of several components, as well as increased applicability to more complicated systems. B. Principles of Multicomponent Kinetic Methods Historical approaches to resolving mixtures of reacting species by kinetic methods have centered on adjusting or taking advantage of differences in the reaction rates of the components”. These approaches are subdivided into two categories, methods for mixtures whose reactions have large rate differences, and methods for mixtures whose reactions have small rate differences. Since methods for mixtures with small rate differences make up the large majority of the work in this field, each of the most popular techniques are discussed in turn. 1. Methods for large Rate Differences If a mixture which is to be determined has large differences in reaction rates, the analysis is quite simple because each species can be treated separately. It is assumed that during any time interval, only a single species is reacting, and that the Other species have already completely reacted, or are reacting so slowly that they do not interfere with the reaction of interest. The time dependence of the concentrations of C1 and C2 in equation 2.2 can be given as: 10 MEL = kit [C1]. 1n_[C2], = k'zt [C2]. (2.3) If both C1 and C2 react during a given time period, a time independent expression for the two concentrations can be obtained by dividing the two equations in equation 2.3 above to yield: 11. I] ln-‘C2 ° [C2]. 1“. IO " RIF, N (2.4) I—A With this equation, Mark, Papa, and Reilley13 were able to determine the ratio of rate constants necessary to obtain a given error in the concentration of [C1]o. They found that if an error of <1% was desired, it was necessary to have a rate constant ratio of >500. If larger errors were tolerable, smaller rate constant ratios could be used; however, the errors always increased greatly as k'1/k'2 got smaller. In general, rate constant ratios of <50 yield errors of greater than 5% by neglecting the slower reaction. The reported errors assume that the measured parameters (e.g., absorbance, diffusion current) are equally as sensitive for C1 and C2. Criteria are also discussed13 whereby one can neglect the reactions of faster reacting components. If the reaction rate of C1 is very large compared to that of C2, then at long times, C1 can be considered to be completely reacted (i.e., [C1]t=0), at which time the 11 change in the signal will be directly proportional to C2. There are two limitations to this method. The first is the time necessary for the determination, and the second is that a reasonable amount of C2 must remain when C1 has completely reacted. This limits the ratio of rate constants under which precise results can be obtained. The actual ratio depends upon conditions, but is generally greater than 10:1. Other methods for simultaneous determination of mixtures of components are discussed13. For example, assume that a binary mixture of C1 and C2 is to be determined and C1 reacts essentially completely in 20 minutes. If k'1/k'2 =SOO, C1 can be easily determined, but a great deal of time must elapse before one can accurately estimate C2. Methods are also discussed for estimating concentrations under such conditions as: changing the reaction temperature after C1 has reacted to completion, changing the concentration of the common reagent, B, after the first reaction is over, or adding a catalyst after the first reaction. Because these methods are rather self-explanatory, they are not discussed further here. 2. Traditional Methods for Small Rate Differences If the ratio of the two rate constants is small, the techniques described above are not applicable because the assumptions necessary to neglect one of the reactions are no longer valid. For this case, special techniques have been developed in order to analyze mixtures. 12 The approaches to solving this problem may behdivided into three general categories: masking methods, methods involving changes in the kinetics of a system, and methods for systems having unchangeable rate constant differenceslz. Each of these traditional techniques is described, and then more specific approaches are discussed. Modern computer-based methods are discussed later. a. Masking methods Masking methods generally involve shifting the equilibrium of an interfering species so that it no longer reacts in the presence of the species of interest13. The most common example of the masking technique is the conversion of all interfering species into complexes of extremely high stability such that they no longer react with the reagentlz. A common example of this process is the titration of metal ions with EDTA. In this case, it is relatively easy to adjust the reaction conditions such that only one metal ion forms a stable complex with EDTA, or to add another complexing agent that selectively complexes an interfering metal ion so that it does not interfere with the ion of interest. b. Methods Involving Changing the Kinetics of a System The second method for determinations of mixtures with small differences in rate constant occurs when it is possible to alter the type of reagent used in order to obtain suitable rate differences”. In other cases, the rates of reaction are shifted into a more suitable 13 region. This could be because the rates of reaction are too rapid to determine the species under normal conditions. An example of this method is again the reaction of metal ions with EDTA. In many cases, the reactions of free metal ions with EDTA are too rapid for common kinetic techniques. If, however, a complexing agent is added to the mixture before the addition of EDTA, the reactions can be made slower. The rates of displacement of the preliminary complexing agent may be different for the ions being studied, in which case the rate constant ratios may be more suitable and analysis more practical. c. Systems With Unchangeable Rate Constant Differences The third method is actually a class of techniques that have been developed in order to obtain analytical information regarding systems with small differences in the rate constants. This is an important area of multicomponent kinetics, because often the rate constants cannot be separated to a sufficient degree by either of the two methods described above. Many of the previous approaches to this problem are well covered in the literature, and so will not be described in detail here. The most popular approaches to multicomponent kinetic determinations historically have been graphical extrapolation methods and the method of proportional equationsZJ3. There have also been several methods developed recently. These computer- based methods are described in section 3. 14 i. Graphical Extrapolation Methods The most common graphical method historically has been the logarithmic extrapolation method13. Linear extrapolation methods have also been used, however, they have not found as wide an application as logarithmic extrapolation methods. The concentration of the products in equation 2.2 can be given by the following if the two products are identical and the initial concentration of the product is zero: [P]! = ([Cllo -[C1]t)+([C2]o -[C2]t) (2'5) The concentrations of the two reactants at any time follow simple first-order kinetics as follows: [C1]. = [Clloe-kh (2 6) [C2].=[C2].°"‘"‘ ' By substituting equation 2.6 into equation 2.5 we obtain the concentration of the product at any time in terms of only the initial reactant concentrations and their respective rate constants. Thus: [p]! = ([c,]o —[c,]°c’*1‘)+([c,L - [odors-‘2') (2.7) Using equation 2.7, we can now obtain the following relationship between the concentrations of the reactants and products: 15 (1C1 11 + [(3211) = ([P]- " [PM = [Czioc-i"zt (2'8) The logarithmic extrapolation method assumes that when the faster reacting component, C1, has reacted to completion, the term [C,]°e""" in equation 2.7 becomes negligible. By taking the logarithm of both sides of equation 2.8 and plotting 1n([C,]‘ +[C2]‘)=ln([P]_° —[P],) versus time, one obtains a straight line with a slope of -k'2 and an extrapolated intercept of ln[C2]o. It is then possible to obtain the value of [C1]o by subtracting [C 2]o from the total initial concentration of reactants. The range of rate constant ratios where this method is applicable depends on a variety of factors, such as concentration ratio and what error is acceptable to the user. A complete error analysis is given by Mark, Papa, and Reilleylz; however, they found that under certain conditions, rate constant ratios of as low as 5 could be tolerated with this method. The linear extrapolation method requires the total amount of the two species be known”. For two reactions following second- order kinetics as described in equation 2.1, and if [R]=[C1]o +[C2]o, the rate of reaction can be expressed by: $=§=klkli¢l+k2[R].[C21. (29) where x is the amount of R consumed at any time, I. If C1 reacts faster than C2, then when the C1 complex has reacted essentially to completion, equation 2.9 gives, after integration: 16 x = k,[C2]‘([R]° - x)t + [c,]° (2.10) The value of [C1]o is then determined from the extrapolated intercept at t=0 of a plot of x vs ([R]o-x)t. ii. Method of Proportional Equations The method of proportional equations is based on the principle of constant fractional life, which applies to first-order reactions and to systems which can be reduced to pseudo-first-order kinetics”. A species is said to have a constant fractional life if, at any given time, a constant fraction of the initial concentrations has reacted irrespective of the initial concentration. This is a well-known property for first-order reactions, where half-lives are often used as a measure of the reaction rate. It can be shown that the concentration of the product at any time is directly proportional to the initial concentration of the reactant. For example, in a one- component system following non-reversible first- or pseudo-first— order kinetics, given by: Cl—lL—aP (2.11) the concentration of the product P1 at any time, t, is given by: [p], =[Cl]°(l-e""") (2.12) 17 01' [P], =x,[c,]° (2.13) where xl =(1—c‘ki‘) (2.14) Therefore, the concentration of P at any time, t, is directly proportional to the initial concentration of C1. This treatment can easily be extended to the case of two reactants in a mixture. If C2 reacts under first-order conditions to form P, the concentration of P at any time, t, will be proportional to the initial concentration of C2 by: [p], = x,[c,]° (2.15) If the two reactions are independent, the concentration of P at time t for a mixture of C1 and C2 can be given by: [P], = x,[c,]o +x,[(:,]o (2.16) and the concentration of P at some later time, t', can be given by: [PL = x,'[c,]a +x,'[c,]o (2.17) 18 In order to use the relationships in equations 2.16 and 2.17 in a real determination, it is first necessary to determine the values of X1 and X2 at times t and t'. These values can be determined experimentally from equation 2.13 by measuring the amount of P produced by known concentrations of C1 only and C2 only during the two time intervals. These values could also be calculated by substituting the known reaction rate constants into equation 2.14. The analysis of a two-component mixture is achieved by measuring the concentration of P at times t and t'. The resulting data are then substituted into equations 2.16 and 2.17 which are solved simultaneously in order to determine the concentrations of C1 and C2. This method can be extended to the case of a larger number of reacting specieslz. However, the errors involved in determining the values of P and of the various X values at different times limits the number of species which can be determined to three, or at the most, four. The accuracy which can be obtained using this method is dependent on a number of variables. A complete error analysis is found elsewherelz; however, this method is generally not applicable to mixtures with rate constant ratios below 5/1. A variation of the method of proportional equations has also been developed. This method is known as two point kinetic simultaneous determination using flow injection analysis (FIA)15. This method is quite similar to the method of proportional equations, except it includes a term for the dispersion in the FIA system at each of the two points measured. Ultimately, however, these approaches, which have seen their popularity decline in recent years, are quite limited in a way which 19 has been overcome in recent years. These basic limitations are discussed for graphical extrapolation methods and for the method of proportional equations below. Graphical extrapolation methods are limited by erroneous assumptions regarding the system model. In other words, this technique makes the assumption that one reaction is complete while the other reaction is only in its initial stages. As the ratio of the rate constants approaches unity, this assumption becomes less and less valid. For this reason, the method of graphical extrapolation becomes limited as the rate constant decreases. The method of proportional equations finds its main limitation from a different source. Since only two data points are used in order to determine concentrations in a two-component mixture, the precision of this method is not very good. In a kinetic method, usually a fairly large number of data points are obtained, often on the order of 50-1000. But since the method of proportional equations only uses a small fraction of the available data, the precision of this method is limited. Kopanica and Staral6 have summarized the conditions in reaction systems under which each of the above techniques is applicable. iii. Single Point Method The single point method can be applied to systems such as those described in equation 2.2. This method requires knowledge of the total concentration of the mixture and the extent of reaction at a 20 single selected time during the course of the reaction”. The method is based on a plot of ([c,]‘ +[c,]‘) / ([c,]o +[c,]°) or ([C]- -[C],)/ [C], at any time, t, versus initial mole fraction of C1 in the mixture. This plot yields a straight line which has a slope of (6“ -c"‘3‘) and intercepts of 6"” and e4" at a mole fraction C1 equal to zero and unity, respectively. The plot is used as a calibration curve, and is easily constructed by measuring the extent of reaction of pure C1 and C2 at the chosen value of t, and then drawing a straight line through the points. This method, like the method of proportional equations, suffers from a lack of precision. Using data from only one time, when much more is readily available, severely limits the precision which can be obtained with this method. It is likely that this lack of precision is one of the main reasons that the single point method has not found wider applicability. 3. Modern Computer-Based Methods The introduction of digital computers in the chemical laboratory has brought with it an increase in the power of kinetic methods of analysis. Since kinetic systems most often follow non- linear relationships, and computers are quite able to model these, many new applications have been discovered. Computers bring several advantages to chemists. Among these are speed of computation, ability to perform non-linear regression rapidly, automation of data acquisition and processing, and graphical applications. Several different approaches involving the use of 21 computers for multicomponent kinetic determinations have been utilized. Several of the more popular methods are discussed in this section. a. General Methods Since kinetic systems are non-linear in nature (except for zero order systems), methods developed before computers became readily available generally tried to transform the data to fit a linear function. However, since non-linear regression methods are able to fit the function without the need for linearization, computers allow the non-transformed model to be used. Pausch and Margerum13 were the first to apply digital computers to kinetic methods of analysis. They used an IBM 7094 to determine mixtures of magnesium, calcium, strontium, and barium by the method of proportional equations, using 30-60 data points instead of just two, as had been used prior to the advent of computers. This method is described by extending equations 2.16 and 2.17 to having 30-60 simultaneous equations, each of the form [PL = x“ [q]o + x,‘ [C,]° (2.18) As long as all values of Xj are known, this becomes a system of 30-60 simultaneous equations with two unknowns. Using this greater number of data points allows for greater precision. A simplified linear least squares method was developed for the analysis of two- and three-component mixtures”. The two- 22 component system can be described by equation 2.2. If P1=P2, the concentration of the product, P, at any time can be expressed by: [p], f[C,]°(l—e’ki‘)+[C2]o(l—e'k3‘)+X (2.19) where X is the background signal. At each time, t, the value of Pi can be expressed by: [P]t = [C,]Oa, +[C,]°b, +X (2.20) where at =(l—e’kl‘) and bt =(l—e""‘). The shape of the response-time curve (formed by [P]1 versus time) is determined by the values of [C1]o and [C2]o and not by X. The goal of this method is to find the values of these initial concentrations which best fit the response-time curve. It is possible to measure [P]; at hundreds of times, and therefore to have hundreds of simultaneous equations of the form given by equation 2.20. The error between the measured value of [P]1 and the predicted value is denoted by ilt and is given by: u, =a,[C,L+b,[C,]°+X-P, (2.21) In order to obtain the least squares error, the sum of the squares of the errors for each data point from t=1 to n is obtained as given by: 2x10“):=‘2:‘1w,(a,[Cl]a+b,[C2]o +x-[P],) (2.22) t- 23 where wt is a wighting factor associated with each [P]1 value. The 20102 values are minimized with respect to each coefficient ([C1]o, [C2]o, and X) by partial differentiation and (omitting wt and the summation bounds) three equations resulting in three unknowns. The equations are: 2(af[C,]° + a,b,[C2]o + a,X) = 231W], 2(a,b,[c,]o + bf[c,]o + bx) = 2b,[P]t (2.23) 2(aic. 1. + big]. + x) = 2m. The linear combination coefficients are resolved by solution of the matrix equation, Ewtatz zwtatbt Ewtat C110 thaIIPL zwab, 2‘”th EWJJ, C21, = 2w,b,[P], (2.24) zwtat thbt 2‘”: X ZWIIPL where the summations for all times, t, are taken. The weighting factor is: —k.t -k',r w, =°—fif]—°— (2.25) to account for differences in weight of the measured data in terms of both time and magnitude. The exponential terms in the numerator deemphasize data taken near the end of each reaction. The product concentration term in the denominator is included to prevent the data from slower reacting components from being weighted too heat com; simu com: pseu cont systt be : regr the 11192 estir athi. dete; 24 heavily because of reaction product resulting from a faster reacting component. Nonlinear regression has also been applied for the treatment of simultaneous multicomponent kinetic data”. Assuming a two component system as described in equation 2.2 undergoes first- or pseudo-first order kinetics, and that the products are identical, the concentration of P at any time is given by equation 2.19. This system represents a nonlinear model. Data from such a system can be analyzed by any of a number of the common non-linear regression methodsZO. All of the non-linear regression methods have the same goal, which is to minimize the differences between the measured value for each point and the estimated value. Good initial estimates of the regression parameters are necessary in order to achieve accurate results. The Kalman filter has also been applied for the kinetic determination of mixtures“. The Kalman filter is a recursive filter and offers the advantages of linear least squares, but is simpler and more efficient. This filter is able to extract parameters from noisy data and to model complex systemszz. This technique is a major subject of this thesis and is described in much greater detail in Chapter III. Schechter7-3 developed a new error-compensating algorithm for kinetic determinations of dual-component mixtures without prior knowledge of reaction orders and rate constants. If a system described by equation 2.1 follows kinetics of unknown order, the rate of change of the two components are: 25 M = _kl[Cl]nt d‘ (2.26) d c ., 131.12.,” where n1 and n2 are the (unknown) reaction orders. The concentrations of C1 and C2 at time, t, are given by: [C] =[k,(n, -1),°'t+[C](1n.)](”“') l (2.27) [C] =[k,(n2 —1)t+[C ] (1- .,)]V( “2) If ax is the detector sensitivity for compound x, and the signal is linearly dependent on the concentrations, then the observed signal at time, t, is: St = ac,[C1], + arc2 [C2]t + ap’ [P1]. + orP2 [P2]t + X =[aCl —q,oz,,llC,]t +[ozC2 —q,oz,,2][C,]t +X (2.28) =fi1[Cl]t +fl2[cl]t +x where [31 is a constant proportionality factor, qi account for the reaction stoichiomeuies and X is the constant background. Equation 2.28 is then fit to the experimental data by using a modified Levenberg-Marquardt algorithm in order to solve the nonlinear least-squares problem. The method was tested using simulated data, and it was found to be applicable for the determination of concentrations, rate constants, and reaction orders in a wide range. However, there are some severe limitations and restrictions to the method, such as sensitivity of the detector to each of the components, 26 range of rate constants that are applicable, and the initial estimates of each of the parameters. Schechter and SchroderZ4 have developed an algorithm based on the Levenberg-Marquardt method for nonlinear least squares, which can be used for kinetic determinations in systems of mixed first- and second-order reactions. The development of the model for use in the algorithm is similar to that described above for systems of unknown order. The system can be described as follows: The algorithm was tested by simulated data, and the influence of the noise level, the range of the rate constants and other parameters were studied. It was found to be suitable for a wide range of parameters, and only a poor estimate of them was needed in order for the algorithm to converge. A method applicable to kinetic methods has been developed which estimates component amplitudes in multiexponential data25. This technique obtains quantitative information about individual component contributions to multiexponential data by a reiterative regression algorithm that employs a linear least-squares determination of component amplitudes within a nonlinear least- squares search for the exponential decay times. This is a unique application of a combination of the linear and nonlinear least squares methods described above. 27 b. Multiwavelength Methods As multiwavelength detectors such as photodiode arrays and charge-coupled devices (CCD) have become more popular, interest has grown in applying the power of these detectors to kinetic methods of analysis. As more of these devices are incorporated into analytical laboratories, it is expected that many new applications will be developed to take advantage of their simultaneous, multiwavlength capabilities. Multiwavelength methods were first developed in the 1970's using vidicon rapid scanning spectrometer326. The use of multiple wavelengths in kinetics experiments allows a greater diversity of systems to be successfully analyzed. Systems which have rate constants that are too close to analyze under normal means may be determined if there are spectral differences present. Factor analysis has also been applied to the resolution of simultaneous kinetic processe527»23»29,30:31. This method has the advantage of not requiring standards of the components to be determined or assumptions on the shape of their spectra. The data to be factor analyzed are required to be arranged in a matrix, D, the dimensions of which are NT x NW, where NT is the number of spectra acquired at different times, and NW is the number of wavelengths acquired for each spectrum. Each row of D corresponds to the spectrum of the mixture at a given point of its kinetics, and each colum matches the kinetic curve at a given wavelength. The algorithm for factor analysis will not be included, as it is well covered in the literature31,32. 28 A new procedure known as the kinetic wavelength-pair method has been recently reported33. This process involves measuring the difference in the rate of change of the absorbance with respect to time at two preset wavelength pairs. Since the rate of change of absorbance at any wavelength is dependent on the wavelength used, by measuring the rate of change of absorbance at two wavelength pairs, the concentrations of two species in a mixture can be determined using a single sample. For a mixture of two components, as in equation 2.1 above, application of the method involves measuring the rate of change of absorbance at two wavelength pairs, (11, 12) and (l3, 2.4), which are chosen so that the difference between the rates of change of the absorbance for the first wavelength pair, (21, 22) is as large as possible for one of the components (e.g., C1) and as small as possible for the other. The reverse should be true for the second wavelength pair. Data were collected with a diode-array spectrophotometer. c. Miscellaneous Methods In addition to the approaches described above, several other methods have been developed which are not easily classified in the categories already discussed. These so-called miscellaneous methods are described here. An H-point standard additions method (HPSAM) has been deveIOped for the ultraviolet-visible spectrosc0pic kinetic analysis of two-component systems3 4. The method is described under two cases: the analysis of two species of which only one evolves with 29 time, and the analysis of two species with overlapped time evolutions. The determination of the concentration of C1 by the HPSAM at equilibrium entails selecting two wavelengths M and 12 lying on each side of the absorption maximum of C2, so that the absorbance of the latter component is the same at both wavelengths. Then, known amounts of C1 are successively added to the mixture and the resulting absorbances are measured at the two wavelengths. The two straight lines thus obtained intersect at the so-called "H-point", (-C11, A11) where -C11=(-Cc1) is the unknown concentration of C1 and AH=(Ac2) is the analytical signal of C2. When one of the species evolves with time, the variables to be fixed are two times, t1 and t2. Species C2 which should not evolve with time over this range, should have constant absorbance. This contrasts with the equilibrium HPSAM, where two wavelengths are chosen. A plot is then prepared of added [C1] on the ordinate, and AA for the two times on the abscissa. The point where AA=O is equal to {Q}. When both species evolve with time, [C1] can be calculated by plotting AA for the two times against the added concentration of C1 at two wavelengths, A1 and 22, provided the absorbances of C2 at these two wavelengths are the same. The continuous addition of reagent technique is also applicable to multicomponent kinetic determinations”. This method is based on the continuous addition of reagent at a constant rate to the species to be determined36. The model for this method is developed here for a single component, but is applicable to multiple reacting species. Consider a reaction given by equation 2.1 where only one analyte is 30 present. If a solution of a reagent, R, with a concentration, [R]o is added at a constant rate, a, to a volume V0 of a solution containing the analyte C1, the overall rate of the process can be given by: —d|C,l=_[d q) _[d c,]] (2.29) dt dt (11 if C1 is the monitored species. The reaction rate can be expressed by: d c,] _ —[ dt 1...... _ k,[c,]R] (2.30) where k1 is the pseudo second-order rate constant. The reagent concentration increases with time according to: [R]: {—2— “IR—L— } (2.31) V +ut The term at, the volume of reagent added in time t, is also a dilution factor for the reagent. The dilution of the analyte C1 takes place according to: [c]: [c ,1, (V. V+ut] (2.32) Taking the derivative of equation 2.32 gives the rate of dilution: 9&1) =[“ )c] (2.33) dt W V+ut 31 By combining equations 2.30, 2.31, and 2.33, the overall rate can be shown to be: J1: =k[“_“‘L}1c,1+[ “ )[q] (2.34) V°+ut V +ut 0 In practice, equation 2.34 forms the basis of the continuous addition of reagent method using the initial-rate technique. Under short reaction times (at << V0), dilution of the analyte can be neglected and equation 2.34 can be rewritten as: d[C1]_ ill—KL ’T‘k( v0 )‘mL (2.35) Equation 2.35 is used for the construction of a linear calibration graph by plotting the initial slopes of the kinetic curves as a function of [C1]o over a fixed time interval in all experiments. C. Applications of Multicomponent Methods There have been several different approaches applied for the simultaneous kinetic resolution of mixtures in recent years. This section discusses the applications where multicomponent kinetic methods have been utilized for methods with small rate constant differences. 32 1. Graphical Extrapolation Methods The most popular graphical technique used for simultaneous kinetic determinations has been the logarithmic extrapolation method, which is described in section B above. This method has been applied to several problems in analytical chemistry. Table 2.1 summarizes some recent applications which have taken advantage of graphical extrapolation methods. As can be seen from Table 2.1, graphical extrapolation methods have performed well as long as there are sufficient differences in rate constants among the species being determined. The smallest ratio yielding accurate results reported in the references in the table was 7.313. Some of the mixtures determined by the reaction of aniline and its derivatives with DMPD37 had rate constant ratios of less than two. However, in these experiments, the total concentration of the two species in the mixture was determined, and not the individual concentrations of each species. These could not be successfully determined. The displacement reaction of the cobalt complex of pyridoxal thiosemicarbazone with DCTA does not occur, which allowed the analysis of mixtures of cobalt and nickel, and cobalt and copper to be successfully determined-i3. Since both nickel and copper complexes undergo the exchange reaction, their binary mixtures with cobalt can be analyzed by graphical methods, as the effective rate constant ratio is infinity. However, the rate constant ratio for the binary mixture of nickel and copper was 4.40, and could not be successfully 33 Table 2.1: Applications utilizing graphical extrapolation methods. Cont. Rate System Range Constant Errors Bet. m Alkaline earth- 106- 1:7.3:109 9.5% 13 DCTAl/Cd(ll) substitution 10‘5 M Mg:Ca:Sr ave. reaction Aniline and 104- * <14% 37 derivatives/DMPDi 10'3 M Fe & Co with pyridoxal 10-5 M 10.5 55% 39 thiosemicarbazone Ni & Cu with pyridoxal 11g mL'l NR 56% 40 thiosemicarbazone Displacement of Co, Ni & Cu. ug mL'1 coll $696 38 with pyridoxal thiosemicarbazone complexes with DCTA Displacement of In & 1.1g mL'l NR <10% 41 Ga/PARfi complexes with EDTAl‘” NR-not reported *-Only total concentration determined l-DCTA=(1,2,diaminocyclohexane, N, N, N ', N'-tetraacetic acid) LDMPD=N,N-dimethyl—p-phenylenediamine Tl-Co does not undergo the exchange reaction (kc(,=0) Il=¢~PAR=4-( 2-pyridylazo)resorcinol ill-EDTA=ethylenediamine tetraacetic acid 34 determined by this method. In order to determine this mixture, the single-point method was used. 2. Proportional Equations Based Methods The method of proportional equations described in section B above has also been used extensively in analytical chemistry. Table 2.2 summarizes some recent applications which have utilized the method of proportional equations. Several of these applications42,43.44’45:46,47,48.49.50 have utlized the technique of flow injection analysis to allow large numbers of samples to be processed in a short amount of time. Sampling rates of 60 h'1 or more have often been reported using this configuration. Several of the systems used also required the use of two procedures47»43:50’51.52. Each of the two procedures was optimized for the response of one of the components. In this manner, it is possible to determine more closely related species. Overall, the method of proportional equations can be quite useful. Its biggest drawback occurs when the species being determined are closely related kinetically. Under these conditions, the results obtained from this technique are not favorable, and it would be best to use a different method. Tal Zr 35 Applications utilizing the method of proportional Table 2.2: equations. Conc. Bate sttem Range Constant Errors Bet. Ratios Dissociation of Mg & Sr 104- 175 <10% 42 complexes of trans-1,2- 103 M diaminocyclohexanete tra- acetic acid Dissociation of 10-4- large’t <10% 43.44.45 metal/cryptand complexes 10‘3 M Dissociation of citrate 11g mL-1 43.4 <10% 15 complexes of Co & Ni Fe(III) & Mn(II) or Ni & Co ng mL'l 27.7(Fe,Mn) <4% 46.51 with ZOH-BATI 10’4 M 24.9(Ni,(b) <10% Z-component mixtures of 104- NR <6% 53 ammonia, hydrazine, and 10‘2 M hydroxylamine with 20H- BAAll Formation of molybdate ng ml:1 130 <5% 47.54 heterome acids of silicate 106- and phosphate 10'5 M Furfural & vanillin with p- ug ml‘l NR <6% 55 amincmhenol 36 Table 2.2 continued 9.0m. Bate sttem Range Constant Errors Ref. Eli—SIS histidine & 1- ug mL-l NR <6% 56 methylhistidine with o- phthaldialdehyde Fe(II) & Fe(III) with 11g mL-l large <5% 57 pyrocatechol violet Fe(ll) & Fe(III) by 1.1g mL'1 NR <10% 48 differential catalysis Co & Ni by ligand exchange ug mL'1 >100 <5% 49 with PSAA=l=‘~t and nitrilotriacetic acid Fe(II), Fe(III), and Ti(lV) 11g mL'1 NR <5% 50 with Tironl‘l“l MO(VI) & W(VI) with 10‘6 M NR <5% 52 DAPm and HLOZ Ethanol & Methanol with 10-6- NR <896 58 alcohol oxidase 10'5 M l-The second dissociation did not occur at room temperature, it was necessary to elevate the temperature after the first dissociation. L20H-BAT=2-hyroxybenzaldehyde thiosemicarbazone TT-ZOH-BAA=2-hydroxybenzaldehyde azine fi-PSAAa-2 -(S-brtrno Z-pyndylamFS-(Npqryl-Nsulfqrppylamim) aniline lfi-Tiron=Sodium 1,2-dihyroxybenzene-3,5 -disulfonate iii-DAP=2 ,4-diaminophenol dihyrochloride ra pr 91" tat DIE iss 4.: Stud med 11111 Used 37 3. Single-Point Method The single-point method, as described in section B above has also been used for several applications. Table 2.3 describes some of the recent studies that have utilized the single-point method. As is apparent from the table, this method is applicable to systems with smaller rate constant ratios than either the logarithmic extrapolation method or the method of proportional equations. However, once this ratio gets below approximately 5, the errors start to increase rapidly. Also, as is mentioned in section B, this method generally has a lower precision than does the logarithmic extrapolation method. The single-point method has been compared to the logarithmic extrapolation method using several of the systems in the tab1e38»39,40»59. Generally, it was found that the graphical method is preferable to the single-point method, unless the rate constant ratio is small. More modern techniques, such as computer methods, have been shown to be superior to either of these methods”. 4. Computer Methods As mentioned previously, the introduction of computers in the analytical laboratory has allowed a broader range of systems to be studied. Table 2.4 summarizes some recent applications of computer methods to kinetic methods of analysis. Also included in the table is which of the computer analysis methods described in section B was used for each application. 38 Table 2.3: Modern applications utilizing the single-point method. Cong. Me Sistem Bans: Constant Errors Rm 53.1115; £29; Fe & Co with pyridoxal 105 M 10.5 <5% 39 thiosemicarbazone Co, Ni, Cu/pyridoxal 11g mL'1 4.4 <5% 38 thiosemicarbazone ligand (Ni:Cu) exchange with DCTA Ni & Cu/pyridoxal 11g mL'1 large <5% 40 thiosemicarbazone Co & Cu and Co & Ni with pg mL'1 large <5% 60 TrADATl Secondary amines with 103- 18-232 <10% 61 carbon disulfide 10'1 M Cortisone & hydrocortisone 11g mL-1 1.8 most 59 with Blue Tetrazolium <30% l—TrADAT=3-( 1 'H-l ',2',4'-triazolyl-3 '-azo)-2,6-diaminotoluene 39 Table 2.4: Applications of computer methods to multicomponent kinetic determinations. hLdrochloride Cont. Bate 5351:2111 Range Constant Errors M21113 Ref. Ratios Alkaline earths with 105- 6.5-1660 <10% MPE& 18.19 CyDTA 10‘4 M (2) L13 <20% (3) Al(III) & Al-citrate 10‘5 M >100,000I <5% NLR 20 with Calcein Blue Metal-Zinconl 106- 12-183 <896 LIS 62 dissociation 10-5 M Epinephrine, roS- 2.4-20.3 <15% us 62 norepinephrine and L- 10‘4 M (> 15% Dopa with ascorbic acid for k= 2.4) H801) & Zn(II) with 10‘5 M 1.6 NR WIS 26 Zincon Amino acids with 105 M 2.5 <1496 KP 21 Tl‘initrobenzenesulfonic acid Alkaline phosphatase 10‘6- 4.6 15% KP 63 isozymes/guanidinium 10'5 M 4O Table 2.4 Continued Cont. Bale 5361:2111 Range Constant Errors Mound Ref. Ratios Cortisone & 1.1g mL'1 1.8 <8% KP 59 hyrocortisone with Blue Tetrazolium Co-EGTA & Ni-EGTA 1.1g mL-1 20 NR PA 31 complex displacement with PAR Co-EGTA & Ni-EGTA 1.1g mL-1 20 <10% NLR 64 complex displacement with PAR Alcohols with alcohol 103- 2.2-3.3 <10% KP 65 dehydrogenase- 10'6 M nicotinamide adenine dinucleotide Simulation studies varied varied varied NLR, 2123.241 L15, ,25,31, K12, 66,67 and PA l—Zincon=( 2-carboxy-2'-hydroxy—5 '-sulfoformazylbenzene) MPE=modified proportional equations NIS=non-linear least squares 'I.I.S=linear least squares Nl.R=non-linear regression WIS-weighted least squares KP=Kalman filter PA=Factor analysis 41 As is mentioned in Table 2.4, simulation studies were reported in several of the references. Simulations allow easy control of the "system" being studied. It is quite easy to change conditions on one, some, or all of the parameters of interest. In a real system, it is often difficult to change just one parameter without changing others. In most of the simulation studies listed in Table 2.4, many factors were varied. Among these are: rate constant ratios, concentration ratios, amount of noise present, and data acquisition rate. The errors found in the simulated systems varied greatly, often from O to >100% depending on the conditions chosen for each "experiment". Some of the computer methods also utlized detection at multiple wavelength526t31. This allows a broader range of systems to be studied. Since using data at multiple wavelengths allows resolution of signal due to both spectral and kinetic differences, these methods are applicable to even more systems. It is expected that more work in this area will appear shortly. S. Multiwavelength Methods Rapid-scanning spectrometers were the first devices used to perform multiwavelength studies for multicomponent kinetic analysis. Since then, diode array spectrometers have become popular tools for multiwavelength kinetic determinations. Table 2.5 summarizes the progress in multiwavelength methods to date. Table 2.5: Applications utilizing multiwavelength detection. 42 Cong. Rate Sststem Range Constant Errors Bet. m Hg(Il) & Zn(II) with Zincon 10—5 M 1.64 NR 26 Hemoglobin and 10'3 M NR <7% 68 Methemoglobin with cyanide and hexacyanoferrate(IH) Formaldehyde and acrolein 11g mL'1 2.9 <10% 33 with MBTHT Co & Ni EGTA complex 11g mL‘1 20 NR 31 displacement with PAR Co & Ni EGTA complex 11g mL'1 20 <10% 64 displacement with PAR l-MBTH=3-methylbenzothiazolin-2-one hydrazone 43 6. Miscellaneous Methods In addition to the methods for multicomponent kinetic determinations categorized above, there have been several applications which have not used one of the common techniques. These methods are discussed in section B above, and the applications which utilize these are summarized in Table 2.6. Although these methods have had few applications to date, they both seem quite promising. It is entirely possible that either or both of these techniques will become widely applied to multicomponent kinetic determinations. D. Conclusions The introduction of computers in the analytical laboratory has had a great impact in kinetic methods. This advance has allowed kinetic methods to become better automated, as well as the changes in the data analysis methods previously discussed. Kinetic methods are well suited to intelligent automationlo, since these techniques require precise timing and careful control over such reaction conditions as temperature, pH, ionic strength, and reagent concentrations. These advances have led to new data processing methods, such as the Kalman filter and non-linear regression methods. These techniques generally lead to more accurate and faster results than those that were used before the advent of computers. 44 Table 2.6: Miscellaneous applications. by the Laffé method F Cont. Rate sttem Range Constant Errors Metlm Bet. Ratios Zineb and maneb with 11g mL'1 5.4 <5% CAR 69 zincon Cu & Pe with pyridoxal 11g mL'1 4.67 f the parameters are used to calculate the estimated absorbance at he time the first measurement is to be made. The measurement is :hen made, and the difference between the actual and predicted absorbance values is calculated. This difference is then used to readjust the parameters and the weighting of the filter, and to provide an estimate of the absorbance of the system when the second measurement is made. The second absorbance measurement is then made, and the process is repeated for each of the remaining measurements. If a valid model of the system being studied is provided, the estimates of the absorbance will converge on the true (measured) absorbance. The example shown in Figure 3.1 is an ideal (noiseless) case. In reality, all measurements have some noise associated with them, and the differences between the actual and measured values should have a random distribution with zero mean. Kalman's original filter required the use of a linear model, however, modified versions of the filter have appeared in recent years. There are three main implementations of the filter in current use. These are: the linear (original) Kalman filter, the extended Kalman filterl, and the adaptive Kalman filter4. All versions of the Kalman filter use a "model" of the system, which is constructed of state variables as described below. Also used is a sequence of weighted measurements made on the system to obtain improved (filtered) estimates for the state variablesl. The original Kalman filter was developed for use with linear models. However, many of its applications have dealt with nonlinear 55 problems. In order to deal with a nonlinear system of state variables, the extended Kalman filter is used. The adaptive Kalman filter is used to compensate for the presence of model errors or outlying data points4. The idea is that data points which are inconsistent with the model are attributed to random noise in order that the data points not be used to corrupt the parameter estimates. Some of the applications of all three forms of the filter are described, followed by the algorithm for the extended Kalman filter, as this is the form of the filter that is most used in the current work. B. The Extended Kalman Filter Algorithm Although the Kalman filter was originally developed for use with linear models, it has been applied to several nonlinear problems, as described above. These problems often require the use of models with nonlinear system dynamics, nonlinear measurement models, or both. Systems of this type require the use of the extended Kalman filter, which is a modified version of the filter that allows the use of nonlinear systems. The algorithm for the extended Kalman filter is well covered in the literature, but is included here for the sake of completeness. First, we consider a time-variant process, called the system, from which noisy measurements are obtained at discrete time 'ntervals, t. Like other filters, the process does not need to be time 'ariant, but must consist of measurements made at discrete 1 tervals. For example, the noisy measurements could be made at iscrete wavelengths. The extended Kalman filter is an algorithm 56 which recursively estimates state parameters from the series of noisy measurements in agreement with a system model. The parameters of this model, which are unknown prior to filtering, are the state parameters, and collectively form the state vector, which is not required to be static. One way in which the Kalman filter differs from ordinary digital filters is that it is a state space methods. This means that the objective of the algorithm is to obtain the best estimate of the state vector, and therefore, of the true signal. In order to do this, we must first define a measurement model for any point in time of the form: Z,=h(X,)+v, (3.1) In this model, Zt is an m x 1 vector which consists of m measurements made at interval t (i.e., for t=1,2,...). The n x 1 state vector is given by Xt, where n is the number of state parameters. The 111 x n observation matrix is given by h(), and provides the connection between Z: and K: under ideal (noiseless) conditions. Finally, V: is an m x 1 vector which describes the noise in the measurement. This term is assumed to be a white noise sequence. In addition to the measurement model, the extended Kalman filter requires a model to define the propagation of the state vector between discrete measurement intervals. This is known as the system dynamics equation, and is given by: X‘. = LX, +G,w, (3.2) 57 where the "U” subscript refers to the predicted values for the state vector, ft is the n x n system dynamics matrix, which describes how the state vector Xt, propagates between measurements. G: is an n x 11 matrix which describes the system noise and wt is an n x 1 white noise vector which accounts for the noise in the state vector (parameter estimates). Equation 3.2 is often simplified when the state vector is time invariant (ft=I, where I is the identity matrix) and considered to be noise free (Gtwt=0). Once the models given by equations 3.1 and 3.2 have been defined, the extended Kalman filter may be applied. However, since the models are non-linear in nature, their form must be adapted for use in the filter. This is done by means of a Taylor series expansion for the observations matrix and the state transition matrix, with the expansion being truncated after the first term. These adjusted models may then be used in the extended Kalman filter as follows. First, the Kalman gain is calculated by: Kt = PVH,1[H,P'.H,T + M" (3.3) In the notation used here, the "T" superscript refers to the transpose of a matrix and the "t" subscript indicates that this is the best estimate of a parameter prior to including the new measurement. K: is the n x m Kalman gain matrix which describes how the new measurement is to be used to update the state vector estimate. The n x 11 matrix, Pt- is the current estimate of the error covariance matrix for the state parameters. The diagonal elements of this matrix contain the variances of each of the elements of the state 58 vector. Usually, since we have no a priori knowledge of the values of the state vector, the diagonal elements of the error covariance matrix are taken to be large, and the off-diagonal elements are taken to be zero. Since this matrix is updated by the extended Kalman filter algorithm, only the initial value needs to be provided. R: is the covariance matrix for the measurement noise element vt, and is of the size m x m. The diagonal elements of this matrix are taken to be the variances of the associated measurements, 2;, and the off- diagonal elements are taken to be zero. One the Kalman gain has been calculated, the next step is to include the measurement at interval t to update the estimate of the state vector. The correction factor between the actual measurement vector Zt and the predicted measurement vector h(Xt) is weighted by the Kalman gain, as shown in equation 3.4: x, = X.— + K,[z, - h(x,)] (3.4) Similarly, the Kalman gain is used to estimate the updated covariance matrix, Pt: P, = (I-K,H,)Pt-(I—K,H,)T +K,R,K,T (3.5) Where I is the n x 11 identity matrix. Ht is the m x n observation matrix. This is the linearized form of h(Xt) obtained from the Taylor series expansion described above. Simpler forms of equation 3.5 are sometimes given, but are rarely used in practice because of poor numerical stability5. pro usl usi “9‘ Dr 59 The next step of the extended Kalman filter algorithm is to project the estimate of Xt ahead to the next measurement interval using the system dynamics matrix ft. x‘, = r,x, (3.6) The final step of the algorithm is to project ahead the estimate of Pt using F(Xt), which is the linearized form of the system dynamics matrix using the Taylor series expansion described above. This update is given by: P‘, = F(X,)P,F(x,)T +0, (3.7) where Q is the covariance matrix for the white noise sequence Gtwt used in the system dynamics descriptor given in equation 3.2. It is often taken to be zero. Equations 3.3 through 3.7 constitute the algorithm for the extended Kalman filter. This is the exact form of the algorithm used for the work presented in this thesis. For each of the systems which are studied in this work, the measurement model, linearized measurement function, and state vector are given in the sections describing the systems under study. In all the work described in this thesis, the state vector is set to be a static entity. This simplifies the algorithm by making equations 3.6 and 3.7 equal to their previous values (since Q; is taken as zero). 60 C. Applications All three versions of the Kalman filter described have found many applications in analytical chemistry. This section describes some of the uses each form has found. 1. Linear Kalman Filter The linear version of the Kalman filter has been used to a much greater extent than either of the other two forms of the filter. Some of the uses of the Kalman filter in chemistry have recently been reviewed1.6»7. The Kalman filter has been applied to the resolution of overlapped chromatographic responses3’9,10,11.111314, and for chromatographic optimizationl 5. Resolving overlapped responses, such as poorly resolved chromatographic techniques, requires the detection of more than one response parameter. This is accomplished either with several different types of detectors, or with multiwavelength array type detectors, such as photodiode arrays or charge-coupled devices (CCD). Another area where the Kalman filter has been applied using data from more than one wavelength is in multicomponent analysis of mixtures at equilibrium. This has been done utilizing both ultraviolet-visible spectrophotometry16»17.13»19:20, and fluorimetryzLZZ. These applications do not require the use of array- type detectors or multiple detectors, however, since the measurements are performed at equilibrium, and the time evolution of a species is not limiting. 61 The linear Kalman filter has also been applied to electrochemical processes. It has been used for anodic stripping voltammetry23, for resolving overlapped electrochemical peaksZ4 and square-wave voltammetric responses-25.26,, for estimating electrochemical charge transfer parameters”, and for evaluating the current-time function in dc polarography2 3. The linear filter has also been applied tp the analysis of enzyme kinetic data29»30»31, and for kinetic determinations of two- component mixtures32.33. The linear Kalman filter has also been used in inductively coupled plasma-atomic emission spectrometry (ICP-AES) for the determination of spectral interferents34,35, for data reduction36, and for the determination of trace elements37. The reliability of Kalman filtering results in ICP-AES has been evaluated”. The filter has also proven well suited for the compensation of drifts in several systems39»40.41»42»43. Other areas in which the linear form of the filter have been applied include: pulsed photoacoustic spectrsocopy for the determination of metal complexation parameters“, setpoint control in continuous titrations45, the planning of sampling“), beam modulation for alpha-particle diagnostics47, for "removal" of interferences in complex sample matrices43, for modeling of chemical response surfaces49, in circular dichroism spectroscopy”, and in x-ray fluorescence spec trometry5 1. As can be seen from the great variety of applications. the linear Kalman filter is quite adaptable and widely used. 62 2 . Extended Kalman Filter The extended Kalman filter has not been as extensively utilized as the linear form. This is probably due to the greater complexity of the modeling required, and the fact that the linear form of the filter is sufficient for many applications. Also, the extended filter is less stable than the linear form. There is more dependence on the initial estimates provided to the filter, and if care is not taken in choosing these estimates, convergence of the estimates is not always obtained. There are two other reasons that this form of the filter has not been as widely applied. The first is that the computations resulting from using this form of the filter are more complex, and therefore more time consuming, and the second is that the filter is not necessarily optimal in the least squares sense. Greater care must be taken using this form of the filter, and this has probably been one of the main reasons it has not been applied more often. However, its demonstrated applicability to several non-linear problems with excellent results will probably cause this method to become more popular- The extended Kalman filter has been applied to optimize simulations for step voltammetry52. The extended Kalman filter has been compared to the simplex and Marquardt procedures in the estimation of parameters from voltammetric curvesS 3. The filter was found to be very favorable to the other two methods. The extended Kalman filter has also been used in kinetic methods of analysis. First-order kinetic parameters of one- component systems have been estimated54.55. Changes in the rate 63 constant during a reaction due to changes in temperature have been modeled and corrected56-57. The extended Kalman filter has also been applied to the estimation of ester hydrolysis parametersSS. 3. Adaptive Kalman Filter The adaptive form of the Kalman filter has been applied to systems where there are errors in the model. When using the adaptive Kalman filter, it is necessary that the model be correct for a portion of the response, and the errors must be attributable to a single component. This method has been applied to compensate for model errors in multicomponent spec trometry5 9. The adaptive filter has also been combined with simplex optimization60»61. In this combined method, the initial guesses for the filter are generated by the simplex algorithm. The adaptive Kalman filter has been employed in thin-layer chromatography (TLC) for background subtraction62»63, and for resolution of severely overlapped responses in TLC64. The adaptive filter has also been compared to other robust regression methods in the analysis of linear calibration data“. The adaptive filter was found to be the best suited for the detection of outliers in small calibration data sets. This method has also been used to find a previously overlooked contribution to gas-liquid partition coefficients from excess-solute-induced solvent reorganization66 haw Qher. it is . [he 1 item the t 64 4. Kalman Filter Networks The most recent implementation of Kalman filtering has been in the area of Kalman filter networks. This techniques employs a set of parallel models through which each data point is passed67. The performance of each of the models is then evaluated, and the best model is selected. This method has the advantage of being a recursive implementation of principal components analysis. To date, there have been only two applications of this method. The first uses this technique for peak purity analysis66. The second uses a set of 40 or more parallel filters for the correction of errors from between-sample variations in the pseudo-first-order rate constant for kinetic methods for a one-component system“. Although there have been very few implementations of this filtering scheme, the technique is quite new, and should be widely applicable to a variety of chemical problems. All three forms of the Kalman filter should work quite well in this implementation as well. D. Conclusions The three forms of the Kalman filter described in this chapter have found an ever growing list of applications in analytical chemistry. Each of the three forms has its own applications to which it is most suited. The extended Kalman filter appears well-suited to the task of multicomponent kinetic determinations. All previous methods of multicomponent kinetic determinations have assumed the rate constant to be invariant from run to run. However, as 65 described in Chapter II, this is seldom the case. Conditions can be held such that the variations in the rate constant are minimized, but these variations will always be present to some extent. By using the extended Kalman filter though, the rate constants can be included as modeled parameters, and the filter can help to compensate for their between run variations. There are a number of advantages in using the Kalman filter and its various forms for the analysis of chemical data1. The first is the possibility of including time in the model, which permits the modeling of drift, reaction kinetics, and of time-dependent inputs (e.g., titrations) in a straightforward manner. A second advantage is that the filters are simple and fast enough to be used on most analytical instruments. Finally, they are applicable to most analytical methods. As long as it is possible to mathematically model the processes occurring, the filter can be applied. Use of adaptive filtering allows unmodeled processes to be compensated for. This chapter discusses the Kalman filter and the modifications to Kalman's original linear filter algorithm that have appeared since its introduction. Applications of the various forms of the filter to chemical problems are discussed. LIST OF REFERENCES LIST OF REFERENCES 1 S.D. Brown, Anal. Chim. Acta, 1986, 181, 1-26. 2 R.E. Kalman, Trans. ASME: J. Basic Eng., 1960, 82, 35-45. 3 P. Eykhoff, system Identification, Wiley, New York, 1974. 4 S.C. Rutan, Anal. Chem., 1991, 63, 1103A-1109A. 5 PD. Wentzell, Ph.D. Dissertation, Michigan State University, 1987. 6 S.D. Brown, T.Q. Barker, R.J. Iarivee, S.L. Monfre, and H.R. Wilk, Anal. Chem., 1988, 60, 252R-273R. 7 S.D. Brown, R.S. Bear, Jr., and TB. Blank, Anal. Chem., 1992, 64, 22R-49R. 8 Y. Hayashi, T. Shibazaki, and M. Uchiyama, Anal. Chim. Acta, 1987, 201, 185-191. 9 Y. Hayashi, T. Shibazaki, R. Matsuda, and M. Uchiyama, Anal. Chim. Acta, 1987, 202, 187-197. 10 Y. Hayashi, R. Matsuda, S. Yoshioka, and Y. Takeda, Anal. Chim. Acta, 1988, 209, 45-55. 11 Y. Hayashi, S. Yoshioka, and Y. Takeda, Anal. Chim. Acta, 1988, 212, 81-94. 12 T.Q, Barker and S.D. Brown, Anal. Chim. Acta, 1989, 225, 53-68. 13 TL. Cecil, R.B. Poe, and S.C. Rutan, Anal. Chim. Acta, 1991, 250, 37- 44. 14 M. Redmond, S.D. Brown, and H.R. Wilk, Anal. Lett., 1989, 22, 963- 979. 66 67 15 Y. Hayashi and R. Matsuda, Anal. Chim. Acta, 1989, 222, 313-322. 16 H.N.J. Poulisse, Anal. Chim. Acta, 1979, 112, 361-374. 17 A. van Loosbroek, H.J.G. Debets, and P.M.J. Coenegracht, Anal. Lett., 1984, 173, 779-792. 13 A. van Loosbroek, H.J.G. Debets, and DA. Doornbos, Anal. Lett., 1984, 17A, 677-688. m 19 Yi-Ming Liu and Ru-an Yu, Talanta, 1988, 35, 707-711. 20 I. Lavagnini, P. Pastore, and F. Magno, Anal. Chim. Acta, 1990, 23 9, 95-106. 21 TL Cecil and S.C. Rutan, Anal. Chem., 1990, 62, 1998-2004. 22 L. Jianjun, Z. Yun'e, and C. Guanquan, Talanta, 1990, 37, 809-813. 23 P.F. Seelig and H.N. Blount, Anal. Chem., 1976, 48, 252-258. 24 T.F. Brown and S.D. Brown, Anal. Chem., 1981, 53, 1410-1417. 25 CA. Scolari and S.D. Brown, Anal. Chim. Acta, 1984, 166, 253-260. 26 CA. Scolari and S.D. Brown, Anal. Chim. Acta, 1985, 1 78, 239-246. 27 T.F. Brown, D.M. Caster, and S.D. Brown, Anal. Chem., 1984, 56, 1214-1221. 23 M. Bos, Talanta, 1986, 33, 583-586. 29 S.C. Rutan and S.D. Brown, Anal. Chim. Acta, 1985, 1 75, 219-229. 30 W.H. Lewis, Jr. and S.C. Rutan, Anal. Chem., 1991, 63, 627-629. 31 E. FOrster, M. Silva, M. Otto, and D. Pérez-Bendito, Anal. Chim. Acta, 1993, 274, 109-116. 32 P.D. Wentzell, M.L Karayannis, and S.R. Crouch, Anal. Chim. Acta, 1989, 224, 263-274. 68 33 R. Xiong, A. Velasco, M. Silva, and D. Pérez-Bendito, Anal. Chim. Acta, 1991, 251, 313-319. 34 E.H. van Veen and M.T.C. de Loos-Vollebregt, Spectrochim. Acta, 1990, 458, 313-328. 35 EB. van Veen, F.J. Oukes, and M.T.C. de Loos-Vollebregt, Spectrochim. Acta, 1990, 458, 1109-1120. 36 EB. van Veen and M.T.C. de Loos-Vollebregt, Anal. Chem., 1991, é 63, 1441-1448. .‘ 37 E.H. van Veen, M.T.C. de Loos-Vollebregt, A.P. Wasskink, and H. Kalter, Anal. Chem., 1992, 64, 1643-1649. 33 J. Yang, Z. Piao, X. Zeng, Z. Zhang, and X. Chen, Spectrochirn. Acta, 1992, 473, 1055-1064. 39 H.N.J. Poulisse and P. Engelen, Anal. Lett., 1980, 13, 1211-1234. 40 P.C. Thijssen, S.M. Wolfrum, and G. Kateman, Anal. Chim. Acta, 1984, 156, 87-101. 41 P.C. Thijssen, Anal. Chim. Acta, 1984, 162, 253-262. 42 P.C. Thijssen, G. Kateman, and H.C. Smit, Anal. Chim. Acta, 1985, 1 73, 265-272. 43 P.C. Thijssen, L.T.M. Prop, G. Kateman, and H.C. Smit, Anal. Chim. Acta, 1985, 174, 27-40. 44 S.C. Rutan and S.D. Brown, Anal. Chem., 1983, 55, 1707-1710. 45 P.C. Thijssen, N.J.M.L. Janssen, G. Kateman, and H.C. Smit, Anal. Chim. Acta, 1985, 1 77, 57-69. 46 F.W. Pijpers, Anal. Chim. Acta, 1986, 190, 79-87. 47 W.S. Cooper, Rev. Sci. Instrum, 1986, 57, 1846-1848. 48 Yi-Ming Liu and Qu-Qin Yu, Analyst, 1987, 112, 1135-1137. 491 50 Act 51 ' 19‘. 52] 53 22. S45 69 49 P.D. Wentzell, A.P. Wade, and S.R. Crouch, Anal. Chem., 1988, 60, 905-91 1. 50 A.F. Guillem, D.D. Shillady, L.F. Jones, and S.C. Rutan, Anal. Chim. Acta, 1991, 246, 1-8. 51 Z. Li, R. Yu, L. Shi, J. Xu, M. Zhang, and Q. Wang, Anal. Chim. Acta, 1991, 248, 257-261. 52 RS. Bear, Jr., and S.D. Brown, Anal. Chem., 1993, 65, 1061-1068. 53 I. Lavagnini, P. Pastore, and F. Magno, Anal. Chim. Acta, 1989, 223, 193-204. 54 S.C. Rutan and S.D. Brown, Anal. Chim. Acta, 1985, 167, 23-37. 55 S.C. Rutan, C.P. Fitzpatrick, J.W. Skoug, W.E. Weiser, and H.L. Pardue, Anal. Chim. Acta, 1989, 224, 243-261. 56 CA. Corcoran and S.C. Rutan, Anal. Chem., 1988, 60, 1146-1153. 57 CA. Corcoran and S.C. Rutan, Anal. Chem., 1988, 60, 2450-2454. 53 S.L. Monfre and S.D. Brown, Anal. Chim. Acta, 1987, 200, 397-410. 59 S.C. Rutan and S.D. Brown, Anal. Chim. Acta, 1984, 160, 99-119. 60 S.C. Rutan and S.D. Brown, Anal. Chim. Acta, 1985, 167, 39-50. 61 H.R. Wilk and S.D. Brown, Anal. Chim. Acta, 1989, 225, 37-52. 62 DD. Gerow and S.C. Rutan, Anal. Chim. Acta, 1986, 184, 53-64. 63 DD. Gerow and S.C. Rutan, Anal. Chem., 1988, 60, 847-852. 64 S.C. Rutan and CB. Motley, Anal. Chem., 1987, 5 9, 2045-2050. 65 S.C. Rutan and P.W. Carr, Anal. Chim. Acta, 1988, 215, 131-142. 66 S.C. Rutan, P.W. Carr, and R.W. Taft, J. Phys. Chem., 1989, 93, 4292-4297. 67- g. 2519 68 p} 181. 70 67 SJ. Vanslyke and PD. Wentzell, Anal. Chem., 1991, 63, 2512- 2519. 68 RD. Wentzell and SJ. Vanslyke, Anal. Chim. Acta, 1992, 25 7, 173- 181. CHAPTER N Simulation Studies of Two-Component Systems A series of simulations was used to test the applicability of the extended Kalman filter for multicomponent kinetic determinations. Simulations offer a number of distinct advantages when testing a data analysis technique. These advantages include the ability to adjust parameters which might affect the results of the filter, the ability to adjust parameters to any value desired without need to search for a chemical system which may be hard to control, and the ability to vary each parameter independently. This allows precise knowledge and control of all of the inputs to the Kalman filter, which can let one develop a greater understanding of the advantages and limitations of the method. A. Kinetic Model 1. Two-Component Mixture Following Irreversible Kinetics A typical system which could be used for multicomponent kinetic determinations is described below. For a two-component mixture following irreversible kinetics, the system can be described as: 71 when denot tone the r simpl there 1'1 =l Unde 111111 equal and [ 110m 72 Cl +R-—"I—9P, C2 + R—"L—->P2 (4° 1) where R is a common reagent reacting with two similar species, denoted by C1 and C2 to form two similar, but not identical products P1 and P2. If the common reagent is held in excess such that its concentration does not change appreciably throughout the course of the reaction and pseudo-first-order kinetics apply, equation 4.1 simplifies to: C,—£i—->Pl (4.2) C2—"1—>P2 where k'1 and k'z are the pseudo-first-order rate constants and: k'r =k1[R]t and k'z =k2[R]t Under these conditions, the concentrations of C1, and C2 are given by: c = C 6““ [ 1]: 1 110 k. (4.3) [C2], =[C2LC- ’1 If [Pl]o-[Pz]o=0, the concentrations of products at time t are given by equation 4.4: P = C — C [I]. [11. [11. (44) [P2]\ = [C210 -[C2]: and [Rh can be obtained by subtracting the concentration of products from the initial reagent concentration as shown in equation 4.5: If\ at the Sim! 73 [R]. =[R].-([C1].-[C1].)-([C2].-[C2].) (4.5) If molecular absorption measurements are used to follow the reaction, and the absorbances of all species are additive (i.e., Beer's law is obeyed), then the absorbance of the system described in equation 4.2 per unit pathlength at any time t, and any wavelength j, is given by: Alt = 81°C. [Cl]: '1' Sic; [C2 1‘ '1' 5511115]: ‘1’ ejP. [P1]! '1' 5,1),[1’21‘ (4-6) If we now take the relations for the concentration of each species as a function of time given in equations 4.3-4.5 above, and substitute them into equation 4.6, followed by collection of terms and simplification, we obtain the measurement function for this system: A), = (55c, + 5m - £11,! )[C,]¢e"‘it + (sic2 + $5,, - sip2 XC2]°e“9‘ (4.7) +(ejp‘ - 85R)[C1]° + (raj,2 — £3,1ch + 85,,[R]o The measurement function described in equation 4.7 is now expressed as a function of five parameters, [C 110, [C 210, [R]o, k'1, and k'z. By using equation 4.7 as a model for a first-order or pseudo- first-order chemical reaction, information about the unknown concentrations can be obtained, even if the spectra of the reactants and products are severely overlapped. Also, because measurements are taken at several wavelengths, a greater amount of information is 74 available to the filter so that the reaction need not be monitored to completion. 2. Use of the Model in the Extended Kalman Filter The extended Kalman filter requires a system to be defined by two equations. The first is the measurement model, and the second is the linearized measurement function. The linearized measurement function relates the quantities in the state vector to measurable quantities. For example, with the two-component non-reversible system following pseudo-first-order kinetics defined by equation 4.1, and the measurement model defined by equation 4.7, then the state vector for this case would be a 5x1 vector and is given by: "[C1]. [02]. (4.8) As each new data point is obtained and passed to the filter, the estimates for each of the elements in the state vector are updated and used as the initial estimates for the next loop of the filter. The kinetic model used to derive equation 4.7 assumes that the reactions behave independently and there are no synergistic effects. The extended Kalman filter also requires a function which relates the quantities in the state vector to a measurable quantity, which in this case would be absorbance. This is done by means of 75 the linearized measurement function, which was described in detail in Chapter III. For the measurement model defined in equation 4.7, and the state vector defined in equation 4.8, the linearized measurement function is a 1x5 vector, and is given in equation 4.9 below: -k'r 11 (55C. + Era " 55p. ' + (£59. ' 53R) 4!: 12 (sic, +5111 75% )e 2 +(3jP. ‘53:!) 1 £51: (4.9) _ -k't 1‘ - —t(£jct + £113 - £191 ){C1loc 1 _ -k': 15 - -t(£lcz + elk ’ 55?; )[CZLC 2 0’ ll I‘VE-'1‘!" The linearized measurement function is simply a Taylor's series expansion about the mean of the measurement model with respect to each component in the state vector and with truncation after the linear term. In the cases described here, this amounts to the first derivative of the measurement function with respect to each component in the state vector. In an ideal reaction, each product would absorb at a unique wavelength, there would be no spectral overlap, and the reactants would not absorb. A more realistic case, however, is where most species absorb to some extent, and there is at least some spectral overlap. 76 B. Experimental 1. Simulated Data Simulated kinetic data were used to test the filter under a variety of conditions which would be difficult to obtain experimentally. Exact knowledge of the reaction conditions also allows the user to change one or more parameters at a time to determine their influence on the filter, and thereby provide a better understanding of its advantages and limitations. Data were generated using a full second-order kinetic model with the common reagent, R, in excess. The filter assumes that the data follow a first- order or a pseudo-first order model. The second order system was used for the data to avoid biasing the filter by using the same model for data analysis as for data generation. Because the filter also estimates the concentration of [R]o, small variations in its concentration are tolerated. Random noise was added to all synthetic data; the algorithm used to generate the noise is discussed below. All studies presented contained 100 data points unless otherwise specified. It was also possible to vary the number of wavelengths that were used for the simulations. In the cases discussed in this chapter, this value was varied from only one wavelength being monitored to 15 wavelengths. Studies were performed on systems containing anywhere from one to eight separate reacting components. The derivation for the measurement model for reactions with n separate components is given in part A above. 77 Since the extended Kalman filter is a recursive algorithm, it is necessary to provide initial estimates of the state parameters. Some knowledge of the rate constant is presumed, and an approximate value is provided to the filter as an initial estimate. The initial values of the rate constants input to the filter were varied by as much as 5096 in the simulations, and the results of this are presented later. Also, the initial concentration of the common reagent, R, should be known, so this value was provided to the algorithm with no error. The concentrations of the reactants are assumed to be unknown, and are the which are desired from the filter. In order to avoid bias, the initial estimates provided to the filter for reactant concentrations were zero in all cases. While providing the expected values of the rate constants and the common reagent concentration to the filter is not strictly necessary, it does help to speed up the filtering process. Acceptable values for all parameters can usually be determined with no prior knowledge of any of the parameters, but multiple passes through the algorithm are necessary. 2. Generation of Random Noise Random noise with a Gaussian distribution was added to the pure (noiseless) absorbance values calculated according to the full second-order kinetic model. In order to add this noise to the pure absorbance values, it is first necessary to generate Gaussian noise. In Microsoft QuickBASIC, a function called RND is available. This function generates random numbers between 0 and 1 with a uniform dish rant ranc‘ tot] dot for I give ZOOt lfldlt pass 3.1] disse 111104 were Dreci hlth mean incluc linear 78 distribution. The method used to obtain normally distributed random variables using the RND function is based on obtaining a random number which is used as the argument in an approximation to the inverse of the standardized cumulative normal distribution function1»2:3. The noise generated has zero mean, and the standard deviation can be defined by the user. Figure 4.1 shows a histogram for the added noise. The theoretical standard deviation of the noise given to the algorithm used to generate the noise was 0.005 and 20000 points were determined. The skewness and kurtosis both indicate that the added noise is Gaussian in nature, and the data also pass the chi-square test for a Gaussian distribution. 3. Data Processing Unless otherwise stated, all computations described in this dissertation were carried out on an IBM-compatible 286-type microcomputer equipped with a math co-processor. All programs were written in Microsoft QuickBASIC (version 4.0) using double precision arithmetic. The functional form of the extended Kalman filter used for these studies is described in Chapter III of this work, with the appropriate changes in the system dynamics and measurement function equations made as necessary. These changes include incorporation of the appropriate measurement function and linearized measurement function as described above. $25300 us Flgul- them 79 2000 1500 1000 # counts 500 0 -0.02 -0.01 0 0.01 0.02 range Figure 4.1: Histogram for added noise. 20000 total points, theoretical standard deviation=0.005. simuh' noise) norma‘ both si Appen numbe of wav COmp( filter siStet We used throu the cc in Tat WEre ( 11an 80 The listing of the program NRNDAT, which was used to simulate the data for irreversible reactions is included in Appendix A. This includes the algorithm for the generation of the pure (no noise) absorbance values, and the algorithm for generation of normally distributed random noise described in part B.2. above. The listing of the program EKFNRN, which was used to process both simulated and real data for irreversible reactions is included in Appendix B. This program is applicable to a system containing any number of irreversible reactions which are monitored at any number of wavelengths. C. Results and Discussion The great majority of simulations were done on the Mo- component system. Using this system is a compromise because the filter is quite capable of accurate results using more complex systems. However, the smaller number of variables makes this system more amenable to rigorous study. The concentrations, rate constants, and molar absorptivities used are given in Table 4.1. These values were held constant throughout all simulations unless otherwise noted. In cases where the concentrations were varied, the starting values were those given in Table 4.1. For all studies discussed in this section, one hundred points were collected over a constant total elapsed time period. Different numbers of wavelengths were used, but data were always generated for each wavelength at each point in time. This was done because 81 Table 4.1: Values of molar absorptivity times pathlength (L mol’1) used for simulations unless otherwise noted. [C1]o=0.0001 M, [c2]o=o.ooor M, [R]o=0.01 M, k1=50 3'1,k2=10-50 s-1 20 # €c1 €c2 8R 8P1 8P2 1 0 0 10 5000 2500 2 O 0 5 2500 5000 la the the (On 11102 11101 am 82 the diode array spectrometer system used for experimental verification of the simulation studies is capable of storing only 100 scans. Therefore, in order to match experimental capabilities, the simulated data were limited to only 100 points. However, the effect of changing the data interval, and the total number of points was investigated, and it was found that the most important factor affecting the accuracy of the filter estimates was not the number of data points, or the data interval (although at extremely short intervals oversampling became problematic), but was the percent completion to which all reactions were sampled. 1. Error vs. Percent Completion Throughout the course of the simulations, it became apparent that one of the main factors affecting the accuracy of the estimates was the extent of the reaction that was monitored. Figure 4.2 shows the error in estimated concentration versus percent completion of the reaction. Although these results were obtained with only one component, the same trend holds true for each component in a multicomponent mixture. As can be seen, the estimates become more accurate as the extent of reaction increases. This increased accuracy is primarily a result of the time required for the Kalman filter to converge as can be seen in Figure 4.2. For analytical purposes it would be desirable to obtain results as rapidly as possible and thus to use a small extent of reaction. In order to obtain accurate results (<596 error), it was generally necessary to follow each reaction to approximately 50—60% completion. This is in agreement Estf] 83 1.310" 1.210'4 1.110'4 Est[] 1010‘4 9.010'5 8.010'5 7.010'5 0 2 0 40 60 8 0 1 00 % completion Figure 4.2: Error in estimated concentration of [C110 vs. 96 completion of the reaction. [C1]o=1x104 M, k=100 sec-1 84 with results reported by Mieling and Pardue4 for the predictive kinetic approach. These workers found that it was necessary to follow a reaction for one to two half-lives (SO-75% complete) in order to predict accurately the final equilibrium position. However, in competing parallel reactions, if the faster reacting species is greater than 50% reacted, while the slower species is not, the estimates from the first reaction can be accurate, while those from the slower reaction are not. (The estimates from the faster reaction are not, of course, as accurate as they would be in the absence of the slower reaction.) The accuracy of the estimated concentrations is fairly independent of other species in the mixture, as long as there are sufficient spectral differences. The effect of varying the concentrations of the reactants was investigated for a two component mixture. This study was done in tandem with another in which the ratio of the two rate constants was varied. For each rate constant ratio and concentration ratio, the simulation was repeated nine times, and the average error was obtained. The results of these studies appear in Figures 4.3-4.6. Figures 4.3 and 4.4 show the average error in the estimated concentrations of C1 and C2 when C2 is held constant, while Figures 4.5 and 4.6 show the average error when C1 is held constant. Note that in Figure 4.3 that for all values tested, the largest error in estimated concentration was still less than 5%. Figures 4.4 and 4.6 show that even though the errors in the estimated concentration of C2 are greater, the largest errors come at rate constant ratios corresponding to a small percent completion of the reaction of C2 and at concentration ratios where C2 is a much small value than C1 . Also, 85 Figure 4.3: Error in estimated C1 as a function of concentration ratio and rate constant ratio. C2 constant for all experiments. Flgu; - ’r,’/¢ , . <39“ ‘-” o‘é’ -: away-€7.- - 1‘" 7 /" ““5593”; ’ ” 0 5 0.1 x‘ . «7;;‘1”,“5':-?‘ ’ o ” 0 9 Figllre 4. 4: Error in estimated C2 as a function of concentration ratio and rate constant ratio. C2 constant for all experiments. 87 4.5 4 3.5 3 E 2.5 "Of "1’ / 1.5 :1. l’ 1 i ’7 ’ ; // . 0.5 - "mm“; A , . - . grieve-stadtllit/x ‘° 0' s’.‘ . . ,Q ) I O ‘ ' f/"A ‘ "’ "4»? 41“?»2'1‘. If A9- I - ‘x’;’, A \ Figure 4.5: Error in estimated C1 as a function of concentration ratio and rate constant ratio. C1 constant for all experiments. F12 88 i 16- E -16 I I '1: 12 a | P 104 ”1“”, 1 "'10 Error 8" ’1‘“, 1 n —8 €2,036) 6- 1111‘ 111-» "6 /‘/ A “1 LI, JgI/9},f\ ’4 4 t t 10‘ l; )1) 1'4 ‘r 7,, "I“ ’1' “ 4 ".1 (13‘1“ (/ ’I/ 7"I/£fi//,,j//’f’/,’ ”2 2 I .1. ‘lt.“’/’///////////,,,/ to 0 .10710‘s“~‘-!es.7:-.'-,-;.;;;:522 r 6 5 I i , 7)" \‘Vli’: :{Cgo 4'5 43.5 “41134;; 0-7 2 32.5 2 0.3 c2 /c1 1.5 1 Figure 4.6: Error in estimated C2 as a function of concentration ratio and rate constant ratio. C1 constant for all experiments. re re de Vat as: We 011 C0: the 89 in the great majority of experiments, the error was less than 5%. It Seems quite likely that if the reaction were followed for a longer time in the cases where errors were the highest, accuracy would be inrproved. For example, under the conditions used for these Simulations, at the end of data collection when the rate constant ratio is 5.0, the second reaction is only 39% complete, while the first reaction is 92% complete at the end of the analysis. Following the second reaction to 50-60% completion would improve the results for both components, as can be seen in Figure 4.1; however, the improved estimates would be greater for the second reaction. These experiments follow the general trend of improving accuracy as the responses of the components become approximately equal, and the reaction is sampled for at least 1 to 1% half-lives. 2° Spec tral Overlap Simulations studies were performed for the kinetic deteI‘Ininations of mixtures in which the products absorb with varying amounts of spectral overlap. The same kinetic behavior was assumed for these studies as that shown in Table 4.1. Three cases were tested: no spectral overlap, medium overlap, and severe overlap. Molar absorptivities for the three cases are given in Table 4.2. Figures 4.7 to 4.10 show the estimated concentration of component C 2 for no overlap, medium overlap, and severe overlap, l’eSpec tively. The response for component C1 was better, since it is the faster reacting species, and therefore, is sampled to a larger eXtEnt of reaction. The solid line represents the actual concentration, 90 Table 4.2: Molar absorptivities for spectral overlap study. Concentrations and rate constants appear in Table 4.1. nomerlan medium totalmerlan oletlan Zn 1.2 11 2a 1.1 2.2 £c1 0 O O 0 O 0 8C2 0 0 0 0 0 0 8R 10 5 10 S 10 5 3P1 5000 0 5000 2500 5000 5000 8P2 0 5000 2500 5000 5000 5000 91 while the dotted lines represent 5% errors in estimated concentration. Each rate constant ratio was repeated 9 times, and all points are included. Note particularly that the error in estimated concentration actually decreases with the rate constant ratio for the no overlap and medium overlap cases. This result was unexpected and has not been observed previously. As a comparison, the medium spectral overlap case was repeated using data from only the first wavelength. Results are presented in Figure 4.9. From this figure, it is apparent that the addition of the second wavelength causes a large increase in the accuracy of the results. Also, there is no improvement with decreasing rate constant ratio. In several cases, at relatively large rate constant ratios, the estimated concentration was in error by a large amount (> 100%). The reason for this is not clear; however, the statistics provided by the filter (e.g., covariance matrix) let the user know that there is a problem with the fit, and another trial can be done. Many values of molar absorptivities for all components were tested, although not as thoroughly as those described here. Examples were tested where one, some, or all components exhibited absorbance at the wavelengths used. The only cases where the filter was not able to obtain accurate estimates were when there was no absorbance change throughout the course of the reaction. Of course, other methods would also fail under these Conditions. Mixtures of species with rate constant ratios as low as unity were determined as long as there were sufficient spectral differences. In these cases, the rate constant ratio did not affect the ability of the filter to obtain accurate concentration estimates. As 92 1.1 1.06 1.02 4 02x10 0.98 0.94 ' 0.9 11.5 2 2.5 3 3.5 4 4.5 5 k1/k 2 Figure 4.7: Estimated concentration of C2 as a function of k1/k2 with go spectral overlap. (Solid line)-True concentration, (---)-5% error ars. 93 1.1 C o 9 : 0' I 1.05 ------------------- .-- .. m... . z. .. ..:::.::: ‘ ..; V .. .0 . E 1 . . 2 0.. .. :'§::':.E;:0':'.. "0: UN ..r 0|": 0....'°o: '0. . :;9 :. .:::s . to. O. 0.95 -------------------- .9- - O- - ...-.. .00 0 0.9 11.5 2 2.5 3 3.5 4 4.5 5 k1/lk2 Figure 4.8: Estimated concentration of C2 as a function of k1/k2 with IIledium spectral overlap. (Solid line)-True concentration, (--)-5% eI‘ror bars. Two wavelengths. 94 2 1.6 . a... ”o 12 0.. 0. . O O. :01 ”:7 '..‘11H1111n‘...111 '7 9. . 7. o o 1 * . '0 0 0.8 3?. '4 ' 3:, " 3" 3"" ' fig. .0 9 o 0.4 ° ' o 7 01 o . C 0 11.5 2 2.5 3 3.5 4 4.5 5 Figure 4.9: Estimated concentration of C2 as a function of k1/k2 with One wavelength. (Solid line)-True concentration, (--)-5% error bars. Fig 101; 95 11.5 2 2.5 3 3.5 4 4.5 5 k1/k2 1:igure 4.10: Estimated concentration of C2 as a function of k1/k2 with tOtal spectral overlap. (Solid line)-True concentration, (--)-5% error ars. 96 long as all reactions were followed for a long enough period, accurate results were obtained. Quite severe spectral overlap still allowed determination of reactant concentrations in most cases. As the amount of spectral overlap increases, it becomes necessary to use more wavelengths in order to obtain accurate concentrations from the filter. The filter has been tested using up to six wavelengths, and it will accept as many as the user desires up to the memory limit of the computer used. The problem with using a greater number of wavelengths is the trade-off in the amount of time required to obtain the estimates. This increase in data analysis time is not linear with increasing number of wavelengths used. This is due to the necessity to invert an m x m matrix, where m is the number of wavelengths used. In general, it is advisable to use the minimum number of wavelengths necessary in order to obtain estimates with an acceptable amount of error to the user. In practice, this may require some initial experimentation, but once the number of wavelengths required for a given system is known, it should remain constant. Of course, as microprocessors become faster, the number of wavelengths that can be analyzed in a given time should increase, and this method should become even more powerful. 3 - Effect of Random Noise The effect of varying the amount of noise added to the simulated data was investigated. The noise added to the true absorbance was Gaussian in nature, and its standard deviation was user-contr‘olled. The standard deviations for this study varied from 0.1 (0] 11h tea C011 9113 1101: C011 [0 y 97 0.001 to 0.9 absorbance units, and five replicate determinations were made with each set of conditions. The signal-to-noise ratio was 848 when both reactions are complete and s=0.001, and was 0.942 when both reactions are complete and s=0.9. Figure 4.1 1 shows the error in the average estimated concentration of C2 versus the standard deviation of the noise and the rate constant ratio, while Figure 4.12 shows the corresponding standard deviation of the mean estimated concentration of C2. Errors for component C1 were generally less severe as this reaction is followed to a greater extent completion in all cases except when the two rate constants are equal. In general, excellent accuracy (<5% error) was maintained as long as the standard deviation of the noise was less than approximately 0.1, which corresponds to a S/N ratio of =8.5 at the completion of both reactions. Figure 4.13 shows the absorbance versus time curve for an example with a standard deviation of the noise of 0.1, and a rate constant ratio of 1.2. The estimated concentrations for this experiment were: C1=9.63 x 10'5 and C2=1.05 x 1074. If an error of <20% is acceptable, the standard deviation of the noise can be as high as 0.5, corresponding to a S/ N of 1.7 at Completion. Overall, the filter is highly tolerant of noise, and is able to yield accurate results at very poor signal to noise ratios. 98 ’4’! l '0101 ‘- ‘\ 11, (iii. ’1 1111 ""1111 111 at... 13 111111111111]! ‘1“‘11/1111‘1"é1' :3; —L "_ # a ‘5‘» .1 1: Error in the average estimated concentration of C2 as a e 4 tion of rate constant ratio and the standard deviation of the 11 [81 99 70* 607 50— 40" b; 1 1 I 30- ll; ,1 1‘11 20 MIMI“ i" ”’3’" Figure 4.12: Standard deviation of the mean of C2 as a function of rate constant ratio and the standard deviation of the noise. 100 Abs Figure 4.13: Absorbance and estimated absorbance for k'1/k'2 =1.20 and standard deviation of the noise of 0.10. (Solid line)-True absorbance, (+)-estimated absorbance. 101 4. Effect of Varying Estimated Rate Constant Since the Kalman filter is recursive in nature, it requires initial estimates for the parameters. As mentioned previously, to avoid biasing the filter, the initial estimates provided to the filter for the unknown concentrations are 0. Also, since the concentration of the common reagent is known, this value is provided to the filter. Some knowledge of the value of the rate constants is necessary in order for the filter to calculate accurate estimates of concentrations. However, rate constants are dependent on temperature, pH, and ionic strength, as well as other factors. Therefore, the value of the rate constant may change somewhat from sample to sample. With carefully controlled reaction conditions, this variation can be minimized, but not totally eliminated. In order to determine the effect of inaccurate estimates of the rate constants on the results of the filter, a study was performed in Which the values of the rate constants given to the filter were varied from 50 to 150% of the actual value. Three cases were studied. In all cases, the value of k1 was kept constant at 50 3‘1 while kg was given values of 10, 25, and SO 5‘1 in the three cases. This corresponds to rate constant ratios of S, 2, and 1, respectively. Plots 0f error in the estimated concentrations of both components versus the input values of the rate constants appear in Figures 4.14 to 4.19. As is apparent from the charts, errors are at a minimum for values of R1 and k2 closest to the true values, and generally get worse as eStimates get farther from the true values. As long as the estimate for the rate constants are within 10%, the error in the estimated 102 1 initial k initial k =5.00. Figure 4.14: Error in estimated concentration of C1 as a function of the initial estimates for the rate constants when lq/kz 103 250- / :I/ \\\\\\\\\ initial k 1 initial k 2 Figure 4.15: Error in estimated concentration of C2 as a function of the initial estimates for the rate constants when k1/k2 =S.OO. Fig? 104 m, ”IIIIIIIIIIII II 1 II II . [Inn/15’ ’IIIIIIIIIIg I h 8 "III ”IIIIIIIll/ , . i$fiflZ“tIifltfifififilflnl‘lljyzza‘llll E; I II:III’IIIIII””II%, 4 \\ s \x ‘ \ I’l’Illllcggcl -\~\\ l77711115,,5 \ \‘ \\ 'I” I’IIIIIII‘III I, s \s {“l. ‘2 III/1 [’3’ ‘ ‘ \\‘\ ’1‘ ‘ \ \ \» I’,’ ’16”, “ “fit/15;; ‘ €512” \ \ - , 35 75 65 initialk 1 initialk 2 38 Figure 4.16: Error in estimated concentration of C1 as a function of the initial estimates for the rate constants when k1/k2 =2.00. 105 23 75 initial k z 38 initial k 1 Figure 4.17: Error in estimated concentration of C2 as a function of the initial estimates for the rate constants when lq/kz =2.00. 106 / 25 25- -20 2°“ IIIIIIIIIIIII ‘5 IIIIII c (96) ”III” _ 1 0 1 . 10- 5 - '5 25 55 65 75 65 Initial k 1 initial k 2 75 Figure 4.18: Error in estimated concentration of C1 as a function of the initial estimates for the rate constants when k1/k2 =1.00. Figl 107 Error initial k 2 75 Figure 4.19: Error in estimated concentration of C2 as a function of the initial estimates for the rate constants when lq/kz =1.00. 108 value of C1 was less than 2.5%, and the error for C2 was less than 16% (this error decreases to 2.5% as k2 approaches k1). Also, in general, it seems to be preferable to overestimate the rate constants, rather than to underestimate them. 5. Other Considerations The filter performed remarkable well over a large range of concentrations, as long as the absorbance change was approximately 0.05-O.l absorbance units, and as long as there was a reasonable signal-to—noise level. The filter was also tested at a variety of molar absorptivity values for all components. The only area where it failed was when the reactants and products had the same molar absorptivity such that there was no absorbance change observable. The filter was able to estimate concentrations from this level up to a level where the accuracy of the absorbance measurements would be instrument dependent. Generally, this allows for linearity of two to four orders of magnitude. The filter is also able to perform other functions. For example, using equation 4.7, the final absorbance for the reaction can be estimated at each point. The accuracy of this estimate is also dependent on the extent of the reactions, and the error becomes quite acceptable when 1 to 1.5 half-lives have been completed. This provides a way of performing predictive kinetic determinations“. Since the absorbance of a system at equilibrium is often less dependent on temperature, pH, and ionic strength, than is a dynamic syst (Off. 1381‘. adv preI SOII‘. 10% con 8th ma rapi long 109 system, such a predictive approach would provide additional error compensation. Use of the extended Kalman filter and multiple wavelengths to perform multicomponent kinetic determinations offers several advantages. First, the nonlinear filter does not depend on knowing precisely the rate constants of the reacting components. Although some knowledge of rate constants is necessary, accurate results can be obtained as long as the value given the filter initially is within 10% of the true value. By using multiple wavelengths, the kinetic constraints of previous approaches have been overcome. As long as spectral differences exist, accurate results are obtained even if the reactions occur at exactly the same rate. The filter processes data rapidly, which allows the possibility of data analysis in real-time as long as the reactions under investigation are sufficiently slow. LIST OF REFERENCES LIST OF REFERENCES 1 M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, 1968, Ch. 26. 2 L. Kaufman, Trends in Anal. Chem., 1983, 2, 244. 3 Peter D. Wentzell, Ph.D. Dissertation, Michigan State University, 1986. 4 G. E. Mieling and H. L. Pardue, Anal. Chem., 1978, 50, 1611-1618. 110 CHAPTER V Simulations of Complex Chemical Reactions The simulations described in Chapter IV were continued using models of other chemical systems. Several other types of reacting systems were studied; some with greater than two components present, and others with reversible chemical reactions present. Although these complex systems were not studied in as great detail as was the two-component irreversible system , it is shown that the extended Kalman filtering method is applicable to complex systems as well. These systems and the results of the simulation studies are discussed in this chapter. A. Kinetic Models 1. General Irreversible Reactions The two-component irreversible model developed in Chapter N Can be generalized to a system with any number of components all reacting with a common reagent. Such a general case can be dEScribed as follows: 111 whi con app assm fum SpeI the Rite Equa [em This 112 C,+R—"1—)P, C2+R—"J—-)P2 (5.1) Cn+R—5—+Pn which simplifies as was shown previously in equation 4.2 if the common reagent, R, is in excess, and pseudo-first-order conditions apply. In the general case, the absorbance at any wavelength, j, and at any time, t, is given by: A], = 85c. [C1]t + rajc2 [C2]‘ + . . . + 81c, [Cal + 51?.[1’1]. (5.2) +8].p2 [P2]t + . . . + em [Pu], + 5:1.R[R]t assuming that the pathlength is unity. If we now substitute the functions for the time-dependent concentrations of each of the Species in equation 5.2, collect the terms and simplify, we arrive at the general absorbance function for any number of components, n, given by: It +(ejc_ + 85R - £1,_)[Cn]°c“"-‘ + (53?. - €51: )[C,]o (5 .3) +(ejpz - ejR )[Cz], + . . . + (£11,. - a), )[Cn], + sz[R]o _ -k': -k't A. -(£jC, +£iR-EiPrXClioe l +(85C2 +£5R-85P2 Xczloc 2 +0“ Equation 5.3 can then be truncated after the appropriate number of terIns to give the correct model for the extended Kalman filter to use. This measurement function is expressed in (2*n)+1 parameters {81' If H firs: Spet The ; descr COHQ “381111 113 where n is the number of unknown components. The parameters that are unknown to the extended Kalman filter are: [C 110, [C2]o, ..., [CnJo, [R]o, k'1, k'2, ..., and k'n. This equation can be used for any number of parallel first— or pseudo-first-order irreversible reactions, and for any number of wavelengths. 2. Reversible Models If conditions are such that the reaction is reversible, yet another model is required for the extended Kalman filter. A reversible chemical system can be described by the following: C,+R€$P, (5.4) If the concentration of R is held in large excess such that pseudo- first-order reversible conditions apply, the concentration of each sPecies with respect to time can be given by: _ [c1 [C1] (19+ +k_)[k -1 _kl ] kIC] 131:“ (k —k—+_,)[l-c [R]: = [R10 - ([CIL " [C1].) —(k,+k_, )1] (55) The absorbance at any time, t, and any wavelength, j, of the system described in equation 5.4 can then be obtained by substituting, C0llecting terms, and simplifying the terms in equation 5.5. The re$ult of this is: filII esu' 1'91' was abo If we the 114 Art =(8iCr+£ J'R(—[C—]_—)[ k H: )‘1-[k R‘s-(WW1) (5.6) +£jPl[(—_——:ll[ Cll‘]___o_1 ll__c)[ c,-(k +k ;)t]]+£jR[R]o This equation is now expressed in four parameters, [C 1]o, [R]o, k1 and k1. These four parameters are those that the extended Kalman filter must estimate. This form of the filter could now be used to estimate these parameters from a one-component system following reversible kinetics. We can again develop a model for n-parallel reacting species as was done for the pseudo-first-order non-reversible reactions case above. This case can be shown by: kI C1 + RfvP, 1:. C2 + Rsz (5.7) kl Cn+RfoPn If we now extend equation 5.5 to include all reactions, and substitute the time-dependence of each of the species into the absorbance function, we obtain the following: 'meq issht 115 Aj,=(ejcl +5 4);“ k Ci—-T°-—-l)[k__ 4,4441% gjpl[(k‘[ Ci—Lol )[H 4.. +4.41}. 44: 4 4 1H [4 (816- +8 {ah—I21]?— ..k'“)[ Iii-“JWUH’ ((1: “Ch 1:1:1-6)[( ”1)+ 8,-an1. As is apparent from the above derivations, it is possible to derive a model to be used in the extended Kalman filter as long as one can describe the system mathematically. More complex systems tend to have more limitations, and are not applied to the filter as easily; these problems are discussed in more detail later. 3. Use of Models with the Kalman Filter Similarly, the state vector for the n-component case described in equations 5.1 to 5.3 can be shown to be an (2*n)+1x1 vector which is shown in equation 5.9 below: The linearized measurement function for the n-component system is a 1x(2*n)+1 vector, and is analogous to the result obtained above for the two-component case extrapolated to 11 components. In the case of reversible reactions, described for the one- component case in equations 5.4 to 5.6, the state vector consists of a 4x1 vector given below: ’[CJJ [R]. kl The linearized measurement function for this system is a 1x4 vector given by the following: L k-1 (5.9) (5.10) the she (01] Vets gen, talc Derf Iepe 117 1 - + -l t k " 1+ 4‘ L“ =£3C1(k +k )(k_l+klc (k, k ))-|-(£jpl -£jk)(k +lk yl—c (k k )) l -l 1 -1 L =8- [C110 —k-l +c‘(kl*k-I)‘ 1-6(1 __ kl ‘3 ’C' k,+k_ k H: k «H: +(£~p -£- a) [‘—C—'—]°1L—c'("”"')t 1- tk + k " k +k_ k +k_, +_k 1‘1 L :3. [C110 1____+kck-1-(k.+k-.): (5.11) 14 ’C‘ kid-lg1 k1+k- +(8jp._£ja)[kl[Cl]° -1 +c‘(‘““‘")‘ 11:11:) k1+ k_1 k1 + l(_1 This system can easily be extended to the general case where there are n parallel reactions. The linearized measurement function for this case is a 1x(3*n)+1 vector, and is analogous to the equations shown above for the one-component case extrapolated to 11 components. B. Experimental Simulated data were generated in order to test the different versions of the extended Kalman filter described above. Data were generated using an appropriate model for the system under study. Random noise was then added to the noiseless absorbance data calculated as is described in Chapter N. Data processing was also Performed as described in Chapter IV, and so its description is not repeated here. 311 C0. 118 C. Results and Discussion 1. Three-Component Irreversible System A three-component system was tested in order to determine the effect that the values of the three rate constants have on the errors in estimated concentration of the components. Such a system is described by equation 5.1 where n=3. In the system studied, R1 was held constant while IQ and k3 were varied. The values of the concentrations, rate constants, and molar absorptivities are given in Table 5.1. For this study, three wavelengths were used, and data were simulated for 100 times. Plots of error versus IQ and k3 appear in Figures 5.1 to 5.3 for each of the three components. From these, it is apparent that the general trend mentioned in Chapter IV for the two component case is again followed here. In order to obtain accurate results, it is necessary for there to be either sufficient spectral or kinetic differences, or a combination of the two. Also, it is necessary to follow the reaction of interest to a sufficient degree. In Chapter IV, it was found that, in general, it was necessary to follow a particular reaction to at least 50-60% completion. Figure 5.2 shows a high error in the estimated concentration of C2 when the rate constant for the second reaction is equal to 10 5'1, and moderate errors (IO-15%) occur when the rate constant is equal to 20 5'1. For these two cases, the reaction has not progressed to an extent Sufficient for accurate results to be obtained. When the rate constant is equal to 10 $4, the reaction is only ~40% complete when the data C011ection is complete. In Figure 5.3, the high error in the estimated Table 5.1: Values of molar absorptivity for the three component SYStem Studied. [C110 =[ C2]0 = [C3]0 =1 X 10'4 M, [R]O = 0.01 M, k1 = 100 8'1, k2 = 10-100 8'1, k3 = 10-100 5'1. 8C2 8C3 8R 8P1 8P2 8P3 11 5 10 15 2 750 2000 2800 12 20 16 9 4 3000 1800 1200 13 1 1 20 5 3 1500 3000 1700 Figure 5.1: Error in estimated value of C1 versus IQ and k3. 121 #120 Figure 5.2: Error in estimated value of C2 versus IQ and k3. 122 mwmwmmmmwo Figure 5.3: Error in estimated value of C3 versus kg and k3. tncentrat' reaction is the estt natal or COMBINE 2. Pour-C 123 concentration of C3 occurs when the rate constant for the third reaction is equal to 10 5'1. Figure 5.1 shows no areas of large error in the estimated concentration of C1. This indicates that sufficient spectral or kinetic differences are present in this example so that the concentrations can be accurately determined. 2. Four-Component Irreversible System A system consisting of eight simultaneously reacting components was simulated. From this system two sets of four components each were chosen with different degrees of spectral overlap. These components are the same as those used for the eight- component study which is described later. A four-component system can be described by equation 5.1 where n=4. For this trial, two separate sets of conditions were used. The first consisted of a moderate amount of spectral overlap, and the second consisted of a greater extent of spectral overlap. Synthetic spectra were generated for eight separate reactions at 15 different wavelengths. Each of the reacting components was given a molar absorptivity of 0, while the common reagent was given a small (but not 0) molar absorptivity. Each of the eight products was given a Gaussian absorbance profile, with an identical peak response. For the case with moderate spectral overlap, the peak-to-peak separation was 0.8 standard deviations from one peak to the next, while for the case with greater spectral overlap, the peak-to-peak separation was 0.4 standard deviations. Figure 5.4 shows the synthetic molar absorptivities for the eight Components and the common reagent, R, 124 5000 if 3000 .z:::.. 1000 1 3 5 7 9 11 13 15 wavelength # Figure 5.4: Synthetic spectra of eight components. (solid line)- products, (dotted line)-common reagent. 125 at all fifteen wavelengths. The moderate overlap case used the first, third, fifth, and seventh spectra, from left to right in Figure 5.4, while the case with the greater spectral overlap used the first four spectra. For both cases, the initial concentrations of all reactants was 5 x 10‘5 M, the concentration of the common reagent was 0.01 M, and the rate constants were: k1=60 5'1, k2=40 5'1, k3=20 5'1, k4=30 S'l, giving a maximum rate constant ratio of three. The number of wavelengths used was varied from three to six, and the wavelengths used corresponded to the first n wavelengths from Figure 5.4. Table 5.2 shows the final estimated concentrations and errors for each of the components for the moderate overlap case, while Table 5.3 shows the final estimated concentrations and errors for each of the components for the greater spectral overlap case. As can be seen from the tables, the errors from the moderate overlap case are smaller than those from the greater overlap case. It is also apparent that the errors decrease with increasing number of wavelengths used. However, this is a trade-off, since the amount of time necessary to complete the data analysis increases with increasing number of wavelengths used. These examples also illustrate the general rule that it is necessary to use at least as many wavelengths as there are components present, and sometimes more. The number of wavelengths necessary for a particular analysis is dependent on several factors, such as: number of components present, amount of noise present, and spectral and kinetic overlap of the species. If there is either total spectral or total kinetic overlap of the species, and no noise is present, the number of wavelengths necessary would be equal to the number of components. This is so that there are at 126 Table 5.2: Final estimated concentration, and error in the estimated concentration for the four-component case with moderate spectral overlap. Numhemfflaxelengths 3 4 s 6 [C1]o/1O'S 5.00 5.01 5.02 5.03 [C2]o/ 10-5 5.8 5.09 5.05 5.03 [C3]0/1O'5 4.79 4.93 5.00 5.05 [C4]o/1O-5 5.27 5.16 5.06 5.03 err C1/% 0.0 0.2 0.4 0.6 err C2/96 2.8 1.8 1.0 0.6 err C3/96 -4.2 -1.4 0.0 1.0 err C4/% 5.4 3.2 1.2 0.6 127 Table 5.3: Final estimated concentration, and error in the estimated concentration for the four-component case with greater spectral overlap. WW 3 4 5 6 [C1]o/ 10’5 4.87 5.20 5.6 4.92 [C2]o/10'S 5.48 4.26 4.83 5.39 [C3]o/10'5 4.49 6.23 5.25 4.58 [C4]0/10'S 5.23 4.5 1 4.94 5.23 err C1/% -2.6 4.0 2.0 -1.6 err C2/% 9.6 -14.8 -3.4 7.8 err C3/% -10.2 24.6 5.0 -8.4 err C4/% 4.6 -9.8 -1.2 4.6 C0: 128 least as many equations as there are unknowns. However, with noise present and varying degrees of overlap present, the actual number of wavelengths required is dependent on the factors above. For example, one wavelength is sufficient as long as enough kinetic difference is present, but as the kinetics of the systems under study exhibit more overlap, a greater number of wavelengths are required. Since the number of wavelengths used in the data analysis is the limiting factor in the amount of time the analysis takes, it is generally advisable to use the smallest number of wavelengths that provides accurate results. This number can be determined from an analysis of a known mixture, and should be constant as long as the experimental conditions do not vary too greatly from the standards. The addition of extra wavelengths after accurate results have been obtained does not provide much more accuracy. This is illustrated in Table 5.2 in which the increase from three to four wavelengths provides a fair improvement in the accuracy of the estimates, while increasing beyond four wavelengths does not provide much more improvement. This can also be seen in Table 5.3 where an increase beyond five wavelengths does not provide much improvement in the accuracy. 3. Eight-Component Irreversible System A study using eight components was performed in order to test the limits of the filter. This system can be given by equation 5.1 where n is equal to eight, and corresponds to a state vector that contains 17 parameters. The molar absorptivities for all reactants in 129 this study was taken to be 0 at all wavelengths used. The common reagent, R, was assigned a small molar absorptivity at all wavelengths, described by a Gaussian distribution. The spectra of all products was given by a Gaussian distribution, and the peak-to-peak separation of each of the components was 0.4 standard deviations. Figure 5.4 shows the spectra for all eight of the components plus the common reagent. This system was tested using both the first ten wavelengths, and all fifteen wavelengths, and it was found that all fifteen wavelengths were necessary in order to obtain accurate estimates of the concentrations of the species. The concentrations of all species were taken to be 5 x 10'5 M, the concentration of the common reagent was taken to be 0.01 M, and the rate constants for each of the reactions were: k1=50 3'1, k2=40 s-l, k3=30 5'1, k4=20 3'1, k5=25 S'l, k6=35 8'1, k7=55 8'1, and k3=60 5'1. This yields a maximum rate constant ratio of three. A summary of the estimated concentrations of each of the components, and their errors is given in Table 5.4. It is apparent that increasing the number of wavelengths from 10 to 15 greatly improves the accuracy of the determination. Unfortunately, this analysis was quite time consuming. Running on a 286-type microcomputer equipped with a math coprocessor running at 10 MHz required 3400 seconds to carry out the analysis using 15 wavelengths. Running the same analysis on a 486-DX type microcomputer running at 33 MHz required approximately 1700 seconds. It is certainly possible to perform determinations of even greater numbers of components. However, if the spectral overlap is large, as it often is in systems of this type, the analysis time can Table 5.4: Estimated concentration of each component and error in components for eight-component case. OO\l0\U1AUJNt—Iu—o 130 10 wavelengths 8.85 77.0 -9.27 -285.4 30.7 5 14.0 -22.3 -546.0 22.7 354.0 —0.67 —1 13.4 4.85 -3.0 5.43 8.6 15 wavelengths 5.15 4.59 6.15 4.15 4.75 6.00 4.42 5.20 3.0 -8.2 23.0 -17.0 -5.0 20.0 -1 1.6 4.0 real We 131 become prohibitive. Of course, as processors become even faster, this will not be as great a problem. Whether or not a system of this complexity should be analyzed using this method depends on a number of factors. These include: suitability of other methods for the analysis, the time other methods would take, and availability of powerful computers. 4. Reversible Systems The extended Kalman filter method has also been tested using systems which follow reversible kinetics such as the reaction given in equation 5.4. Such a system is sometimes seen in ligand displacement reactions. For example, the displacement of metal/EGTA complexes with 4-(2-pyridylazo)resorcinol (PAR) shows reversible behavior for some metals. This version of the filter was tested to see which ratios of the forward and reverse rate constants are tolerated. For this study, a one component model was used, and the value of the forward rate constant was held constant at 100 5'1, while the reverse rate constant was varied from 0.1 to 500 $4, giving ratios of k1/k-1 of 1000 to 0.2. A plot of the error in the estimated concentration of C1 is plotted against this ratio in Figure 5.5. As can be seen, all ratios tested yielded very accurate results, and the errors only started to become appreciable as k1 became much larger than k.1. Of course, under these conditions, the reverse reaction is negligible, and the irreversible form of the filter discussed Previously is applicable. 132 2.5 2 E 1.5 LIJ 1 0.5 OJ 1 10 100 1000 k1/k _1 Figure 5.5: Error in estimated C1 versus ratio of forward to reverse rate constants. 133 A two-component version of this system was also studied. This system can be described by extending equation 5.4 to two parallel reactions. This system was tested for the effect of the ratio of the rate constants of the two forward reactions. In this case, the rate constant for the first reaction, k1, was held constant at 200 8’1, and the forward rate constant for the second reaction was varied from 20 to 200 s*1, which gives values of the ratio of the forward rate constants of 10 to 1. For all cases, the rate constants of the two reverse reactions were held constant at 20 5'1. Figure 5.6 shows the error in the estimated concentrations of the two components as a function of the ratio of the forward rate constants. As can be seen, the error in the first component is relatively constant, and is less than 1.0 96 for all cases, while the error in the second component becomes smaller as the rate constant ratio approaches unity. This is again due to the filter seeing more of the second reaction. These systems were also tested to see the effect of the initial estimates of the rate constants and of added noise. Results were similar to those found for the irreversible reactions discussed in Chapter IV. These systems are tolerant of noise to quite a high degree, and are also capable of converging on accurate results even if the initial values for the rate constants are not completely accurate. As long as the estimates are within 10 to 20% of the true values, the estimates from the filter are quite accurate. Fig 1‘l 11;; 134 Error (%) 1 2.8 4.6 6.4 8.2 10 k1"‘2 Figure 5.6: Error in estimated concentration of C 1(0) and C2(+) versus kl/kz. 135 D. Conclusions The extended Kalman filter has been successfully applied to a number of systems of varying complexity and differing kinetic parameters. This method has been shown to be quite versatile and tolerant of the initial estimates of the kinetic parameters supplied. The filter has also been shown to be able to resolve systems containing a larger number of components with highly overlapped spectral and kinetic characteristics. The method of the extended Kalman filter is applicable to any kinetic system as long as the behavior of the system is known and sufficient kinetic and spectral differences are present. These differences do not have to be very large in order for accurate results to be obtained. This can be useful for providing information on systems with several similar species present in a sample, all of which have highly overlapped responses to the detection method and highly overlapped kinetic responses. An example of such a system is given in Chapter V1 with the complexation reactions of a series of metals. This technique could also be applied to functional group determinations in organic mixtures with a reagent specific to a certain functional group. Other areas of application include resolution of fluorescent and phosphorescent species. Since the decay of luminescence is exponential, and different species have different luminescent lifetimes, the extended Kalman filter technique described here can be used to determine multiple species in these methods as well. 136 The extended Kalman filter is the main form of the filter used in the current work. An extension of this work could be the use of a combined adaptive Kalman filter-extended Kalman filter. The adaptive Kalman filter is not applicable to nonlinear systems. However, by using the adaptive filter to identify problems with the model used, and the extended Kalman filter to analyze the data points with no problems, interferents could be detected more readily. Another extension of the current work involves the use of Kalman filter networks for principal components analysis. This would be a very interesting application. By using networks of extended Kalman filters, the number of reacting species could be determined, as well as their concentrations. Chapters IV and V describe simulation studies of multicomponent kinetic systems using the extended Kalman filter. Reactions of differing complexity are modeled and tested for their suitability for use with the filter. Another area where the work could be continued is in the use of mixed models. This would allow analysis of systems following mixed-order kinetic processes. Another area where the work could be furthered would be the use of a Kalman filter on a chip. It is possible to produce integrated circuit chips to perform many of the calculations used in the Kalman filter. If these chips become readily available, they could greatly increase the speed of the calculations. The advantages of this are obvious: faster processing of the data would allow fast reactions to be analyzed in real time, and more complex models containing a greater number of reacting species could be utilized. CHAPTER VI Simultaneous Kinetic Determination of Lanthanides with the Extended Kalman Filter A. Introduction Many metal complexing agents are in common use in analytical chemistry. One of these is 4-(2-pyridylazo)resorcinol (PAR) whose structure is shown in Figure 6.1. PAR has proven to be very versatile for the determination of metalsl. PAR reacts with over 50 elements to form either colored complexes and/ or precipitates from aqueous solutions. The resulting complexes which are formed when OH Figure 6.1: Structure of 4-(2-pyridylazo)resorcinol (PAR) a metal and PAR are reacted are often highly colored, with molar absorptivities ranging from 104 to 105. These highly colored 1‘37 138 solutions offer the ability to determine trace amounts of metal species in aqueous solutions. The versatility of PAR as a reagent for many metals is also one of its shortcomings. This is because the complexes resulting from the reaction of a metal and PAR are quite similar spectrally. Therefore, determination of most metals with PAR is interfered with by the presence of other metals. This can be compensated for by maximizing the formation of one of the metal-PAR complexes while trying to minimize the formation of all other complexes. Selective formation can be accomplished in some cases by pH control. The spectra of the resulting complexes are quite overlapped, and do not offer enough differences by themselves to lead to a separation of signals for many species. However, if these spectral differences are combined with kinetic differences, the combination can be enough for resolution of mixtures. The approach of combining the spectral and kinetic differences for the simultaneous determination of mixtures using PAR as a common reagent has been evaluated. For all of the previous studies, an exchange reaction has been used such as that described by equation 6.1: M-EDTA+PAR4=M-PAR+EDTA (6.1) where EDTA is ethylenediamintetraacetic acid, and M stands for any metal. Using this exchange reaction serves two purposes. First, most complexation reactions of a metal with PAR are rapid and occur on the stopped-flow time scale, while the exchange reaction occurs on a longer time scale, with reactions often being monitored for several 139 minutes. The second purpose of the exchange reaction is that the differences in the rate constants of the exchange reaction are usually greater than for the metal-PAR complexation reaction. This approach has been used for the determination of gallium and indiumz. The exchange reaction described in equation 6.1 has also been applied with a substitution of EGTA (ethyleneglycol bis(2- aminoethylether)N,N,N',N'-tetraacetic acid) for EDTA. With this reaction, mixtures of cobalt and nickel have been analyzed by logarithmic extrapolation3, and more recently by non-linear regression and factor analysis with multiple wavelengths4»5. If no exchange reaction is used, the sample throughput could be much greater. Since the complexation is quite fast, monitoring of the direct reaction would be necessary for only several seconds at most. This would allow many samples to be tested in a short time. The current work focuses on the reaction of PAR with a group of lanthanides. The straight complexation reaction is used, and the amount of time that the reactions were monitored was less than two seconds. Two and three-component mixtures of lanthanum, praseodymium, and neodymium were determined based on their spectral and kinetic differences. The reaction can be described by equation 5.1, and the measurement function used for the data analysis can be given by equation 5.3, when both equations include the appropriate number of terms. 140 B. Experimental Section 1. Apparatus A Tracor Northern (Model TN-6123) SlZ-element intensified linear photodiode array was used to acquire spectra in the 300-700 nm wavelength range. Software to control the data acquisition was written in house. A stopped-flow reagent mixing system, constructed in house6, was used to mix the reagents. 2. Reagents 4-(2-pyridylazo)resorcinol (Sigma Chemical Co. #P-8019) was used as the common reagent. Salts of all metals used were dried and used to prepare stock solutions. No buffering or ionic strength control was used. 3. Procedure Solutions of 4-(2-pyridylazo)resorcinol (PAR) were mixed with solutions of the appropriate mixture of metals in the stopped-flow mixing system. Progress of the reaction was followed spectrophotometrically with the linear photodiode array (LPDA). Using the LPDA, it is possible to take one complete spectrum every 3.5 ms, and up to 100 spectra may be acquired in a single run. Actual length of time between acquisition of spectra was varied as necessary. 141 4. Computational Aspects The extended Kalman filter programs used were written in QuickBASIC (Microsoft Corp.) version 4.0; the algorithm and the measurement model described in Chapter III were used. Calculations were carried out on a 286 type computer equipped with a math coprocessor. C. Results and Discussion In all situations described below, certain information needs to be supplied to the filter for it to be able to obtain accurate estimates of concentrations. First of all, it is necessary for the algorithm to have the pure component spectra available to it for all species present. This is done in practice by having separate data files containing molar absorptivities for each component in a reaction. In order to provide estimates of the rate constants to the filter, they were first determined by fitting a nonlinear function to a single run of each metal. The relative rate constants for lanthanum, praseodymium, and neodymium are 1:1.7:1.9. 1. Single Component Studies In order to determine how well the filter works under the conditions used in this work, it was first tested in a one-component system. 142 Varying amounts of each metal solution and PAR were mixed in the stopped flow under the conditions given in Table 6.1. The extended Kalman filter algorithm described above was then applied to the data from each diode individually in order to determine which wavelength regions were best for this study. Figure 6.2 shows the molar absorptivities for each of the products of the complexation reaction. As can be seen, there is a great deal of spectral overlap among the species. This spectral overlap, combined with the closeness of the reaction kinetics for the three species makes this system an excellent test of the extended Kalman filter method for multicomponent kinetics. Figures 6.3 to 6.5 show the estimated concentration of each of the three lanthanides versus the wavelengths used for analysis for an example run. Each of the three lanthanides was run at several concentrations from 7 x 10‘6 M to 3 x 10'5 M. The error bars represent the standard deviation in the filter estimates for triplicate determinations. Relative standard deviations were approximately 10% or less for each of the three metals in the region of the optimum wavelength. If a greater number of wavelengths is used by the filter for the analysis, the precision improves. However, the greater the number of wavelengths used, the longer the filter takes to analyze the data. It is apparent from Figures 6.3 to 6.5 that the area from approximately 450 nm to 550 nm is analytically useful for all metals while the rest of the spectral region studied is less useful. The spectral region below 450 nm is not useful due to the strong absorbance of PAR. Since PAR is in excess to maintain pseudo-first- order conditions, the overall absorbance in this region changes little 143 Table 6.1: Conditions for spectral acquisition of single component data Metal Time between Delay before Total data scans (ms) first scan (ms) collection time (ms) La(III) 8.5 or 10.5 13.0 or 15.0 863 or 1065 Pr(III) 8.5 or 10.5 13.0 or 15.0 863 or 1065 Nd(III) 8.5 or 10.5 13.0 or 15.0 863 or 1065 2.00 x 10‘4 M PAR used for all reactions 144 300 350 400 450 500 550 600 650 700 nm Figure 6.2: Molar absorptivities of lanthanides and 4-(2- pyridylazo)resorcinol (PAR). (o)-PAR, (+)-La(PAR) complex, (__)- Pr(PAR) complex, (--)-Nd(PAR) complex. 145 5.5 10‘5 5.0 10'5 4.5 10'5 M 4.0 10'5 est [La] 3.5 10'5 3.0 10'5 2.510’5 425 468.7 512.5 556.2 600 nm Figure 6.3: Estimated initial concentration of 2.88 x 10'5 M lanthanum with respect to wavelength. True concentration shown by solid line. 8.5 ms scans. 146 425 468.7 512.5 556.2 600 nm Figure 6.4: Estimated initial concentration of 1.42 x 10‘5 M praseodymium with respect to wavelength. True concentration shown by solid line. 8.5 ms scans. 147 3.0 10'5 2.5 10'5 2.0 10'5 est [Nd], M _L '01 —L O O! 1.010‘5 5.0 10'6 0 . 425 468.7 512.5 556.2 600 nm Figure 6.5: Estimated initial concentration of 1.39 x 10'5 M neodymium with respect to wavelength. True concentration shown by solid line. 8.5 ms scans. 148 throughout the course of the reaction. Therefore, the Kalman filter has little information to work with in this region, and is not able to obtain accurate estimates of the metal concentration. Above 550 nm, the metal/PAR complexes exhibit little absorbance, and so again, there is little absorbance change throughout the course of the reaction, and little information for the filter to work with. Of course, in these wavelength regions, other techniques would also fail to obtain accurate results, since there is little useful analytical information. Figures 6.6 to 6.8 show the estimated concentrations of samples of lanthanum, praseodymium, and neodymium with respect to time that result from the filter. These figures represent single component data processed at only a single wavelength (495 nm). For consistency, the wavelength was kept the same for all three determinations. Note that this wavelength does not necessarily provide the Optimal estimates for any of the metals. However, this wavelength does provide quite accurate results for all three metals, while other wavelengths may be better for some of the metals. It is apparent that the filter settles onto an approximate value early in the reaction, and after that, the estimated concentration values vary little from point to point. Error bars are included on the plots, and these represent square roots of values taken from the diagonal of the covariance matrix, which yields the standard deviation of the estimated values. As the reaction progresses and more data are available to the filter, the precision of the estimated values improves while the accuracy is essentially unchanged once the filter settles onto a stable estimate. Note that the size of the error bars is quite 149 2.8 10'5 2.4 10'5 2.0 10'5 21.610‘5 .--------------------------------------. "1.210‘5 8.010'6 4.010‘6 0 0.2 0.4 0.6 0.8 1 1.2 Time(s) Figure 6.6: Estimated concentration using 1 pixel at 495 nm of 1.44 x 10‘5 M lanthanum with respect to time. Error bars :1 standard deviation. (_)-True concentration, (--)-10% error in concentration. 10.5 ms scans. 150 2.0 10'5 1.610'5 ...1.210'5 ’ est 8.0 10'6 4.0 10'6 0 0.2 0.4 0.6 0.8 1 1.2 Time(s) Figure 6.7: Estimated concentration using 1 pixel at 495 nm of 1.42 x 10-5 M praseodymium with respect to time. Error bars 21:1 standard deviation. (__)-True concentration, (--)-10% error in concentration. 10.5 ms scans. 151 2.0 10'5 1.610'5 1.210’5 ’ est[] 8.0 10"6 4.0 10'6 0 0.2 0.4 0.6 0.8 1 1.2 Time (s) Figure 6.8: Estimated concentration using 1 pixel at 495 nm of 1.39 x 10'5 M neodymium with respect to time. Error bars :tl standard deviation. (__)-True concentration, (--)-10% error in concentration. 10.5 ms scans. 152 small. This provides an indication that the filter has enough information to give accurate estimates of concentrations using only a single wavelength. This is not surprising, as most multicomponent kinetic methods only use data from a single detector. However, as the number of components increases, and the spectral and kinetic overlap increases, the errors can increase if only a single wavelength is utilized. Figure 6.9 shows the residual absorbance values for a reaction of neodymium. This figure is from the same run as the concentration estimate shown in Figure 6.8. Note that the residual values are randomly distributed about zero. This shows that the model used to fit the data in the extended Kalman filter algorithm is adequate. The only residual that is large and is not due to fitting error is the very first point. This residual is due to the initial values provided to the filter. The values of the rate constant provided to the filter were also varied from the calculated value. These values were varied from 50 to 150% of the true value. Although the best errors were usually observed when the calculated value of the rate constant was provided to the filter, even values that were off by as much as 50% provided excellent results. The maximum increase in error for any of the trials tested was less than 3%. This seeming insensitivity to the value of the rate constant is due to the nonlinear form of the Kalman filter. Since the rate constant value is one of the parameters that is adjusted at each iteration, as long as the rate constant provided to the filter initially is approximately correct, the filter is able to calculate accurate results. 153 0.06 0.05 0.04 0.03 0.02 0.01 Residual 0 . 1‘1“.“ 11"l'u‘1 "1‘" 'i’i' V -0.01 ' 1 ' 1 ' -0.02 0 0.2 0.4 0.6 0.8 1 1.2 Time (s) Figure 6.9: Residual absorbance values using 1 pixel at 495 nm for a reaction of 1.39 x 10'5 M neodymium with respect to time. 154 2. Two- and Three-Component Studies All three possible two-component mixtures of the three lanthanides were prepared, as well as three component mixtures. These mixtures were then reacted with the PAR reagent in the stopped flow apparatus and the reactions followed spectrophotometrically as described above. Conditions for the data collection are described in Table 6.2 for each of the mixtures studied. The number of wavelengths necessary to achieve accurate estimates of species concentrations was tested. For greater than one wavelength used, the choice of which wavelengths to utilize for the analysis was made by a simplex optimization procedure7»8. Once these wavelengths were determined, they were kept constant for all remaining analyses of that mixture. Tables 6.3 to 6.5 summarize the final estimated concentrations from the extended Kalman filter for each of the possible two-component mixtures using from one to three wavelengths. Table 6.6 summarizes the final estimated concentrations for a three component mixture using from one to four wavelengths. Note that in general, the more wavelengths used for the analysis, the more accurate the results from the extended Kalman filter. One notable exception to this is in the mixture of praseodymium and neodymium where the increase from two to three wavelengths actually decreases the accuracy of the determination. However, the precision of the determination always increases with more wavelengths being utilized, and the differences for this determination using two or three wavelengths are within the relative standard deviations of the estimates. Another interesting Tabi com llit‘ 1&ij PI] 155 Table 6.2: Conditions for spectral acquisition for two and three component mixtures Mixture Time between Delay before Total data scans (ms) first scan (ms) collection time (m5) Ia/Pr 8.5 13.0 863 Ia/Nd 8.5 13.0 863 Pr/ Nd 6.5 1 1 .0 66 1 La/Pr/Nd 6.5 1 1.0 66 1 2.00 x 10‘4 M PAR used for all reactions iii for 156 Table 6.3: Final extended Kalman filter estimates of concentration for mixture of lanthanum and praseodymium. La E_I r Iaken Found Error Taken Bound Error ME. LIQSM LIQSM 4%) LIQSM [J.Q'SM Mi 1 1.44 0.82 —43.0 1.42 1.49 4.9 2 1.44 1.21 -16.0 1.42 1.27 -10.6 3 1.44 1.29 -10.4 1.42 1.42 0.0 157 Table 6.4: Final extended Kalman filter estimates of concentration for mixture of lanthanum and neodymium. 44% SQ 8 18km Eormd Error Taken Eormd Error M 21051! LlQ‘SM 21%) LID-5M 2105M M 1 1.44 0.74 -48.4 1.39 1.45 4.3 2 1.44 1.24 -13.7 1.39 1.35 ~2.8 3 1.44 1.26 -12.5 1.39 1.37 -1.4 158 BI Table 6.5: Final extended Kalman filter estimates of concentration for mixture of praseodymium and neodymium. IakenEoundError 2105M LLQ'SM M 1.06 -0.39 -137 1.06 1.06 0.0 1.06 1.17 10.2 NJ Taken Found Error 2105M 2105M 2126) 1.04 1.67 60.6 1.04 1.08 4.0 1.04 1.16 11.8 159 Table 6.6: Final extended Kalman filter estimates of concentration for mixture of lanthanum, praseodymium, and neodymium. All; fl N4— 1 Iaken Eound Error Iaken mung Error Taken Eounsi Error M 210511 LIQ‘SM 2061 2105M 2105M 21261 2105M 0.0-SM 21416; 1 1.44 0.67 -53.3 1.06 23.6 2126 1.04 -4.69 -551 2 1.44 -0.49 -134 1.06 3.05 188 1.04 1.41 35.6 3 1.44 0.81 -43.5 1.06 1.99 87.7 1.04 0.98 -6.3 4 1.44 1.24 -13.9 1.06 0.87 -17.9 1.04 0.97 -6.7 160 result is seen in Table 6.4 for the analysis of lanthanum and neodymium. The estimated concentration of neodymium using only one wavelength appears to be quite accurate. However, the relative standard deviation for this number from the covariance matrix is approximately 80%. It is generally advisable to use three wavelengths to determine two-component mixtures, and four wavelengths to determine three component mixtures. Under these conditions, the relative standard deviations from the covariance matrix were generally less than 5%. Increasing the number of wavelengths used beyond these three or four will yield an improvement in the precision of the determination, but will also take longer to analyze. Relative standard deviations from multiple runs of samples were generally less than 10% for two-component mixtures and less than 15% for three-component mixtures. Different concentrations of each of the species were tested to see the range over which the method is applicable to this system. Table 6.7 summarizes the results of these studies. The praseodymium/ neodymium system was tested the most thoroughly, and excellent results were obtained for concentration ratios from 5:1 to 1:2. At ratios higher than these, excellent results were obtained for the more concentrated species, while the errors for the less concentrated species increased. This is most likely due to the limits of detection for the chemical system, and not a flaw in the method. Simulations presented in Chapter IV show excellent results for broader concentration ratios. As a comparison to other methods for multicomponent kinetic determinations, the concentrations of praseodymium and 161 Table 6.7: Final extended Kalman filter estimates of concentration for mixtures of lanthanum, praseodymium, and neodymium. La if; M taken £61m Error Taken Bound Error Iaken Eound Error LIQ‘SM 2105M 2061 210% LID‘SM 21261 LIQ'SM LLQ’SM 21261 1.44 1.29 -10.4 1.42 1.42 0 - - - 1.50 1.09 -27.3 3.00 4.22 40.7 - - - 1.44 1.26 -12.5 - - - 1.39 1.37 -1.40 - - - 1.06 1.06 0 1.04 1.08 4.0 - - - 1.50 1.72 14.7 0.75 0.72 -3.5 - - - 1.50 1.40 -6.7 0.30 0.29 -4.0 - - - 0.750 0.756 0.8 1.50 1.35 -10.0 - - - 1.50 1.31 -12.7 0.15 0.25 67.0 1.44 1.24 -13.9 1.06 0.87 -17.9 1.04 0.97 -6.7 1.44 1.58 9.9 1.42 1.61 13.7 1.39 1.51 8.9 162 neodymium in their two-component mixture was determined by the method of proportional equations. This method was modified to be used at the two wavelengths used for the extended Kalman filter analysis shown in Table 6.5. Table 6.8 shows a comparison of the two methods. As can be seen in the table, the extended Kalman filter shows a vast improvement over the method of proportional equations. The rate constants provided to the filter for two—component mixtures were varied from 50% to 150% of the true value for each of the components. The worst estimates of concentration occurred when both rate constants were below the true values. The maximum change in the errors in concentration for the mixture of praseodymium and neodymium described in Table 6.5 using two wavelengths was 6% when both rate constants were underestimated significantly. The maximum error for either component under these conditions was 9.6%, compared to 4.0% when accurate rate constants for both reactions are provided to the filter. Residual absorbance values at each of the wavelengths in multicomponent mixtures are similar to those shown in Figure 6.9. The residuals are random and centered about zero. The number of iterations that it takes the filter to settle onto accurate estimates is longer in multicomponent mixtures than it is in single component samples. Figure 6.10 shows the progress of the estimated concentrations of praseodymium and neodymium with time. Again, however, once the estimated concentrations have settled, they remain quite stable. 163 Table 6.8: Comparison of the method of proportional equations and the extended Kalman filter. 1?: N0 Method Taken Found Error IakenEound Error 2105M 21058 426) 2105M [.LQ‘SM 4g) Prop. Eqn. 1.06 2.25 112.3 1.04 -0.07 -106.7 Extended 1.06 1.06 0.0 1.04 1.08 4.0 164 2.510“5 0 2.010‘5 31.5105 0 1.010‘5 5.010’6 . 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Time (s) Figure 6.10: Estimated concentrations of praseodymium (o) and neodymium (+) versus time using two wavelengths. True concentrations: [Pr]=l.06 x 10'5 M, [Nd]=1.04 x 10-5 M, [PAR]=2.00 x 10‘4 M. 165 D. Conclusions The extended Kalman filter has been shown to be a powerful technique for multicomponent kinetic determinations. There are several advantages to this method over most previous approaches. Determinations are based on differences in both spectral and kinetic behavior, which allow differentiation of more closely related species. The recursive nature of the Kalman filter provides fast computations. Finally, the use of the nonlinear, extended form of the Kalman filter avoids the most common assumption used in other methods, which is that the rate constant has been required to be invariant from run to run. This method does not require that this assumption be valid. Previous multicomponent kinetic methods require strict control of reaction conditions, such as pH, ionic strength, and temperature. These controls were not used in the determinations described in this chapter. ' This chapter describes the analysis of mixtures of lanthanides reacting with a common complexing agent by the extended Kalman filter. The filter was shown to be applicable to systems with highly overlapped kinetic and spectral responses. A mixture of three reacting metals was successfully analyzed even though the ratio of their pseudo-first-order rate constants was less than 2:1. A further application of this work would be to use a greater number of reacting components, and to use parallel Kalman filter networks for the data analysis. By using this method to analyze data, both the number of reacting metals, and their identities could be determined. Unknown 166 samples could be qualitatively and quantitatively analyzed with very little effort on the end-users part. This area would also be very amenable to the use of faster processors and to the so-called Kalman filter on a chip. Since the amount of data processing necessary increases greatly with an increasing number of parallel Kalman filters, these systems are numerically very intensive. With the increase in computational power that is sure to materialize, the amount of time necessary for data analysis using these complex models will not be as great as is currently the case. Another area where this method could be extended is to the variation of a rate constant within a run due to changes in temperature or pH during the reaction. This has already been done for single component systems by monitoring temperature as well as absorbance, and using an Arrhenius-type temperature dependence in the model9. This could easily be extended to a multiple components case. If the rate constant for a reaction is dependent on pH then simultaneously monitoring the pH of the solution as well as its absorbance and including a term in the model for the pH dependence could allow this correction to be included. The extended Kalman filter with data from multiple wavelengths has been shown to be a powerful data analysis tool for multicomponent kinetic studies. This is the first method for multicomponent kinetics which has not required that the rate constants be invariant from run-to-run. This has allowed a greater flexibility in the conditions used to perform the reactions. It is also not necessary to use such rigid controls on the reaction conditions as has previously been the case. This method has been shown be an 167 effective data analysis tool even if the kinetics or spectra of a system are too close to be determined by other methods. The technique has been shown to be more flexible, more powerful, and more accurate than previous methods for multicomponent kinetic determinations. LIST OF REFERENCES LIST OF REFERENCES 1 S. Shibata, Z-Pyn'dylazo Compounds in Analytical Chemistry, H.A. Flaschka and AJ. Barnard, Jr., eds., Chelates in Analytical Chemistry, Vol. IV, Marcel Dekker, New York, 1972, pp. 3-13, pp. 116-165. 2 O. Abollino, E. Mentasti, C. Sarzanini, and V. Porta, Analyst, 1991, 116, 1167-1170. 3 M. Tanaka, S. Funahashi, and K. Shirai, Anal. Chim. Acta, 1967, 39, 43 7-445. 4 A. Cladera, E. Gomez, J .M. Estela, V. Cerda, and J .L. Cerda, Anal. Chim. Acta, 1993, 272, 339-344. 5 A. Cladera, E. Gomez, J .M. Estela, and V. Gerda, Anal. Chem., 1993, 65, 707-715. 6 PM. Beckwith and S.R. Crouch, Anal. Chem., 1972, 44, 221-227. 7 D. Betteridge, A.P. Wade, and A.G. Howard, Talanta, 1985, 32, 709- 722. 8 D. Betteridge, A.P. Wade, and A.G. Howard, Talanta, 1985, 32, 723- 734. 9 CA. Corcoran and S.C. Rutan, Anal. Chem., 1988, 60, 1146-1153. 168 APPENDDI I sir. res COI Appendix I Program NRNDAT The sciences do not try to explain, they hardly even try to interpret, they mainly make models. By a model is meant a mathematical construct which, with the addition of certain verbal interpretations, describes observed phenomena. The justification of such a mathematical construct is solely and precisely that it is expected to work. John Von Neumann This appendix describes the program NRNDAT, which was written to simulate multicomponent kinetics data. The system that this program is modeled after is described by equation 5.1, the simultaneous, irreversible reaction of n components with a common reagent to form 11 different products. The actual number of components and number of wavelengths to be used for a particular run is input by the user. The user also inputs the rate constants for each of the reactions, and the names of the files containing the molar absorptivities for each of the components. The program then correctly dimensions the necessary arrays to hold the data. A full second-order irreversible kinetic model is used to generate the concentration of each species at each point in time. Since this model contains an interdependency between the concentrations of the common reagent and each of the reactants, the program iterates the calculation of each of these concentrations at each point in time until 169 it [0 wk filt 170 the difference in the calculated concentration of each of the components is less than 0.001% between two successive iterations. The program then uses the molar absorptivity files to obtain the true absorbance for each component for each time interval. The absorbance due to all of the components is then summed, and the total true absorbance for each wavelength at each time is obtained. As it is written, the program simulates absorbances for 100 time intervals each separated by 0.05 time units. The true absorbances are then passed to the subroutine ADDNS2 which adds randomly distributed noise to the true absorbance values1’2»3. The standard deviation of the noise to be added is set by the variable STD, and has the units of absorbance units. The resulting simulated absorbances for each point in time and for each wavelength are then printed out to a file specified by the user. These files are of the form: TIME, ABS(1), ABS(Z), ..., ABS(N) where ABS(n) is the absorbance at the nth wavelength. These files are then in the form necessary to be read by the extended Kalman filter program described in Appendix II. 171 DECLARE SUB ADDNS2 (IDIMll, IDZI, SIGMAl, BCKGRDI, ISEEDl, IMODE!) PROGRAM NRNDAT.BAS VERSION 1.0 This program was written in QuickBASIC 4.0 This version uses the full second order model to determine []'s and absorbance values. This program generates kinetic data for the system: C1 + R -> Pl rate constant k(l) C2 + R -> P2 rate constant k(2) Cu + R --> Pn rate constant k(n) at a user specified number of wavelengths. The program prompts the user for values for the rate constant, the initial absorbance of each reactant, and the molar absorptivities of each reactant and each product at each wavelength. The program outputs 100 time, absorbance at wavelength 1, absorbance at wavelength 2, etc. sets. A random noise term is added to each absorbance to simulate data that would be obtained in a real experiment. Written by: Brett Quencer Dept. of Chemistry Michigan State University East Lansing, MI 48824 Date completed: Feb. 23, 1993 Modified: Subroutine ADDNSZ adapted from a (very) similar one written in FORTRAN-77 by Pete Wentzell. His version of the subroutine may be found on p. 170 of his Ph.D. dissertation (MSU, 1987). DEFDBL T CIS INPUT "# OF COMPONENTS"; NUMCOMP NUMSTATES = (2 * NUMCOMP) + 1 INPUT "# OF WAVELENGTHS"; NWAVE 'NWAVE = 3 'INPUT "ENTER INTERVAL TO USE"; INTER INTER = .05 'INPUT "ENTER # POINTS"; PTS .I’I‘S = 100 REDIM BACK(NWAVE), AFIN(NWAVE, I’I‘S) AS DOUBLE REDIM EPS(NUMSTATES, NWAVE) AS DOUBLE REDIM K(NUMCOMP) AS DOUBLE, CONC(NUMCOMP + 1) AS DOUBLE REDIM TCONC(NUMSTATIB, P'IS) AS DOUBLE 172 REDIM RCONC(NUMSTATES) AS DOUBLE REDIM CONC1(NUMCOMP, PTS) AS DOUBLE, T(NUMCOMP * 2) AS DOUBLE REDIM XX(NUMCOMP + 1) AS DOUBLE REDIM TEMP(NUMSTATI~S) AS DOUBLE REDIM ABSV(NUMSTATES) AS STRING Prompt user for wavelengths, initial concentrations of components, rate constants for each component, molar absorptivities of each component, and the background absorbance at each wavelength. CLS FOR I = 1 TO NUMCOMP PRINT "INITIAL CONCENTRATION OF COMPONENT C("; l; " "; INPUT CONC(I) NEXT 1 INPUT "INITIAL CONCENTRATION OF COMPONENT R "; CONC(NUMCOMP + 1) FOR I = 1 TO NUMCOMP PRINT "RATE CONSTANT FOR REACTION "; 1; INPUT 1((1) NEXT I FOR I = I TO NUMCOMP PRINT "Enter molar absorptivity file for C("; I; ")"; INPUT ABSV(I) NEXT I PRINT "Enter molar absorptivity file for R"; INPUT ABSV(NUMCOMP + 1) FOR I = 1 TO NUMCOMP PRINT "Enter molar absorptivity file for P("; I; ")"; INPUT ABSV(NUMCOMP + 1+ 1) NEXT 1 FOR I = 1 TO NUMSTATB FILNUM = 30 + 1 OPEN ABSV(I) FOR INPUT AS #FILNUM NEXT I FOR I = 1 TO NWAVE FOR] = 1 TO NUMSTATES FILNUM = 30 +1 INPUT #FILNUM, EPSU, 1) NEXT 1 NEXT I FOR I = 1 TO NUMSTATES FILNUM = 30 + I CLOSE #FILNUM NEXT I ' Determine concentration of each component. Then determine the absorbance at each wavelength and time. (Time is in arbitrary units, equally spaced), and find the absorbances at each wavelength. Need to loop the calculation of concentration at each time, since each component (A, B, and C) is dependent on the concentration of R at time T, which is in turn dependent on the concentration of A, B, & C. The calculations are iterative until there is less than .01% change in ' the concentration of each component at that time. FORI= l TONUMCOMP+1 RCONC(I) = CONC(I) 'Starting concentrations 173 'are set to original NEXT I 'concentrations. FOR I = 1 TO NUMCOMP T(I) = CONC(I) / CONC(NUMCOMP + 1) 'Simplifying equations NEXTI FOR I = 1 TO NUMCOMP T(NUMCOMP + I) = CONC(I) - CONC(NUMCOMP + 1) NEXT I FORJ=1TOPTS 'Loop forall times L = INTER * J 'Set time spacing to 0.1 time units FORI=1TONUMCOMP+1 42 : XX(I) = RCONC(I) 'Set for previous concentration NEXT I 'for iterations RCONC(NUMCOMP + l) = CONC(NUMCOMP + l) ' Determine concentrations K = NUMCOMP + 1 FOR I = I TO NUMCOMP RCONC(I) = (T(l) * RCONC(K)) * EXP(T(NUMCOMP + I) * 1((1) * L) NEXTI FOR I = I TO NUMCOMP RCONC(NUMCOMP + I + 1) = CONC(I) - RCONC(I) RCONC(K) = RCONC(K) - RCONC(NUMCOMP + I + 1) NEXT I FORI: ITO NUMCOMP+1 XX(I) = XX(I) / RCONC(I) 'Calculates ratios IF XX(I) > 1.00001 THEN 42 'Convergence test NEXT I FOR I = 1 TO NWAVE FOR K = 1 TO NUMSTATB TEMP(I() = EPS(I(, I) * RCONC(K) AFIN(I, J) = AFIN(I, J) + TEMP(K) NEXT K NEXTI FORI = 1 TO NUMCOMP CONC1(I, J) = CONC(I) - RCONC(I) 'Conc. of component C(I) at time t NEXTI NEXT J Now we have the pure (no noise) absorbance values. Call the subroutine ADDNSZ to generate random noise, and add to or subtract from the pure absorbance values. 'INPUT "ENTER STD. DEV."; STD STD = .001 'INPUT "ENTER SEED"; SE SE = RND(250) * 1000 CALL ADDNSZ(NWAVE, PIS, STD, 0, SEE, 1) ' Print results to output file CIS INPUT "NAME OF OUTPUT FILE FOR KALMAN FILTER"; GS OPEN GS FOR OUTPUT AS #1 FOR I = 1 TO P'IS PRINT #1, USING "##.### "; 1* INTER, FOR J = 1 TO NWAVE 174 PRINT #1, USING "#.##### "; AFIN(J, I), NEXT] PRINT #1, NEXT I CLOSE #1 END DEFSNG T SUB ADDNSZ (IDIMl, ID2, SIGMA, BCKGRD, ISEED, IMODE) [This subroutine was originally written in FORTRAN-77 by Pete Wentzell. It was adapted to QuickBASIC by Brett Quencer, and contains all of the comments of the original. New comments are delimited by square brackets] This subroutine adds random noise with a Gaussian distribution to the data in an array. The parameters involved are: DATA- is the array containing the data to which the error is to be added or the array which will contain the errors generated, depending on the setting of IMODE. [NOTE DATA is named AFIN here and is a shared array.] IDIM#- is the size of the data array [in the dimension in #]. SIGMA- is the desired standard deviation of the errors to be added (value is absolute). BCKGRD- is a constant background level to be added to all of the errors generated. ISEED- is the random number seed to be used in generating the random errors. IMODE- is the mode of error generation: IMODE=1 indicates errors generated are to be added to the data already in the array IMODE=2 indicates errors generated are to be put directly into the array (not added). ' To generate the normally distributed random errors, this subroutine ' uses the random number generator, which gives a rectangular distribution, ' in conjunction with the cumulative Gaussian distribution. The random ' number generated is used as the argument in an approximation to the ' inverse of the standardized cumulative normal distribution function (see ' M. Abramowitz and LA. Stegun, 'Handbook of Mathematical Functions' (1968), ' Equation 26.2.23). This method is briefly described by L Kaufman in ' 'Trends in Anal. Chem.', vol. 2, p. 244 ' Set up house SHARED AFIN() AS DOUBLE C0 = 2.515517 C1 = .802853 C2 = .010328 D1 = 1.432788 D2 = .189269 D3 = .001308 175 ' Zero data array if mode 2 is being used IF IMODE = 2 THEN FOR I = 1 TO IDIMl FOR J = 0 TO ID2 AFIN(I, J) = 0! NEXT J NEXT I END IF ' Generate errors and add to array ' Generate random number. Needs to be in range 0 0. Compute error, scale, and add background FORI=1TOIDIM1 FORJ = 1 TO ID2 P = RND(ISEED) IF P > .5 THEN SIGN = -II P = P - .5 ELSE SIGN = 1! IF P = 0! THEN P = .00001 END IF T = SQR(-2! * LOG(P)) XP=T-(C0+Cl*T+C2*T*T)/(l +D1*T+D2*T*T+D3*T*T) AFIN(I, J) = AFIN(I, J) + SIGN * SIGMA * XP + BCKGRD NEXT J NEXT I ' All finished, now go home END SUB LIST OF REFERENCES LIST OF REFERENCES 1 M. Abramowitz and LA. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, 1968, Ch. 26. 2 L. Kaufman, Trends in Anal. Chem., 1983, 2, 244. 3 RD. Wentzell, Ph.D. Dissertation, Michigan State University, 1986. 176 APPENDIX II WaVG Appendix II Program EKFNRN You think you know when you learn, are more sure when you can write, even more when you can teach, but certain when you can program. Alan J. Perlis This appendix contains the program listing for EKFNRN. The program was written to process multicomponent kinetics data using the extended Kalman filter. This particular version was written for the analysis of n simultaneously reacting components with a common reagent to form n different products following pseudo-first-order kinetics as described in equation 5.1. The program utilizes the extended Kalman filter algorithm described in equations 3.3 to 3.7. The measurement function for this system is given by equation 5.3. The user enters the number of components to determine as well as the estimated rate constants for their corresponding reactions. The user also enters the initial concentration of the common reagent, as this would be known, and the filenames for the files containing the molar absorptivities of each of the components. The program then dimensions all arrays according to the number of components and the number of wavelengths used. The filter then proceeds through the extended Kalman filter algorithm for each point in time. The absorbance at each of the wavelengths used is printed out to a file (name specified by user) for 177 SI fi pr VE c0 ex to star the 0f 1 Star the Star Out- file. Silm 178 each cycle of the filter. The estimated values and their corresponding standard deviations for all of the parameters in the state vector are printed to a second file for each loop through the filter. Finally, after all points have been processed, the program prints out the final estimates for each of the parameters in the state vector along with their standard deviations. The standard deviations are taken from the covariance matrix. The diagonal of this matrix contains the variances of the parameters of the state vector. For example, the value of the first row and the first column in the covariance matrix contains the variance in the first parameter in the state vector. The square root of this value is computed to achieve the standard deviation in that parameter. This is repeated for each of the elements in the state vector. Also printed to the screen at the end of the run are the standard error of the estimate for each of the wavelengths used, and the sum for all wavelengths used. Also, the inverse of the sum of the standard errors is also printed to the screen. Finally, the program outputs the time required to process all of the data from the input file. I This program has been utilized for systems containing anywhere from a single component to systems containing eight simultaneous reacting components. This latter case is described in Chapter V. ' ' ' - ' ' ' ' ' ' ' ' U ' ' ' F ' ' ' ' ' 5 ‘ ' ' ' ' ' 179 DECLARE SUB INV3 (M11, AINV#(), A#(), IK!(), JK!()) I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I PROGRAM EKFNRNBAS VERSION 1.0 This program was written in QuickBASIC v. 4.0 This version handles any number of wavelengths. This version takes into account the change in the pseudo—first-order rate constant because of the decrease in R(t). (K(1)=k(1)*R(t) This program employs the extended Kalman filter algorithm and the Potter-Schmidt square root algorithm to determine the starting concentrations of reactants and the rate constant in the kinetic system defined as: C1 + R -> Pl , rate constant k(l) C2 + R -> P2 , rate constant k(2) Cu + R -—> Pn , rate constant k(n) where all components absorb. It is assumed that the system is pseudo- first—order in C(1), C(2), ..., C(n). Parameters required to be input by the user are: -name of input file -name of output file -number of component (NUMCOMP) -number of wavelengths (Ml) -molar absorptivity of each component at each wavelength(EPS(2*NUMCOMP+ l , M1)) —initial guess for the concentration of C(1), C(2),..., C(n), & R. (X(1.1).X(2.l).---.X(n.l).X(n+1.1)) -initial guess for the rate constants (X(n+2,1),X(n+3,1),...,X(2n+l,l) -initial guess for the error covariance matrix (P(2n+1,2n+l)) -variance in data at each wavelength (R(M1,M1)) The input file is expected to be in the form: TIME, ABS(1), ABS(Z), ..., ABS(M) where the number in parentheses is the number of the wavelength. The initial set up for the program will handle data provided at two wavelengths. To increase the number of wavelengths, a new array needs to be defined for each additional wavelength (Zx(COUNT)) where x is the wavelength number. Also, the INPUT #1 and INPUT #2 statements need to be changed to read in the correct number of variables. There will be two output files. The first will be in the form: TIMB(N). est [C(1)](0). est [C(2)](0).--.. est [C(n)](0). est k(l), est k(2),..., est k(n) The second will be in the form: 180 TIME(N), est Abs(wavelengthl )(N), est Abs(wavelength2)(N) The name of each output file and the input file can be specified by the user. The input files can contain any number of points, but must be in the form specified above. Written by: Dept. of Chemistry Michigan State University E. Lansing, MI 48824 ' Brett Quencer ' Date completed: Feb. 18, 1993 ' Modified: ' Find out what we need to know from the user: CLS INPUT "How many components do you want to test for"; NUMCOMP INPUT "What is the name of the input file "; AS 'INPUT "What do you want to name the output file for [] & k "; CS CS = "CONIC.OU'I" 'INPUT "What do you want to name the output file for abs.'s "; DS DS = "CON 1 OUT“ INPUT "How many wavelengths were used "; M1 'M1 = 3 NUMSTATES = (NUMCOMP * 2) + 1 ST: CIS ' Set up the arrays REDIM X(NUMSTATES + 1, 1) AS DOUBLE, P(NUMSTATES, NUMSTATES) AS DOUBLE REDIM H(Ml, NUMSTATES) AS DOUBLE REDIM EPS(NUMSTATI:S, M1) AS DOUBLE, PHALF(NUMSTATES, NUMSTATES) AS DOUBLE REDIM IDENT(NUMSTATES, NUMSTATES) AS DOUBLE, R(Ml, M1) AS DOUBLE REDIM GAIN(NUMSTATES, M1) AS DOUBLE REDIM YA(NUMSTATES), YB(NUMSTATES) REDIM ABSV(NUMSTATES) AS STRING DEFSNG Y CIS FOR I = 1 TO NUMCOMP PRINT "Enter molar absorptivity file for C("; I; ")"; INPUT ABSV(I) NEXT I PRINT "Enter molar absorptivity file for R"; INPUT ABSV(NUMCOMP + 1) FOR I = 1 TO NUMCOMP PRINT "Enter molar absorptivity file for P("; I; ")"; INPUT ABSV(NUMCOMP + I + 1) NEXT I INPUT "Enter filename for wavelength variances"; VARS FOR I = 1 TO NUMSTATFS FILNUM = 30 + I OPEN ABSV(I) FOR INPUT AS #FILNUM NEXT 1 \l.£rr. EFF. DP RCF Nmrru Pl 1 F. N U ' ' - ' V ' 181 OPEN VARS FOR INPUT AS #21 FOR I = 1 TO Ml FORJ = 1 TO NUMSTATES FILNUM = 30 + J INPUT #FILNUM, EPS(J, I) NEXT J INPUT #21, R(I, 1) NEXT 1 FORI=1TONUMSTATFB FILNUM = 30 + I CLOSE #FILNUM NEXT I CLOSE #21 FOR I = 1 TO NUMCOMP X(I, l) = 0 NEXT I 20 : INPUT "What is the initial estimate of R0"; X(NUMCOMP + l, 1) FOR I = I TO NUMCOMP PRINT "What is the initial estimate of k("; I; ")"; INPUT X(NUMCOMP + I + l, 1) NEXT I PINIT = 50 FOR] = 1 TO NUMSTATB PRINT X(J, 1) NEXT J PRINT PINIT PRINT "Are the above values correct? (Y/N)" INPUT W5 IF W$ = "N" OR WS = "n" THEN 20 R0 = X(NUMCOMP + l, 1) 'Set R0 to be initial conc of reagent CLS 'for use in later calculations FORI=1TONUMSTATFB PRINT X(I, 1) NEXT I PRINT: PRINT LOCATE 12, 25, 0 PRINT "CALCULATING, PLEASE WAIT . . ." ' Set up covariance array. Also, define the identity matrix. FOR I = I TO NUMSTATES FOR] = 1 TO NUMSTATES IF I = J THEN P(I, J) = PINIT ELSE P(I, J) = 0! IF I = J THEN IDENT(I, J) = I! ELSE IDENT(I, J) = 0! NEXT J NEXT I ' Calculate the upper triangular square root matrix of the covariance ' matrix. This will be used for propagation of the covariance ' throughout the filtering. This is the main change for the Potter- ' Schmidt square root algorithm. The formulas for calculating the ' upper triangular square root matrix are found in the book: ' V.N. Faddeeva, "Computational Methods of Linear Algebra", Dover ' Publications, NY, 1959, pp. 81-4. FORI= ITO NUMSTATES PHALFU. 1) = SQR(P(I. 1)) NEXT I (538%? .5. nggan‘ég< 9.2.. 182 ' Find out number of data points REDIM DIFF(M1) AS DOUBLE, TSUM(M1) AS DOUBLE, RAVE(M1) AS DOUBLE P = M1 + 1 REDIM DUM(P) COUNT = 0 OPEN AS FOR INPUT AS #1 WHILE NOT EOF( 1) FOR I = 1 TO P INPUT #1, DUM(I) NEXT I COUNT = COUNT + 1 WEND CLOSE #1 ' Determine constant sums that will be used in calculating the linearized ' measurement function vector. REDIM TEMPSUM(Z * NUMCOMP, M1) AS DOUBLE FOR I = 1 TO M1 FOR] = 1 TO NUMCOMP TEMPSUM(J, I) = EPS(J, I) + EPS(NUMCOMP + 1, I) - EPS(NUMCOMP + J + l, I) TEMPSUM(NUMCOMP + J, I) = EPS(NUMCOMP + J + l, I) - EPS(NUMCOMP + 1, 1) NEXT J DIFF(I) = 01 NEXT I ' Set up all the arrays REDIM TIM(COUNT) AS DOUBLE, Z(Ml, COUNT) AS DOUBLE, RESID(M1, COUNT) DEFDBL T ' Open the data file and start filtering OPEN AS FOR INPUT AS #2 OPEN CS FOR OUTPUT AS #5 OPEN DS FOR OUTPUT AS #4 PRINT #5, " Time "; FOR I = 1 TO NUMCOMP PRINT #5, ” est C("; I; ") "; " Std Dev"; NEXT I PRINT #5, " est Ro "; " Std Dev "; FOR I = I TO NUMCOMP PRINT #5, " EST k"; I; " Std Dev "; NEXT I PRINT #5, IT] = TIMER FOR F = 1 T O COUNT INPUT #2, TIM(F) FOR I = 1 TO Ml INPUT #2, Z(I, F) NEXT I m********************************************************** “*‘k** ' Calculate the linearized measurement function vector (H(M1,NUMSTATES)) ' But first calculate R(t) (X(NUMSTATFS+1,1)) REDIM TMP(NUMCOMP) AS DOUBLE, T(NUMSTATES) AS DOUBLE X(NUMSTATFS + 1, 1) = R0 FOR I = 1 TO NUMCOMP K = NUMCOMP + I + l 183 L = NUMSTATES + 1 TMPU) = (X(I, 1) -X(I. 1) * EXP(-X(K. 1)* X(L, 1) * TIM(F))) X(L. 1) = X0. 1) - TMP(I) NEXT I FOR I = 1 TO Ml FOR] = 1 TO NUMCOMP T0) = (TEMPSUM(J, I) * (EXP(-X(NUMSTATFS + l, l) * X(l + NUMCOMP + J, l) * TIM(F)))) 110, J) = T(J) + TEMPSUM(NUMCOMP + J, I) T(NUMCOMP + J + 1) = (-X(NUMSTATFS + 1, 1) * TIM(F)) * (TEMPSUM(J, 1)) H0, NUMCOMP + J + 1) = T(NUMCOMP + J + l) * X(J, 1) * (EXP(-X(NUMSTATES + l, l) * X(J + NUMCOMP +1,1)* TIM(F))) NEXT J H(l, NUMCOMP + l) = EPS(NUMCOMP + l, 1) NEXT 1 From the algorithm, the state estimate extrapolation is equal to the last update. Also, the error covariance extrapolation is equal to the last update plus the system covariance matrix. But, since this is a static system, the system covariance matrix is equal to 0, so the error covar. extrapolation is just equal to the last update. Therefore, there are no updates performed here. *****m**m********************** ************************** *‘k**** ' Now calculate the Kalman gain ' Now zero the temporary arrays REDIM TRANSP(NUMSTATB, M1) AS DOUBLE, TEMP1(NUMSTATES, M1) AS DOUBLE REDIM TEMP2(M1, Ml) AS DOUBLE REDIM T EMP3(M1, M1) AS DOUBLE, TEMP4(NUMSTATES, M1) AS DOUBLE REDIM IK(M1), ]K(Ml) FOR I = 1 TO Ml FORJ = 1 TO Ml TEMP2(I, J) = 0! NEXT J NEXT I ' Now find the transpose of the linearized measurement function vector ' H(M1,NUMSTATB), and zero TEMP1(), & GAIN() FOR I = 1 TO NUMSTATB FORJ = 1 TO Ml TRANSP(I. l) = HO. 1) TEMP1(I, J) = 0# GAIN(I, J) = 0# NEXT] NEXT I ' Multiply PHALF() by TRANSPO, and store in TEMP1(). FOR I = 1 TO NUMSTATES FOR] = 1 TO Ml FOR K = I TO NUMSTATES TEMP1(I, J) = TEMP1(I, J) + (PHALF(I, K) * TRANSP(K, J)) NEXT K NEXT] 184 NEXT I ' Next, multiply the linearized measurement function vector H() by ' TEMP1(), and store in TEMP2(). FOR I = 1 TO Ml FORJ = 1 TO M1 FOR K = 1 TO NUMSTATES TEMP2(I, J) = TEMP2(I, J) + (H(I, K) * TEMP1(K, J)) NEXT K NEXT J NEXT I ' Now add R() to TEMP2() and put the result back into TEMP2() FOR I = 1 TO Ml FORJ = 1 TO M1 TEMP2(I, J) = TEMP2(I, J) + R(I, J) NEXT J NEXT I ' Now invert TEMP2() and put into TEMP3(). CALL INV3(M1, TEMP3(), TEMP2(), IK(), JK()) ' Now multiply TEMP1() by TEMP3(), and this is GAIN(). FOR I = 1 TO NUMSTATES FORJ = I TO M1 FOR K = 1 TO Ml GAIN(I, J) = GAIN(I, J) + (TEMP1(I, K) * TEMP3(K, J)) NEXT K NEXT] NEXT 1 ************************************************************* ****** ' Now update the state estimate. ' Zero temporary arrays. REDIM TEMP21(M1, 1)AS DOUBLE, TEMP22(M1, 1) AS DOUBLE REDIM TEMP23(NUMSTATES, 1) AS DOUBLE, TMl(NUMSTATFB, M1) AS DOUBLE FOR I = 1 TO M1 TEMP21(I, 1) = 0! TEMP22(I, l) = 01 NEXT I FOR I = 1 TO NUMSTATES TEMP23(I, l) = 0! NEXT 1 ' Determine the estimated absorbance at each wavelength. FOR I = 1 TO Ml FORJ = 1 TO NUMCOMP K = l + NUMCOMP + J L = NUMSTATES + 1 TMl(J, 1) - TEMPSUM(J. 1) * X(J. 1) * (EXP(-X(L, 1) * X(K. 1) * TIM(F))) TMl(K - l, I) = TEMPSUM(K - l, I) * X(J,1) NEXT J . TMl(NUMSTATFB, I) = EPS(NUMCOMP + l, I) * X(NUMCOMP + 1, 1) FOR J = 1 TO NUMSTATES TEMP21(I, l) = TEMP21(I, 1) + TMl(J, I) NEXT J NEXT I ' Subtract TEMP21() from Z(), and store in TEMP22(). 185 FOR I = 1 TO M1 TEMP22(I, l) = Z(I, F) - TEMP21(I, 1) RFSID(I, F) = TEMP22(I, 1) TSUM(I) = TSUMU) + TEMP22(I, 1) NEXT I ' Multiply the Kalman gain GAIN() by TEMP22() and store in TEMP23(). FOR I = ITO NUMSTATES FOR K = 1 TO Ml TEMP23(I, l) = TEMP23(I, l) + (GAIN(I, K) * TEMP22(K, 1)) NEXT K NEXTI ' Add TEMP23() to the previous state estimate to get the state estimate ' update X(). FOR I = 1 TO NUMSTATFB X(I, l) = X(I, 1)+ TEMP23(I, 1) NEXTI '************************************************************* ******* ' Now update the covariance matrix ' First, set up and clear temporary arrays. REDIM TEMP01(NUMSTATES, NUMSTATES) AS DOUBLE, TEMP02(NUMSTATES, NUMSTATES) AS DOUBLE REDIM TEMP03(NUMSTATES, NUMSTATES) AS DOUBLE, TEMP04(NUMSTATES, NUMSTATIS) AS DOUBLE REDIM TEMP05(NUMSTATES, NUMSTATES) AS DOUBLE, TEMP10(NUMSTATES, Ml) AS DOUBLE REDIM TEMP11(M1, NUMSTATES) AS DOUBLE, TEMP12(NUMSTATES, NUMSTATES) AS DOUBLE FOR I = 1 TO NUMSTATES FORJ: ITO NUMSTATES TEMP01(I, J) = 0! TEMP04(I, J) = 0! TEMP05(I, J) = 0! NEXT J FOR K = 1 TO Ml TEMP10(I, K) = 01 TEMPl l(K, I) = 01 NEXT K NEXT I ' Multiply the gain by the linearized measurement function. Store in ' TEMP01(). FOR I = 1 TO NUMSTATES FORJ = 1 TO NUMSTATES FOR K = I TO Ml TEMP01(I. J) = TEMP01(I. l) + (GAIN(I, K) * H(K, J)) NEXT K NEXT J NEXT I ' Subtract this from the identity matrix, and store in TEMP02(). ' Also, transpose TEMP02() and store in TEMP03(). FOR I = 1 TO NUMSTATES FORJ = 1 TO NUMSTATES TEMP02(I, J) = IDENT(I, J) - TEMP01(I, J) TEMP03(J, I) = TEMP02(I, J) 186 NEXT J NEXT I ' Multiply TEMP02() by the previous covariance matrix and store in ' TEMP04() FOR I = 1 TO NUMSTATES FORJ = I TO NUMSTATES FOR K = 1 TO NUMSTATES TEMP04(I, J) = TEMP04(I, J) + (TEMP02(I, K) * PHALF(K, J)) NEXT K NEXT J NEXT I ' Multiply TEMP04() by TEMP03(), and put the result into TEMP05(). FORI= ITO NUMSTATB FORJ = 1 TO NUMSTATES FOR K = 1 TO NUMSTATES TEMP05(I, J) = TEMP05(I, J) + (TEMP04(I, K) * TEMP03(K, J)) NEXT K NEXT J NEXT I ' Multiply the Kalman gain GAIN () by the measurement variance R(), ' and store in TEMP10(). Transpose the gain, and store in TEMP11(). FORI=1TONUMSTATES FORJ = 1 TO Ml FOR K = l T 0 M1 TEMP10(I, J) = TEMP10(I, J) + (GAIN(I, K) * R(K, J)) NEXT K TEMP11(J, I) = GAIN(I, J) NEXT J NEXT I ' Multiply TEMPIOO by TEMP11() and store in TEMP12() FOR I = 1 TO NUMSTATES FORJ = 1 TO NUMSTATES FOR K = 1 TO Ml TEMP12(I, J) = TEMP12(I, J) + (TEMP10(I, K) * TEMP11(K, J)) NEXT K NEXT J NEXT I ' Add TEMP05() and TEMPIZO together to get the new covariance update. FOR I = 1 TO NUMSTATES FORJ = 1 TO NUMSTATES PHALF(I, J) = TEMP05(I, J) + TEMP12(I, J) NEXT J NEXT I Imm******************************************************* m**** ' Now print out the results. PRINT #5, USING "##.### "; TIM(F); PRINT #4, USING "##.### "; TIM(F); FOR I = 1 TO NUMSTATES YA(I) = X(I, l) YB(1) = SQR(PHALF(I. 1)) NEXT 1 FOR I = I TO NUMSTATES PRINT #5, USING "#.####NVV\ "; YA(I); 187 PRINT #5, USING "##.###AAAA "; YB(I); NEXT I PRINT #5, ' TEMP21( ) gives the predicted absorbance at each wavlength. FOR I = 1 TO Ml PRINT #4, USING "##.###### "; TEMP21(I, 1), NEXT 1 PRINT #4, NEXT F F = F - 1 ' Because F=101, we need to subtract 1. 1001 : 1T2 = TIMER CLOSE #2 CLOSE #5 CLOSE #4 ' Determine standard deviations of residuals. REDIM SD(M1) AS DOUBLE FOR I = 1 TO Ml RAVE(I) = TSUM(1) / COUNT s0(1)= 01 NEXTI FOR1= 1 TO COUNT FORJ = ITO M1 80(1) = SDU) + ((RESIDU. I) - RAVE(J)) A 2) NEXT] NEXTI FORI: 1 TO M1 SD“) = SQR(SD(I) / (F- 1)) SSD = SSD + SD(I) PRINT "ST D. ERROR OF THE ESTIMATE FOR WAVELENGTH "; I; " = "; SD(I) NEXT 1 PRINT ' Print final estimates. PRINT "NUMBER OF POINTS: "; F FOR I = 1 TO NUMCOMP PRINT "C("; I; ")o=";X(1, l); " +/— "; SQR(PHALF(I, 1)) NEXT I PRINT "Ro="; X(NUMCOMP + 1, 1); " +/- "; SQR(PHALF(NUMCOMP + 1, NUMCOMP + 1)) FOR I = 1 TO NUMCOMP J=NUMCOMP+I+1 PRINT "k("; I; ")=": X(J. 1): " +/- "; SQR(PHALF(J. 1)) NEXT 1 PRINT PRINT "SUM OF STD. ERRORS = "; SSD PRINT "INVERSE = "; 1 / SSD PRINT PRINT "TIME ELAPSED= ”; TTZ - 'ITl END DEFDBL A, D, S SUB INV3 (M. AINV#(). A#(). IK(), JK()) 188 ' Based on the subroutine MATINV (originally written in FORTRAN), found ' on pp 302-3 of "Data Reduction and Error Analysis for the Physical ' Sciences", by Philip R. Bevington, McGraw-Hill, 1969 ' DET=determinant ' A#()=array to be inverted ' AINV#()=inverted array ' M=order of array ' N=number of parameters ' Set AINV#() to be equal to A() FOR I = 1 TO M FORJ = 1 TO M AINV#(I. l) = M0. 1) NEXT J NEXT 1 10 : DET = 1# ll:FORK= l TOM ' Find largest element in rest of AINV#(I,J) AMAX# = 0# 21 : FOR 1= K TO M FORJ = K TO M 23 : IF (ABS(AINV#(I, J)) >= ABS(AMAX#)) THEN 24 ELSE 30 24: AMAX# = AINV#(I, J) IK(K) = I JK(K) = l 30: NEXT] NEXT 1 ' Interchange rows and columns to put AMAX# in AINV#(K,K). 31 : IF (AMAX#) <> OTHEN 41 ELSE 32 32 : DET = 0# 41 : I = IK(K) 42:1F(1-K) <0THEN 21 IF(1-K)=0THEN51 1F(I-K)>0THEN 43 43:FORJ= l TOM SAV = AINV#(K, J) AINV#(K. J) = AINV#(I. l) 50: AINV#(I, J) = -SAV NEXT J 51 : J = JK(K) IF (J- K) <0THEN 21 IF(J—K)=0THEN61 IF(J-K)>0THEN S3 S3 : FOR I = 1 TO M SAV = AINV#(I, K) AINV#(I, K) = AINV#(I, J) 60: AINV#(I, J) = -SAV NEXT 1 ' Accumulate elements of inverse matrix 61: FORI: 1 TOM 1F(I- K) <> OTHEN 63 ELSE70 63 : AINV#(I, K) = -AINV#(I, K) / AMAX# 70: NEXT 1 71: FORI=1TOM 189 FORJ = 1 TO M IF(I-K) <>0THEN74ELSE80 74: IF(J-K)<>0THEN75EISE80 75 : AINV#(I, J) = AINV#(I, J) + AINV#(I, K) * AINV#(K, J) 80 : NEXTJ NEXTI 81 : FORJ = 1 TO M IF (J - K) <> 0 THEN 83 ELSE 90 83 : AINV#(K, J) = AINV#(K, J) / AMAX# 90: NEXT] 91: AINV#(K, K) = 1# / AMAX# 100: DET = DET * AMAX# NEXT K ' Restore ordering of matrix 101 : FORL= l TOM K = M - L + 1 J= IK(K) IF (J-K) <= OTHEN 111 ELSE 105 105: FORI: l TOM SAV = AINV#(I, K) AINV#(I, K) = —AINV#(I, J) 110: AINV#(I, J) = SAV NEXTI 111 : I= JK(K) IF (I- K) <= 0 THEN 130 EISE113 113: FOR]: 1 TOM SAV = AINV#(K, J) AINV#(K, l) = -AINV#(I. l) 120: AINV#(I, J) = SAV 130 : NEXT] NEXT L END SUB MICHIGAN S 1i111111111111