Optimization with sequential simplex of variable size
by E. G. Romero-Blanco <eromero@costarricense.cr> and J. F. Ogilvie <ogilvie@cecm.sfu.ca>
Universidad de Costa Rica, Escuela de Quimica, Ciudad Universitaria Rodrigo Facio, San Jose, Costa Rica.
@ 2002 July
abstract
We present an algorithm for unconstrained optimisation of experimental data in chemical kinetics using as method a sequential simplex of variable size. An animated plot of progress of a triangular simplex descending across a surface of to converge to an absolute minimum illustrates an implementation of this approach.
introduction
A common problem in treating chemical and physical systems is a necessity to evaluate an optimum response -- a maximum or minimum, which is explicitly or implicitly a function of several factors. Although many methods exist for this purpose, conceptually the most simple is the sequential simplex method or its modifications. This approach has as its objective to locate efficiently the region or point of optimum response in a space defined by the factors and to find simultaneously the value of this extremum and to evaluate the magnitudes of factors that produce it [1, 2].
A simplex is a geometric figure defined on a response surface, or hypersurface in multidimensional space, of which the number of vertices or corners exceeds by one the number of factors, so . Such a simplex as a geometric figure is described as belonging to a factor space of dimensions [1,2]. The basic sequential simplex or algorithm for a sequential simplex of fixed size, introduced by Spendley, Hext and Himsworth [3], provides rules to force such a simplex to move by reflexions of vertices of an initially selected simplex to a region of optimum response.
Although originally developed for experimental optimization, the simplex method and its modifications have been applied to mathematical optimisation. For a purpose of applying the simplex method to a least-square fit of data to nonlinear models, Nelder and Mead [4] reported a sequential simplex with variable sizeon introducing two modifications into the original algorithm of Spendley et alii . These two modifications lead the simplex to expand in directions that are favourable and to contract in directions that are unfavourable, so as to provide a means according to which this figure can accelerate its progress toward a region of an optimum and can then diminish its search region until the extremum is located within desired accuracy.
Advantages of using a sequential simplex in fitting nonlinear models over other methods include its compact size, its relative independence from initial values of factors or parameters, and no necessity to calculate derivatives with respect to parameters. Like other approaches to fit a model nonlinear in parameters, the sequential simplex requires initial estimates of those parameters. Whereas the sequential simplex method converges from almost any initial values, in sets, approaches based on derivatives require initial estimates typically near ultimate values of parameters; this property becomes important during fitting of models with many parameters.
The purpose of this work is to employ Maple 's capabilities to demonstrate the applicability of a sequential simplex method, based on a least sum of squares of residuals and with variable size of simplex, to fit, to previously published experimental data [2] , a model well known in chemical kinetics of first order for the variation of concentration of a reaction product as a function of time,
in which is absorbance of a product of a reaction of first kinetic order and is taken to be directly proportional to the concentration of that product, t is time, is absorbance at infinite time, and k is the rate coefficient; such a relation implies that progress of a chemical reaction is being followed by means of the extent of absorption of a coloured product as a function of time. The particular objective of solution of our optimisation task is to find numerical values of factors and k that minimize the response function , or chi squared, which is defined as a sum of squared residuals, with each residual as a difference between an observed value of absorbance and the corresponding expectation for experimental data of number n .
With respect to the chemical system, the response to variation of independent variable time, denoted , is the absorbance of the product of the chemical reaction in the, for instance, visible spectrum, whereas with respect to the motion of thesimplex the response is as a result of variation of parameters and in the course of fitting the kinetic model to measured values of and .
basis of the calculation
To begin, we specify, as an arrow function, a model or objective function to fit, in terms of an independent variable time and two unknown factors -- the rate coefficient k and the ultimate absorbance of the coloured product at completion of reaction, represented as A_inf .
Experimental data to be fitted correspond to nine duplicated measurements of absorbance of a chemical system reacting according to first kinetic order during a defined period; we name lists of values of time and for values of absorbance .
A plot of these data shows the behaviour of this system of first kinetic order.
Through a Maple command we directly evaluate parameters and k that fit the model to these data on minimizing a sum of squared residuals, , over all experimental data.
For future reference we take as exact values these results for parameters and a corresponding magnitude of , denoted respectively as A_min , k_min and chisq_min , which we extract from the output above.
We view the topography of the response surface of for the particular data by means of a plot in three dimensions, one for each parameter and k and the third for , over values of parameters in a convenient range. The surface of presents a curved valley with parabolic cross sections near the minimum; that minimum response is readily located at the deepest point of the valley. On this plot in three dimensions we superimpose contours of for selected values,
and plot separately these contours.
In the upper three-dimensional plot, thin black lines indicate contours of at constant values increasing ina geometric progression, 0.008, 0.008 , 0.008 and so forth, the same values as for orange lines in the lower contour plot; a small circle within the innermost contour indicates the minimum in both cases.
algorithm of a sequential simplex
To implement the method of sequential simplex with variable size we begin by constructing an initial simplex. The efficiency of any sequential simplex method depends to some extent on the size, orientation and location of the initial simplex. In this case, because two factors, , namely and k , define explicitly the response that we apply as a sum of squared residuals, , to fix an initial simplex we must define coordinates of three vertices, , that form a triangle, in a factor space having two dimensions. For vertices of an initial simplex we select coordinates of points on a surface of that enable us to observethe behaviour of successive simplices in a sequence with varied sizes, orientations and locations during progress toward a minimum; each initial vertex is defined as a value of a subscripted variable, so that the number of coordinates equals the number of factors. For ourtest case we take initial vertices to have these coordinates [ , k ]:
[0.600, 0.800], [0.638, 0.826] and [0.613, 0.897],
which define a small simplex on the side of the valley representing a surface of remote from its minimum; generating further simplices from this initial condition, we can view motion of a simplex across the surface.
We show this initial simplex on the contour plot.
For this initial simplex each vertex of the three has acorresponding value of response .
We must rank, according to its magnitude, the response for each vertex as best, B, intermediate, N, and worst, W, such that the best vertex has the minimum response and the worst vertex has the maximum response for these three cases. For a case in which factors number more than two, the vertex next to worst would be selected instead of that here called intermediate.
The initial simplex is so characterized with its best values of two parameters and and of response
To make the first movement of the triangular initial simplex we reject its worst vertex W and replace it with a new vertex R, generated on reflexion of vertex W through midpoint P of the opposite edge between best B and intermediate N vertices:
After that initial reflexion, further movements of the sequential simplex proceed according to modifications by Nelder and Mead of the original algorithm; this modified algorithm contains rules necessary to effect the operations reflexion, expansion and contraction of an initial simplex until it reaches an optimum response, or corresponding minimum value of . The conditions " > " , " > " , " < " and " < " imply avertex that has a response " better than ", " better than or equal to ", " worse than " and " worse than or equal to " [1, 2] according to the following rules.
We apply exactly these rules in the following programmed loop of commands until convergence according to a defined criterion. Before the loop begins we supply information about the initial simplex and the first reflexion .
Following this preparation we execute a loop to construct a sequential simplex with variable size.
generation of sequentialsimplices
After a criterion for convergence is satisfied, we obtain the following results for the final simplex.
These final simplex results of A_inf_final and k_final correspond to values of parameters that fit the first-order kinetic model to experimental data with equal to chisq_final . We plot the fitted curve with the experimental points.
We compare these final simplex results with exact values obtained above.
We calculate a fractional error as a convenient comparison of results in the two sets; this bias depends on a criterion for convergence fixed in the loop above.
These errors can be decreased on diminishing a criterion of convergence in the above loop, but such a diminution has no discernible effect on movement of the simplex on a scale of contours shown above, and likely involve fitting data within the range of their uncertainties.
movement of simplex towards convergence
The following plot shows the movement of the simplex from its selected initial position on the surface of , according to its contours, to its ultimate destination at the location of a minimum as evaluated above. Play this animation!
The next plot tracks the sequence of simplices from an initial location until convergence, and shows how the size of each simplex expands or contracts from the preceding figure according to the slope of the surface: the simplex rapidly expands down the slope into the valley, contracts when reaching the valley, accelerates along the valley toward the minimum and finally collapses at the minimum after contractions in a continuous sequence.
Here is a corresponding plot of values of in a sequence towards convergence; the behaviour shown for this plot is consistent with the sequential movement of the simplex over the response surface. The variations of response produced are initially large as the simplex moves down the slope into the valley, but thereafter variations between consecutive figures are small as the simplex moves along the valley as far as the optimum.
The following plots show how values of parameters and vary during the progress of the simplex from its initial location until convergence. The behaviour that is observed for both parameters reflects movements of the simplex across the surface: rapid expansions and contractions at the beginning of the sequential simplex produce important -- even erratic -- alterations in the magnitudeof each parameter. When the simplex approaches the optimum by mean of sucessive contractions, there result smaller alterations of parameters until criterion for convergence criterion is satisfied.
table of numerical simplex progress
Here is a table of numerical values of parameters and response in progress of the simplex during optimization.
other initial vertices
Varying size and location of an initial simplex results in varied behaviour across the response surface until the minimum. To explore this behaviour we suggest for an alternative initial vertex these coordinates [ , k ] :
conclusion
Although simplex algorithms generally applied in mathematical optimisation workwith simplices of fixed or variable size, important modifications are reported in the literature to involve boundary conditions for constrained optimisation, weighted-centroid algorithms, and fitting of multiparametric models with estimates of errors of parameters [1]. Important published modifications include a super-modified simplex [5] or a composite-modified simplex [6] that seems to effect increased efficiency in the simplex. Here we have shown the applicability of Maple 's environment to study optimisation through fitting of a model with two parameters according to a criterion of least sum of squared residuals using an algorithm for a sequential simplex of variable size; introduction and analysis of any other simplex extension is straightforward, reflecting Maple 's capabilities for graphics and mathematical tools, easy handling of variables and powerful programming.
references
[1] Walters, F. H.; Parker, L. R.; Morgan, S. L.; Deming, S. N. Sequential Simplex Optimisation , CRC Press LLC, Florida, USA, 1991 .
[2] Deming, S. N.; Morgan, S. L. "Simplex optimisation of variables in analytical chemistry", Analytical Chemistry , 1973 , 45, 278A--283A.
[3] Spendly, W.; Hext, G. R.; Himsworth, F. R. "Sequential application of simplex designs in optimisation and evolutionary operation", Technometrics , 1962 , 4, 441--461.
[4] Nelder, A; Mead, R. "A simplex method for function minimisation", Computer Journal , 1965 , 7, 308--313.
[5] Routh, M. W.; Swartz, R. A.; Denton, M. B. "Performance of the super-modified simplex", Analytical Chemistry , 1977 , 49, 1422--1428.
[6] Betteridge, D.; Wade, R. A.; Howard, A. G. "Reflections on the modified simplex - I, II", Talanta , 1985 , 32, 709--734.
