* Bug 16365 fixed: median(m,'r'|'c') was wrong after 5dc990
[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 // Copyright (C) 2009-2010 - DIGITEO - Michael Baudin
4 // Copyright (C) 2010 - DIGITEO - Allan CORNET
5 // Copyright (C) 2012 - Scilab Enterprises - Adeline CARNIS
6 //
7 // Copyright (C) 2012 - 2016 - Scilab Enterprises
8 //
9 // This file is hereby licensed under the terms of the GNU GPL v2.0,
10 // pursuant to article 5.3.4 of the CeCILL v.2.1.
11 // This file was originally licensed under the terms of the CeCILL v2.1,
12 // and continues to be available under such terms.
13 // For more information, see the COPYING file which you should have received
14 // along with this program.
15
16 //
17 // nmplot_boxproblemA.sce --
18 //   Show that the Box algorithm is able to reproduce the
19 //   numerical experiment presented in Box's paper.
20 //
21
22 function demo_boxproblemB()
23
24     filename = "neldermead_boxproblemB.sce";
25     dname = get_absolute_file_path(filename);
26
27     mprintf(_("Illustrates Box'' algorithm on Box problem B.\n"));
28
29     mprintf("M.J. Box, \n");
30     mprintf(_("""A new method of constrained optimization \n"));
31     mprintf(_("and a comparison with other methods"".\n"));
32     mprintf("The Computer Journal, Volume 8, Number 1, 1965, 42--52\n");
33     mprintf(_("Problem B\n"));
34
35     //
36     //  Reference:
37     //
38     //   M.J.Box,
39     //   "A new method of constrained optimization
40     //   and a comparison with other methods".
41     //   The Computer Journal, Volume 8, Number 1, 1965, 42--52
42     //   Problem A
43     //
44     //   Algorithm 454: the complex method for constrained optimization [E4]
45     //   Communications of the ACM, Volume 16 ,  Issue 8  (August 1973)
46     //   Pages: 487 - 489
47     //
48     //   Richardson and Kuester Results :
49     //   F=1.0000
50     //   X1 = 3.0000
51     //   X2 = 1.7320
52     //   Iterations : 68
53     //
54     //
55     // Note
56     //    The maximum bound for the parameter x1
57     //    is not indicated in Box's paper, but indicated in "Algo 454".
58     //    The maximum bound for x2 is set to 100/sqrt(3) and satisfies the constraint on x2.
59     //    The original problem was to maximize the cost function.
60     //    Here, the problem is to minimize the cost function.
61
62     //
63     // boxproblemB --
64     //   Computes the Box problem B cost function and
65     //   inequality constraints.
66     //
67     // Arguments
68     //   x: the point where to compute the function
69     //   index : the stuff to compute
70     //   data : the parameters of Box cost function
71     //
72     function [ f , c , index ] = boxproblemB ( x , index )
73         f = []
74         c = []
75         x3 = x(1) + sqrt(3.0) * x(2)
76
77         if ( index==2 | index==6 ) then
78             f = -(9.0 - (x(1) - 3.0) ^ 2) * x(2) ^ 3 / 27.0 / sqrt(3.0)
79         end
80
81         if ( index==5 | index==6 ) then
82             c1 = x(1) / sqrt(3.0) - x(2)
83             c2 = x3
84             c3 = 6.0 - x3
85             c = [c1 c2 c3]
86         end
87
88     endfunction
89
90
91     //
92     // Initialize the random number generator, so that the results are always the
93     // same.
94     //
95     rand("seed" , 0)
96
97     x0 = [1.0 0.5].';
98     // Compute f(x0) : should be close to -0.0133645895646
99     fx0 = boxproblemB ( x0 , 2 );
100     mprintf("Computed fx0 = %e (expected = %e)\n",fx0 , -0.0133645895646 );
101     [fx0 , cx0 , index ] = boxproblemB ( x0 , 6 );
102     mprintf("Computed Constraints(x0) = [%e %e %e]\n", ...
103     cx0(1), cx0(2), cx0(3) );
104     mprintf("Expected Constraints(x0) = [%e %e %e]\n", ...
105     0.0773503 , 1.8660254 , 4.1339746 );
106
107
108     xopt = [3.0 1.7320508075688774].';
109     // Compute f(xopt) : should be -1.0
110     fopt = boxproblemB ( xopt , 2 );
111     mprintf("Computed fopt = %e (expected = %e)\n", fopt , -1.0 );
112
113     nm = neldermead_new ();
114     nm = neldermead_configure(nm,"-numberofvariables",2);
115     nm = neldermead_configure(nm,"-function",boxproblemB);
116     nm = neldermead_configure(nm,"-x0",x0);
117     nm = neldermead_configure(nm,"-maxiter",300);
118     nm = neldermead_configure(nm,"-maxfunevals",300);
119     nm = neldermead_configure(nm,"-method","box");
120     nm = neldermead_configure(nm,"-boundsmin",[0.0 0.0]);
121     nm = neldermead_configure(nm,"-boundsmax",[100.0 57.735026918962582]);
122     // Configure like Box
123     nm = neldermead_configure(nm,"-simplex0method","randbounds");
124     nm = neldermead_configure(nm,"-nbineqconst",3);
125     nm = neldermead_configure(nm,"-tolxmethod" , %f );
126     nm = neldermead_configure(nm,"-tolsimplexizemethod",%f);
127     nm = neldermead_configure(nm,"-boxtermination" , %t );
128     nm = neldermead_configure(nm,"-boxtolf" , 0.001 );
129     nm = neldermead_configure(nm,"-boxboundsalpha" , 0.0001 );
130
131     //
132     // Check that the cost function is correctly connected.
133     //
134     [ nm , f ] = neldermead_function ( nm , x0 );
135
136     //
137     // Perform optimization
138     //
139     mprintf(_("Searching (please wait) ...\n"));
140     nm = neldermead_search(nm);
141     //
142     // Print a summary
143     //
144     exec(fullfile(dname,"neldermead_summary.sci"),-1);
145     neldermead_summary(nm)
146     mprintf("==========================\n");
147     xcomp = neldermead_get(nm,"-xopt");
148     mprintf("x expected = [%s]\n", strcat(string(xopt), " "));
149     shift = norm(xcomp-xopt)/norm(xopt);
150     mprintf(_("Relative error on x = %e\n"), shift);
151     fcomp = neldermead_get(nm,"-fopt");
152     mprintf(_("f expected = %f\n"), fopt);
153     shift = abs(fcomp-fopt)/abs(fopt);
154     mprintf(_("Relative error on f = %e\n"), shift);
155     nm = neldermead_destroy(nm);
156     mprintf(_("End of demo.\n"));
157
158     //
159     // Load this script into the editor
160     //
161     m = messagebox(_("View Code?"), "Question", "question", [_("Yes") _("No")], "modal")
162     if(m == 1)
163         editor ( dname + filename, "readonly" );
164     end
165 endfunction
166
167 demo_boxproblemB();
168 clear demo_boxproblemB;
169
170
171
172
173
174
175
176