Avoid problem with localization tools
[scilab.git] / scilab / modules / optimization / demos / neldermead / neldermead_boxproblemA.sce
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab\r
2 // Copyright (C) 2008-2009 - INRIA - Michael Baudin\r
3 // Copyright (C) 2009-2010 - DIGITEO - Michael Baudin\r
4 // Copyright (C) 2010 - DIGITEO - Allan CORNET\r
5 // Copyright (C) 2012 - Scilab Enterprises - Adeline CARNIS\r
6 //\r
7 // This file must be used under the terms of the CeCILL.\r
8 // This source file is licensed as described in the file COPYING, which\r
9 // you should have received as part of this distribution.  The terms\r
10 // are also available at\r
11 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt\r
12 \r
13 //\r
14 // nmplot_boxproblemA.sce --\r
15 //   Show that the Box algorithm is able to reproduce the \r
16 //   numerical experiment presented in Box's paper.\r
17 //\r
18 \r
19 function demo_boxproblemA()\r
20 \r
21     filename = 'neldermead_boxproblemA.sce';\r
22     dname = get_absolute_file_path(filename);\r
23 \r
24     mprintf(_("Illustrates Box'' algorithm on Box problem A.\n"));\r
25 \r
26     mprintf("M.J. Box, \n");\r
27     mprintf(_("""A new method of constrained optimization \n"));\r
28     mprintf(_("and a comparison with other methods"".\n"));\r
29     mprintf("The Computer Journal, Volume 8, Number 1, 1965, 42--52\n");\r
30     mprintf(_("Problem A\n"));\r
31 \r
32 \r
33     //\r
34     //  Reference:\r
35     //\r
36     //   M.J. Box, \r
37     //   "A new method of constrained optimization \r
38     //   and a comparison with other methods".\r
39     //   The Computer Journal, Volume 8, Number 1, 1965, 42--52\r
40     //   Problem A\r
41     //\r
42     //   Algorithm 454: the complex method for constrained optimization [E4]\r
43     //   Communications of the ACM, Volume 16 ,  Issue 8  (August 1973)\r
44     //   Pages: 487 - 489   \r
45     //\r
46 \r
47     //\r
48     // boxproblemA --\r
49     //   Computes the Box problem A cost function and \r
50     //   inequality constraints.\r
51     //\r
52     // Arguments\r
53     //   x: the point where to compute the function\r
54     //   index : the stuff to compute\r
55     //   data : the parameters of Box cost function\r
56     //\r
57     function [ f , c , index ] = boxproblemA ( x , index , data )\r
58         f = []\r
59         c = []\r
60         b = x(2) + 0.01 * x(3)\r
61         x6 = (data.k1 + data.k2 * x(2) ...\r
62         + data.k3 * x(3) + data.k4 * x(4) + data.k5 * x(5)) * x(1)\r
63         y1 = data.k6 + data.k7 * x(2) + data.k8 * x(3) ...\r
64         + data.k9 * x(4) + data.k10 * x(5)\r
65         y2 = data.k11 + data.k12 * x(2) + data.k13 * x(3) ...\r
66         + data.k14 * x(4) + data.k15 * x(5)\r
67         y3 = data.k16 + data.k17 * x(2) + data.k18 * x(3) ...\r
68         + data.k19 * x(4) + data.k20 * x(5)\r
69         y4 = data.k21 + data.k22 * x(2) + data.k23 * x(3) ...\r
70         + data.k24 * x(4) + data.k25 * x(5)\r
71         x7 = ( y1 + y2 + y3 ) * x(1)\r
72         x8 = (data.k26 + data.k27 * x(2) + data.k28 * x(3) ...\r
73         + data.k29 * x(4) ...\r
74         + data.k30 * x(5) ) * x(1) + x6 + x7\r
75 \r
76         if ( index == 2 | index == 6 ) then\r
77             f = (data.a2 * y1 + data.a3 * y2 + data.a4 * y3 + data.a5 * y4 ...\r
78             + 7840 * data.a6 - 100000 * data.a0 ...\r
79             - 50800 * b * data.a7 + data.k31 + data.k32 * x(2) + data.k33 * x(3) ...\r
80             + data.k34 * x(4) + data.k35 * x(5)) * x(1) ...\r
81             - 24345 + data.a1 * x6\r
82             f = -f\r
83         end\r
84 \r
85         if ( index == 5 | index == 6 ) then\r
86             c1 = x6\r
87             c2 = 294000 - x6\r
88             c3 = x7\r
89             c4 = 294000 - x7\r
90             c5 = x8\r
91             c6 = 277200 - x8\r
92             c = [c1 c2 c3 c4 c5 c6]\r
93         end\r
94     endfunction\r
95 \r
96     boxparams = struct();\r
97     boxparams.a0 = 9;\r
98     boxparams.a1 = 15;\r
99     boxparams.a2 = 50;\r
100     boxparams.a3 = 9.583;\r
101     boxparams.a4 = 20;\r
102     boxparams.a5 = 15;\r
103     boxparams.a6 = 6;\r
104     boxparams.a7 = 0.75;\r
105     boxparams.k1 = -145421.402;\r
106     boxparams.k2 = 2931.1506;\r
107     boxparams.k3 = -40.427932;\r
108     boxparams.k4 = 5106.192;\r
109     boxparams.k5 = 15711.36;\r
110     boxparams.k6 = -161622.577;\r
111     boxparams.k7 = 4176.15328;\r
112     boxparams.k8 = 2.8260078;\r
113     boxparams.k9 = 9200.476;\r
114     boxparams.k10 = 13160.295;\r
115     boxparams.k11 = -21686.9194;\r
116     boxparams.k12 = 123.56928;\r
117     boxparams.k13 = -21.1188894;\r
118     boxparams.k14 = 706.834;\r
119     boxparams.k15 = 2898.573;\r
120     boxparams.k16 = 28298.388;\r
121     boxparams.k17 = 60.81096;\r
122     boxparams.k18 = 31.242116;\r
123     boxparams.k19 = 329.574;\r
124     boxparams.k20 = -2882.082;\r
125     boxparams.k21 = 74095.3845;\r
126     boxparams.k22 = -306.262544;\r
127     boxparams.k23 = 16.243649;\r
128     boxparams.k24 = -3094.252;\r
129     boxparams.k25 = -5566.2628;\r
130     boxparams.k26 = -26237.0;\r
131     boxparams.k27 = 99.0;\r
132     boxparams.k28 = -0.42;\r
133     boxparams.k29 = 1300.0;\r
134     boxparams.k30 = 2100.0;\r
135     boxparams.k31 = 925548.252;\r
136     boxparams.k32 = -61968.8432;\r
137     boxparams.k33 = 23.3088196;\r
138     boxparams.k34 = -27097.648;\r
139     boxparams.k35 = -50843.766;\r
140 \r
141     //\r
142     // Initialize the random number generator, so that the results are always the\r
143     // same.\r
144     //\r
145     rand("seed" , 0)\r
146 \r
147     x0 = [2.52 2.0 37.5 9.25 6.8].';\r
148     // Compute f(x0) : should be close to -2351244.0\r
149     [ fx0 , c , index ] = boxproblemA ( x0 , 2 , boxparams );\r
150     mprintf("Computed fx0 = %e (expected = %e)\n",fx0 , -2351244. );\r
151 \r
152     xopt = [4.53743 2.4 60.0 9.3 7.0].';\r
153     // Compute f(xopt) : should be -5280334.0\r
154     [ fopt , c , index ] = boxproblemA ( xopt , 2 , boxparams );\r
155     mprintf("Computed fopt = %e (expected = %e)\n", fopt , -5280334.0 );\r
156 \r
157     nm = neldermead_new ();\r
158     nm = neldermead_configure(nm,"-numberofvariables",5);\r
159     nm = neldermead_configure(nm,"-function",list(boxproblemA,boxparams));\r
160     nm = neldermead_configure(nm,"-x0",x0);\r
161     nm = neldermead_configure(nm,"-maxiter",300);\r
162     nm = neldermead_configure(nm,"-maxfunevals",1000);\r
163     nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);\r
164     nm = neldermead_configure(nm,"-method","box");\r
165 \r
166     // Configure like Box\r
167     nm = neldermead_configure(nm,"-boundsmin",[0.0 1.2 20.0 9.0 6.0]);\r
168     nm = neldermead_configure(nm,"-boundsmax",[5.0 2.5 60.0 9.3 7.0]);\r
169     nm = neldermead_configure(nm,"-simplex0method","randbounds");\r
170     nm = neldermead_configure(nm,"-nbineqconst",6);\r
171 \r
172     //\r
173     // Check that the cost function is correctly connected.\r
174     //\r
175     [ nm , f ] = neldermead_function ( nm , x0 );\r
176 \r
177     //\r
178     // Perform optimization\r
179     //\r
180     mprintf(_("Searching (please wait)...\n"));\r
181     nm = neldermead_search(nm);\r
182     mprintf(_("...Done\n"));\r
183     //\r
184     // Print a summary\r
185     //\r
186     exec(fullfile(dname,"neldermead_summary.sci"),-1);\r
187     neldermead_summary(nm)\r
188     mprintf("==========================\n");\r
189     xcomp = neldermead_get(nm,"-xopt");\r
190     mprintf("x expected = [%s]\n",strcat(string(xopt)," "));\r
191     shift = norm(xcomp-xopt)/norm(xopt);\r
192     mprintf(_("Relative error on x = %e\n"),shift);\r
193     fcomp = neldermead_get(nm,"-fopt");\r
194     mprintf("f expected = %f\n",fopt);\r
195     shift = abs(fcomp-fopt)/abs(fopt);\r
196     mprintf("Relative error on f = %e\n",shift);\r
197     nm = neldermead_destroy(nm);\r
198     mprintf(_("End of demo.\n"));\r
199 \r
200     //\r
201     // Load this script into the editor\r
202     //\r
203     m = messagebox(_("View Code?"), "Question", "question", [_("Yes") _("No")], "modal")\r
204     if(m == 1)\r
205         editor ( dname + filename, "readonly" );\r
206     end\r
207 endfunction \r
208 \r
209 demo_boxproblemA();\r
210 clear demo_boxproblemA;\r
211 \r
212 \r
213 \r
214 \r
215 \r
216 \r
217 \r\r