fminsearch: added (please wait) messages in the demos
[scilab.git] / scilab / modules / optimization / demos / neldermead / neldermead_boxproblemA.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_boxproblemA.sce --
12 //   Show that the Box algorithm is able to reproduce the 
13 //   numerical experiment presented in Box's paper.
14 //
15
16 mprintf("Illustrates Box'' algorithm on Box problem A.\n");
17 mprintf("Defining Box Problem A function...\n");
18
19 //
20 //  Reference:
21 //
22 //   M.J.Box, 
23 //   "A new method of constrained optimization 
24 //   and a comparison with other methods".
25 //   The Computer Journal, Volume 8, Number 1, 1965, 42--52
26 //   Problem A
27 //
28 //   Algorithm 454: the complex method for constrained optimization [E4]
29 //   Communications of the ACM, Volume 16 ,  Issue 8  (August 1973)
30 //   Pages: 487 - 489   
31 //
32
33 //
34 // boxproblemA --
35 //   Computes the Box problem A cost function and 
36 //   inequality constraints.
37 //
38 // Arguments
39 //   x: the point where to compute the function
40 //   index : the stuff to compute
41 //   data : the parameters of Box cost function
42 //
43 function [ f , c , index , data ] = boxproblemA ( x , index , data )
44   f = []
45   c = []
46   b = x(2) + 0.01 * x(3)
47   x6 = (data.k1 + data.k2 * x(2) ...
48             + data.k3 * x(3) + data.k4 * x(4) + data.k5 * x(5)) * x(1)
49   y1 = data.k6 + data.k7 * x(2) + data.k8 * x(3) ...
50             + data.k9 * x(4) + data.k10 * x(5)
51   y2 = data.k11 + data.k12 * x(2) + data.k13 * x(3) ...
52             + data.k14 * x(4) + data.k15 * x(5)
53   y3 = data.k16 + data.k17 * x(2) + data.k18 * x(3) ...
54             + data.k19 * x(4) + data.k20 * x(5)
55   y4 = data.k21 + data.k22 * x(2) + data.k23 * x(3) ...
56             + data.k24 * x(4) + data.k25 * x(5)
57   x7 = ( y1 + y2 + y3 ) * x(1)
58   x8 = (data.k26 + data.k27 * x(2) + data.k28 * x(3) ...
59             + data.k29 * x(4) ...
60             + data.k30 * x(5) ) * x(1) + x6 + x7
61   if ( index == 2 | index == 6 ) then
62     f = (data.a2 * y1 + data.a3 * y2 + data.a4 * y3 + data.a5 * y4 ...
63              + 7840 * data.a6 - 100000 * data.a0 ...
64              - 50800 * b * data.a7 + data.k31 + data.k32 * x(2) + data.k33 * x(3) ...
65              + data.k34 * x(4) + data.k35 * x(5)) * x(1) ...
66              - 24345 + data.a1 * x6
67     f = -f
68   end
69   if ( index == 5 | index == 6 ) then
70       c1 = x6
71       c2 = 294000 - x6
72       c3 = x7
73       c4 = 294000 - x7
74       c5 = x8
75       c6 = 277200 - x8
76       c = [c1 c2 c3 c4 c5 c6]
77   end
78 endfunction
79
80 boxparams = struct();
81 boxparams.a0 = 9;
82 boxparams.a1 = 15;
83 boxparams.a2 = 50;
84 boxparams.a3 = 9.583;
85 boxparams.a4 = 20;
86 boxparams.a5 = 15;
87 boxparams.a6 = 6;
88 boxparams.a7 = 0.75;
89 boxparams.k1 = -145421.402;
90 boxparams.k2 = 2931.1506;
91 boxparams.k3 = -40.427932;
92 boxparams.k4 = 5106.192;
93 boxparams.k5 = 15711.36;
94 boxparams.k6 = -161622.577;
95 boxparams.k7 = 4176.15328;
96 boxparams.k8 = 2.8260078;
97 boxparams.k9 = 9200.476;
98 boxparams.k10 = 13160.295;
99 boxparams.k11 = -21686.9194;
100 boxparams.k12 = 123.56928;
101 boxparams.k13 = -21.1188894;
102 boxparams.k14 = 706.834;
103 boxparams.k15 = 2898.573;
104 boxparams.k16 = 28298.388;
105 boxparams.k17 = 60.81096;
106 boxparams.k18 = 31.242116;
107 boxparams.k19 = 329.574;
108 boxparams.k20 = -2882.082;
109 boxparams.k21 = 74095.3845;
110 boxparams.k22 = -306.262544;
111 boxparams.k23 = 16.243649;
112 boxparams.k24 = -3094.252;
113 boxparams.k25 = -5566.2628;
114 boxparams.k26 = -26237.0;
115 boxparams.k27 = 99.0;
116 boxparams.k28 = -0.42;
117 boxparams.k29 = 1300.0;
118 boxparams.k30 = 2100.0;
119 boxparams.k31 = 925548.252;
120 boxparams.k32 = -61968.8432;
121 boxparams.k33 = 23.3088196;
122 boxparams.k34 = -27097.648;
123 boxparams.k35 = -50843.766;
124
125 //
126 // Initialize the random number generator, so that the results are always the
127 // same.
128 //
129 rand("seed" , 0)
130
131 x0 = [2.52 2.0 37.5 9.25 6.8].';
132 // Compute f(x0) : should be close to -2351244.0
133 [ fx0 , c , index , data ] = boxproblemA ( x0 , 2 , boxparams );
134 mprintf("Computed fx0 = %e (expected = %e)\n",fx0 , -2351244. );
135
136 xopt = [4.53743 2.4 60.0 9.3 7.0].';
137 // Compute f(xopt) : should be -5280334.0
138 [ fopt , c , index , data ] = boxproblemA ( xopt , 2 , boxparams );
139 mprintf("Computed fopt = %e (expected = %e)\n", fopt , -5280334.0 );
140
141 nm = neldermead_new ();
142 nm = neldermead_configure(nm,"-numberofvariables",5);
143 nm = neldermead_configure(nm,"-function",boxproblemA);
144 nm = neldermead_configure(nm,"-costfargument",boxparams);
145 nm = neldermead_configure(nm,"-x0",x0);
146 nm = neldermead_configure(nm,"-maxiter",300);
147 nm = neldermead_configure(nm,"-maxfunevals",1000);
148 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
149 nm = neldermead_configure(nm,"-method","box");
150 nm = neldermead_configure(nm,"-verbose",1);
151 nm = neldermead_configure(nm,"-logfile" , "boxproblemA.txt" );
152 nm = neldermead_configure(nm,"-verbosetermination",1);
153 // Configure like Box
154 nm = neldermead_configure(nm,"-boundsmin",[0.0 1.2 20.0 9.0 6.0]);
155 nm = neldermead_configure(nm,"-boundsmax",[5.0 2.5 60.0 9.3 7.0]);
156 nm = neldermead_configure(nm,"-simplex0method","randbounds");
157 nm = neldermead_configure(nm,"-nbineqconst",6);
158 //
159 // Check that the cost function is correctly connected.
160 //
161 [ nm , f ] = neldermead_function ( nm , x0 );
162 //
163 // Perform optimization
164 //
165 mprintf("Searching (please wait)...\n");
166 nm = neldermead_search(nm);
167 mprintf("...Done\n");
168 neldermead_display(nm);
169 xcomp = neldermead_get(nm,"-xopt");
170 mprintf("x computed=%s\n",strcat(string(xcomp)," "));
171 mprintf("x expected=%s\n",strcat(string(xopt)," "));
172 shift = norm(xcomp-xopt)/norm(xopt);
173 mprintf("Shift =%f\n",shift);
174 fcomp = neldermead_get(nm,"-fopt");
175 mprintf("f computed=%f\n",fcomp);
176 mprintf("f expected=%f\n",fopt);
177 shift = abs(fcomp-fopt)/abs(fopt);
178 mprintf("Shift =%f\n",shift);
179 nm = neldermead_destroy(nm);
180 deletefile ( "boxproblemA.txt" )
181 mprintf("End of demo.\n");
182
183 //
184 // Load this script into the editor
185 //
186 filename = 'neldermead_boxproblemA.sce';
187 dname = get_absolute_file_path(filename);
188 editor ( dname + filename );
189