Anal. Chem. 1997, 69, 1909-1918
Potentiometric Nonlinear Multivariate Calibration with Genetic Algorithm and Simplex Optimization Margaret Hartnett† and Dermot Diamond*
School of Chemical Sciences, Dublin City University, Glasnevin, Dublin 9, Ireland
In this contribution, a genetic algorithm (GA) and a modified Simplex technique are investigated as a means of developing nonlinear multivariate calibration models for an array of ion-selective electrodes. The responses of an array of ammonium-, sodium-, potassium-, and calciumselective electrodes employed in a flow injection analysis system were modeled over the concentration range of 1 × 10-4 to 1 × 10-2 M using the GA and simplex techniques to optimize the cell potentials, slopes, and selectivity coefficient parameters of the Nikolskii-Eisenman equation for each electrode. Correlations between activities predicted from the calibration model and the actual activities of the solutions presented to the array ranged from 0.98 to 0.88 for the four ions.
sample, E°j is the standard cell potential of the electrode, and Sj is the slope of the electrode or the change in potential of the electrode per decade change in activity of the primary ion in the absence of interferents. The activity of the primary ion k in the ith sample is aik and the interfering ions l are given by ail. The term kpot jkl refers to the selectivity coefficient of the electrode against the lth interferent with respect to the primary ion, and 3k/3i the summation ∑ kpot represents the error arising from the jkl ail l*k
net contribution of all the interferents. A simplified form of the Nikolskii-Eisenman expression (described in expression 1) was used for modeling the responses of the ISEs in these studies. This simplified expression is described in
E ) E° + S log(ai + Methods for the determination of ammonium are becoming increasingly important as the understanding of this ion’s role in the environment grows. Among its many functions, it is known that ammonium acts as an important nutrient in estuaries and that increased concentrations of the ion in these waters can lead to problems associated with algal bloom. Consequently, methods for the determination of ammonium are finding increasing application in a range of areas, including the monitoring of water quality. A range of techniques are available for the determination of ammonium including gas-sensing ISEs,1 flow injection amperometric methods,2 and optical methods with evanescent waveguides.3 In this contribution, we consider the application of an array of sensors for determining ammonium and some other important cations, namely, sodium, potassium, and calcium. A nonlinear multivariate regression technique is developed which would in principle enable the dynamic recalibration of the sensors if necessitated by changing operating conditions (a potential problem which might occur if the sensor were used in the field). Specifically, we consider the application of genetic algorithms (GAs) to modeling the response of an array of ion-selective electrodes (ISEs) used in a flow injection analysis (FIA) system. The response of an ISE is described by the well-known Nikolskii-Eisenman expression
Eij ) E°j + Sj log(aik +
∑k
pot 3k/3i ) jkl ail
(1)
l*k
where Eij is the measured potential of the jth electrode to the ith † Current address: The Advanced Control Engineering Research Group, Department of Electrical and Electronic Engineering, The Queen’s University of Belfast, Belfast, U.K. (1) Arnold, M. A. Anal. Chim. Acta 1983, 154, 33. (2) Liu, R.; Sun, B.; Johns, I. Analyst 1995, 120, 2845. (3) Brandenburg, A.; Edelha¨user, R.; Werner, T.; He, H.; Wolfbeis, O. S. Mikrochim. Acta 1995, 121, 95.
S0003-2700(96)00489-1 CCC: $14.00
© 1997 American Chemical Society
∑k* a ) ij j
(2)
j
In this expression k*ij is a “conditional selectivity constant” which is related to the selectivity coefficient defined in eq 1 as follows: zi/zj k*ijaj ) kpot ij aj
(3)
Hence k*ij ) kijpot when zi ) zj. As discussed in ref 4, k*ij is not a global constant but is rather a constant within the constraints of the conditions of the experiment and calibration design. An example of the simplified Nikolskii-Eisenman expression for an ammonium ISE is depicted in eq 4.
EiNH4+ ) E°NH4+ + SNH4+ log(aiNH4+ + k* NH4+Na+aiNa+ + k* NH4+K+aiK + k* NH4+Ca2+aiCa2+) (4) In this expression, EiNH4+ refers to the potential (mV) response of the ISE to the ith sample. E°NH4+ refers to the standard cell potential of the ammonium ISE and SNH4+ refers to the slope of the ISE. aiNH4+, aiNa+, aiK+, and aiCa2+ refer to the activities of ammonium, sodium, potassium and calcium in the ith sample. * +Na+, kNH * +K+, and k*NH +Ca2+ refer to the conditional selectivity kNH 4 4 4 constant of the ammonium ISE against sodium, potassium, and calcium, respectively. The purpose of the calibration procedure is to determine the values of the cell potentials, slopes, and conditional selectivity coefficients of each electrode (e.g., for the ammonium ISE, * +Na+, kNH * +K+, and kNH * +Ca2+) that best describe E°iNH4+, SiNH4+, kNH 4 4 4 (4) Sa´ez deViteri, F. J.; Diamond, D. Analyst 1994, 119, 749. (5) Forster, R. J.; Regan, F.; Diamond, D. Anal. Chem. 1991, 63, 876.
Analytical Chemistry, Vol. 69, No. 10, May 15, 1997 1909
the primary ion/matrix response of the electrode over the calibration range chosen. This is achieved by using the parameters of and the Nikolskii-Eisenman expression for each electrode to predict the potential (mV) responses of the electrodes to a series of calibration solutions. The predicted potentials are then compared with the experimentally measured values for the potential (mV) responses of the electrodes to the same solutions. The extent to which the predicted and measured potentials differ is described by a sum squared error (SSE) term defined in eq 5. In this expression N, refers to the total number of samples N
SSE )
∑ s)1
( ) E-E ˜
2
E
(5)
used for calibration, E refers to the measured potential (mV) of an ISE to sample s, and E˜ refers to the predicted potential of the electrode based on the modified Nikolskii-Eisenman expression (see eq 4 for an example of this expression for the ammonium ISE). The objective of the optimization technique is thus to determine the values of the cell potential, slope, and conditional selectivity coefficient that minimize the SSE for the calibration solutions. Previous work of this kind had already been performed using simplex techniques.4-6 However, Betteridge et al.7 noted that the final result of a simplex is dependent on the estimates of the function variables with which the simplex is initialized. It was also clear that the size of the simplex could have quite a dramatic effect on the final result obtained as the complexity of the surfaces studied grows, due to increased possibility of local minima existing. With this problem in mind, it was decided to investigate the use of the GA since according to literature sources8 it is useful in high-dimensional search spaces and makes less assumptions about the search space than strong optimization methods like the simplex technique (for further details see the following discussion of genetic algorithms). Genetic Algorithms. While GAs have been used in a range of disciplines such as engineering9 and finance10 for some time, it is only recently that they have been applied to problems in analytical chemistry.11-13 GAs are a group of mathematical techniques that were initially designed to simulate the behavior of biologically based adaptive systems. These techniques were developed to study what kind of emergent behavior arose from a set of simple rules and how changes in the algorithm would affect this behavior. While the GA was not initially developed as a function optimizer itself, modifying the GA can produce powerful GA-based optimizers.14 (6) Forster, R. J.; Diamond, D. Anal. Chem. 1992, 64, 1721. (7) Betteridge, D.; Wade, A. P.; Howard, A. G. Talanta 1985, 32, 709. (8) Lucasius, C. B. Towards Genetic Algorithm Methodology in Chemometrics. Ph.D. Thesis, Katholieke Universiteit Nijmegen, Netherlands, 1993. (9) Orero, S. O.; Irving, M. A. Int. J. Elec. Power Energy Syst. 1996, 18, 19. (10) Allen, F.; Karjalainen, R. J. Finance 1995, 5013, 945. (11) Li, T.-H.; Lucasius, C. B.; Kateman, G. Anal. Chim. Acta 1992, 268, 123. (12) Wienke, D.; Lucasius, C. B.; Buydens, L; Kateman, G. Anal. Chim. Acta 1993, 271, 253. (13) Wehrens, R.; Lucasius, C. B.; Buydens, L; Kateman, G. Anal. Chim. Acta 1993, 277, 313. (14) De Jong, K. A. Are Genetic Algorithms Function Optimizers? In Parallel Problem Solving From Nature 2; Ma¨nner, R., Manderick, B., Ed.; Elsevier Publishers B. V.: Amsterdam, 1992.
1910
Analytical Chemistry, Vol. 69, No. 10, May 15, 1997
Genetic algorithms behave in a fashion similar to the class of moderate search techniques. These techniques make fewer assumptions about the response landscape than strong methods and are less computationally demanding than weak methods. Strong methods15 such as gradient search and simplex techniques use heuristics which concentrate on local areas of the response surface, by searching for an optimum, in the neighborhood of a particular point. These methods make assumptions about the response surface (e.g., the existence of a derivative for gradient search) to aid the exploration. These two characteristics of local searching and assumptions about the response surface can limit the applicability of strong methods on surfaces that are rough, discontinuous, or multimodal. However, in cases when the surface is smooth and the search technique is close to an optimum, these methods will converge very rapidly and efficiently to the optimum position on the surface. Weak methods, which include grid searching, make very few assumptions about the surface and tend to search in an enumerative or a random manner. While these methods are very robust, they are also computationally inefficient and become quickly limited as the dimensionality of a problem increases. Moderate techniques such as genetic algorithms and simulated annealing are another class of optimization techniques which bridge the division between the weak and strong methods. GAs use random events, which are directed by information about the previously unknown surface, accumulated during the search for the optimum, in order to identify and focus on regions of the search space likely to contain the required optimum.8,16 Moderate techniques tend to be particularly useful in situations where a problem is very complex (as the assumptions on which strong methods are based will make them fail) and also in cases where the dimensionality of the problem is high (as the computational time required for weak techniques to find an optimum becomes impractical). GAs use concepts gleaned from Darwin’s theory of evolution by natural selection and also from the mechanisms involved in the alteration and transfer of genetic information to individuals in a population.17 These processes enable a population to both survive in the surrounding environment and continuously adapt to any changes that may occur in it. A GA that mimics these evolutionary processes is implemented as an iterative procedure which maintains a constant size population P of possible solutions to a particular optimization problem. The P possible solutions to the optimization problem contained within the population are known as chromosomes, and the M variables of the search space are encoded on the chromosomes as M genes. The numeric values of the variables which might represent concentration, potential, wavelength, or any other analytical parameter are encoded in a uniform alphabet on the genes. Quite commonly the alphabet used is binary, in which case the number of bits used to represent the variable in the gene determine the resolution with which the variables are optimized. The basic configuration of the GA used in this study was composed of an initial population setup followed by a cyclic (15) Nash, J. C. Compact Numerical Methods for Computers, 2nd Ed.; Adam Hilgen: Bristol, U.K., 1990. (16) Goldberg, D. E. Genetic Algorithms in Search, Optimization and Machine Learning; Addison-Wesley: Reading, MA, 1989. (17) Becker, J. D., Eisele, I., Mundemann F. W., Eds. Parallelism, Learning, Evolution, Workshop on evolutionary models and strategies, Germany 1989. Workshop on parallel processing: logic, organisation and technology, WOPPLOT 89, Germany 1989. Springer Verlag: New York, 1991.
Figure 1. Basic configuration of a genetic algorithm.
repetition of evaluation, scaled reproduction, single-point crossover and single-point mutation (as depicted in Figure 1). Stages in the Operation of a GA. (1) Population Initialization. In this study, each parameter to be optimized is encoded on a chromosome in a binary form as a 16-bit word. This was achieved by establishing a range for each parameter within which the SSE (see expression 5) of the calibration model would be expected to be minimized. The range for each parameter would thus define a continuous but finite area of the search space for each parameter within which the search would proceed. The maximum and minimum values of these ranges were then mapped into the numeric ranges 0-216 - 1 using the following transform.
(216 - 1)[U - Umin] rescale ) Umax - Umin
(6)
where [Umin, Umax] represents the interval within which the optimization is to proceed for a particular parameter U. The ranges for these variables were set with different widths depending on the significance of the contribution of each variable to the error of the calibration model. The individual genes corresponding to the parameters of the Nikolskii-Eisenman expression were then concatenated to form a single bit string or chromosome. Variablesized populations were filled using a random number generator to produce numbers in the range 0-216 - 1 for each gene. In view of the difficulties that arise with system-supplied random number generators,18a code for a portable random number generation based on Knuth’s subtractive method16,18b was developed and seeded from the PC clock. (18) Press, W. H.; Teukolsky, S. A.; Flannery, B. P.; Vetterling, W. T. Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press: New York, 1991. (a) System supplied random number generators p 206. (b) Knuth’s subtractive method p 212.
(2) Evaluation. Each candidate model described by a particular chromosome was evaluated in order to determine how closely it described the experimental data. This was achieved by decoding the binary genes on the chromosomes into values of the model parameters. The resulting model was then used to predict the values of the experimental data. A SSE expression (as described by eq 5) was thus used as an objective function to quantitatively express how well the candidate model described the experimental data. Since a genetic algorithm is normally concerned with maximizing the performance or fitness of a population of chromosomes, the reciprocal of the error (i.e., 1/SSE) was used to describe the fitness of each candidate chromosome. (3) Scaled Reproduction. The existing calibration models in a particular population are exploited to improve the future performance of the population by performing biased reproduction, ensuring that particularly good calibration models have greater probability of being reproduced in the next generation than less good models. This involves sorting the models according to their SSE by means of a quick-sort procedure. (4) Crossover. In this study, single-point crossover is implemented by randomly selecting two chromosomes from the current population. The position at which crossover was to occur on both bit strings was chosen at random. The genes on both chromosomes in which the cross point occurred were converted from integers in the range 0-216 - 1 into their binary equivalents. This was achieved by running a bit mask across the integer representations of the genes and storing the resulting binary representations in character arrays. The elements of the two character arrays (corresponding to the two genes within which crossover was to occur) to the right of the crosspoint were swapped as depicted in Figure 2. The two resulting binary character arrays were reconverted into genes of an integer form according to the following relation. L
x)
∑a 2
i-1
i
(7)
i)1
Where x represents the parameter (integer value) and A ) aL, aL-1, ..., a2, a1 represents the L element bit string. (5) Mutation. Mutation is an operator that inverts the value of a bit which is chosen at random on a chromosome (Figure 3). If it were to be used as the main way of generating new individuals in the population, a GA would become a randomized search which, as described earlier, is inefficient for problems of high dimensionality. Instead, mutation acts to prevent genes that may occur in bad combinations with other genes from being completely lost to the population (by selective reproduction), when they may become more useful later in the GA in a combination with better genes. In this study, chromosomes were mutated at randomly chosen points on the bit strings of chromosomes by inverting the binary value of the bit at the mutation point through a logical XOR operation with a bit of value 1. In summary, the simple genetic algorithm (SGA) involves encoding a problem of M variables onto a population of N chromosomes with M genes. The chromosomes are allowed to reproduce according to their fitness determined by the objective function and are randomly subjected to genetic operators such as crossover and mutation. This iterative process is continued until some termination criterion is met. It can be seen that the Analytical Chemistry, Vol. 69, No. 10, May 15, 1997
1911
Figure 2. Processes occurring during crossover at the bit level. Chromosome 1 (dotted) and Chromosome 2 (white) are chosen to cross at bit position 9 in gene 2. Bits to the right of position 9 in gene 2 are switched between the two chromosomes as are the genes to the right of gene 2 in each chromosome.
Figure 3. Processes occurring during mutation. Bit 3 is has a value of 0 before mutation, after mutation bit 3 is toggled to a value of 1, (all the other bits in the gene remain the same).
GA is a highly parallel search technique whose parallelism is achieved by the simultaneous movement of the P search points through the M dimensional multivariate search space.19 For further discussion of the theory and application of GAs references 8, 16, 19, and 20 are recommended. AIMS AND NOMENCLATURE The aims of this contribution are to investigate the use of the GA and simplex optimization methods for the determination of the coefficients in the Nikolskii-Eisenman expression (eq 1) that best describe the responses of the ISEs over the calibration range used. To clarify the discussion, a superscript will be used to depict (19) Holland J. H. Adaptation in Natural and Artificial Systems; Bradford Books Edition, USA, MIT Press: Cambridge, MA, 1992. (20) De Weijer, A. P.; Lucasius, C. B.; Buydens, L.; Kateman, G.; Heuvel, H. M.; Mannee, H. Anal. Chem. 1994, 66, 23.
1912
Analytical Chemistry, Vol. 69, No. 10, May 15, 1997
the optimization technique used for the determination of a coefficient. For example, E°j (simp) indicates a standard cell potential of the electrode j, determined by a simplex method, similarly Ej°(GA) indicates a standard cell potential of electrode j, determined by the GA method. De Weijer et al.20 described the potential for improvement of the GA approach by hybridizing it with other numerical methods. In an attempt to improve on the GA method being used in the study, it was hybridized with a simplex technique and a Nikolskii-Eisenman coefficient determined in this manner will be described by the superscript (GA-simp) (e.g., k*ij (Ga-simp) is the conditional selectivity coefficient of the electrode i against species j, determined by the GA-simplex hybrid). Another note concerning the nomenclature used in this study refers to the use of a standard deviation on the NikolskiiEisenman coefficient estimates. The GA and simplex methods were repeated a number of times to allow for different initializations of each method. These repetitions also provided the opportunity to study the spread of the determined NikolskiiEisenman coefficients. The standard deviation determined in this manner is denoted by the subscript (rep). When a model was finally determined by the GA-simplex hybrid, the standard errors of the parameters were determined by the method given in Appendix A. Standard errors determined in this manner are denoted by SE(fin). EXPERIMENTAL DETAILS The data used for this study were provided courtesy of Dr. F. J. Sa´ez de Viteri, School of Chemical Sciences, Dublin City University. The data were acquired from an array of electrodes
selective for ammonium, sodium, potassium, and calcium which were employed in a flow injection system. The electrodes were based on PVC membranes containing the following ionophores: (i) ammonium ionophore I for the ammonium-selective electrode, (ii) p-tert-butylcalix[4]arene methoxyethyl ester21 for the sodiumselective electrode, (iii) valinomycin for the potassium-selective electrode, and (iv) calcium ionophore II (ETH 129) for the calciumselective electrode. Details of the construction of the electrodes and the flow injection analysis system can be found in ref 4. The calibration solutions were prepared according to a twolevel, four-factor, partial factorial experimental design, developed in a fashion to evoke a significant response from the electrodes to their interferents in comparison to that of the corresponding primary ions. This approach resulted in the use of eight solutions specifically designed for the calibration of a particular electrode. Since this system involved four electrodes, a total of 32 solutions were required to calibrate the electrodes. The concentration range for the calibration solutions was 10-410-2 M for the primary ion of each of the 4 electrodes and 10-410-3 M for their respective interferent ions. The pH of the solutions was controlled via the addition of 10-2 M Tris buffer. The actual composition of the solutions is described in more detail in ref 4. The activity coefficients (γi) for ions in these solutions were calculated by means of the Davies eq22 described in expression 8. In this expression, I refers to the ionic strength of
log γi ) -0.5z2i
(
xI
1 + xI
)
- 0.2I
(8)
the solution and zi refers to the charge of the ion in question. The resulting activities of the calibration solutions and the potential (mV) responses of the ISEs to these solutions are also described in ref 4. All 32 calibration solutions were injected into the FIA system four times and the resulting traces smoothed with a second order Savitzky-Golay filter. The heights of the FIA peaks were measured; the purpose of this study is to relate these measurements with the activities of the ions in the corresponding solutions. Since only eight specifically designed solutions are required to build a model for a given ISE (according to the two-level, fourfactor, partial factorial experimental design), the remaining 24 solutions could be used to test the performance of the electrode model. This approach was used in the study, whereby the 8 solutions specifically designed for the calibration of the electrode were used to build the model and the remaining 24 solutions were used to test it. Future discussions of residuals will refer to residuals obtained from the model prediction of the ion activities in the test set for a given electrode. The experimental data were acquired via an Analog devices RTI-815 data acquisition card fitted inside an IBM 386 compatible PC. Data acquisition and processing software were written in Microsoft QuickBasic and Visual Basic. The GA and simplex (21) Mc Kervey, M. A.; Seward, E. M.; Ferguson, G.; Ruhl, B.; Harris, S. J. J. Chem. Soc., Chem. Commun. 1985, 388. (22) Horvath, A. L. Handbook of Aqueous Electrolyte Solutions: Physical Properties, Estimation and Correlation Methods; Ellis Horwood Ltd.: Chichester, W. Sussex, England, 1985; p 213. (23) Baker, J. E. In Adaptive Selection Methods for Genetic Algorithms; Grefenstette, J. J., Ed.; Proceedings of the First International Conference on Genetic Algorithms. Lawrence Erlbaum, NJ, 1985; p 101.
Table 1. General Configuration of Genetic Algorithms Used for This Study population scaling rank scaling upper limit selection crossover crossover restriction
cross rate mutation mutation rate
300 chromosomes rank scaling19 1.1 stochiastic remainder12 single point12 chromosome mates must be nonidentical and must not have already crossed in the same iteration of the GA 0.65 single point 1.0/(population size)
software were written in Borland Turbo C++ operating on a 386 PC. RESULTS Application to an Array of ISEs Used in the FIA Regime. The configuration of the GA used for this study is depicted in Table 1. The GA was posthybridized with a simplex method to refine the solution determined by the GA method. This was achieved by repeating the GA six times; the calibration parameters obtained from the chromosomes that most closely modeled the response of the ISE from each repetition of the GA were then employed as vertices of a simplex optimization procedure used to refine the eventual calibration model. The activities of the ions in the calibration solutions were predicted by linearization of the simplified Nikolskii-Eisenman expressions for each electrode (using the values of the cell potentials, slopes, and selectivity coefficients determined by the GA-simplex hybrid), followed by Gauss-Jordan elimination of the selectivity coefficient matrix. Each calibration parameter for each electrode was initially encoded over a broad range (e.g., for the potassium ISE 100 < + E°(GA)
