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
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_boxproblemB()\r
20 \r
22     dname = get_absolute_file_path(filename);\r
23 \r
24     mprintf(_("Illustrates Box'' algorithm on Box problem B.\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 B\n"));\r
31 \r
32     //\r
33     //  Reference:\r
34     //\r
35     //   M.J.Box, \r
36     //   "A new method of constrained optimization \r
37     //   and a comparison with other methods".\r
38     //   The Computer Journal, Volume 8, Number 1, 1965, 42--52\r
39     //   Problem A\r
40     //\r
41     //   Algorithm 454: the complex method for constrained optimization [E4]\r
42     //   Communications of the ACM, Volume 16 ,  Issue 8  (August 1973)\r
43     //   Pages: 487 - 489   \r
44     //\r
45     //   Richardson and Kuester Results :\r
46     //   F=1.0000\r
47     //   X1 = 3.0000\r
48     //   X2 = 1.7320\r
49     //   Iterations : 68\r
50     //\r
51     //\r
52     // Note \r
53     //    The maximum bound for the parameter x1\r
54     //    is not indicated in Box's paper, but indicated in "Algo 454".\r
55     //    The maximum bound for x2 is set to 100/sqrt(3) and satisfies the constraint on x2.\r
56     //    The original problem was to maximize the cost function.\r
57     //    Here, the problem is to minimize the cost function.\r
58 \r
59     //\r
60     // boxproblemB --\r
61     //   Computes the Box problem B cost function and \r
62     //   inequality constraints.\r
63     //\r
64     // Arguments\r
65     //   x: the point where to compute the function\r
66     //   index : the stuff to compute\r
67     //   data : the parameters of Box cost function\r
68     //\r
69     function [ f , c , index ] = boxproblemB ( x , index )\r
70         f = []\r
71         c = []\r
72         x3 = x(1) + sqrt(3.0) * x(2)\r
73 \r
74         if ( index==2 | index==6 ) then\r
75             f = -(9.0 - (x(1) - 3.0) ^ 2) * x(2) ^ 3 / 27.0 / sqrt(3.0)\r
76         end\r
77 \r
78         if ( index==5 | index==6 ) then\r
79             c1 = x(1) / sqrt(3.0) - x(2)\r
80             c2 = x3\r
81             c3 = 6.0 - x3\r
82             c = [c1 c2 c3]\r
83         end\r
84 \r
85     endfunction\r
86 \r
87 \r
88     //\r
89     // Initialize the random number generator, so that the results are always the\r
90     // same.\r
91     //\r
92     rand("seed" , 0)\r
93 \r
94     x0 = [1.0 0.5].';\r
95     // Compute f(x0) : should be close to -0.0133645895646\r
96     fx0 = boxproblemB ( x0 , 2 );\r
97     mprintf("Computed fx0 = %e (expected = %e)\n",fx0 , -0.0133645895646 );\r
98     [fx0 , cx0 , index ] = boxproblemB ( x0 , 6 );\r
99     mprintf("Computed Constraints(x0) = [%e %e %e]\n", ...\r
100     cx0(1), cx0(2), cx0(3) );\r
101     mprintf("Expected Constraints(x0) = [%e %e %e]\n", ...\r
102     0.0773503 , 1.8660254 , 4.1339746 );\r
103 \r
104 \r
105     xopt = [3.0 1.7320508075688774].';\r
106     // Compute f(xopt) : should be -1.0\r
107     fopt = boxproblemB ( xopt , 2 );\r
108     mprintf("Computed fopt = %e (expected = %e)\n", fopt , -1.0 );\r
109 \r
119     // Configure like Box\r
122     nm = neldermead_configure(nm,"-tolxmethod" , %f );\r
124     nm = neldermead_configure(nm,"-boxtermination" , %t );\r
125     nm = neldermead_configure(nm,"-boxtolf" , 0.001 );\r
126     nm = neldermead_configure(nm,"-boxboundsalpha" , 0.0001 );\r
127 \r
128     //\r
129     // Check that the cost function is correctly connected.\r
130     //\r
131     [ nm , f ] = neldermead_function ( nm , x0 );\r
132 \r
133     //\r
134     // Perform optimization\r
135     //\r
138     //\r
139     // Print a summary\r
140     //\r
143     mprintf("==========================\n");\r
145     mprintf("x expected = [%s]\n", strcat(string(xopt), " "));\r
146     shift = norm(xcomp-xopt)/norm(xopt);\r
147     mprintf(_("Relative error on x = %e\n"), shift);\r
149     mprintf(_("f expected = %f\n"), fopt);\r
150     shift = abs(fcomp-fopt)/abs(fopt);\r
151     mprintf(_("Relative error on f = %e\n"), shift);\r
153     mprintf(_("End of demo.\n"));\r
154 \r
155     //\r
156     // Load this script into the editor\r
157     //\r
158     m = messagebox(_("View Code?"), "Question", "question", [_("Yes") _("No")], "modal")\r
159     if(m == 1)\r
160         editor ( dname + filename, "readonly" );\r
161     end\r
162 endfunction\r
163 \r
164 demo_boxproblemB();\r
165 clear demo_boxproblemB;\r
166 \r
167 \r
168 \r
169 \r
170 \r
171 \r
172 \r\r