scilab/modules/optimization/demos/neldermead/neldermead_boxproblemB.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_boxproblemB()\r
  20 \r
  21     filename = 'neldermead_boxproblemB.sce';\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
 110     nm = neldermead_new ();\r
 111     nm = neldermead_configure(nm,"-numberofvariables",2);\r
 112     nm = neldermead_configure(nm,"-function",boxproblemB);\r
 113     nm = neldermead_configure(nm,"-x0",x0);\r
 114     nm = neldermead_configure(nm,"-maxiter",300);\r
 115     nm = neldermead_configure(nm,"-maxfunevals",300);\r
 116     nm = neldermead_configure(nm,"-method","box");\r
 117     nm = neldermead_configure(nm,"-boundsmin",[0.0 0.0]);\r
 118     nm = neldermead_configure(nm,"-boundsmax",[100.0 57.735026918962582]);\r
 119     // Configure like Box\r
 120     nm = neldermead_configure(nm,"-simplex0method","randbounds");\r
 121     nm = neldermead_configure(nm,"-nbineqconst",3);\r
 122     nm = neldermead_configure(nm,"-tolxmethod" , %f );\r
 123     nm = neldermead_configure(nm,"-tolsimplexizemethod",%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
 136     mprintf(_("Searching (please wait) ...\n"));\r
 137     nm = neldermead_search(nm);\r
 138     //\r
 139     // Print a summary\r
 140     //\r
 141     exec(fullfile(dname,"neldermead_summary.sci"),-1);\r
 142     neldermead_summary(nm)\r
 143     mprintf("==========================\n");\r
 144     xcomp = neldermead_get(nm,"-xopt");\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
 148     fcomp = neldermead_get(nm,"-fopt");\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
 152     nm = neldermead_destroy(nm);\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