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
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 //
16 mprintf("Illustrates Box'' algorithm on Box problem B.\n");
17 mprintf("Defining Box Problem A function...\n");
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.
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
72 //
73 // Initialize the random number generator, so that the results are always the
74 // same.
75 //
76 rand("seed" , 0)
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 );
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 );
102 nm = neldermead_configure(nm,"-logfile" , "boxproblemB.txt" );
106 // Configure like Box
109 nm = neldermead_configure(nm,"-tolxmethod" , %f );
111 nm = neldermead_configure(nm,"-boxtermination" , %t );
112 nm = neldermead_configure(nm,"-boxtolf" , 0.001 );
113 nm = neldermead_configure(nm,"-boxboundsalpha" , 0.0001 );
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 //
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);