b4acff36205caad1b6911dc71fa4ccb46b787f74
[scilab.git] / scilab / modules / optimization / demos / neldermead / neldermead_boxproblemB.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 B.\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 //   Richardson and Kuester Results :
33 //   F=1.0000
34 //   X1 = 3.0000
35 //   X2 = 1.7320
36 //   Iterations : 68
37 //
38 //
39 // Note 
40 //    The maximum bound for the parameter x1
41 //    is not indicated in Box's paper, but indicated in "Algo 454".
42 //    The maximum bound for x2 is set to 100/sqrt(3) and satisfies the constraint on x2.
43 //    The original problem was to maximize the cost function.
44 //    Here, the problem is to minimize the cost function.
45
46 //
47 // boxproblemB --
48 //   Computes the Box problem B cost function and 
49 //   inequality constraints.
50 //
51 // Arguments
52 //   x: the point where to compute the function
53 //   index : the stuff to compute
54 //   data : the parameters of Box cost function
55 //
56 function [ f , c , index ] = boxproblemB ( x , index )
57   f = []
58   c = []
59   x3 = x(1) + sqrt(3.0) * x(2)
60   if ( index==1 | index==3 ) then
61     f = -(9.0 - (x(1) - 3.0) ^ 2) * x(2) ^ 3 / 27.0 / sqrt(3.0)
62   end
63   if ( index==2 | index==3 ) then
64       c1 = x(1) / sqrt(3.0) - x(2)
65       c2 = x3
66       c3 = 6.0 - x3
67       c = [c1 c2 c3]
68   end
69 endfunction
70
71
72 //
73 // Initialize the random number generator, so that the results are always the
74 // same.
75 //
76 rand("seed" , 0)
77
78 x0 = [1.0 0.5].';
79 // Compute f(x0) : should be close to -0.0133645895646
80 fx0 = boxproblemB ( x0 );
81 mprintf("Computed fx0 = %e (expected = %e)\n",fx0 , -0.0133645895646 );
82 result = boxproblemB ( x0 , 2 );
83 mprintf("Computed Constraints(x0) = [%e %e %e]\n", ...
84   result(1), result(2), result(3) );
85 mprintf("Expected Constraints(x0) = [%e %e %e]\n", ...
86   0.0773503 , 1.8660254 , 4.1339746 );
87
88
89 xopt = [3.0 1.7320508075688774].';
90 // Compute f(xopt) : should be -1.0
91 fopt = boxproblemB ( xopt );
92 mprintf("Computed fopt = %e (expected = %e)\n", fopt , -1.0 );
93
94 nm = neldermead_new ();
95 nm = neldermead_configure(nm,"-numberofvariables",2);
96 nm = neldermead_configure(nm,"-function",boxproblemB);
97 nm = neldermead_configure(nm,"-x0",x0);
98 nm = neldermead_configure(nm,"-maxiter",300);
99 nm = neldermead_configure(nm,"-maxfunevals",300);
100 nm = neldermead_configure(nm,"-method","box");
101 nm = neldermead_configure(nm,"-verbose",1);
102 nm = neldermead_configure(nm,"-logfile" , "boxproblemB.txt" );
103 nm = neldermead_configure(nm,"-verbosetermination",1);
104 nm = neldermead_configure(nm,"-boundsmin",[0.0 0.0]);
105 nm = neldermead_configure(nm,"-boundsmax",[100.0 57.735026918962582]);
106 // Configure like Box
107 nm = neldermead_configure(nm,"-simplex0method","randbounds");
108 nm = neldermead_configure(nm,"-nbineqconst",3);
109 nm = neldermead_configure(nm,"-tolxmethod" , %f );
110 nm = neldermead_configure(nm,"-tolsimplexizemethod",%f);
111 nm = neldermead_configure(nm,"-boxtermination" , %t );
112 nm = neldermead_configure(nm,"-boxtolf" , 0.001 );
113 nm = neldermead_configure(nm,"-boxboundsalpha" , 0.0001 );
114
115 //
116 // Check that the cost function is correctly connected.
117 // The index must be provided, because the additionnal argument "data"
118 // comes after.
119 //
120 [ nm , result ] = neldermead_function ( nm , x0 , 1 );
121 //
122 // Perform optimization
123 //
124 nm = neldermead_search(nm);
125 neldermead_display(nm);
126 xcomp = neldermead_get(nm,"-xopt");
127 mprintf("x computed=%s\n",strcat(string(xcomp)," "));
128 mprintf("x expected=%s\n",strcat(string(xopt)," "));
129 shift = norm(xcomp-xopt)/norm(xopt);
130 mprintf("Shift =%f\n",shift);
131 fcomp = neldermead_get(nm,"-fopt");
132 mprintf("f computed=%f\n",fcomp);
133 mprintf("f expected=%f\n",fopt);
134 shift = abs(fcomp-fopt)/abs(fopt);
135 mprintf("Shift =%f\n",shift);
136 nm = neldermead_destroy(nm);
137 deletefile ( "boxproblemB.txt" )
138 //
139 // Load this script into the editor
140 //
141 filename = 'neldermead_boxproblemB.sce';
142 dname = get_absolute_file_path(filename);
143 editor ( dname + filename );
144