fminsearch: added (please wait) messages in the demos
[scilab.git] / scilab / modules / optimization / demos / neldermead / neldermead_boxpost.sce
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2008-2009 - INRIA - Michael Baudin
3 //
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution.  The terms
7 // are also available at
8 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9
10 //
11 // nmplot_boxpost.sce --
12 //   Show that the Box algorithm is able to reproduce the 
13 //   numerical experiment presented in Richardson and Kuester's paper.
14 //   Rosenbrock's Post Office
15 //
16
17 mprintf("Illustrates Box'' algorithm on Box problem A.\n");
18 mprintf("Defining Box Problem A function...\n");
19
20 //
21 //  Reference:
22 //
23 //   M.J.Box, 
24 //   "A new method of constrained optimization 
25 //   and a comparison with other methods".
26 //   The Computer Journal, Volume 8, Number 1, 1965, 42--52
27 //   Problem A
28 //
29 //   Algorithm 454: the complex method for constrained optimization [E4]
30 //   Communications of the ACM, Volume 16 ,  Issue 8  (August 1973)
31 //   Pages: 487 - 489   
32 //
33 //   Richardson and Kuester Results :
34 //   F=1.0000
35 //   X1 = 3.0000
36 //   X2 = 1.7320
37 //   Iterations : 68
38 //
39 //
40 // Note 
41 //    The maximum bound for the parameter x1
42 //    is not indicated in Box's paper, but indicated in "Algo 454".
43 //    The maximum bound for x2 is set to 100/sqrt(3) and satisfies the constraint on x2.
44 //    The original problem was to maximize the cost function.
45 //    Here, the problem is to minimize the cost function.
46
47 //
48 // boxproblemB --
49 //   Computes the Box problem B cost function and 
50 //   inequality constraints.
51 //
52 // Arguments
53 //   x: the point where to compute the function
54 //   index : the stuff to compute
55 //   data : the parameters of Box cost function
56 // Note
57 //  The following protocol is used
58 //  * if index=1, or no index, returns the value of the cost 
59 //    function (default case)
60 //  * if index=2, returns the value of the nonlinear inequality 
61 //    constraints, as a row array
62 //  * if index=3, returns an array which contains
63 //    at index #0, the value of the cost function  
64 //    at index #1 to the end is the list of the values of the nonlinear 
65 //    constraints
66 //  The inequality constraints are expected to be positive.
67 //
68 function [ f , c , index ] = boxproblemB ( x , index )
69   if (~isdef('index','local')) then
70     index = 1
71   end
72   x3 = x(1) + sqrt(3.0) * x(2)
73   if ( index==1 | index==3 ) then
74     f = -(9.0 - (x(1) - 3.0) ^ 2) * x(2) ^ 3 / 27.0 / sqrt(3.0)
75   end
76   if ( index==2 | index==3 ) then
77       c1 = x(1) / sqrt(3.0) - x(2)
78       c2 = x3
79       c3 = 6.0 - x3
80       c = [c1 c2 c3]
81   end
82 endfunction
83
84
85 //
86 // Initialize the random number generator, so that the results are always the
87 // same.
88 //
89 rand("seed" , 0)
90
91 x0 = [1.0 0.5].';
92 // Compute f(x0) : should be close to -0.0133645895646
93 [ fx0 , c, index ] = boxproblemB ( x0 , 1 );
94 mprintf("Computed fx0 = %e (expected = %e)\n",fx0 , -0.0133645895646 );
95 result = boxproblemB ( x0 , 2 );
96 mprintf("Computed Constraints(x0) = [%e %e %e]\n", ...
97   result(1), result(2), result(3) );
98 mprintf("Expected Constraints(x0) = [%e %e %e]\n", ...
99   0.0773503 , 1.8660254 , 4.1339746 );
100
101
102 xopt = [3.0 1.7320508075688774].';
103 // Compute f(xopt) : should be -1.0
104 fopt = boxproblemB ( xopt );
105 mprintf("Computed fopt = %e (expected = %e)\n", fopt , -1.0 );
106
107 nm = neldermead_new ();
108 nm = neldermead_configure(nm,"-numberofvariables",2);
109 nm = neldermead_configure(nm,"-function",boxproblemB);
110 nm = neldermead_configure(nm,"-x0",x0);
111 nm = neldermead_configure(nm,"-maxiter",300);
112 nm = neldermead_configure(nm,"-maxfunevals",300);
113 nm = neldermead_configure(nm,"-method","box");
114 nm = neldermead_configure(nm,"-verbose",1);
115 nm = neldermead_configure(nm,"-logfile" , "postoffice.txt" );
116 nm = neldermead_configure(nm,"-verbosetermination",1);
117 nm = neldermead_configure(nm,"-boundsmin",[0.0 0.0]);
118 nm = neldermead_configure(nm,"-boundsmax",[100.0 57.735026918962582]);
119 // Configure like Box
120 nm = neldermead_configure(nm,"-simplex0method","randbounds");
121 nm = neldermead_configure(nm,"-nbineqconst",3);
122 nm = neldermead_configure(nm,"-tolxmethod" , %f );
123 nm = neldermead_configure(nm,"-tolsimplexizemethod",%f);
124 nm = neldermead_configure(nm,"-boxtermination" , %t );
125 nm = neldermead_configure(nm,"-boxtolf" , 0.001 );
126 nm = neldermead_configure(nm,"-boxboundsalpha" , 0.0001 );
127
128 //
129 // Check that the cost function is correctly connected.
130 // The index must be provided, because the additionnal argument "data"
131 // comes after.
132 //
133 [ nm , result ] = neldermead_function ( nm , x0 , 1 );
134 //
135 // Perform optimization
136 //
137 nm = neldermead_search(nm);
138 neldermead_display(nm);
139 xcomp = neldermead_get(nm,"-xopt");
140 mprintf("x computed=%s\n",strcat(string(xcomp)," "));
141 mprintf("x expected=%s\n",strcat(string(xopt)," "));
142 shift = norm(xcomp-xopt)/norm(xopt);
143 mprintf("Shift =%f\n",shift);
144 fcomp = neldermead_get(nm,"-fopt");
145 mprintf("f computed=%f\n",fcomp);
146 mprintf("f expected=%f\n",fopt);
147 shift = abs(fcomp-fopt)/abs(fopt);
148 mprintf("Shift =%f\n",shift);
149 nm = neldermead_destroy(nm);
150 deletefile ( "postoffice.txt" )
151 mprintf("End of demo.\n");
152
153 //
154 // Load this script into the editor
155 //
156 filename = 'neldermead_boxpost.sce';
157 dname = get_absolute_file_path(filename);
158 editor ( dname + filename );
159