1 <?xml version="1.0" encoding="ISO-8859-1"?>
3 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4 * Copyright (C) 2008 - 2009 - INRIA - Michael Baudin
5 * Copyright (C) 2009 - 2011 - DIGITEO - Michael Baudin
7 * This file must be used under the terms of the CeCILL.
8 * This source file is licensed as described in the file COPYING, which
9 * you should have received as part of this distribution. The terms
10 * are also available at
11 * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
14 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="neldermead" xml:lang="fr">
16 <refname>neldermead</refname>
18 Provides direct search optimization algorithms.
22 <title>SYNOPSIS</title>
24 newobj = neldermead_new ()
25 this = neldermead_destroy (this)
26 this = neldermead_configure (this,key,value)
27 value = neldermead_cget (this,key)
28 value = neldermead_get ( this , key )
29 this = neldermead_search ( this )
30 this = neldermead_restart ( this )
31 [ this , result ] = neldermead_function ( this , x )
32 stop = neldermead_defaultoutput(state, data)
36 <title>Description</title>
38 This class provides several direct search optimization algorithms
39 based on the simplex method.
42 The optimization problem to solve is the minimization of a cost
43 function, with bounds and nonlinear constraints
47 l_i <= x_i <= h_i, i = 1,n
48 g_i(x) >= 0, i = 1,nbineq
55 <para>number of variables</para>
61 <para>number of inequality constraints</para>
66 The provided algorithms are direct search algorithms, i.e.
67 algorithms which do not use the derivative of the cost function. They are
68 based on the update of a simplex, which is a set of k>=n+1 vertices,
69 where each vertex is associated with one point and one function
72 <para>The following algorithms are available :</para>
75 <term>Spendley, Hext and Himsworth fixed shape simplex method</term>
78 This algorithm solves an unconstrained optimization problem
79 with a fixed shape simplex made of k=n+1 vertices.
84 <term>Nelder and Mead variable shape simplex method</term>
87 This algorithm solves an unconstrained optimization problem
88 with a variable shape simplex made of k=n+1 vertices.
93 <term>Box complex method</term>
96 This algorithm solves an constrained optimization problem with
97 a variable shape simplex made of an arbitrary k number of vertices
98 (k=2n is recommended by Box).
104 See the demonstrations, in the Optimization section, for an overview
108 See the "Nelder-Mead User's Manual" on Scilab's wiki and on the
109 Scilab forge for further information.
113 <title>Design</title>
115 The neldermead component is built on top of the <link linkend="optimbase_overview">optimbase</link> and <link linkend="optimsimplex_overview">optimsimplex</link> components.
119 <title>Functions</title>
120 <para>The following functions are available.</para>
123 <term xml:id="neldermead_new">newobj = neldermead_new ()</term>
125 <para>Creates a new neldermead object.</para>
130 <para>The new object.</para>
137 <term xml:id="neldermead_destroy">this = neldermead_destroy (this)</term>
139 <para>Destroy the given object.</para>
144 <para>The current object.</para>
151 <term xml:id="neldermead_configure">this = neldermead_configure (this,key,value)</term>
154 Configure the current object with the given value for the
161 <para>The current object.</para>
168 the key to configure. The following keys are
172 <emphasis>Basic.</emphasis>
176 <term>-numberofvariables</term>
179 a 1-by-1 matrix of doubles, positive, integer value,
180 the number of variables to optimize (default numberofvariables = 0).
183 This does not need to be set if <literal>-x0</literal> key has been set before,
184 because the <literal>-x0</literal> key sets <literal>-numberofvariables</literal>
185 to the size of the column vector <literal>x0</literal>.
190 <term>-function</term>
193 a function or a list, the objective function.
194 This function computes the value
195 of the cost and the non linear constraints, if
199 There is no default value, i.e. the user must
200 provide <literal>f</literal>.
203 See below for the details of the communication
204 between the optimization system and the cost
213 a n-by-1 matrix of doubles, where
214 n is the number of variables, the initial guess.
217 There is no default value, i.e. the user must
218 provide <literal>x0</literal>.
224 <emphasis>Output.</emphasis>
228 <term>-outputcommand</term>
231 a function or a list, a function which is called back for output.
234 The default output function is empty, meaning that there is
238 See below for the details of the communication
239 between the optimization system and the output command
245 <term>-storehistory</term>
248 a 1-by-1 matrix of booleans, set to %t to enable the history storing (default
254 <term>-verbose</term>
257 a 1-by-1 matrix of doubles, positive, integer value, set to 1 to
258 enable verbose logging (default verbose = 0).
263 <term>-verbosetermination</term>
266 a 1-by-1 matrix of doubles, positive, integer value,
267 set to 1 to enable verbose termination logging (default verbosetermination = 0).
273 <emphasis>Bounds and constraints.</emphasis>
277 <term>-boundsmin</term>
280 a n-by-1 matrix of doubles, the minimum bounds for the parameters
281 where n is the number of variables (default boundsmin = [], i.e. there are no
287 <term>-boundsmax</term>
290 a n-by-1 matrix of doubles, the maximum bounds for the parameters
291 where n is the number of variables (default boundsmax = [], i.e. there are no
297 <term>-nbineqconst</term>
300 a 1-by-1 matrix of doubles, positive, integer value,
301 the number of inequality constraints (default nbineqconst = 0).
307 <emphasis>Initialization.</emphasis>
311 <term>-simplex0method</term>
314 a 1-by-1 matrix of strings, the method to use to compute the initial
315 simplex (default simplex0method = "axes").
316 The first vertex in the simplex is always the initial
317 guess associated with the -x0 option. The following
318 methods are available :
325 the coordinates associated with the -coords0
326 option are used to compute the initial simplex, with
327 arbitrary number of vertices.
330 This allow the user to setup the initial
331 simplex by a specific method which is not provided
332 by the current component (for example with a simplex
333 computed from a design of experiments). This allows
334 also to configure the initial simplex so that a
335 specific behaviour of the algorithm an be reproduced
336 (for example the Mac Kinnon test case).
339 The given matrix is expected to have n rows
340 and k columns, where n is the dimension of the
341 problem and k is the number of vertices.
349 the simplex is computed from the coordinate
350 axes and the length associated with the
351 -simplex0length option.
356 <term>"spendley"</term>
359 the simplex is computed so that it is regular
360 with the length associated with the -simplex0length
361 option (i.e. all the edges have the same
367 <term>"pfeffer"</term>
370 the simplex is computed from a heuristic, in
371 the neighborhood of the initial guess. This initial
372 simplex depends on the -simplex0deltausual and
378 <term>"randbounds"</term>
381 the simplex is computed from the bounds and a
382 random number. This option is available only if
383 bounds are available : if bounds are not available,
384 an error is generated. This method is usually
385 associated with Box's algorithm. The number of
386 vertices in the simplex is taken from the
395 <term>-coords0</term>
398 a nbve-by-n matrix of doubles, where nbve is the number of vertices and n is
399 the number of variables, the coordinates of the vertices of the initial
400 simplex (default coords0=[]).
401 If the -simplex0method option is set to
402 "given", these coordinates are used to compute the
408 <term>-simplex0length</term>
411 a 1-by-1 matrix of doubles, the length to use when the initial simplex is
412 computed with the "axes" or "spendley" methods (default simplex0length = 1).
413 If the initial simplex is computed from "spendley" method, the
414 length is expected to be a scalar value. If the initial
415 simplex is computed from "axes" method, it may be either
416 a scalar value or a vector of values, with rank n, where
417 n is the number of variables.
422 <term>-simplex0deltausual</term>
425 a 1-by-1 matrix of doubles, the relative delta for non-zero parameters in
426 "pfeffer" method (default simplex0deltausual = 0.05).
431 <term>-simplex0deltazero</term>
434 a 1-by-1 matrix of doubles, the absolute delta for non-zero parameters in
435 "pfeffer" method (default simplex0deltazero = 0.0075).
441 <emphasis>Termination.</emphasis>
445 <term>-maxfunevals</term>
448 a 1-by-1 matrix of doubles, positive, integer value,
449 the maximum number of function evaluations
450 (default maxfunevals = 100).
451 If this criteria is triggered, the status of the optimization is set to
457 <term>-maxiter</term>
460 a 1-by-1 matrix of doubles, positive, integer value,
461 the maximum number of iterations (default maxiter = 100).
462 If this criteria is triggered, the status of the
463 optimization is set to "maxiter".
468 <term>-tolfunabsolute</term>
471 a 1-by-1 matrix of doubles, positive, the absolute tolerance for the function value
472 (default tolfunabsolute = 0).
477 <term>-tolfunrelative</term>
480 a 1-by-1 matrix of doubles, positive, the relative tolerance for the function value
481 (default tolfunrelative = %eps).
486 <term>-tolfunmethod</term>
489 a 1-by-1 matrix of booleans, set to %t to enable termination with
490 tolerance on function value (default tolfunmethod = %f).
493 If this criteria is triggered, the
494 status of the optimization is set to "tolf".
499 <term>-tolxabsolute</term>
502 a 1-by-1 matrix of doubles, positive, the absolute tolerance
503 on x (default tolxabsolute = 0).
508 <term>-tolxrelative</term>
511 a 1-by-1 matrix of doubles, positive, the relative
512 tolerance on x (default tolxrelative = sqrt(%eps)).
517 <term>-tolxmethod</term>
520 a 1-by-1 matrix of booleans, set to %t to enable the tolerance on x in the
521 termination criteria (default tolxmethod = %t).
524 If this criteria is triggered, the
525 status of the optimization is set to "tolx".
530 <term>-tolsimplexizemethod</term>
533 a 1-by-1 matrix of booleans, set to %f to disable the tolerance on the simplex
534 size (default tolsimplexizemethod = %t). If this criteria is
535 triggered, the status of the optimization is set to
539 When this criteria is enabled, the values of the
540 options -tolsimplexizeabsolute and
541 -tolsimplexizerelative are used in the termination
542 criteria. The method to compute the size is the
548 <term>-tolsimplexizeabsolute</term>
551 a 1-by-1 matrix of doubles, positive, the absolute tolerance on the
552 simplex size (default tolsimplexizeabsolute = 0).
557 <term>-tolsimplexizerelative</term>
560 a 1-by-1 matrix of doubles, positive, the relative tolerance on the
561 simplex size (default tolsimplexizerelative = %eps).
566 <term>-tolssizedeltafvmethod</term>
569 a 1-by-1 matrix of booleans, set to %t to enable the termination criteria based
570 on the size of the simplex and the difference of
571 function value in the simplex (default tolssizedeltafvmethod = %f).
572 If this criteria is triggered, the status of the
573 optimization is set to "tolsizedeltafv".
576 This termination criteria uses the values of the
577 options -tolsimplexizeabsolute and -toldeltafv. This
578 criteria is identical to Matlab's fminsearch.
583 <term>-toldeltafv</term>
586 a 1-by-1 matrix of doubles, positive, the absolute tolerance on the difference between
587 the highest and the lowest function values (default toldeltafv = %eps).
593 <emphasis>Algorithm.</emphasis>
600 a 1-by-1 matrix of strings, the name of the algorithm to use (default method = "variable").
601 The following methods are available :
608 the Spendley et al. fixed simplex shape
609 algorithm. This algorithm is for unconstrained
610 problems (i.e. bounds and non linear constraints are
611 not taken into account)
616 <term>"variable"</term>
619 the Nelder-Mead variable simplex shape
620 algorithm. This algorithm is for unconstrained
621 problems (i.e. bounds and non linear constraints are
622 not taken into account)
630 the Box complex algorithm. This algorithm
631 takes into account bounds and nonlinear inequality
640 the user-defined algorithm, associated with the
641 <literal>-mymethod</literal>option. See below for details.
649 <term>-mymethod</term>
652 a function, a user-derined simplex algorithm. See below for
653 details (default is empty).
659 <emphasis>Options of the "variable" algorithm.</emphasis>
666 a 1-by-1 matrix of doubles, the reflection coefficient. This parameter is used
667 when the -method option is set to "fixed" or "variable" (default rho = 1).
675 a 1-by-1 matrix of doubles, the expansion coefficient. This parameter is used
676 when the -method option is set to "variable" (default chi = 2).
684 a 1-by-1 matrix of doubles, the contraction coefficient. This parameter is
685 used when the -method option is set to "variable" (default gamma = 0.5).
693 a 1-by-1 matrix of doubles, the shrinkage coefficient. This parameter is used
694 when the -method option is set to "fixed" or "variable" (default sigma = 0.5).
700 <emphasis>Option of "box" algorithm.</emphasis>
704 <term>-scalingsimplex0</term>
707 a 1-by-1 matrix of strings, the algorithm used to scale the initial simplex into
708 the nonlinear constraints (default scalingsimplex0 = "tox0").
709 The following two algorithms are provided :
714 "tox0": scales the vertices toward the initial guess.
719 "tocenter": scales the vertices toward the centroid, as recommended by Box.
724 If the centroid happens to be unfeasible, because the constraints are
725 not convex, the scaling of the initial simplex toward the centroid
726 may fail. Since the initial guess is always feasible, scaling toward
727 the initial guess cannot fail.
728 The default value is "tox0".
733 <term>-boxnbpoints</term>
736 a 1-by-1 matrix of doubles, positive, integer value, the number
737 of points in the initial simplex, when
738 the -simplex0method is set to <literal>"randbounds"</literal>
739 (default boxnbpoints = 2*n, where n is the number of variables of the problem).
740 The value of this option is also use to update the simplex when a
741 restart is performed and the <literal>-restartsimplexmethod</literal>
742 option is set to <literal>"randbounds"</literal>.
747 <term>-boxineqscaling</term>
750 a 1-by-1 matrix of doubles, in the interval [0,1], the scaling coefficient used to scale the trial
751 point for function improvement or into the constraints of Box's
752 algorithm (default boxineqscaling = 0.5).
757 <term>-guinalphamin</term>
760 a 1-by-1 matrix of doubles, positive, the minimum value of alpha when scaling the vertices of the
761 simplex into nonlinear constraints in Box's algorithm (default guinalphamin = 1.e-5).
766 <term>-boxreflect</term>
769 a 1-by-1 matrix of doubles, positive, the reflection factor in
770 Box's algorithm (default = 1.3).
775 <term>-boxtermination</term>
778 a 1-by-1 matrix of booleans, set to %t to enable Box termination criteria (default boxtermination = %f).
783 <term>-boxtolf</term>
786 a 1-by-1 matrix of doubles, positive, the absolute tolerance on difference of function values in the
787 simplex, suggested by Box (default boxtolf = 1.e-5). This tolerance is used if
788 the -boxtermination option is set to %t.
794 <term>-boxnbmatch</term>
797 a 1-by-1 matrix of doubles, positive, integer value,
798 the number of consecutive match of Box termination criteria (default boxnbmatch = 5).
803 <term>-boxboundsalpha</term>
806 a 1-by-1 matrix of doubles, positive, the parameter used to project the vertices into the
807 bounds in Box's algorithm (default boxboundsalpha = 1.e-6).
813 <emphasis>Auto-Restart.</emphasis>
817 <term>-kelleystagnationflag</term>
820 a 1-by-1 matrix of booleans, set to %t to enable the termination criteria using
821 Kelley's stagnation detection, based on sufficient
822 decrease condition (default kelleystagnationflag = %f). If this
823 criteria is triggered, the status of the optimization is
824 set to "kelleystagnation".
829 <term>-kelleynormalizationflag</term>
832 a 1-by-1 matrix of booleans, set to %f to disable the normalization of the
833 alpha coefficient in Kelley's stagnation detection, i.e.
834 use the value of the option -kelleystagnationalpha0 as
835 is (default kelleynormalizationflag = %t, i.e. the simplex gradient of
836 the initial simplex is taken into account in the
837 stagnation detection).
842 <term>-kelleystagnationalpha0</term>
845 a 1-by-1 matrix of doubles, the parameter used in Kelley's stagnation
846 detection (default kelleystagnationalpha0 = 1.e-4).
851 <term>-restartflag</term>
854 a 1-by-1 matrix of booleans, set to %t to enable the automatic restart of the
855 algorithm (default restartflag = %f).
860 <term>-restartdetection</term>
863 a 1-by-1 matrix of strings, the method to detect if the automatic restart must
864 be performed (default restartdetection = "oneil").
865 The following methods are available:
869 <term>"oneill"</term>
872 the factorial local optimality test by O'Neill
873 is used. If the test finds a local point which is
874 better than the computed optimum, a restart is
880 <term>"kelley"</term>
883 the sufficient decrease condition by O'Neill
884 is used. If the test finds that the status of the
885 optimization is "kelleystagnation", a restart is
886 performed. This status may be generated if the
887 -kelleystagnationflag option is set to %t.
892 <para>The default method is "oneill".</para>
896 <term>-restartmax</term>
899 a 1-by-1 matrix of doubles, the maximum number of restarts, when automatic
900 restart is enabled via the -restartflag option (default
906 <term>-restarteps</term>
909 a 1-by-1 matrix of doubles, the relative epsilon value used to check for
910 optimality in the factorial O'Neill restart detection (default restarteps = %eps).
915 <term>-restartstep</term>
918 a 1-by-1 or a n-by-1 matrix of doubles, positive, where n is the number of
919 variables in the problem, the absolute step length used to check for
920 optimality in the factorial O'Neill restart detection (default restartstep = 1).
925 <term>-restartsimplexmethod</term>
928 a 1-by-1 matrix of strings, the method to compute the initial simplex after a
929 restart (default restartsimplexmethod = "oriented"). The following methods are available.
936 the coordinates associated with the -coords0
937 option are used to compute the initial simplex, with
938 arbitrary number of vertices.
941 This allow the user to setup the initial
942 simplex by a specific method which is not provided
943 by the current component (for example with a simplex
944 computed from a design of experiments). This allows
945 also to configure the initial simplex so that a
946 specific behaviour of the algorithm an be reproduced
947 (for example the Mc Kinnon test case).
950 The given matrix is expected to have n rows
951 and k columns, where n is the dimension of the
952 problem and k is the number of vertices.
960 the simplex is computed from the coordinate
961 axes and the length associated with the
962 -simplex0length option.
967 <term>"spendley"</term>
970 the simplex is computed so that it is regular
971 with the length associated with the -simplex0length
972 option (i.e. all the edges have the same
978 <term>"pfeffer"</term>
981 the simplex is computed from a heuristic, in
982 the neighborhood of the initial guess. This initial
983 simplex depends on the -simplex0deltausual and
989 <term>"randbounds"</term>
992 the simplex is computed from the bounds and a
993 random number. This option is available only if
994 bounds are available : if bounds are not available,
995 an error is generated. This method is usually
996 associated with Box's algorithm. The number of
997 vertices in the simplex is taken from the
1003 <term>"oriented"</term>
1006 the simplex is computed so that it is
1007 oriented, as suggested by C.T. Kelley.
1012 <para>The default method is "oriented".</para>
1021 <para>the value.</para>
1028 <term xml:id="neldermead_cget">value = neldermead_cget (this,key)</term>
1031 Get the value for the given key. If the key is unknown,
1038 <para>The current object.</para>
1045 the name of the key to quiery. The list of available
1046 keys is the same as for the neldermead_configure
1055 <term xml:id="neldermead_get">value = neldermead_get ( this , key )</term>
1058 Get the value for the given key. If the key is unknown,
1062 Most fields are available only after an optimization has
1063 been performed with one call to the neldermead_search
1070 <para>The current object.</para>
1076 <para>the key to get.</para>
1077 <para>The following keys are available :</para>
1080 <term>-funevals</term>
1082 <para>the number of function evaluations</para>
1086 <term>-iterations</term>
1088 <para>the number of iterations</para>
1095 the x optimum, as a n x 1 column vector, where n
1096 is the number of variables.
1103 <para>the optimum cost function value</para>
1107 <term>-historyxopt</term>
1110 an array, with nbiter values, containing the
1111 history of x during the iterations.
1114 This array is available after optimization if the
1115 history storing was enabled with the -storehistory
1121 <term>-historyfopt</term>
1124 an array, with nbiter values, containing the
1125 history of the function value during the
1129 This array is available after optimization if the
1130 history storing was enabled with the -storehistory
1138 <para>the function value for the initial guess</para>
1142 <term>-status</term>
1145 a string containing the status of the
1146 optimization. See below for details about the
1147 optimization status.
1152 <term>-historysimplex</term>
1155 a matrix containing the history of the simplex
1156 during the iterations. This matrix has rank nbiter x
1157 nbve x n, where nbiter is the number of iterations, nbve
1158 is the number of vertices in the simplex and n is the
1159 number of variables.
1164 <term>-simplex0</term>
1167 the initial simplex. This is a simplex object,
1168 which is suitable for processing with the optimsimplex
1174 <term>-simplexopt</term>
1177 the optimum simplex. This is a simplex object,
1178 which is suitable for processing with the optimsimplex
1184 <term>-restartnb</term>
1186 <para>the number of actual restarts performed.</para>
1196 <term xml:id="neldermead_search">this = neldermead_search ( this )</term>
1199 Performs the optimization associated with the method
1200 associated with the -method option and find the optimum.
1206 <para>The current object.</para>
1211 If the -restartflag option is enabled, automatic restarts are
1212 performed, based on the -restartdetection option.
1217 <term xml:id="neldermead_restart">this = neldermead_restart ( this )</term>
1220 Restarts the optimization by updating the simplex and
1221 performing a new search.
1227 <para>The current object.</para>
1234 <term xml:id="neldermead_function">[ this , result ] = neldermead_function ( this , x )</term>
1236 <para>Call the cost function and return the value.</para>
1241 <para>The current object.</para>
1247 <para>the point where the function is to be evaluated</para>
1254 optional, a flag to pass to the cost function (default
1255 = 1). See the section "The cost function" for available values
1264 <term xml:id="neldermead_defaultoutput">stop = neldermead_defaultoutput(state, data)</term>
1267 Prints messages at an iteration.
1270 This function provides a default implementation for the output function.
1271 There is one line by iteration, presenting the number of iterations, the
1272 number of function evaluations, the current function value and the
1273 current algorithm step.
1276 See "The output function" section below for a description of the input and
1278 See in the Examples section below for examples of this function.
1285 <title>The cost function</title>
1287 The option <literal>-function</literal> allows to configure the cost function. The cost
1288 function is used to compute the objective function value <literal>f</literal>.
1289 If the <literal>-nbineqconst</literal> option is configured to a non-zero value, the cost function
1290 must also compute the value of the nonlinear, positive, inequality constraints <literal>c</literal>.
1291 Depending of these options, the cost function can have one of the following headers :
1294 [ f , index ] = costf ( x , index )
1295 [ f , c , index ] = costf ( x , index )
1302 <para>the current point, as a column vector</para>
1308 <para>optional, an integer representing the value to compute</para>
1314 <para>the value of the cost function</para>
1320 <para>the value of the non-linear, positive, inequality constraints</para>
1325 The index input parameter tells to the cost function what is expected
1326 in the output arguments. It has the following meaning
1330 <term>index = 2</term>
1333 compute <literal>f</literal>
1338 <term>index = 5</term>
1341 compute <literal>c</literal>
1346 <term>index = 6</term>
1349 compute <literal>f</literal> and <literal>c</literal>
1355 In the most simplex case, there is no additionnal cost function argument and no nonlinear constraints.
1356 In this case, the cost function is expected to have the following header
1359 [ f , index ]= myfunction ( x , index )
1362 It might happen that the function requires additionnal arguments to be evaluated.
1363 In this case, we can use the following feature.
1364 The argument <literal>fun</literal> can also be the list <literal>(myfun,a1,a2,...)</literal>.
1365 In this case <literal>myfun</literal>, the first element in the list, must be a function and must
1368 [ f , index ] = myfun ( x , index , a1, a2, ... )
1369 [ f , c , index ] = myfun ( x , index , a1, a2, ...)
1371 where the input arguments <literal>a1, a2, ...</literal>
1372 are automatically appended at the end of the calling sequence.
1376 <title>The output function</title>
1378 The option -outputcommand allows to configure a command which is
1379 called back at the start of the optimization, at each iteration and at the
1380 end of the optimization.
1382 <para>The output function must have the following header</para>
1384 stop = outputcmd(state, data)
1392 a string representing the current state of the algorithm.
1393 Available values are "init", "iter", "done".
1400 <para>a data structure containing at least the following entries</para>
1405 <para>the current optimum</para>
1411 <para>the current function value</para>
1415 <term>iteration</term>
1417 <para>the current iteration index</para>
1421 <term>funccount</term>
1423 <para>the number of function evaluations</para>
1427 <term>simplex</term>
1429 <para>the current simplex</para>
1436 the previous step in the algorithm. The following values
1437 are available : "init", "done", "reflection", "expansion",
1438 "insidecontraction", "outsidecontraction", "reflectionnext",
1450 a 1-by-1 matrix of booleans, set to true to stop the
1451 algorithm, set to false to continue the optimization.
1457 It might happen that the output function requires additionnal arguments to be evaluated.
1458 In this case, we can use the following feature.
1459 The argument <literal>outputcmd</literal> can also be the list <literal>(outf,a1,a2,...)</literal>.
1460 In this case <literal>outf</literal>, the first element in the list, must be a function and must
1463 stop = outf ( state, data, a1, a2, ... )
1465 where the input arguments <literal>a1, a2, ...</literal>
1466 are automatically appended at the end of the calling sequence.
1469 If the output function sets the <literal>stop</literal> variable to true,
1470 then the optimization alorithm stops and the status of the optimization is
1471 set to <literal>"userstop"</literal>.
1475 <title>Termination</title>
1477 The current component takes into account for several generic
1478 termination criterias.
1480 <para>The following termination criterias are enabled by default :</para>
1483 <para>-maxiter,</para>
1486 <para>-maxfunevals,</para>
1489 <para>-tolxmethod.</para>
1492 <para>-tolsimplexizemethod.</para>
1496 The optimization_terminate function uses a set of rules to compute
1497 if the termination occurs, which leads to an optimization status which is
1498 equal to one of the following : "continue", "maxiter", "maxfunevals",
1499 "tolf", "tolx", "tolsize", "tolsizedeltafv",
1500 "kelleystagnation", "tolboxf", "tolvariance". The value of the status may also
1501 be a user-defined string, in the case where a user-defined termination
1502 function has been set.
1504 <para>The following set of rules is examined in this order.</para>
1508 By default, the status is <literal>"continue"</literal> and the <literal>terminate</literal> flag is
1514 The number of iterations is examined and compared to the
1515 <literal>-maxiter</literal> option : if the following condition
1518 iterations >= maxiter
1521 is true, then the status is set to "maxiter" and terminate is
1527 The number of function evaluations and compared to the
1528 <literal>-maxfunevals</literal> option is examined : if the following condition
1531 funevals >= maxfunevals
1534 is true, then the status is set to <literal>"maxfuneval"</literal> and <literal>terminate</literal> is
1540 The tolerance on function value is examined depending on the
1541 value of the <literal>-tolfunmethod</literal>.
1547 <para>then the criteria is just ignored.</para>
1553 <para>if the following condition</para>
1555 abs(currentfopt) < tolfunrelative * abs(previousfopt) + tolfunabsolute
1558 is true, then the status is set to "tolf" and terminate is
1565 The relative termination criteria on the function value works
1566 well if the function value at optimum is near zero. In that case, the
1567 function value at initial guess fx0 may be used as
1568 <literal>previousfopt</literal>.
1571 This criteria is sensitive to the <literal>-tolfunrelative</literal>
1572 and <literal>-tolfunabsolute</literal> options.
1575 The absolute termination criteria on the function value works if
1576 the user has an accurate idea of the optimum function value.
1581 The tolerance on x is examined depending on the value of the
1588 <para>then the criteria is just ignored.</para>
1594 <para>if the following condition</para>
1596 norm(xopt - previousxopt) < tolxrelative * norm(xopt) + tolxabsolute
1599 is true, then the status is set to <literal>"tolx"</literal> and <literal>terminate</literal> is
1606 This criteria is sensitive to the <literal>-tolxrelative</literal>
1607 and <literal>-tolxabsolute</literal> options.
1610 The relative termination criteria on x works well if x at
1611 optimum is different from zero. In that case, the condition measures
1612 the distance between two iterates.
1615 The absolute termination criteria on x works if the user has an
1616 accurate idea of the scale of the optimum x. If the optimum x is near
1617 0, the relative tolerance will not work and the absolute tolerance is
1623 The tolerance on simplex size is examined depending on
1624 the value of the <literal>-tolsimplexizemethod</literal> option.
1630 <para>then the criteria is just ignored.</para>
1636 <para>if the following condition</para>
1638 ssize < tolsimplexizerelative * simplexsize0 + tolsimplexizeabsolute
1641 is true where <literal>simplexsize0</literal> is the size of the simplex at
1642 iteration 0, then the <literal>status</literal> is set to <literal>"tolsize"</literal>
1643 and <literal>terminate</literal> is set to %t.
1646 The size of the simplex is computed from the "sigmaplus" method of the
1647 <literal>optimsimplex</literal> component.
1648 This criteria is sensitive to the <literal>-tolsimplexizeabsolute</literal> and
1649 the <literal>-tolsimplexizerelative</literal> options.
1657 The absolute tolerance on simplex size and absolute difference
1658 of function value is examined depending on the value of the
1659 -tolssizedeltafvmethod option.
1665 <para>then the criteria is just ignored.</para>
1671 <para>if both the following conditions</para>
1673 ssize < tolsimplexizeabsolute
1676 shiftfv < toldeltafv
1679 is true where <literal>ssize</literal> is the current simplex size and
1680 <literal>shiftfv</literal> is the absolute value of the difference of function
1681 value between the highest and lowest vertices, then the status
1682 is set to <literal>"tolsizedeltafv"</literal> and <literal>terminate</literal> is set to %t.
1690 The stagnation condition based on Kelley sufficient decrease
1691 condition is examined depending on the value of the
1692 <literal>-kelleystagnationflag</literal> option.
1698 <para>then the criteria is just ignored.</para>
1704 <para>if the following condition</para>
1706 newfvmean <= oldfvmean - alpha * sg' * sg
1709 is true where <literal>newfvmean</literal> (resp. <literal>oldfvmean</literal>) is the function
1710 value average in the current iteration (resp. in the previous
1711 iteration), then the status is set to "kelleystagnation" and
1712 terminate is set to %t. Here, <literal>alpha</literal> is a non-dimensional
1713 coefficient and <literal>sg</literal> is the simplex gradient.
1721 The termination condition suggested by Box
1722 is examined depending on the value of the
1723 -boxtermination option.
1729 <para>then the criteria is just ignored.</para>
1735 <para>if both the following conditions</para>
1737 shiftfv < boxtolf
1740 boxkount == boxnbmatch
1743 is true where <literal>shiftfv </literal>is the difference of function value
1744 between the best and worst vertices, and <literal>boxkount</literal> is the number of
1745 consecutive iterations where this criteria is met,
1746 then the status is set to "tolboxf" and
1747 terminate is set to %t.
1748 Here, the <literal>boxtolf</literal> parameter is the value associated with
1749 the <literal>-boxtolf</literal> option
1750 and is a user-defined absolute tolerance on the function value.
1751 The <literal>boxnbmatch</literal> parameter is the value associated with
1752 the <literal>-boxnbmatch</literal> option
1753 and is the user-defined number of consecutive match.
1761 The termination condition based on the variance of the function values in the simplex
1762 is examined depending on the value of the
1763 <literal>-tolvarianceflag</literal> option.
1769 <para>then the criteria is just ignored.</para>
1775 <para>if the following condition</para>
1777 var < tolrelativevariance * variancesimplex0 + tolabsolutevariance
1780 is true where <literal>var </literal>is the variance of the function values
1782 then the status is set to "tolvariance" and
1783 terminate is set to %t.
1784 Here, the <literal>tolrelativevariance</literal> parameter is the value associated with
1785 the <literal>-tolrelativevariance</literal> option
1786 and is a user-defined relative tolerance on the variance of the function values.
1787 The <literal>tolabsolutevariance</literal> parameter is the value associated with
1788 the <literal>-tolabsolutevariance</literal> option
1789 and is the user-defined absolute tolerance of the variance of the function values.
1798 <title>Kelley's stagnation detection</title>
1800 The stagnation detection criteria suggested by Kelley is based on a
1801 sufficient decrease condition, which requires a parameter alpha > 0 to
1802 be defined. The -kelleynormalizationflag option allows to configure the
1803 method to use to compute this alpha parameter : two methods are available,
1804 where each method corresponds to a different paper by Kelley :
1808 <term>constant</term>
1811 In "Detection and Remediation of Stagnation in the
1812 Nelder--Mead Algorithm Using a Sufficient Decrease Condition",
1813 Kelley uses a constant alpha, with the suggested value 1.e-4, which
1814 is is typical choice for line search method.
1819 <term>normalized</term>
1822 in "Iterative Methods for Optimization", Kelley uses a
1823 normalized alpha, computed from the following formula
1826 alpha = alpha0 * sigma0 / nsg
1829 where sigma0 is the size of the initial simplex and nsg is the
1830 norm of the simplex gradient for the initial guess point.
1837 <title>O'Neill factorial optimality test</title>
1839 In "Algorithm AS47 - Function minimization using a simplex
1840 procedure", R. O'Neill presents a fortran 77 implementation of the simplex
1841 method. A factorial test is used to check if the computed optimum point is
1842 a local minimum. If the -restartdetection option is set to "oneill", that
1843 factorial test is used to see if a restart should be performed.
1844 O'Neill's factorial test requires <literal>2n</literal> function evaluations, where <literal>n</literal>
1845 is the number of variables.
1849 <title>User-defined algorithm</title>
1851 The <literal>-mymethod</literal> option allows to configure a user-defined
1852 simplex-based algorithm. The reason for this option is that many simplex-based
1853 variants of Nelder-Mead's algorithm have been developed over the years,
1854 with specific goals. While it is not possible to provide them all, it is very convenient
1855 to use the current structure without being forced to make many developments.
1858 The value of the <literal>-mymethod</literal> option is expected to be
1859 a Scilab function with the following header
1862 this = myalgorithm ( this )
1865 where <literal>this</literal> is the current object.
1868 In order to use the user-defined algorithm, the <literal>-method</literal> option must
1869 be set to "mine". In this case, the component performs the optimization
1870 exactly as if the user-defined algorithm was provided by the component.
1873 The user interested in that feature may use the internal scripts provided in the
1874 distribution as templates and tune his own algorithm from that point.
1875 There is of course no warranty that the
1876 user-defined algorithm improves on the standard algorithm, so that users
1877 use this feature at their own risks.
1881 <title>Example #1: basic use</title>
1883 In the following example, we solve a simple quadratic test case. We
1884 begin by defining the cost function, which takes 2 input arguments and
1885 returns the objective. The classical starting point [-1.2 1.0] is used.
1886 The neldermead_new creates a new neldermead object. Then we use the
1887 neldermead_configure method to configure the parameters of the problem. We
1888 use all default settings and perform the search for the optimum. The
1889 neldermead_display function is used to display the state of the
1890 optimization and the neldermead_get is used to retrieve the optimum
1893 <programlisting role="example"><![CDATA[
1894 function [ f , index ] = quadratic ( x , index )
1895 f = x(1)^2 + x(2)^2;
1898 nm = neldermead_new ();
1899 nm = neldermead_configure(nm,"-numberofvariables",2);
1900 nm = neldermead_configure(nm,"-function",quadratic);
1901 nm = neldermead_configure(nm,"-x0",x0);
1902 nm = neldermead_search(nm);
1904 xopt = neldermead_get(nm,"-xopt");
1905 nm = neldermead_destroy(nm);
1906 ]]></programlisting>
1909 <title>Example #2: customized use</title>
1911 In the following example, we solve the Rosenbrock test case. We
1912 begin by defining the Rosenbrock function, which takes 2 input arguments
1913 and returns the objective. The classical starting point [-1.2 1.0] is
1914 used. The neldermead_new creates a new neldermead object. Then we use the
1915 neldermead_configure method to configure the parameters of the problem.
1916 The initial simplex is computed from the axes and the single length 1.0
1917 (this is the default, but is explicitly written here as an example). The
1918 variable simplex algorithm by Nelder and Mead is used, which corresponds
1919 to the -method "variable" option. The neldermead_search function performs
1920 the search for the minimum. Once the minimum is found, the
1921 neldermead_contour allows to compute the data required by the contour
1922 function. This is possible since our problem involves only 2 parameters.
1923 This function uses the cost function previously configured to compute the
1924 required data. The contour plot is directly drawn from the data provided
1925 by neldermead_contour. Then we plot the initial guess on the contour plot
1926 as a blue dot. The neldermead_get function is used to get the optimum,
1927 which is associated with the -xopt option. The optimum is plot on the
1928 contour plot as a red dot.
1930 <programlisting role="example"><![CDATA[
1931 function [ f , index ] = rosenbrock ( x , index )
1932 f = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;
1935 nm = neldermead_new ();
1936 nm = neldermead_configure(nm,"-numberofvariables",2);
1937 nm = neldermead_configure(nm,"-function",rosenbrock);
1938 nm = neldermead_configure(nm,"-x0",x0);
1939 nm = neldermead_configure(nm,"-maxiter",200);
1940 nm = neldermead_configure(nm,"-maxfunevals",300);
1941 nm = neldermead_configure(nm,"-tolfunrelative",10*%eps);
1942 nm = neldermead_configure(nm,"-tolxrelative",10*%eps);
1943 nm = neldermead_configure(nm,"-simplex0method","axes");
1944 nm = neldermead_configure(nm,"-simplex0length",1.0);
1945 nm = neldermead_configure(nm,"-method","variable");
1946 nm = neldermead_configure(nm,"-verbose",1);
1947 nm = neldermead_configure(nm,"-verbosetermination",1);
1948 nm = neldermead_search(nm);
1949 xopt = neldermead_get(nm,"-xopt")
1950 nm = neldermead_destroy(nm);
1951 // Contouring the function.
1952 function f = rosenbrockC ( x1 , x2 )
1954 [ f , index ] = rosenbrock ( [x1,x2]' , index )
1956 xdata = linspace ( -2 , 2 , 100 );
1957 ydata = linspace ( -1 , 2 , 100 );
1959 contour ( xdata , ydata , rosenbrockC , [2 10 100 500 1000 2000] )
1960 // Plot starting point: x0 : blue dot
1961 plot(x0(1),x0(2),"bo");
1963 plot(xopt(1),xopt(2),"r*");
1964 ]]></programlisting>
1966 The -verbose option allows to get detailed information about the
1967 current optimization process. The following is a sample output for an
1968 optimization based on the Nelder and Mead variable-shape simplex
1969 algorithm. Only the output corresponding to the iteration #156 is
1970 displayed. In order to display specific outputs (or to create specific
1971 output files and graphics), the -outputcommand option should be
1976 Iteration #156 (total = 156)
1983 Optim Simplex Object:
1984 =====================
1990 > iterations=156 >= maxiter=200
1991 > funevals=299 >= maxfunevals=300
1992 > e(x)=8.798D-15 < 2.220D-15 * 1.4142136 + 0
1993 > Terminate = F, status = continue
1994 > simplex size=2.549D-13 < 0 + 2.220D-16 * 1
1995 > Terminate = F, status = continue
1998 Function Evaluation #300, index=2, x= [1 1]
1999 xr=[1 1], f(xr)=0.000000
2001 Function Evaluation #301, index=2, x= [1 1]
2002 xc=1 1, f(xc)=0.000000
2003 > Perform Inside Contraction
2009 <title>Example #3: use output function</title>
2011 There are several ways to print intermediate messages or plots during the
2012 optimization process.
2013 The first method is to set the "-verbose" option to 1,
2014 which prints a lot of detailed information.
2015 The other method is to use the <literal>"-outputcommand"</literal> option.
2016 We can either set it to the <literal>neldermead_defaultoutput</literal> or
2017 define our own function.
2018 In this section, we present the methods based on the <literal>"-outputcommand"</literal> option.
2021 In the following example, we use the <literal>"-outputcommand"</literal>option and
2022 set it to the <literal>neldermead_defaultoutput</literal>
2023 default output function.
2024 This function prints one line by iteration, with the main optimization information.
2026 <programlisting role="example"><![CDATA[
2027 function [ f , index ] = quadratic ( x , index )
2028 f = x(1)^2 + x(2)^2;
2031 nm = neldermead_new ();
2032 nm = neldermead_configure(nm,"-numberofvariables",2);
2033 nm = neldermead_configure(nm,"-function",quadratic);
2034 nm = neldermead_configure(nm,"-x0",x0);
2035 nm = neldermead_configure(nm,"-outputcommand",neldermead_defaultoutput);
2036 nm = neldermead_search(nm);
2037 nm = neldermead_destroy(nm);
2038 ]]></programlisting>
2040 The previous script produces the following output.
2044 Iter. #0, Feval #5, Fval = 2 -- init
2045 Iter. #1, Feval #5, Fval = 2 -- init
2046 Iter. #2, Feval #6, Fval = 2 -- reflection
2047 Iter. #3, Feval #8, Fval = 0.5 -- expansion
2048 Iter. #4, Feval #9, Fval = 0.5 -- reflection
2050 Iter. #48, Feval #92, Fval = 8.557D-13 -- reflection
2051 Iter. #49, Feval #94, Fval = 7.893D-13 -- insidecontraction
2052 Iter. #50, Feval #96, Fval = 1.601D-13 -- insidecontraction
2053 Iter. #51, Feval #98, Fval = 1.291D-13 -- insidecontraction
2054 Iter. #52, Feval #100, Fval = 3.139D-14 -- outsidecontraction
2055 =================================
2057 Iter. #52, Feval #100, Fval = 3.139D-14 -- done
2060 In the following example, we define our own output function "myoutputcmd", which takes the current state
2061 as its first argument. The state is a string which can contain "init",
2062 "iter" or "done", depending on the status of the optimization. The data
2063 input argument is a tlist, which contains the data associated with the
2064 current iteration. In this case, we use the fields to print a message in
2067 <programlisting role="example"><![CDATA[
2068 function [ f , index ] = rosenbrock ( x , index )
2069 f = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
2072 function stop = myoutputcmd ( state , data )
2073 iter = data.iteration
2074 if ( state == "init" ) then
2075 mprintf ( "=================================\n");
2076 mprintf ( "Initialization\n");
2077 elseif ( state == "done" ) then
2078 mprintf ( "=================================\n");
2079 mprintf ( "End of Optimization\n");
2084 simplex = data.simplex
2085 // Simplex is a data structure, which can be managed
2086 // by the optimsimplex class.
2087 ssize = optimsimplex_size ( simplex )
2088 mprintf ( "Iter. #%3d, Feval #%3d, Fval = %.1e, x = %s, S = %.1e\n", ..
2089 iter, fc, fval, strcat(string(x)," "), ssize);
2093 nm = neldermead_new ();
2094 nm = neldermead_configure(nm,"-numberofvariables",2);
2095 nm = neldermead_configure(nm,"-function",rosenbrock);
2096 nm = neldermead_configure(nm,"-x0",[-1.2 1.0]');
2097 nm = neldermead_configure(nm,"-maxiter",200);
2098 nm = neldermead_configure(nm,"-maxfunevals",300);
2099 nm = neldermead_configure(nm,"-tolfunrelative",10*%eps);
2100 nm = neldermead_configure(nm,"-tolxrelative",10*%eps);
2101 nm = neldermead_configure(nm,"-simplex0method","axes");
2102 nm = neldermead_configure(nm,"-simplex0length",1.0);
2103 nm = neldermead_configure(nm,"-method","variable");
2104 nm = neldermead_configure(nm,"-verbose",0);
2105 nm = neldermead_configure(nm,"-verbosetermination",0);
2106 nm = neldermead_configure(nm,"-outputcommand",myoutputcmd);
2107 nm = neldermead_search(nm);
2108 nm = neldermead_destroy(nm);
2109 ]]></programlisting>
2111 The previous script produces the following output.
2114 =================================
2116 Iter. # 0, Feval # 5, Fval = 2.4e+001, x = -1.2 1, S = 1.0e+000
2117 Iter. # 1, Feval # 5, Fval = 2.4e+001, x = -1.2 1, S = 1.0e+000
2118 Iter. # 2, Feval # 7, Fval = 2.4e+001, x = -1.2 1, S = 1.0e+000
2119 Iter. # 3, Feval # 9, Fval = 2.4e+001, x = -1.2 1, S = 1.0e+000
2120 Iter. # 4, Feval # 11, Fval = 1.0e+001, x = -1.0125 0.78125, S = 6.0e-001
2121 Iter. # 5, Feval # 13, Fval = 4.7e+000, x = -1.028125 1.1328125, S = 3.5e-001
2123 Iter. #155, Feval #297, Fval = 2.0e-026, x = 1 1, S = 4.6e-013
2124 Iter. #156, Feval #299, Fval = 6.9e-027, x = 1 1, S = 2.5e-013
2125 Iter. #157, Feval #301, Fval = 6.0e-027, x = 1 1, S = 2.8e-013
2126 =================================
2128 Iter. #157, Feval #301, Fval = 6.0e-027, x = 1 1, S = 2.8e-013
2131 As another example of use, we could format the message so
2132 that it uses LaTeX formatting rules, which may allow the user to directly
2133 copy and paste the output into a LaTeX report.
2137 <title>Example #4: Optimization with bounds</title>
2139 The <literal>neldermead</literal> solver can optimize problems with
2141 To do this, we can use Box's algorithm, which projects the simplex
2142 into the bounds during the optimization.
2143 In this case, the initial guess must be located within the bounds.
2146 In the following example, we find the minimum of a quadratic function
2147 within given bounds.
2148 In order to compute the initial simplex, we use randomized bounds, that is,
2149 we compute k random vertices uniformly distributed within the bounds.
2150 The default value is so that the number of points is twice the number of
2151 variables of the problem.
2152 In this particular case, we have n=2 variables and k=4 vertices.
2154 <programlisting role="example"><![CDATA[
2155 function [ f , index ] = myquad ( x , index )
2160 nm = neldermead_new ();
2161 nm = neldermead_configure(nm,"-numberofvariables",2);
2162 nm = neldermead_configure(nm,"-function",myquad);
2163 nm = neldermead_configure(nm,"-x0",x0);
2164 nm = neldermead_configure(nm,"-method","box");
2165 nm = neldermead_configure(nm,"-boundsmin",[1 1]);
2166 nm = neldermead_configure(nm,"-boundsmax",[2 2]);
2167 nm = neldermead_configure(nm,"-simplex0method","randbounds");
2168 nm = neldermead_search(nm);
2169 xopt = neldermead_get(nm,"-xopt") // Should be [1 1]
2170 fopt = neldermead_get(nm,"-fopt") // Should be 2
2171 nm = neldermead_destroy(nm);
2172 ]]></programlisting>
2175 <title>Example #5: Optimization with nonlinear constraints</title>
2177 The <literal>neldermead</literal> solver can optimize problems with
2178 nonlinear constraints.
2179 In the following example, we solve Rosenbrock's Post Office problem, which
2180 has both bounds and linear constraints.
2181 In our example, we will manage the linear constraints as general non-linear
2182 constraints (i.e. the solver does not make a difference if the constraints are
2183 linear or non-linear).
2184 This example was first
2185 presented in "An automatic method for finding the greatest or least value of a function",
2187 This example was first used with the complex method of Box in
2188 "Algorithm 454: The complex method for constrained optimization" by
2189 Richardson, Kuester, 1971.
2190 Richardson and Kuester found the minimum function value
2191 F=-3456, with X1 = 24.01, X2 = 12.00, X3 = 12.00 and
2192 72 Iterations were necessary for them to get this result.
2195 In the following function, we define the function <literal>fpostoffice</literal>,
2196 which returns both the objective function <literal>f</literal> and the
2197 constraint value <literal>c</literal>.
2198 The original constraint is the "double" inequality constraint
2199 <literal>0<=x(1) + 2 * x(2) + 2 * x(3) <=72</literal>.
2200 To take this constraint into account, we turn it into two separate, positive,
2201 constraints and set <literal>c</literal> as a 1-by-2 matrix of doubles.
2203 <programlisting role="example"><![CDATA[
2204 function [ f , c , index ] = fpostoffice ( x , index )
2207 if ( index==2 | index==6 ) then
2208 f = -x(1) * x(2) * x(3)
2211 if ( index==5 | index==6 ) then
2212 c1 = x(1) + 2 * x(2) + 2 * x(3)
2217 ]]></programlisting>
2219 In the following script, we solve Rosenbrock's Post Office problem.
2220 First, we initialize the random number generator, so that the results are always the
2222 Then, we check that the cost function is correctly defined and that the
2223 constraints are satisfied at the initial guess.
2224 Then we configure the algorithm so that Box's algorithm is used and
2225 setup the bounds of the problem.
2226 We configure the parameters of the algorithm as suggested by Box.
2228 <programlisting role="example"><![CDATA[
2230 x0 = [1.0 1.0 1.0].';
2231 // Compute f(x0) : should be close to -1
2232 fx0 = fpostoffice ( x0 , 2 )
2233 // Compute the constraints: cx0 should be [5 67]
2234 [ fx0 , cx0, index ] = fpostoffice ( x0 , 6 )
2235 // Compute f(xopt) : fopt should be -3456
2236 xopt = [24 12 12].';
2237 fopt = fpostoffice ( xopt );
2238 // Setup optimization
2239 nm = neldermead_new ();
2240 nm = neldermead_configure(nm,"-numberofvariables",3);
2241 nm = neldermead_configure(nm,"-function",fpostoffice);
2242 nm = neldermead_configure(nm,"-x0",x0);
2243 nm = neldermead_configure(nm,"-maxiter",300);
2244 nm = neldermead_configure(nm,"-maxfunevals",300);
2245 nm = neldermead_configure(nm,"-method","box");
2246 nm = neldermead_configure(nm,"-boundsmin",[0.0 0.0 0.0]);
2247 nm = neldermead_configure(nm,"-boundsmax",[42.0 42.0 42.0]);
2248 // Configure like Box
2249 nm = neldermead_configure(nm,"-simplex0method","randbounds");
2250 nm = neldermead_configure(nm,"-nbineqconst",2);
2251 nm = neldermead_configure(nm,"-tolxmethod" , %f );
2252 nm = neldermead_configure(nm,"-tolsimplexizemethod",%f);
2253 nm = neldermead_configure(nm,"-boxtermination" , %t );
2254 nm = neldermead_configure(nm,"-boxtolf" , 0.001 );
2255 nm = neldermead_configure(nm,"-boxboundsalpha" , 0.0001 );
2257 // Check that the cost function is correctly connected.
2258 [ nm , result ] = neldermead_function ( nm , x0 );
2260 // Perform optimization
2261 nm = neldermead_search(nm);
2262 xcomp = neldermead_get(nm,"-xopt")
2263 // Compare with the exact optimum:
2265 fcomp = neldermead_get(nm,"-fopt")
2266 // Compare with the exact function value:
2268 nm = neldermead_destroy(nm);
2269 ]]></programlisting>
2271 In general, we should not expect too much from this algorithm with
2272 nonlinear constraints.
2273 Indeed, some cases require thousands of iterations to converge to
2274 an optimum, because the nonlinear constraints leave a too small
2275 space for the simplex to evolve.
2279 <title>Example #6: Passing extra parameters</title>
2281 In the following example, we solve a simple quadratic test case.
2282 Notice that the objective function has two extra parameters
2283 <literal>a</literal> and <literal>b</literal>.
2284 This is why the "-function" option is set as a list,
2285 where the first element is the function and the
2286 remaining elements are the extra parameters.
2288 <programlisting role="example"><![CDATA[
2289 function [ f , index ] = quadratic_ab ( x , index , a , b )
2290 f = a * x(1)^2 + b * x(2)^2;
2293 nm = neldermead_new ();
2294 nm = neldermead_configure(nm,"-numberofvariables",2);
2297 nm = neldermead_configure(nm,"-function",list(quadratic_ab,a,b));
2298 nm = neldermead_configure(nm,"-x0",x0);
2299 nm = neldermead_search(nm);
2300 xopt = neldermead_get(nm,"-xopt")
2301 nm = neldermead_destroy(nm);
2302 ]]></programlisting>
2305 <title>Example #7: Restarting without bounds</title>
2307 In the following example, we reproduce the experiment
2308 published by Ken McKinnon in 1998.
2309 For this particular function and this particular initial simplex,
2310 the Nelder-Mead algorithm converges to a nonstationnary point.
2313 We first define the objective function, the initial simplex
2314 and the expected solution of this unconstrained optimization problem.
2316 <programlisting role="example"><![CDATA[
2317 function [ f , index ] = mckinnon ( x , index )
2322 f = theta*phi*abs(x(1))^tau+x(2)*(1+x(2))
2324 f = theta*x(1)^tau+x(2)*(1+x(2))
2328 // The initial simplex
2329 lambda1 = (1.0 + sqrt(33))/8;
2330 lambda2 = (1.0 - sqrt(33))/8;
2337 // The expected solution
2340 ]]></programlisting>
2342 Then we run the algorithm two times in sequence.
2343 At the end of the first optimization process, the algorithm has
2344 converged to the point [0,0] which is nonstationnary.
2345 This is why we restart the algorithm and get the correct minimum.
2347 <programlisting role="example"><![CDATA[
2348 nm = neldermead_new ();
2349 nm = neldermead_configure(nm,"-numberofvariables",2);
2350 nm = neldermead_configure(nm,"-function",mckinnon);
2351 nm = neldermead_configure(nm,"-x0",[1.0 1.0]');
2352 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
2353 nm = neldermead_configure(nm, "-maxiter",200);
2354 nm = neldermead_configure(nm, "-maxfunevals",500);
2355 nm = neldermead_configure(nm,"-simplex0method","given");
2356 nm = neldermead_configure(nm,"-coords0",coords0);
2357 nm = neldermead_configure(nm,"-method","variable");
2359 nm = neldermead_search(nm);
2360 xopt = neldermead_get(nm,"-xopt")
2361 fopt = neldermead_get(nm,"-fopt")
2362 iterations = neldermead_get(nm,"-iterations")
2363 status = neldermead_get(nm,"-status")
2364 // Search #2: succeeds
2365 nm = neldermead_restart ( nm );
2366 xopt = neldermead_get(nm,"-xopt")
2367 fopt = neldermead_get(nm,"-fopt")
2368 iterations = neldermead_get(nm,"-iterations")
2369 status = neldermead_get(nm,"-status")
2370 nm = neldermead_destroy(nm);
2371 ]]></programlisting>
2373 We can also use the automatic stagnation detection method
2374 created by Kelley, so that the algorithm automatically restart
2375 the algorithm when needed.
2377 <programlisting role="example"><![CDATA[
2378 nm = neldermead_new ();
2379 nm = neldermead_configure(nm,"-numberofvariables",2);
2380 nm = neldermead_configure(nm,"-function",mckinnon);
2381 nm = neldermead_configure(nm,"-x0",[1.0 1.0]');
2382 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
2383 nm = neldermead_configure(nm, "-maxiter",200);
2384 nm = neldermead_configure(nm, "-maxfunevals",500);
2385 nm = neldermead_configure(nm,"-simplex0method","given");
2386 nm = neldermead_configure(nm,"-coords0",coords0);
2387 nm = neldermead_configure(nm,"-method","variable");
2388 nm = neldermead_configure(nm,"-kelleystagnationflag",%t);
2389 nm = neldermead_configure(nm,"-restartflag",%t);
2390 nm = neldermead_configure(nm,"-restartdetection","kelley");
2391 nm = neldermead_search(nm);
2392 xopt = neldermead_get(nm,"-xopt")
2393 fopt = neldermead_get(nm,"-fopt")
2394 iterations = neldermead_get(nm,"-iterations")
2395 restartnb = neldermead_get ( nm , "-restartnb" )
2396 status = neldermead_get(nm,"-status")
2397 nm = neldermead_destroy(nm);
2398 ]]></programlisting>
2400 See the demonstrations to get a graphical plot of the
2401 intermediate simplices in Mc Kinnon's experiment.
2405 <title>Example #8: Restarting with bounds</title>
2407 In the following experimeant, we solve an optimization problem
2409 We use Box's algorithm, which is the only algorithm which manages bounds.
2410 We use the randomized bounds simplex both for the initial simplex and
2411 for the restart simplex.
2413 <programlisting role="example"><![CDATA[
2414 function [ f , index ] = myquad ( x , index )
2415 f = x(1)^2 + x(2)^2 + x(3)^2
2417 x0 = [1.2 1.9,1.5].';
2422 nm = neldermead_new ();
2423 nm = neldermead_configure(nm,"-numberofvariables",3);
2424 nm = neldermead_configure(nm,"-function",myquad);
2425 nm = neldermead_configure(nm,"-x0",x0);
2426 nm = neldermead_configure(nm,"-method","box");
2427 nm = neldermead_configure(nm,"-boundsmin",[1 1 1]);
2428 nm = neldermead_configure(nm,"-boundsmax",[2 2 2]);
2429 nm = neldermead_configure(nm,"-simplex0method","randbounds");
2430 nm = neldermead_search(nm);
2431 nm = neldermead_configure(nm,"-maxiter",200);
2432 nm = neldermead_configure(nm,"-maxfunevals",200);
2433 nm = neldermead_configure(nm,"-restartsimplexmethod","randbounds");
2434 nm = neldermead_restart(nm);
2435 xopt = neldermead_get(nm,"-xopt")
2436 fopt = neldermead_get(nm,"-fopt")
2437 status = neldermead_get(nm,"-status")
2438 nm = neldermead_destroy(nm);
2439 ]]></programlisting>
2442 <title>Changes in Scilab 5.4</title>
2444 Many changes have been done in Scilab 5.4, which simplify the use of the
2445 neldermead component.
2448 Tagged -costfargument option of optimbase as obsolete: will be
2449 maintained for backward compatibility until 5.4.1.
2450 The -fun option can now be a list, where the element #1 is a
2451 function, and the elements #2 to the end are automatically appended to
2452 the calling sequence.
2453 To update your code, replace:
2456 nm = neldermead_configure(nm,"-function",myfun);
2457 nm = neldermead_configure(nm,"-costfargument",mystuff);
2463 nm = neldermead_configure(nm,"-function",list(myfun,mystuff));
2466 Tagged -outputcommandarg option of optimbase as obsolete: will be
2467 maintained for backward compatibility until 5.4.1.
2468 The -outputcommand option can now be a list, where the element #1 is
2469 a function, and the elements #2 to the end are automatically appended
2470 to the calling sequence.
2471 To update your code, replace:
2474 nm = neldermead_configure(nm,"-outputcommand",myoutputfun);
2475 nm = neldermead_configure(nm,"-outputcommandarg",mystuff);
2481 nm = neldermead_configure(nm,"-outputcommand",list(myoutputfun,mystuff));
2484 Tagged "outputfun(x,optimValues,state)" calling sequence of fminsearch
2485 as obsolete: will be maintained for backward compatibility until
2487 The new calling sequence is "stop=outputfun(x,optimValues,state)"
2488 To update your code, replace:
2491 function outfun ( x , optimValues , state )
2499 function stop = outfun ( x , optimValues , state )
2505 Tagged "myoutputfun(state,data)" calling sequence of neldermead
2506 as obsolete: will be maintained for backward compatibility until
2508 The new calling sequence is "stop=myoutputfun(state,data)"
2509 To update your code, replace:
2512 function myoutputfun ( state , data )
2520 function stop = myoutputfun ( state , data )
2526 Tagged "-myterminateflag" and "-myterminate" options as obsolete:
2527 will be maintained for backward compatibility until 5.4.1.
2528 To update your code, replace:
2531 function [ this , terminate , status ] = myoldterminate ( this , simplex )
2532 ssize = optimsimplex_size ( simplex , "sigmaplus" );
2533 if ( ssize < 1.e-2 ) then
2545 function stop = myoutputcmd ( state , data )
2546 simplex = data.simplex
2547 ssize = optimsimplex_size ( simplex , "sigmaplus" );
2548 if ( ssize < 1.e-2 ) then
2556 and replace the configuration:
2559 nm = neldermead_configure(nm,"-myterminateflag",%t);
2560 nm = neldermead_configure(nm,"-myterminate",myoldterminate);
2566 nm = neldermead_configure(nm,"-outputcommand",myoutputcmd);
2569 Tagged "-tolvarianceflag", "-tolabsolutevariance", and
2570 "-tolrelativevariance" options as obsolete:
2571 will be maintained for backward compatibility until 5.4.1.
2572 To update your code, create an output function:
2575 function stop = myoutputcmd ( state, data, tolrelativevariance, tolabsolutevariance, variancesimplex0 )
2576 simplex = data.simplex
2578 if ( state == "iter") then
2579 var = optimsimplex_fvvariance ( simplex )
2580 if ( var < tolrelativevariance * variancesimplex0 + tolabsolutevariance ) then
2587 Create the initial simplex and compute the variance
2588 of the function values:
2592 simplex0 = optimsimplex_new ( "axes" , x0.' );
2593 coords0 = optimsimplex_getallx(simplex0);
2594 variancesimplex0 = optimsimplex_fvvariance ( simplex0 );
2597 Finally, replace the configuration:
2600 nm = neldermead_configure(nm,"-tolvarianceflag",%t);
2601 nm = neldermead_configure(nm,"-tolabsolutevariance",1.e-4);
2602 nm = neldermead_configure(nm,"-tolrelativevariance",1.e-4);
2608 tolabsolutevariance = 1.e-4;
2609 tolrelativevariance = 1.e-4;
2610 stopfun = list(myoutputcmd, tolrelativevariance, tolabsolutevariance, variancesimplex0);
2611 nm = neldermead_configure(nm,"-outputcommand",stopfun);
2615 <title>Spendley et al. implementation notes</title>
2617 The original paper may be implemented with several variations, which
2618 might lead to different results. This section defines what algorithmic
2619 choices have been used.
2621 <para>The paper states the following rules.</para>
2625 "Rule 1. Ascertain the lowest reading y, of yi ... yk+1 Complete
2626 a new simplex Sp by excluding the point Vp corresponding to y, and
2627 replacing it by V* defined as above."
2632 "Rule 2. If a result has occurred in (k + 1) successive
2633 simplexes, and is not then eliminated by application of Rule 1, do not
2634 move in the direction indicated by Rule 1, or at all, but discard the
2635 result and replace it by a new observation at the same point."
2640 "Rule 3. If y is the lowest reading in So , and if the next
2641 observation made, y* , is the lowest reading in the new simplex S , do
2642 not apply Rule 1 and return to So from Sp . Move out of S, by
2643 rejecting the second lowest reading (which is also the second lowest
2649 We implement the following "rules" of the Spendley et al.
2655 Rule 1 is strictly applied, but the reflection is done by
2656 reflection the high point, since we minimize a function instead of
2657 maximizing it, like Spendley.
2662 Rule 2 is NOT implemented, as we expect that the function
2663 evaluation is not subject to errors.
2668 Rule 3 is applied, ie reflection with respect to next to high
2674 The original paper does not mention any shrink step. When the
2675 original algorithm cannot improve the function value with reflection
2676 steps, the basic algorithm stops. In order to make the current
2677 implementation of practical value, a shrink step is included, with
2678 shrinkage factor sigma. This perfectly fits into to the spirit of the
2679 original paper. Notice that the shrink step make the rule #3 (reflection
2680 with respect to next-to-worst vertex) unnecessary. Indeed, the minimum
2681 required steps are the reflection and shrinkage. Never the less, the rule
2682 #3 has been kept in order to make the algorithm as close as it can be to
2687 <title>Nelder-Mead implementation notes</title>
2689 The purpose of this section is to analyse the current implementation of Nelder-Mead's algorithm.
2692 The algorithm that we use is described in "Iterative Methods for Optimization" by C. T. Kelley.
2695 The original paper uses a "greedy" expansion, in which the expansion point
2696 is accepted whatever its function value. The current implementation,
2697 as most implementations, uses the expansion point only if it improves
2698 over the reflection point, that is,
2703 if fe<fr, then the expansion point is accepted,
2708 if not, the reflection point is accepted.
2713 The termination criteria suggested by Nelder and Mead is based on an
2714 absolute tolerance on the standard deviation of the function values in the simplex.
2715 We provide this original termination criteria with the <literal>-tolvarianceflag</literal>
2716 option, which is disabled by default.
2720 <title>Box's complex algorithm implementation notes</title>
2722 In this section, we analyse the current implementation of Box's complex method.
2725 The initial simplex can be computed as in Box's paper, but this may not
2726 be safe. In his paper, Box suggest that if a vertex of the initial simplex
2727 does not satisfy the non linear constraints, then it should be "moved halfway
2728 toward the centroid of those points already selected". This behaviour
2729 is available when the <literal>-scalingsimplex0</literal> option is set to
2730 <literal>"tocenter"</literal>. It may happen, as suggested by Guin, that
2731 the centroid is not feasible. This may happen if the constraints are not
2732 convex. In this case, the initial simplex cannot be computed. This is why
2733 we provide the <literal>"tox0"</literal> option, which allows to compute the
2734 initial simplex by scaling toward the initial guess, which is always feasible.
2737 In Box's paper, the scaling into the non linear constraints is performed
2738 "toward" the centroid, that is, by using a scaling factor equal to 0.5.
2739 This default scaling factor might be sub-optimal in certain situations.
2740 This is why we provide the <literal>-boxineqscaling</literal> option,
2741 which allows to configure the scaling factor.
2744 In Box's paper, whether we are concerned with the initial simplex or with the
2745 simplex at a given iteration, the scaling for the non linear constraints is performed
2746 without end. This is because Box's hypothesis is that "ultimately, a satisfactory
2747 point will be found". As suggested by Guin, if the process fails, the algorithm
2748 goes into an infinite loop. In order to avoid this, we perform the scaling until
2749 a minimum scaling value is reached, as defined by the <literal>-guinalphamin</literal>
2753 We have taken into account for the comments by Guin, but it should be emphasized
2754 that the current implementation is still as close as possible to Box's
2755 algorithm and is not Guin's algorithm. More precisely, during the iterations,
2756 the scaling for the non linear constraints is still performed toward the centroid,
2757 be it feasible or not.
2761 <title>Bibliography</title>
2763 "Sequential Application of Simplex Designs in Optimisation and
2764 Evolutionary Operation", Spendley, W. and Hext, G. R. and Himsworth, F.
2765 R., American Statistical Association and American Society for Quality,
2769 "A Simplex Method for Function Minimization", Nelder, J. A. and
2770 Mead, R., The Computer Journal, 1965
2773 "A New Method of Constrained Optimization and a Comparison With
2774 Other Methods", M. J. Box, The Computer Journal 1965 8(1):42-52, 1965 by
2775 British Computer Society
2778 "Discussion and correspondence: modification of the complex method
2779 of constrained optimization", J. A. Guin, The Computer Journal,
2783 "Detection and Remediation of Stagnation in the Nelder--Mead
2784 Algorithm Using a Sufficient Decrease Condition", Kelley C. T., SIAM J. on
2788 "Iterative Methods for Optimization", C. T. Kelley, SIAM Frontiers
2789 in Applied Mathematics, 1999
2792 "Algorithm AS47 - Function minimization using a simplex procedure",
2793 O'Neill, R., Applied Statistics, 1971
2796 "Nelder Mead's User Manual", Consortium Scilab - Digiteo, Michael
2800 Ken McKinnon, Convergence of the Nelder-Mead simplex method to a nonstationary point,
2801 SIAM Journal on Optimization, Volume 9, Number 1, 1998, pages 148-158.
2804 <refsection role="see also">
2805 <title>See Also</title>
2806 <simplelist type="inline">
2808 <link linkend="optimbase_new">optimbase_new</link>
2811 <link linkend="optimsimplex_new">optimsimplex_new</link>
2814 <link linkend="optimbase_overview">optimbase</link>
2817 <link linkend="optimsimplex_overview">optimsimplex</link>
2820 <link linkend="nmplot">nmplot</link>
projects / scilab.git / blob