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-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).
185 <term>-function</term>
188 a function or a list, the objective function.
189 This function computes the value
190 of the cost and the non linear constraints, if
194 There is no default value, i.e. the user must
195 provide <literal>f</literal>.
198 See below for the details of the communication
199 between the optimization system and the cost
208 a n-by-1 matrix of doubles, where
209 n is the number of variables, the initial guess.
212 There is no default value, i.e. the user must
213 provide <literal>x0</literal>.
219 <emphasis>Output.</emphasis>
223 <term>-outputcommand</term>
226 a function or a list, a function which is called back for output.
229 The default output function is empty, meaning that there is
233 See below for the details of the communication
234 between the optimization system and the output command
240 <term>-storehistory</term>
243 a 1-by-1 matrix of booleans, set to %t to enable the history storing (default
249 <term>-verbose</term>
252 a 1-by-1 matrix of doubles, positive, integer value, set to 1 to
253 enable verbose logging (default verbose = 0).
258 <term>-verbosetermination</term>
261 a 1-by-1 matrix of doubles, positive, integer value,
262 set to 1 to enable verbose termination logging (default verbosetermination = 0).
268 <emphasis>Bounds and constraints.</emphasis>
272 <term>-boundsmin</term>
275 a n-by-1 matrix of doubles, the minimum bounds for the parameters
276 where n is the number of variables (default boundsmin = [], i.e. there are no
282 <term>-boundsmax</term>
285 a n-by-1 matrix of doubles, the maximum bounds for the parameters
286 where n is the number of variables (default boundsmax = [], i.e. there are no
292 <term>-nbineqconst</term>
295 a 1-by-1 matrix of doubles, positive, integer value,
296 the number of inequality constraints (default nbineqconst = 0).
302 <emphasis>Initialization.</emphasis>
306 <term>-simplex0method</term>
309 a 1-by-1 matrix of strings, the method to use to compute the initial
310 simplex (default simplex0method = "axes").
311 The first vertex in the simplex is always the initial
312 guess associated with the -x0 option. The following
313 methods are available :
320 the coordinates associated with the -coords0
321 option are used to compute the initial simplex, with
322 arbitrary number of vertices.
325 This allow the user to setup the initial
326 simplex by a specific method which is not provided
327 by the current component (for example with a simplex
328 computed from a design of experiments). This allows
329 also to configure the initial simplex so that a
330 specific behaviour of the algorithm an be reproduced
331 (for example the Mac Kinnon test case).
334 The given matrix is expected to have n rows
335 and k columns, where n is the dimension of the
336 problem and k is the number of vertices.
344 the simplex is computed from the coordinate
345 axes and the length associated with the
346 -simplex0length option.
351 <term>"spendley"</term>
354 the simplex is computed so that it is regular
355 with the length associated with the -simplex0length
356 option (i.e. all the edges have the same
362 <term>"pfeffer"</term>
365 the simplex is computed from a heuristic, in
366 the neighborhood of the initial guess. This initial
367 simplex depends on the -simplex0deltausual and
373 <term>"randbounds"</term>
376 the simplex is computed from the bounds and a
377 random number. This option is available only if
378 bounds are available : if bounds are not available,
379 an error is generated. This method is usually
380 associated with Box's algorithm. The number of
381 vertices in the simplex is taken from the
390 <term>-coords0</term>
393 a nbve-by-n matrix of doubles, where nbve is the number of vertices and n is
394 the number of variables, the coordinates of the vertices of the initial
395 simplex (default coords0=[]).
396 If the -simplex0method option is set to
397 "given", these coordinates are used to compute the
403 <term>-simplex0length</term>
406 a 1-by-1 matrix of doubles, the length to use when the initial simplex is
407 computed with the "axes" or "spendley" methods (default simplex0length = 1).
408 If the initial simplex is computed from "spendley" method, the
409 length is expected to be a scalar value. If the initial
410 simplex is computed from "axes" method, it may be either
411 a scalar value or a vector of values, with rank n, where
412 n is the number of variables.
417 <term>-simplex0deltausual</term>
420 a 1-by-1 matrix of doubles, the relative delta for non-zero parameters in
421 "pfeffer" method (default simplex0deltausual = 0.05).
426 <term>-simplex0deltazero</term>
429 a 1-by-1 matrix of doubles, the absolute delta for non-zero parameters in
430 "pfeffer" method (default simplex0deltazero = 0.0075).
436 <emphasis>Termination.</emphasis>
440 <term>-maxfunevals</term>
443 a 1-by-1 matrix of doubles, positive, integer value,
444 the maximum number of function evaluations
445 (default maxfunevals = 100).
446 If this criteria is triggered, the status of the optimization is set to
452 <term>-maxiter</term>
455 a 1-by-1 matrix of doubles, positive, integer value,
456 the maximum number of iterations (default maxiter = 100).
457 If this criteria is triggered, the status of the
458 optimization is set to "maxiter".
463 <term>-tolfunabsolute</term>
466 a 1-by-1 matrix of doubles, positive, the absolute tolerance for the function value
467 (default tolfunabsolute = 0).
472 <term>-tolfunrelative</term>
475 a 1-by-1 matrix of doubles, positive, the relative tolerance for the function value
476 (default tolfunrelative = %eps).
481 <term>-tolfunmethod</term>
484 a 1-by-1 matrix of booleans, set to %t to enable termination with
485 tolerance on function value (default tolfunmethod = %f).
488 If this criteria is triggered, the
489 status of the optimization is set to "tolf".
494 <term>-tolxabsolute</term>
497 a 1-by-1 matrix of doubles, positive, the absolute tolerance
498 on x (default tolxabsolute = 0).
503 <term>-tolxrelative</term>
506 a 1-by-1 matrix of doubles, positive, the relative
507 tolerance on x (default tolxrelative = sqrt(%eps)).
512 <term>-tolxmethod</term>
515 a 1-by-1 matrix of booleans, set to %t to enable the tolerance on x in the
516 termination criteria (default tolxmethod = %t).
519 If this criteria is triggered, the
520 status of the optimization is set to "tolx".
525 <term>-tolsimplexizemethod</term>
528 a 1-by-1 matrix of booleans, set to %f to disable the tolerance on the simplex
529 size (default tolsimplexizemethod = %t). If this criteria is
530 triggered, the status of the optimization is set to
534 When this criteria is enabled, the values of the
535 options -tolsimplexizeabsolute and
536 -tolsimplexizerelative are used in the termination
537 criteria. The method to compute the size is the
543 <term>-tolsimplexizeabsolute</term>
546 a 1-by-1 matrix of doubles, positive, the absolute tolerance on the
547 simplex size (default tolsimplexizeabsolute = 0).
552 <term>-tolsimplexizerelative</term>
555 a 1-by-1 matrix of doubles, positive, the relative tolerance on the
556 simplex size (default tolsimplexizerelative = %eps).
561 <term>-tolssizedeltafvmethod</term>
564 a 1-by-1 matrix of booleans, set to %t to enable the termination criteria based
565 on the size of the simplex and the difference of
566 function value in the simplex (default tolssizedeltafvmethod = %f).
567 If this criteria is triggered, the status of the
568 optimization is set to "tolsizedeltafv".
571 This termination criteria uses the values of the
572 options -tolsimplexizeabsolute and -toldeltafv. This
573 criteria is identical to Matlab's fminsearch.
578 <term>-toldeltafv</term>
581 a 1-by-1 matrix of doubles, positive, the absolute tolerance on the difference between
582 the highest and the lowest function values (default toldeltafv = %eps).
588 <emphasis>Algorithm.</emphasis>
595 a 1-by-1 matrix of strings, the name of the algorithm to use (default method = "variable").
596 The following methods are available :
603 the Spendley et al. fixed simplex shape
604 algorithm. This algorithm is for unconstrained
605 problems (i.e. bounds and non linear constraints are
606 not taken into account)
611 <term>"variable"</term>
614 the Nelder-Mead variable simplex shape
615 algorithm. This algorithm is for unconstrained
616 problems (i.e. bounds and non linear constraints are
617 not taken into account)
625 the Box complex algorithm. This algorithm
626 takes into account bounds and nonlinear inequality
635 the user-defined algorithm, associated with the
636 <literal>-mymethod</literal>option. See below for details.
644 <term>-mymethod</term>
647 a function, a user-derined simplex algorithm. See below for
648 details (default is empty).
654 <emphasis>Options of the "variable" algorithm.</emphasis>
661 a 1-by-1 matrix of doubles, the reflection coefficient. This parameter is used
662 when the -method option is set to "fixed" or "variable" (default rho = 1).
670 a 1-by-1 matrix of doubles, the expansion coefficient. This parameter is used
671 when the -method option is set to "variable" (default chi = 2).
679 a 1-by-1 matrix of doubles, the contraction coefficient. This parameter is
680 used when the -method option is set to "variable" (default gamma = 0.5).
688 a 1-by-1 matrix of doubles, the shrinkage coefficient. This parameter is used
689 when the -method option is set to "fixed" or "variable" (default sigma = 0.5).
695 <emphasis>Option of "box" algorithm.</emphasis>
699 <term>-scalingsimplex0</term>
702 a 1-by-1 matrix of strings, the algorithm used to scale the initial simplex into
703 the nonlinear constraints (default scalingsimplex0 = "tox0").
704 The following two algorithms are provided :
709 "tox0": scales the vertices toward the initial guess.
714 "tocenter": scales the vertices toward the centroid, as recommended by Box.
719 If the centroid happens to be unfeasible, because the constraints are
720 not convex, the scaling of the initial simplex toward the centroid
721 may fail. Since the initial guess is always feasible, scaling toward
722 the initial guess cannot fail.
723 The default value is "tox0".
728 <term>-boxnbpoints</term>
731 a 1-by-1 matrix of doubles, positive, integer value, the number
732 of points in the initial simplex, when
733 the -simplex0method is set to <literal>"randbounds"</literal>
734 (default boxnbpoints = 2*n, where n is the number of variables of the problem).
735 The value of this option is also use to update the simplex when a
736 restart is performed and the <literal>-restartsimplexmethod</literal>
737 option is set to <literal>"randbounds"</literal>.
742 <term>-boxineqscaling</term>
745 a 1-by-1 matrix of doubles, in the interval [0,1], the scaling coefficient used to scale the trial
746 point for function improvement or into the constraints of Box's
747 algorithm (default boxineqscaling = 0.5).
752 <term>-guinalphamin</term>
755 a 1-by-1 matrix of doubles, positive, the minimum value of alpha when scaling the vertices of the
756 simplex into nonlinear constraints in Box's algorithm (default guinalphamin = 1.e-5).
761 <term>-boxreflect</term>
764 a 1-by-1 matrix of doubles, positive, the reflection factor in
765 Box's algorithm (default = 1.3).
770 <term>-boxtermination</term>
773 a 1-by-1 matrix of booleans, set to %t to enable Box termination criteria (default boxtermination = %f).
778 <term>-boxtolf</term>
781 a 1-by-1 matrix of doubles, positive, the absolute tolerance on difference of function values in the
782 simplex, suggested by Box (default boxtolf = 1.e-5). This tolerance is used if
783 the -boxtermination option is set to %t.
789 <term>-boxnbmatch</term>
792 a 1-by-1 matrix of doubles, positive, integer value,
793 the number of consecutive match of Box termination criteria (default boxnbmatch = 5).
798 <term>-boxboundsalpha</term>
801 a 1-by-1 matrix of doubles, positive, the parameter used to project the vertices into the
802 bounds in Box's algorithm (default boxboundsalpha = 1.e-6).
808 <emphasis>Auto-Restart.</emphasis>
812 <term>-kelleystagnationflag</term>
815 a 1-by-1 matrix of booleans, set to %t to enable the termination criteria using
816 Kelley's stagnation detection, based on sufficient
817 decrease condition (default kelleystagnationflag = %f). If this
818 criteria is triggered, the status of the optimization is
819 set to "kelleystagnation".
824 <term>-kelleynormalizationflag</term>
827 a 1-by-1 matrix of booleans, set to %f to disable the normalization of the
828 alpha coefficient in Kelley's stagnation detection, i.e.
829 use the value of the option -kelleystagnationalpha0 as
830 is (default kelleynormalizationflag = %t, i.e. the simplex gradient of
831 the initial simplex is taken into account in the
832 stagnation detection).
837 <term>-kelleystagnationalpha0</term>
840 a 1-by-1 matrix of doubles, the parameter used in Kelley's stagnation
841 detection (default kelleystagnationalpha0 = 1.e-4).
846 <term>-restartflag</term>
849 a 1-by-1 matrix of booleans, set to %t to enable the automatic restart of the
850 algorithm (default restartflag = %f).
855 <term>-restartdetection</term>
858 a 1-by-1 matrix of strings, the method to detect if the automatic restart must
859 be performed (default restartdetection = "oneil").
860 The following methods are available:
864 <term>"oneill"</term>
867 the factorial local optimality test by O'Neill
868 is used. If the test finds a local point which is
869 better than the computed optimum, a restart is
875 <term>"kelley"</term>
878 the sufficient decrease condition by O'Neill
879 is used. If the test finds that the status of the
880 optimization is "kelleystagnation", a restart is
881 performed. This status may be generated if the
882 -kelleystagnationflag option is set to %t.
887 <para>The default method is "oneill".</para>
891 <term>-restartmax</term>
894 a 1-by-1 matrix of doubles, the maximum number of restarts, when automatic
895 restart is enabled via the -restartflag option (default
901 <term>-restarteps</term>
904 a 1-by-1 matrix of doubles, the relative epsilon value used to check for
905 optimality in the factorial O'Neill restart detection (default restarteps = %eps).
910 <term>-restartstep</term>
913 a 1-by-1 or a n-by-1 matrix of doubles, positive, where n is the number of
914 variables in the problem, the absolute step length used to check for
915 optimality in the factorial O'Neill restart detection (default restartstep = 1).
920 <term>-restartsimplexmethod</term>
923 a 1-by-1 matrix of strings, the method to compute the initial simplex after a
924 restart (default restartsimplexmethod = "oriented"). The following methods are available.
931 the coordinates associated with the -coords0
932 option are used to compute the initial simplex, with
933 arbitrary number of vertices.
936 This allow the user to setup the initial
937 simplex by a specific method which is not provided
938 by the current component (for example with a simplex
939 computed from a design of experiments). This allows
940 also to configure the initial simplex so that a
941 specific behaviour of the algorithm an be reproduced
942 (for example the Mc Kinnon test case).
945 The given matrix is expected to have n rows
946 and k columns, where n is the dimension of the
947 problem and k is the number of vertices.
955 the simplex is computed from the coordinate
956 axes and the length associated with the
957 -simplex0length option.
962 <term>"spendley"</term>
965 the simplex is computed so that it is regular
966 with the length associated with the -simplex0length
967 option (i.e. all the edges have the same
973 <term>"pfeffer"</term>
976 the simplex is computed from a heuristic, in
977 the neighborhood of the initial guess. This initial
978 simplex depends on the -simplex0deltausual and
984 <term>"randbounds"</term>
987 the simplex is computed from the bounds and a
988 random number. This option is available only if
989 bounds are available : if bounds are not available,
990 an error is generated. This method is usually
991 associated with Box's algorithm. The number of
992 vertices in the simplex is taken from the
998 <term>"oriented"</term>
1001 the simplex is computed so that it is
1002 oriented, as suggested by C.T. Kelley.
1007 <para>The default method is "oriented".</para>
1016 <para>the value.</para>
1023 <term xml:id="neldermead_cget">value = neldermead_cget (this,key)</term>
1026 Get the value for the given key. If the key is unknown,
1033 <para>The current object.</para>
1040 the name of the key to quiery. The list of available
1041 keys is the same as for the neldermead_configure
1050 <term xml:id="neldermead_get">value = neldermead_get ( this , key )</term>
1053 Get the value for the given key. If the key is unknown,
1057 Most fields are available only after an optimization has
1058 been performed with one call to the neldermead_search
1065 <para>The current object.</para>
1071 <para>the key to get.</para>
1072 <para>The following keys are available :</para>
1075 <term>-funevals</term>
1077 <para>the number of function evaluations</para>
1081 <term>-iterations</term>
1083 <para>the number of iterations</para>
1090 the x optimum, as a n x 1 column vector, where n
1091 is the number of variables.
1098 <para>the optimum cost function value</para>
1102 <term>-historyxopt</term>
1105 an array, with nbiter values, containing the
1106 history of x during the iterations.
1109 This array is available after optimization if the
1110 history storing was enabled with the -storehistory
1116 <term>-historyfopt</term>
1119 an array, with nbiter values, containing the
1120 history of the function value during the
1124 This array is available after optimization if the
1125 history storing was enabled with the -storehistory
1133 <para>the function value for the initial guess</para>
1137 <term>-status</term>
1140 a string containing the status of the
1141 optimization. See below for details about the
1142 optimization status.
1147 <term>-historysimplex</term>
1150 a matrix containing the history of the simplex
1151 during the iterations. This matrix has rank nbiter x
1152 nbve x n, where nbiter is the number of iterations, nbve
1153 is the number of vertices in the simplex and n is the
1154 number of variables.
1159 <term>-simplex0</term>
1162 the initial simplex. This is a simplex object,
1163 which is suitable for processing with the optimsimplex
1169 <term>-simplexopt</term>
1172 the optimum simplex. This is a simplex object,
1173 which is suitable for processing with the optimsimplex
1179 <term>-restartnb</term>
1181 <para>the number of actual restarts performed.</para>
1191 <term xml:id="neldermead_search">this = neldermead_search ( this )</term>
1194 Performs the optimization associated with the method
1195 associated with the -method option and find the optimum.
1201 <para>The current object.</para>
1206 If the -restartflag option is enabled, automatic restarts are
1207 performed, based on the -restartdetection option.
1212 <term xml:id="neldermead_restart">this = neldermead_restart ( this )</term>
1215 Restarts the optimization by updating the simplex and
1216 performing a new search.
1222 <para>The current object.</para>
1229 <term xml:id="neldermead_function">[ this , result ] = neldermead_function ( this , x )</term>
1231 <para>Call the cost function and return the value.</para>
1236 <para>The current object.</para>
1242 <para>the point where the function is to be evaluated</para>
1249 optional, a flag to pass to the cost function (default
1250 = 1). See the section "The cost function" for available values
1259 <term xml:id="neldermead_defaultoutput">stop = neldermead_defaultoutput(state, data)</term>
1262 Prints messages at an iteration.
1265 This function provides a default implementation for the output function.
1266 There is one line by iteration, presenting the number of iterations, the
1267 number of function evaluations, the current function value and the
1268 current algorithm step.
1271 See "The output function" section below for a description of the input and
1273 See in the Examples section below for examples of this function.
1280 <title>The cost function</title>
1282 The option <literal>-function</literal> allows to configure the cost function. The cost
1283 function is used to compute the objective function value <literal>f</literal>.
1284 If the <literal>-nbineqconst</literal> option is configured to a non-zero value, the cost function
1285 must also compute the value of the nonlinear, positive, inequality constraints <literal>c</literal>.
1286 Depending of these options, the cost function can have one of the following headers :
1289 [ f , index ] = costf ( x , index )
1290 [ f , c , index ] = costf ( x , index )
1297 <para>the current point, as a column vector</para>
1303 <para>optional, an integer representing the value to compute</para>
1309 <para>the value of the cost function</para>
1315 <para>the value of the non-linear, positive, inequality constraints</para>
1320 The index input parameter tells to the cost function what is expected
1321 in the output arguments. It has the following meaning
1325 <term>index = 2</term>
1328 compute <literal>f</literal>
1333 <term>index = 5</term>
1336 compute <literal>c</literal>
1341 <term>index = 6</term>
1344 compute <literal>f</literal> and <literal>c</literal>
1350 In the most simplex case, there is no additionnal cost function argument and no nonlinear constraints.
1351 In this case, the cost function is expected to have the following header
1354 [ f , index ]= myfunction ( x , index )
1357 It might happen that the function requires additionnal arguments to be evaluated.
1358 In this case, we can use the following feature.
1359 The argument <literal>fun</literal> can also be the list <literal>(myfun,a1,a2,...)</literal>.
1360 In this case <literal>myfun</literal>, the first element in the list, must be a function and must
1363 [ f , index ] = myfun ( x , index , a1, a2, ... )
1364 [ f , c , index ] = myfun ( x , index , a1, a2, ...)
1366 where the input arguments <literal>a1, a2, ...</literal>
1367 are automatically appended at the end of the calling sequence.
1371 <title>The output function</title>
1373 The option -outputcommand allows to configure a command which is
1374 called back at the start of the optimization, at each iteration and at the
1375 end of the optimization.
1377 <para>The output function must have the following header</para>
1379 stop = outputcmd(state, data)
1387 a string representing the current state of the algorithm.
1388 Available values are "init", "iter", "done".
1395 <para>a data structure containing at least the following entries</para>
1400 <para>the current optimum</para>
1406 <para>the current function value</para>
1410 <term>iteration</term>
1412 <para>the current iteration index</para>
1416 <term>funccount</term>
1418 <para>the number of function evaluations</para>
1422 <term>simplex</term>
1424 <para>the current simplex</para>
1431 the previous step in the algorithm. The following values
1432 are available : "init", "done", "reflection", "expansion",
1433 "insidecontraction", "outsidecontraction", "reflectionnext",
1445 a 1-by-1 matrix of booleans, set to true to stop the
1446 algorithm, set to false to continue the optimization.
1452 It might happen that the output function requires additionnal arguments to be evaluated.
1453 In this case, we can use the following feature.
1454 The argument <literal>outputcmd</literal> can also be the list <literal>(outf,a1,a2,...)</literal>.
1455 In this case <literal>outf</literal>, the first element in the list, must be a function and must
1458 stop = outf ( state, data, a1, a2, ... )
1460 where the input arguments <literal>a1, a2, ...</literal>
1461 are automatically appended at the end of the calling sequence.
1464 If the output function sets the <literal>stop</literal> variable to true,
1465 then the optimization alorithm stops and the status of the optimization is
1466 set to <literal>"userstop"</literal>.
1470 <title>Termination</title>
1472 The current component takes into account for several generic
1473 termination criterias.
1475 <para>The following termination criterias are enabled by default :</para>
1478 <para>-maxiter,</para>
1481 <para>-maxfunevals,</para>
1484 <para>-tolxmethod.</para>
1487 <para>-tolsimplexizemethod.</para>
1491 The optimization_terminate function uses a set of rules to compute
1492 if the termination occurs, which leads to an optimization status which is
1493 equal to one of the following : "continue", "maxiter", "maxfunevals",
1494 "tolf", "tolx", "tolsize", "tolsizedeltafv",
1495 "kelleystagnation", "tolboxf", "tolvariance". The value of the status may also
1496 be a user-defined string, in the case where a user-defined termination
1497 function has been set.
1499 <para>The following set of rules is examined in this order.</para>
1503 By default, the status is <literal>"continue"</literal> and the <literal>terminate</literal> flag is
1509 The number of iterations is examined and compared to the
1510 <literal>-maxiter</literal> option : if the following condition
1513 iterations >= maxiter
1516 is true, then the status is set to "maxiter" and terminate is
1522 The number of function evaluations and compared to the
1523 <literal>-maxfunevals</literal> option is examined : if the following condition
1526 funevals >= maxfunevals
1529 is true, then the status is set to <literal>"maxfuneval"</literal> and <literal>terminate</literal> is
1535 The tolerance on function value is examined depending on the
1536 value of the <literal>-tolfunmethod</literal>.
1542 <para>then the criteria is just ignored.</para>
1548 <para>if the following condition</para>
1550 abs(currentfopt) < tolfunrelative * abs(previousfopt) + tolfunabsolute
1553 is true, then the status is set to "tolf" and terminate is
1560 The relative termination criteria on the function value works
1561 well if the function value at optimum is near zero. In that case, the
1562 function value at initial guess fx0 may be used as
1563 <literal>previousfopt</literal>.
1566 This criteria is sensitive to the <literal>-tolfunrelative</literal>
1567 and <literal>-tolfunabsolute</literal> options.
1570 The absolute termination criteria on the function value works if
1571 the user has an accurate idea of the optimum function value.
1576 The tolerance on x is examined depending on the value of the
1583 <para>then the criteria is just ignored.</para>
1589 <para>if the following condition</para>
1591 norm(xopt - previousxopt) < tolxrelative * norm(xopt) + tolxabsolute
1594 is true, then the status is set to <literal>"tolx"</literal> and <literal>terminate</literal> is
1601 This criteria is sensitive to the <literal>-tolxrelative</literal>
1602 and <literal>-tolxabsolute</literal> options.
1605 The relative termination criteria on x works well if x at
1606 optimum is different from zero. In that case, the condition measures
1607 the distance between two iterates.
1610 The absolute termination criteria on x works if the user has an
1611 accurate idea of the scale of the optimum x. If the optimum x is near
1612 0, the relative tolerance will not work and the absolute tolerance is
1618 The tolerance on simplex size is examined depending on
1619 the value of the <literal>-tolsimplexizemethod</literal> option.
1625 <para>then the criteria is just ignored.</para>
1631 <para>if the following condition</para>
1633 ssize < tolsimplexizerelative * simplexsize0 + tolsimplexizeabsolute
1636 is true where <literal>simplexsize0</literal> is the size of the simplex at
1637 iteration 0, then the <literal>status</literal> is set to <literal>"tolsize"</literal>
1638 and <literal>terminate</literal> is set to %t.
1641 The size of the simplex is computed from the "sigmaplus" method of the
1642 <literal>optimsimplex</literal> component.
1643 This criteria is sensitive to the <literal>-tolsimplexizeabsolute</literal> and
1644 the <literal>-tolsimplexizerelative</literal> options.
1652 The absolute tolerance on simplex size and absolute difference
1653 of function value is examined depending on the value of the
1654 -tolssizedeltafvmethod option.
1660 <para>then the criteria is just ignored.</para>
1666 <para>if both the following conditions</para>
1668 ssize < tolsimplexizeabsolute
1671 shiftfv < toldeltafv
1674 is true where <literal>ssize</literal> is the current simplex size and
1675 <literal>shiftfv</literal> is the absolute value of the difference of function
1676 value between the highest and lowest vertices, then the status
1677 is set to <literal>"tolsizedeltafv"</literal> and <literal>terminate</literal> is set to %t.
1685 The stagnation condition based on Kelley sufficient decrease
1686 condition is examined depending on the value of the
1687 <literal>-kelleystagnationflag</literal> option.
1693 <para>then the criteria is just ignored.</para>
1699 <para>if the following condition</para>
1701 newfvmean <= oldfvmean - alpha * sg' * sg
1704 is true where <literal>newfvmean</literal> (resp. <literal>oldfvmean</literal>) is the function
1705 value average in the current iteration (resp. in the previous
1706 iteration), then the status is set to "kelleystagnation" and
1707 terminate is set to %t. Here, <literal>alpha</literal> is a non-dimensional
1708 coefficient and <literal>sg</literal> is the simplex gradient.
1716 The termination condition suggested by Box
1717 is examined depending on the value of the
1718 -boxtermination option.
1724 <para>then the criteria is just ignored.</para>
1730 <para>if both the following conditions</para>
1732 shiftfv < boxtolf
1735 boxkount == boxnbmatch
1738 is true where <literal>shiftfv </literal>is the difference of function value
1739 between the best and worst vertices, and <literal>boxkount</literal> is the number of
1740 consecutive iterations where this criteria is met,
1741 then the status is set to "tolboxf" and
1742 terminate is set to %t.
1743 Here, the <literal>boxtolf</literal> parameter is the value associated with
1744 the <literal>-boxtolf</literal> option
1745 and is a user-defined absolute tolerance on the function value.
1746 The <literal>boxnbmatch</literal> parameter is the value associated with
1747 the <literal>-boxnbmatch</literal> option
1748 and is the user-defined number of consecutive match.
1756 The termination condition based on the variance of the function values in the simplex
1757 is examined depending on the value of the
1758 <literal>-tolvarianceflag</literal> option.
1764 <para>then the criteria is just ignored.</para>
1770 <para>if the following condition</para>
1772 var < tolrelativevariance * variancesimplex0 + tolabsolutevariance
1775 is true where <literal>var </literal>is the variance of the function values
1777 then the status is set to "tolvariance" and
1778 terminate is set to %t.
1779 Here, the <literal>tolrelativevariance</literal> parameter is the value associated with
1780 the <literal>-tolrelativevariance</literal> option
1781 and is a user-defined relative tolerance on the variance of the function values.
1782 The <literal>tolabsolutevariance</literal> parameter is the value associated with
1783 the <literal>-tolabsolutevariance</literal> option
1784 and is the user-defined absolute tolerance of the variance of the function values.
1793 <title>Kelley's stagnation detection</title>
1795 The stagnation detection criteria suggested by Kelley is based on a
1796 sufficient decrease condition, which requires a parameter alpha > 0 to
1797 be defined. The -kelleynormalizationflag option allows to configure the
1798 method to use to compute this alpha parameter : two methods are available,
1799 where each method corresponds to a different paper by Kelley :
1803 <term>constant</term>
1806 In "Detection and Remediation of Stagnation in the
1807 Nelder--Mead Algorithm Using a Sufficient Decrease Condition",
1808 Kelley uses a constant alpha, with the suggested value 1.e-4, which
1809 is is typical choice for line search method.
1814 <term>normalized</term>
1817 in "Iterative Methods for Optimization", Kelley uses a
1818 normalized alpha, computed from the following formula
1821 alpha = alpha0 * sigma0 / nsg
1824 where sigma0 is the size of the initial simplex and nsg is the
1825 norm of the simplex gradient for the initial guess point.
1832 <title>O'Neill factorial optimality test</title>
1834 In "Algorithm AS47 - Function minimization using a simplex
1835 procedure", R. O'Neill presents a fortran 77 implementation of the simplex
1836 method. A factorial test is used to check if the computed optimum point is
1837 a local minimum. If the -restartdetection option is set to "oneill", that
1838 factorial test is used to see if a restart should be performed.
1839 O'Neill's factorial test requires <literal>2n</literal> function evaluations, where <literal>n</literal>
1840 is the number of variables.
1844 <title>User-defined algorithm</title>
1846 The <literal>-mymethod</literal> option allows to configure a user-defined
1847 simplex-based algorithm. The reason for this option is that many simplex-based
1848 variants of Nelder-Mead's algorithm have been developed over the years,
1849 with specific goals. While it is not possible to provide them all, it is very convenient
1850 to use the current structure without being forced to make many developments.
1853 The value of the <literal>-mymethod</literal> option is expected to be
1854 a Scilab function with the following header
1857 this = myalgorithm ( this )
1860 where <literal>this</literal> is the current object.
1863 In order to use the user-defined algorithm, the <literal>-method</literal> option must
1864 be set to "mine". In this case, the component performs the optimization
1865 exactly as if the user-defined algorithm was provided by the component.
1868 The user interested in that feature may use the internal scripts provided in the
1869 distribution as templates and tune his own algorithm from that point.
1870 There is of course no warranty that the
1871 user-defined algorithm improves on the standard algorithm, so that users
1872 use this feature at their own risks.
1876 <title>Example #1: basic use</title>
1878 In the following example, we solve a simple quadratic test case. We
1879 begin by defining the cost function, which takes 2 input arguments and
1880 returns the objective. The classical starting point [-1.2 1.0] is used.
1881 The neldermead_new creates a new neldermead object. Then we use the
1882 neldermead_configure method to configure the parameters of the problem. We
1883 use all default settings and perform the search for the optimum. The
1884 neldermead_display function is used to display the state of the
1885 optimization and the neldermead_get is used to retrieve the optimum
1888 <programlisting role="example"><![CDATA[
1889 function [ f , index ] = quadratic ( x , index )
1890 f = x(1)^2 + x(2)^2;
1893 nm = neldermead_new ();
1894 nm = neldermead_configure(nm,"-numberofvariables",2);
1895 nm = neldermead_configure(nm,"-function",quadratic);
1896 nm = neldermead_configure(nm,"-x0",x0);
1897 nm = neldermead_search(nm);
1899 xopt = neldermead_get(nm,"-xopt");
1900 nm = neldermead_destroy(nm);
1901 ]]></programlisting>
1904 <title>Example #2: customized use</title>
1906 In the following example, we solve the Rosenbrock test case. We
1907 begin by defining the Rosenbrock function, which takes 2 input arguments
1908 and returns the objective. The classical starting point [-1.2 1.0] is
1909 used. The neldermead_new creates a new neldermead object. Then we use the
1910 neldermead_configure method to configure the parameters of the problem.
1911 The initial simplex is computed from the axes and the single length 1.0
1912 (this is the default, but is explicitly written here as an example). The
1913 variable simplex algorithm by Nelder and Mead is used, which corresponds
1914 to the -method "variable" option. The neldermead_search function performs
1915 the search for the minimum. Once the minimum is found, the
1916 neldermead_contour allows to compute the data required by the contour
1917 function. This is possible since our problem involves only 2 parameters.
1918 This function uses the cost function previously configured to compute the
1919 required data. The contour plot is directly drawn from the data provided
1920 by neldermead_contour. Then we plot the initial guess on the contour plot
1921 as a blue dot. The neldermead_get function is used to get the optimum,
1922 which is associated with the -xopt option. The optimum is plot on the
1923 contour plot as a red dot.
1925 <programlisting role="example"><![CDATA[
1926 function [ f , index ] = rosenbrock ( x , index )
1927 f = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;
1930 nm = neldermead_new ();
1931 nm = neldermead_configure(nm,"-numberofvariables",2);
1932 nm = neldermead_configure(nm,"-function",rosenbrock);
1933 nm = neldermead_configure(nm,"-x0",x0);
1934 nm = neldermead_configure(nm,"-maxiter",200);
1935 nm = neldermead_configure(nm,"-maxfunevals",300);
1936 nm = neldermead_configure(nm,"-tolfunrelative",10*%eps);
1937 nm = neldermead_configure(nm,"-tolxrelative",10*%eps);
1938 nm = neldermead_configure(nm,"-simplex0method","axes");
1939 nm = neldermead_configure(nm,"-simplex0length",1.0);
1940 nm = neldermead_configure(nm,"-method","variable");
1941 nm = neldermead_configure(nm,"-verbose",1);
1942 nm = neldermead_configure(nm,"-verbosetermination",1);
1943 nm = neldermead_search(nm);
1944 xopt = neldermead_get(nm,"-xopt")
1945 nm = neldermead_destroy(nm);
1946 // Contouring the function.
1947 function f = rosenbrockC ( x1 , x2 )
1949 [ f , index ] = rosenbrock ( [x1,x2]' , index )
1951 xdata = linspace ( -2 , 2 , 100 );
1952 ydata = linspace ( -1 , 2 , 100 );
1954 contour ( xdata , ydata , rosenbrockC , [2 10 100 500 1000 2000] )
1955 // Plot starting point: x0 : blue dot
1956 plot(x0(1),x0(2),"bo");
1958 plot(xopt(1),xopt(2),"r*");
1959 ]]></programlisting>
1961 The -verbose option allows to get detailed information about the
1962 current optimization process. The following is a sample output for an
1963 optimization based on the Nelder and Mead variable-shape simplex
1964 algorithm. Only the output corresponding to the iteration #156 is
1965 displayed. In order to display specific outputs (or to create specific
1966 output files and graphics), the -outputcommand option should be
1971 Iteration #156 (total = 156)
1978 Optim Simplex Object:
1979 =====================
1985 > iterations=156 >= maxiter=200
1986 > funevals=299 >= maxfunevals=300
1987 > e(x)=8.798D-15 < 2.220D-15 * 1.4142136 + 0
1988 > Terminate = F, status = continue
1989 > simplex size=2.549D-13 < 0 + 2.220D-16 * 1
1990 > Terminate = F, status = continue
1993 Function Evaluation #300, index=2, x= [1 1]
1994 xr=[1 1], f(xr)=0.000000
1996 Function Evaluation #301, index=2, x= [1 1]
1997 xc=1 1, f(xc)=0.000000
1998 > Perform Inside Contraction
2004 <title>Example #3: use output function</title>
2006 There are several ways to print intermediate messages or plots during the
2007 optimization process.
2008 The first method is to set the "-verbose" option to 1,
2009 which prints a lot of detailed information.
2010 The other method is to use the <literal>"-outputcommand"</literal> option.
2011 We can either set it to the <literal>neldermead_defaultoutput</literal> or
2012 define our own function.
2013 In this section, we present the methods based on the <literal>"-outputcommand"</literal> option.
2016 In the following example, we use the <literal>"-outputcommand"</literal>option and
2017 set it to the <literal>neldermead_defaultoutput</literal>
2018 default output function.
2019 This function prints one line by iteration, with the main optimization information.
2021 <programlisting role="example"><![CDATA[
2022 function [ f , index ] = quadratic ( x , index )
2023 f = x(1)^2 + x(2)^2;
2026 nm = neldermead_new ();
2027 nm = neldermead_configure(nm,"-numberofvariables",2);
2028 nm = neldermead_configure(nm,"-function",quadratic);
2029 nm = neldermead_configure(nm,"-x0",x0);
2030 nm = neldermead_configure(nm,"-outputcommand",neldermead_defaultoutput);
2031 nm = neldermead_search(nm);
2032 nm = neldermead_destroy(nm);
2033 ]]></programlisting>
2035 The previous script produces the following output.
2039 Iter. #0, Feval #5, Fval = 2 -- init
2040 Iter. #1, Feval #5, Fval = 2 -- init
2041 Iter. #2, Feval #6, Fval = 2 -- reflection
2042 Iter. #3, Feval #8, Fval = 0.5 -- expansion
2043 Iter. #4, Feval #9, Fval = 0.5 -- reflection
2045 Iter. #48, Feval #92, Fval = 8.557D-13 -- reflection
2046 Iter. #49, Feval #94, Fval = 7.893D-13 -- insidecontraction
2047 Iter. #50, Feval #96, Fval = 1.601D-13 -- insidecontraction
2048 Iter. #51, Feval #98, Fval = 1.291D-13 -- insidecontraction
2049 Iter. #52, Feval #100, Fval = 3.139D-14 -- outsidecontraction
2050 =================================
2052 Iter. #52, Feval #100, Fval = 3.139D-14 -- done
2055 In the following example, we define our own output function "myoutputcmd", which takes the current state
2056 as its first argument. The state is a string which can contain "init",
2057 "iter" or "done", depending on the status of the optimization. The data
2058 input argument is a tlist, which contains the data associated with the
2059 current iteration. In this case, we use the fields to print a message in
2062 <programlisting role="example"><![CDATA[
2063 function [ f , index ] = rosenbrock ( x , index )
2064 f = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
2067 function stop = myoutputcmd ( state , data )
2068 iter = data.iteration
2069 if ( state == "init" ) then
2070 mprintf ( "=================================\n");
2071 mprintf ( "Initialization\n");
2072 elseif ( state == "done" ) then
2073 mprintf ( "=================================\n");
2074 mprintf ( "End of Optimization\n");
2079 simplex = data.simplex
2080 // Simplex is a data structure, which can be managed
2081 // by the optimsimplex class.
2082 ssize = optimsimplex_size ( simplex )
2083 mprintf ( "Iter. #%3d, Feval #%3d, Fval = %.1e, x = %s, S = %.1e\n", ..
2084 iter, fc, fval, strcat(string(x)," "), ssize);
2088 nm = neldermead_new ();
2089 nm = neldermead_configure(nm,"-numberofvariables",2);
2090 nm = neldermead_configure(nm,"-function",rosenbrock);
2091 nm = neldermead_configure(nm,"-x0",[-1.2 1.0]');
2092 nm = neldermead_configure(nm,"-maxiter",200);
2093 nm = neldermead_configure(nm,"-maxfunevals",300);
2094 nm = neldermead_configure(nm,"-tolfunrelative",10*%eps);
2095 nm = neldermead_configure(nm,"-tolxrelative",10*%eps);
2096 nm = neldermead_configure(nm,"-simplex0method","axes");
2097 nm = neldermead_configure(nm,"-simplex0length",1.0);
2098 nm = neldermead_configure(nm,"-method","variable");
2099 nm = neldermead_configure(nm,"-verbose",0);
2100 nm = neldermead_configure(nm,"-verbosetermination",0);
2101 nm = neldermead_configure(nm,"-outputcommand",myoutputcmd);
2102 nm = neldermead_search(nm);
2103 nm = neldermead_destroy(nm);
2104 ]]></programlisting>
2106 The previous script produces the following output.
2109 =================================
2111 Iter. # 0, Feval # 5, Fval = 2.4e+001, x = -1.2 1, S = 1.0e+000
2112 Iter. # 1, Feval # 5, Fval = 2.4e+001, x = -1.2 1, S = 1.0e+000
2113 Iter. # 2, Feval # 7, Fval = 2.4e+001, x = -1.2 1, S = 1.0e+000
2114 Iter. # 3, Feval # 9, Fval = 2.4e+001, x = -1.2 1, S = 1.0e+000
2115 Iter. # 4, Feval # 11, Fval = 1.0e+001, x = -1.0125 0.78125, S = 6.0e-001
2116 Iter. # 5, Feval # 13, Fval = 4.7e+000, x = -1.028125 1.1328125, S = 3.5e-001
2118 Iter. #155, Feval #297, Fval = 2.0e-026, x = 1 1, S = 4.6e-013
2119 Iter. #156, Feval #299, Fval = 6.9e-027, x = 1 1, S = 2.5e-013
2120 Iter. #157, Feval #301, Fval = 6.0e-027, x = 1 1, S = 2.8e-013
2121 =================================
2123 Iter. #157, Feval #301, Fval = 6.0e-027, x = 1 1, S = 2.8e-013
2126 As another example of use, we could format the message so
2127 that it uses LaTeX formatting rules, which may allow the user to directly
2128 copy and paste the output into a LaTeX report.
2132 <title>Example #4: Optimization with bounds</title>
2134 The <literal>neldermead</literal> solver can optimize problems with
2136 To do this, we can use Box's algorithm, which projects the simplex
2137 into the bounds during the optimization.
2138 In this case, the initial guess must be located within the bounds.
2141 In the following example, we find the minimum of a quadratic function
2142 within given bounds.
2143 In order to compute the initial simplex, we use randomized bounds, that is,
2144 we compute k random vertices uniformly distributed within the bounds.
2145 The default value is so that the number of points is twice the number of
2146 variables of the problem.
2147 In this particular case, we have n=2 variables and k=4 vertices.
2149 <programlisting role="example"><![CDATA[
2150 function [ f , index ] = myquad ( x , index )
2155 nm = neldermead_new ();
2156 nm = neldermead_configure(nm,"-numberofvariables",2);
2157 nm = neldermead_configure(nm,"-function",myquad);
2158 nm = neldermead_configure(nm,"-x0",x0);
2159 nm = neldermead_configure(nm,"-method","box");
2160 nm = neldermead_configure(nm,"-boundsmin",[1 1]);
2161 nm = neldermead_configure(nm,"-boundsmax",[2 2]);
2162 nm = neldermead_configure(nm,"-simplex0method","randbounds");
2163 nm = neldermead_search(nm);
2164 xopt = neldermead_get(nm,"-xopt") // Should be [1 1]
2165 fopt = neldermead_get(nm,"-fopt") // Should be 2
2166 nm = neldermead_destroy(nm);
2167 ]]></programlisting>
2170 <title>Example #5: Optimization with nonlinear constraints</title>
2172 The <literal>neldermead</literal> solver can optimize problems with
2173 nonlinear constraints.
2174 In the following example, we solve Rosenbrock's Post Office problem, which
2175 has both bounds and linear constraints.
2176 In our example, we will manage the linear constraints as general non-linear
2177 constraints (i.e. the solver does not make a difference if the constraints are
2178 linear or non-linear).
2179 This example was first
2180 presented in "An automatic method for finding the greatest or least value of a function",
2182 This example was first used with the complex method of Box in
2183 "Algorithm 454: The complex method for constrained optimization" by
2184 Richardson, Kuester, 1971.
2185 Richardson and Kuester found the minimum function value
2186 F=-3456, with X1 = 24.01, X2 = 12.00, X3 = 12.00 and
2187 72 Iterations were necessary for them to get this result.
2190 In the following function, we define the function <literal>fpostoffice</literal>,
2191 which returns both the objective function <literal>f</literal> and the
2192 constraint value <literal>c</literal>.
2193 The original constraint is the "double" inequality constraint
2194 <literal>0<=x(1) + 2 * x(2) + 2 * x(3) <=72</literal>.
2195 To take this constraint into account, we turn it into two separate, positive,
2196 constraints and set <literal>c</literal> as a 1-by-2 matrix of doubles.
2198 <programlisting role="example"><![CDATA[
2199 function [ f , c , index ] = fpostoffice ( x , index )
2202 if ( index==2 | index==6 ) then
2203 f = -x(1) * x(2) * x(3)
2206 if ( index==5 | index==6 ) then
2207 c1 = x(1) + 2 * x(2) + 2 * x(3)
2212 ]]></programlisting>
2214 In the following script, we solve Rosenbrock's Post Office problem.
2215 First, we initialize the random number generator, so that the results are always the
2217 Then, we check that the cost function is correctly defined and that the
2218 constraints are satisfied at the initial guess.
2219 Then we configure the algorithm so that Box's algorithm is used and
2220 setup the bounds of the problem.
2221 We configure the parameters of the algorithm as suggested by Box.
2223 <programlisting role="example"><![CDATA[
2225 x0 = [1.0 1.0 1.0].';
2226 // Compute f(x0) : should be close to -1
2227 fx0 = fpostoffice ( x0 , 2 )
2228 // Compute the constraints: cx0 should be [5 67]
2229 [ fx0 , cx0, index ] = fpostoffice ( x0 , 6 )
2230 // Compute f(xopt) : fopt should be -3456
2231 xopt = [24 12 12].';
2232 fopt = fpostoffice ( xopt );
2233 // Setup optimization
2234 nm = neldermead_new ();
2235 nm = neldermead_configure(nm,"-numberofvariables",3);
2236 nm = neldermead_configure(nm,"-function",fpostoffice);
2237 nm = neldermead_configure(nm,"-x0",x0);
2238 nm = neldermead_configure(nm,"-maxiter",300);
2239 nm = neldermead_configure(nm,"-maxfunevals",300);
2240 nm = neldermead_configure(nm,"-method","box");
2241 nm = neldermead_configure(nm,"-boundsmin",[0.0 0.0 0.0]);
2242 nm = neldermead_configure(nm,"-boundsmax",[42.0 42.0 42.0]);
2243 // Configure like Box
2244 nm = neldermead_configure(nm,"-simplex0method","randbounds");
2245 nm = neldermead_configure(nm,"-nbineqconst",2);
2246 nm = neldermead_configure(nm,"-tolxmethod" , %f );
2247 nm = neldermead_configure(nm,"-tolsimplexizemethod",%f);
2248 nm = neldermead_configure(nm,"-boxtermination" , %t );
2249 nm = neldermead_configure(nm,"-boxtolf" , 0.001 );
2250 nm = neldermead_configure(nm,"-boxboundsalpha" , 0.0001 );
2252 // Check that the cost function is correctly connected.
2253 [ nm , result ] = neldermead_function ( nm , x0 );
2255 // Perform optimization
2256 nm = neldermead_search(nm);
2257 xcomp = neldermead_get(nm,"-xopt")
2258 // Compare with the exact optimum:
2260 fcomp = neldermead_get(nm,"-fopt")
2261 // Compare with the exact function value:
2263 nm = neldermead_destroy(nm);
2264 ]]></programlisting>
2266 In general, we should not expect too much from this algorithm with
2267 nonlinear constraints.
2268 Indeed, some cases require thousands of iterations to converge to
2269 an optimum, because the nonlinear constraints leave a too small
2270 space for the simplex to evolve.
2274 <title>Example #6: Passing extra parameters</title>
2276 In the following example, we solve a simple quadratic test case.
2277 Notice that the objective function has two extra parameters
2278 <literal>a</literal> and <literal>b</literal>.
2279 This is why the "-function" option is set as a list,
2280 where the first element is the function and the
2281 remaining elements are the extra parameters.
2283 <programlisting role="example"><![CDATA[
2284 function [ f , index ] = quadratic_ab ( x , index , a , b )
2285 f = a * x(1)^2 + b * x(2)^2;
2288 nm = neldermead_new ();
2289 nm = neldermead_configure(nm,"-numberofvariables",2);
2292 nm = neldermead_configure(nm,"-function",list(quadratic_ab,a,b));
2293 nm = neldermead_configure(nm,"-x0",x0);
2294 nm = neldermead_search(nm);
2295 xopt = neldermead_get(nm,"-xopt")
2296 nm = neldermead_destroy(nm);
2297 ]]></programlisting>
2300 <title>Example #7: Restarting without bounds</title>
2302 In the following example, we reproduce the experiment
2303 published by Ken McKinnon in 1998.
2304 For this particular function and this particular initial simplex,
2305 the Nelder-Mead algorithm converges to a nonstationnary point.
2308 We first define the objective function, the initial simplex
2309 and the expected solution of this unconstrained optimization problem.
2311 <programlisting role="example"><![CDATA[
2312 function [ f , index ] = mckinnon ( x , index )
2317 f = theta*phi*abs(x(1))^tau+x(2)*(1+x(2))
2319 f = theta*x(1)^tau+x(2)*(1+x(2))
2323 // The initial simplex
2324 lambda1 = (1.0 + sqrt(33))/8;
2325 lambda2 = (1.0 - sqrt(33))/8;
2332 // The expected solution
2335 ]]></programlisting>
2337 Then we run the algorithm two times in sequence.
2338 At the end of the first optimization process, the algorithm has
2339 converged to the point [0,0] which is nonstationnary.
2340 This is why we restart the algorithm and get the correct minimum.
2342 <programlisting role="example"><![CDATA[
2343 nm = neldermead_new ();
2344 nm = neldermead_configure(nm,"-numberofvariables",2);
2345 nm = neldermead_configure(nm,"-function",mckinnon);
2346 nm = neldermead_configure(nm,"-x0",[1.0 1.0]');
2347 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
2348 nm = neldermead_configure(nm, "-maxiter",200);
2349 nm = neldermead_configure(nm, "-maxfunevals",500);
2350 nm = neldermead_configure(nm,"-simplex0method","given");
2351 nm = neldermead_configure(nm,"-coords0",coords0);
2352 nm = neldermead_configure(nm,"-method","variable");
2354 nm = neldermead_search(nm);
2355 xopt = neldermead_get(nm,"-xopt")
2356 fopt = neldermead_get(nm,"-fopt")
2357 iterations = neldermead_get(nm,"-iterations")
2358 status = neldermead_get(nm,"-status")
2359 // Search #2: succeeds
2360 nm = neldermead_restart ( nm );
2361 xopt = neldermead_get(nm,"-xopt")
2362 fopt = neldermead_get(nm,"-fopt")
2363 iterations = neldermead_get(nm,"-iterations")
2364 status = neldermead_get(nm,"-status")
2365 nm = neldermead_destroy(nm);
2366 ]]></programlisting>
2368 We can also use the automatic stagnation detection method
2369 created by Kelley, so that the algorithm automatically restart
2370 the algorithm when needed.
2372 <programlisting role="example"><![CDATA[
2373 nm = neldermead_new ();
2374 nm = neldermead_configure(nm,"-numberofvariables",2);
2375 nm = neldermead_configure(nm,"-function",mckinnon);
2376 nm = neldermead_configure(nm,"-x0",[1.0 1.0]');
2377 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
2378 nm = neldermead_configure(nm, "-maxiter",200);
2379 nm = neldermead_configure(nm, "-maxfunevals",500);
2380 nm = neldermead_configure(nm,"-simplex0method","given");
2381 nm = neldermead_configure(nm,"-coords0",coords0);
2382 nm = neldermead_configure(nm,"-method","variable");
2383 nm = neldermead_configure(nm,"-kelleystagnationflag",%t);
2384 nm = neldermead_configure(nm,"-restartflag",%t);
2385 nm = neldermead_configure(nm,"-restartdetection","kelley");
2386 nm = neldermead_search(nm);
2387 xopt = neldermead_get(nm,"-xopt")
2388 fopt = neldermead_get(nm,"-fopt")
2389 iterations = neldermead_get(nm,"-iterations")
2390 restartnb = neldermead_get ( nm , "-restartnb" )
2391 status = neldermead_get(nm,"-status")
2392 nm = neldermead_destroy(nm);
2393 ]]></programlisting>
2395 See the demonstrations to get a graphical plot of the
2396 intermediate simplices in Mc Kinnon's experiment.
2400 <title>Example #8: Restarting with bounds</title>
2402 In the following experimeant, we solve an optimization problem
2404 We use Box's algorithm, which is the only algorithm which manages bounds.
2405 We use the randomized bounds simplex both for the initial simplex and
2406 for the restart simplex.
2408 <programlisting role="example"><![CDATA[
2409 function [ f , index ] = myquad ( x , index )
2410 f = x(1)^2 + x(2)^2 + x(3)^2
2412 x0 = [1.2 1.9,1.5].';
2417 nm = neldermead_new ();
2418 nm = neldermead_configure(nm,"-numberofvariables",3);
2419 nm = neldermead_configure(nm,"-function",myquad);
2420 nm = neldermead_configure(nm,"-x0",x0);
2421 nm = neldermead_configure(nm,"-method","box");
2422 nm = neldermead_configure(nm,"-boundsmin",[1 1 1]);
2423 nm = neldermead_configure(nm,"-boundsmax",[2 2 2]);
2424 nm = neldermead_configure(nm,"-simplex0method","randbounds");
2425 nm = neldermead_search(nm);
2426 nm = neldermead_configure(nm,"-maxiter",200);
2427 nm = neldermead_configure(nm,"-maxfunevals",200);
2428 nm = neldermead_configure(nm,"-restartsimplexmethod","randbounds");
2429 nm = neldermead_restart(nm);
2430 xopt = neldermead_get(nm,"-xopt")
2431 fopt = neldermead_get(nm,"-fopt")
2432 status = neldermead_get(nm,"-status")
2433 nm = neldermead_destroy(nm);
2434 ]]></programlisting>
2437 <title>Changes in Scilab 5.4</title>
2439 Many changes have been done in Scilab 5.4, which simplify the use of the
2440 neldermead component.
2443 Tagged -costfargument option of optimbase as obsolete: will be
2444 maintained for backward compatibility until 5.4.1.
2445 The -fun option can now be a list, where the element #1 is a
2446 function, and the elements #2 to the end are automatically appended to
2447 the calling sequence.
2448 To update your code, replace:
2451 nm = neldermead_configure(nm,"-function",myfun);
2452 nm = neldermead_configure(nm,"-costfargument",mystuff);
2458 nm = neldermead_configure(nm,"-function",list(myfun,mystuff));
2461 Tagged -outputcommandarg option of optimbase as obsolete: will be
2462 maintained for backward compatibility until 5.4.1.
2463 The -outputcommand option can now be a list, where the element #1 is
2464 a function, and the elements #2 to the end are automatically appended
2465 to the calling sequence.
2466 To update your code, replace:
2469 nm = neldermead_configure(nm,"-outputcommand",myoutputfun);
2470 nm = neldermead_configure(nm,"-outputcommandarg",mystuff);
2476 nm = neldermead_configure(nm,"-outputcommand",list(myoutputfun,mystuff));
2479 Tagged "outputfun(x,optimValues,state)" calling sequence of fminsearch
2480 as obsolete: will be maintained for backward compatibility until
2482 The new calling sequence is "stop=outputfun(x,optimValues,state)"
2483 To update your code, replace:
2486 function outfun ( x , optimValues , state )
2494 function stop = outfun ( x , optimValues , state )
2500 Tagged "myoutputfun(state,data)" calling sequence of neldermead
2501 as obsolete: will be maintained for backward compatibility until
2503 The new calling sequence is "stop=myoutputfun(state,data)"
2504 To update your code, replace:
2507 function myoutputfun ( state , data )
2515 function stop = myoutputfun ( state , data )
2521 Tagged "-myterminateflag" and "-myterminate" options as obsolete:
2522 will be maintained for backward compatibility until 5.4.1.
2523 To update your code, replace:
2526 function [ this , terminate , status ] = myoldterminate ( this , simplex )
2527 ssize = optimsimplex_size ( simplex , "sigmaplus" );
2528 if ( ssize < 1.e-2 ) then
2540 function stop = myoutputcmd ( state , data )
2541 simplex = data.simplex
2542 ssize = optimsimplex_size ( simplex , "sigmaplus" );
2543 if ( ssize < 1.e-2 ) then
2551 and replace the configuration:
2554 nm = neldermead_configure(nm,"-myterminateflag",%t);
2555 nm = neldermead_configure(nm,"-myterminate",myoldterminate);
2561 nm = neldermead_configure(nm,"-outputcommand",myoutputcmd);
2564 Tagged "-tolvarianceflag", "-tolabsolutevariance", and
2565 "-tolrelativevariance" options as obsolete:
2566 will be maintained for backward compatibility until 5.4.1.
2567 To update your code, create an output function:
2570 function stop = myoutputcmd ( state, data, tolrelativevariance, tolabsolutevariance, variancesimplex0 )
2571 simplex = data.simplex
2573 if ( state == "iter") then
2574 var = optimsimplex_fvvariance ( simplex )
2575 if ( var < tolrelativevariance * variancesimplex0 + tolabsolutevariance ) then
2582 Create the initial simplex and compute the variance
2583 of the function values:
2587 simplex0 = optimsimplex_new ( "axes" , x0.' );
2588 coords0 = optimsimplex_getallx(simplex0);
2589 variancesimplex0 = optimsimplex_fvvariance ( simplex0 );
2592 Finally, replace the configuration:
2595 nm = neldermead_configure(nm,"-tolvarianceflag",%t);
2596 nm = neldermead_configure(nm,"-tolabsolutevariance",1.e-4);
2597 nm = neldermead_configure(nm,"-tolrelativevariance",1.e-4);
2603 tolabsolutevariance = 1.e-4;
2604 tolrelativevariance = 1.e-4;
2605 stopfun = list(myoutputcmd, tolrelativevariance, tolabsolutevariance, variancesimplex0);
2606 nm = neldermead_configure(nm,"-outputcommand",stopfun);
2610 <title>Spendley et al. implementation notes</title>
2612 The original paper may be implemented with several variations, which
2613 might lead to different results. This section defines what algorithmic
2614 choices have been used.
2616 <para>The paper states the following rules.</para>
2620 "Rule 1. Ascertain the lowest reading y, of yi ... yk+1 Complete
2621 a new simplex Sp by excluding the point Vp corresponding to y, and
2622 replacing it by V* defined as above."
2627 "Rule 2. If a result has occurred in (k + 1) successive
2628 simplexes, and is not then eliminated by application of Rule 1, do not
2629 move in the direction indicated by Rule 1, or at all, but discard the
2630 result and replace it by a new observation at the same point."
2635 "Rule 3. If y is the lowest reading in So , and if the next
2636 observation made, y* , is the lowest reading in the new simplex S , do
2637 not apply Rule 1 and return to So from Sp . Move out of S, by
2638 rejecting the second lowest reading (which is also the second lowest
2644 We implement the following "rules" of the Spendley et al.
2650 Rule 1 is strictly applied, but the reflection is done by
2651 reflection the high point, since we minimize a function instead of
2652 maximizing it, like Spendley.
2657 Rule 2 is NOT implemented, as we expect that the function
2658 evaluation is not subject to errors.
2663 Rule 3 is applied, ie reflection with respect to next to high
2669 The original paper does not mention any shrink step. When the
2670 original algorithm cannot improve the function value with reflection
2671 steps, the basic algorithm stops. In order to make the current
2672 implementation of practical value, a shrink step is included, with
2673 shrinkage factor sigma. This perfectly fits into to the spirit of the
2674 original paper. Notice that the shrink step make the rule #3 (reflection
2675 with respect to next-to-worst vertex) unnecessary. Indeed, the minimum
2676 required steps are the reflection and shrinkage. Never the less, the rule
2677 #3 has been kept in order to make the algorithm as close as it can be to
2682 <title>Nelder-Mead implementation notes</title>
2684 The purpose of this section is to analyse the current implementation of Nelder-Mead's algorithm.
2687 The algorithm that we use is described in "Iterative Methods for Optimization" by C. T. Kelley.
2690 The original paper uses a "greedy" expansion, in which the expansion point
2691 is accepted whatever its function value. The current implementation,
2692 as most implementations, uses the expansion point only if it improves
2693 over the reflection point, that is,
2698 if fe<fr, then the expansion point is accepted,
2703 if not, the reflection point is accepted.
2708 The termination criteria suggested by Nelder and Mead is based on an
2709 absolute tolerance on the standard deviation of the function values in the simplex.
2710 We provide this original termination criteria with the <literal>-tolvarianceflag</literal>
2711 option, which is disabled by default.
2715 <title>Box's complex algorithm implementation notes</title>
2717 In this section, we analyse the current implementation of Box's complex method.
2720 The initial simplex can be computed as in Box's paper, but this may not
2721 be safe. In his paper, Box suggest that if a vertex of the initial simplex
2722 does not satisfy the non linear constraints, then it should be "moved halfway
2723 toward the centroid of those points already selected". This behaviour
2724 is available when the <literal>-scalingsimplex0</literal> option is set to
2725 <literal>"tocenter"</literal>. It may happen, as suggested by Guin, that
2726 the centroid is not feasible. This may happen if the constraints are not
2727 convex. In this case, the initial simplex cannot be computed. This is why
2728 we provide the <literal>"tox0"</literal> option, which allows to compute the
2729 initial simplex by scaling toward the initial guess, which is always feasible.
2732 In Box's paper, the scaling into the non linear constraints is performed
2733 "toward" the centroid, that is, by using a scaling factor equal to 0.5.
2734 This default scaling factor might be sub-optimal in certain situations.
2735 This is why we provide the <literal>-boxineqscaling</literal> option,
2736 which allows to configure the scaling factor.
2739 In Box's paper, whether we are concerned with the initial simplex or with the
2740 simplex at a given iteration, the scaling for the non linear constraints is performed
2741 without end. This is because Box's hypothesis is that "ultimately, a satisfactory
2742 point will be found". As suggested by Guin, if the process fails, the algorithm
2743 goes into an infinite loop. In order to avoid this, we perform the scaling until
2744 a minimum scaling value is reached, as defined by the <literal>-guinalphamin</literal>
2748 We have taken into account for the comments by Guin, but it should be emphasized
2749 that the current implementation is still as close as possible to Box's
2750 algorithm and is not Guin's algorithm. More precisely, during the iterations,
2751 the scaling for the non linear constraints is still performed toward the centroid,
2752 be it feasible or not.
2756 <title>Bibliography</title>
2758 "Sequential Application of Simplex Designs in Optimisation and
2759 Evolutionary Operation", Spendley, W. and Hext, G. R. and Himsworth, F.
2760 R., American Statistical Association and American Society for Quality,
2764 "A Simplex Method for Function Minimization", Nelder, J. A. and
2765 Mead, R., The Computer Journal, 1965
2768 "A New Method of Constrained Optimization and a Comparison With
2769 Other Methods", M. J. Box, The Computer Journal 1965 8(1):42-52, 1965 by
2770 British Computer Society
2773 "Discussion and correspondence: modification of the complex method
2774 of constrained optimization", J. A. Guin, The Computer Journal,
2778 "Detection and Remediation of Stagnation in the Nelder--Mead
2779 Algorithm Using a Sufficient Decrease Condition", Kelley C. T., SIAM J. on
2783 "Iterative Methods for Optimization", C. T. Kelley, SIAM Frontiers
2784 in Applied Mathematics, 1999
2787 "Algorithm AS47 - Function minimization using a simplex procedure",
2788 O'Neill, R., Applied Statistics, 1971
2791 "Nelder Mead's User Manual", Consortium Scilab - Digiteo, Michael
2795 Ken McKinnon, Convergence of the Nelder-Mead simplex method to a nonstationary point,
2796 SIAM Journal on Optimization, Volume 9, Number 1, 1998, pages 148-158.
2799 <refsection role="see also">
2800 <title>See Also</title>
2801 <simplelist type="inline">
2803 <link linkend="optimbase_new">optimbase_new</link>
2806 <link linkend="optimsimplex_new">optimsimplex_new</link>
2809 <link linkend="optimbase_overview">optimbase</link>
2812 <link linkend="optimsimplex_overview">optimsimplex</link>
2815 <link linkend="nmplot">nmplot</link>
projects / scilab.git / blob