Deleted vectorized computation feature. Deleted neldermead_contour. Fixed 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 // boxproblemA --
21 //   Computes the Box problem A cost function and 
22 //   inequality constraints.
23 //
24 // Arguments
25 //   x: the point where to compute the function
26 //   index : the stuff to compute
27 //   data : the parameters of Box cost function
28 // Note
29 //  The following protocol is used
30 //  * if index=1, or no index, returns the value of the cost 
31 //    function (default case)
32 //  * if index=2, returns the value of the nonlinear inequality 
33 //    constraints, as a row array
34 //  * if index=3, returns an array which contains
35 //    at index #0, the value of the cost function  
36 //    at index #1 to the end is the list of the values of the nonlinear 
37 //    constraints
38 //  The inequality constraints are expected to be positive.
39 //
40 function result = boxproblemA ( x , index , data )
41   if (~isdef('index','local')) then
42     index = 1
43   end
44   b = x(2) + 0.01 * x(3)
45   x6 = (data.k1 + data.k2 * x(2) ...
46             + data.k3 * x(3) + data.k4 * x(4) + data.k5 * x(5)) * x(1)
47   y1 = data.k6 + data.k7 * x(2) + data.k8 * x(3) ...
48             + data.k9 * x(4) + data.k10 * x(5)
49   y2 = data.k11 + data.k12 * x(2) + data.k13 * x(3) ...
50             + data.k14 * x(4) + data.k15 * x(5)
51   y3 = data.k16 + data.k17 * x(2) + data.k18 * x(3) ...
52             + data.k19 * x(4) + data.k20 * x(5)
53   y4 = data.k21 + data.k22 * x(2) + data.k23 * x(3) ...
54             + data.k24 * x(4) + data.k25 * x(5)
55   x7 = ( y1 + y2 + y3 ) * x(1)
56   x8 = (data.k26 + data.k27 * x(2) + data.k28 * x(3) ...
57             + data.k29 * x(4) ...
58             + data.k30 * x(5) ) * x(1) + x6 + x7
59   if ( index==1 | index==3 ) then
60     f = (data.a2 * y1 + data.a3 * y2 + data.a4 * y3 + data.a5 * y4 ...
61              + 7840 * data.a6 - 100000 * data.a0 ...
62              - 50800 * b * data.a7 + data.k31 + data.k32 * x(2) + data.k33 * x(3) ...
63              + data.k34 * x(4) + data.k35 * x(5)) * x(1) ...
64              - 24345 + data.a1 * x6
65     f = -f
66   end
67   if ( index==2 | index==3 ) then
68       c1 = x6
69       c2 = 294000 - x6
70       c3 = x7
71       c4 = 294000 - x7
72       c5 = x8
73       c6 = 277200 - x8
74   end
75   select index
76   case 1 then
77     result = f
78   case 2 then
79     result = [c1 c2 c3 c4 c5 c6]
80   case 3 then
81     result = [f c1 c2 c3 c4 c5 c6]
82   else
83     errmsg = sprintf("Unexpected index %d" , index);
84     error(errmsg);
85   end
86 endfunction
87
88 boxparams = struct();
89 boxparams.a0 = 9;
90 boxparams.a1 = 15;
91 boxparams.a2 = 50;
92 boxparams.a3 = 9.583;
93 boxparams.a4 = 20;
94 boxparams.a5 = 15;
95 boxparams.a6 = 6;
96 boxparams.a7 = 0.75;
97 boxparams.k1 = -145421.402;
98 boxparams.k2 = 2931.1506;
99 boxparams.k3 = -40.427932;
100 boxparams.k4 = 5106.192;
101 boxparams.k5 = 15711.36;
102 boxparams.k6 = -161622.577;
103 boxparams.k7 = 4176.15328;
104 boxparams.k8 = 2.8260078;
105 boxparams.k9 = 9200.476;
106 boxparams.k10 = 13160.295;
107 boxparams.k11 = -21686.9194;
108 boxparams.k12 = 123.56928;
109 boxparams.k13 = -21.1188894;
110 boxparams.k14 = 706.834;
111 boxparams.k15 = 2898.573;
112 boxparams.k16 = 28298.388;
113 boxparams.k17 = 60.81096;
114 boxparams.k18 = 31.242116;
115 boxparams.k19 = 329.574;
116 boxparams.k20 = -2882.082;
117 boxparams.k21 = 74095.3845;
118 boxparams.k22 = -306.262544;
119 boxparams.k23 = 16.243649;
120 boxparams.k24 = -3094.252;
121 boxparams.k25 = -5566.2628;
122 boxparams.k26 = -26237.0;
123 boxparams.k27 = 99.0;
124 boxparams.k28 = -0.42;
125 boxparams.k29 = 1300.0;
126 boxparams.k30 = 2100.0;
127 boxparams.k31 = 925548.252;
128 boxparams.k32 = -61968.8432;
129 boxparams.k33 = 23.3088196;
130 boxparams.k34 = -27097.648;
131 boxparams.k35 = -50843.766;
132
133     //
134     // Initialize the random number generator, so that the results are always the
135     // same.
136     //
137     rand("seed" , 0)
138
139 x0 = [2.52 2.0 37.5 9.25 6.8].';
140 // Compute f(x0) : should be close to -2351244.0
141 fx0 = boxproblemA ( x0 , data = boxparams );
142 xopt = [4.53743 2.4 60.0 9.3 7.0].';
143 // Compute f(xopt) : should be -5280334.0
144 fopt = boxproblemA ( xopt , data = boxparams );
145
146 nm = neldermead_new ();
147 nm = neldermead_configure(nm,"-numberofvariables",5);
148 nm = neldermead_configure(nm,"-function",boxproblemA);
149 nm = neldermead_configure(nm,"-costfargument",boxparams);
150 nm = neldermead_configure(nm,"-x0",x0);
151 nm = neldermead_configure(nm,"-maxiter",300);
152 nm = neldermead_configure(nm,"-maxfunevals",1000);
153 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
154 nm = neldermead_configure(nm,"-method","box");
155 //nm = neldermead_configure(nm,"-verbose",1);
156 nm = neldermead_configure(nm,"-verbosetermination",1);
157 // Configure like Box
158 nm = neldermead_configure(nm,"-boundsmin",[0.0 1.2 20.0 9.0 6.0]);
159 nm = neldermead_configure(nm,"-boundsmax",[5.0 2.5 60.0 9.3 7.0]);
160 nm = neldermead_configure(nm,"-simplex0method","randbounds");
161 nm = neldermead_configure(nm,"-nbineqconst",6);
162 //
163 // Check that the cost function is correctly connected.
164 // The index must be provided, because the additionnal argument "data"
165 // comes after.
166 //
167 [ nm , result ] = neldermead_function ( nm , x0 , 1 );
168 //
169 // Perform optimization
170 //
171 nm = neldermead_search(nm);
172 neldermead_display(nm);
173 xcomp = neldermead_get(nm,"-xopt");
174 mprintf("x computed=%s\n",strcat(string(xcomp)," "));
175 mprintf("x expected=%s\n",strcat(string(xopt)," "));
176 shift = norm(xcomp-xopt)/norm(xopt);
177 mprintf("Shift =%f\n",shift);
178 fcomp = neldermead_get(nm,"-fopt");
179 mprintf("f computed=%f\n",fcomp);
180 mprintf("f expected=%f\n",fopt);
181 shift = abs(fcomp-fopt)/abs(fopt);
182 mprintf("Shift =%f\n",shift);
183 nm = neldermead_destroy(nm);
184
185