1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 2// Copyright (C) 2008-2009 - INRIA - Michael Baudin 3// Copyright (C) 2011 - DIGITEO - Michael Baudin 4// 5// Copyright (C) 2012 - 2016 - Scilab Enterprises 6// 7// This file is hereby licensed under the terms of the GNU GPL v2.0, 8// pursuant to article 5.3.4 of the CeCILL v.2.1. 9// This file was originally licensed under the terms of the CeCILL v2.1, 10// and continues to be available under such terms. 11// For more information, see the COPYING file which you should have received 12// along with this program. 13// <-- CLI SHELL MODE --> 14// <-- ENGLISH IMPOSED --> 15// 16// Reference: 17// 18// An extension of the simplex method to constrained 19// nonlinear optimization 20// M.B. Subrahmanyam 21// Journal of optimization theory and applications 22// Vol. 62, August 1989 23// 24// Gould F.J. 25// Nonlinear Tolerance Programming 26// Numerical methods for Nonlinear optimization 27// Edited by F.A. Lootsma, pp 349-366, 1972 28// 29// optimtestcase -- 30// Non linear inequality constraints are positive. 31// 32// Arguments 33// x: the point where to compute the function 34// index : the stuff to compute 35// 36function [ f , c , index ] = optimtestcase ( x , index ) 37 f = [] 38 c = [] 39 if ( ( index == 2 ) | ( index == 6 ) ) then 40 f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ... 41 - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4) 42 end 43 if ( ( index == 5 ) | ( index == 6 ) ) then 44 c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ... 45 - x(1) + x(2) - x(3) + x(4) + 8 46 c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ... 47 + x(1) + x(4) + 10.0 48 c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ... 49 + x(2) + x(4) + 5.0 50 c = [c1 c2 c3] 51 end 52endfunction 53// 54// neldermead_constraints -- 55// The Nelder-Mead algorithm, with variable-size simplex 56// and modifications for bounds and 57// inequality constraints. 58// 59function this = neldermead_constraints ( this ) 60 // Check settings correspond to algo 61 [ this.optbase , hascons ] = optimbase_hasconstraints ( this.optbase ); 62 if ( ~hascons ) then 63 errmsg = msprintf(gettext("%s: Problem has no constraints, but variable algorithm is designed for them."), "neldermead_constraints") 64 error(errmsg) 65 end 66 verbose = optimbase_cget ( this.optbase , "-verbose" ) 67 // 68 // Order the vertices for the first time 69 // 70 simplex = this.simplex0; 71 n = optimbase_cget ( this.optbase , "-numberofvariables" ); 72 fvinitial = optimbase_get ( this.optbase , "-fx0" ); 73 // Sort function values and x points by increasing function value order 74 this = neldermead_log (this,"Step #1 : order"); 75 simplex = optimsimplex_sort ( simplex ); 76 currentcenter = optimsimplex_center ( simplex ); 77 currentxopt = optimbase_cget ( this.optbase , "-x0" ); 78 newfvmean = optimsimplex_fvmean ( simplex ); 79 nbve = optimsimplex_getnbve ( simplex ); 80 ihigh = nbve; 81 inext = ihigh - 1 82 ilow = 1 83 [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase ); 84 nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" ) 85 // 86 // Initialize 87 // 88 terminate = %f; 89 iter = 0; 90 step = "init"; 91 // 92 // Nelder-Mead Loop 93 // 94 while ( ~terminate ) 95 this.optbase = optimbase_incriter ( this.optbase ); 96 iter = iter + 1; 97 xlow = optimsimplex_getx ( simplex , ilow ) 98 flow = optimsimplex_getfv ( simplex , ilow ) 99 xhigh = optimsimplex_getx ( simplex , ihigh ) 100 fhigh = optimsimplex_getfv ( simplex , ihigh ) 101 xn = optimsimplex_getx ( simplex , inext ) 102 fn = optimsimplex_getfv ( simplex , inext ) 103 // 104 // Store history 105 // 106 xcoords = optimsimplex_getallx ( simplex ) 107 this = neldermead_storehistory ( this , n , flow , xlow , xcoords ); 108 currentfopt = flow; 109 previousxopt = currentxopt; 110 currentxopt = xlow; 111 previouscenter = currentcenter; 112 currentcenter = optimsimplex_center ( simplex ); 113 oldfvmean = newfvmean; 114 newfvmean = optimsimplex_fvmean ( simplex ); 115 if ( verbose == 1 ) then 116 deltafv = abs(optimsimplex_deltafvmax ( simplex )); 117 totaliter = optimbase_get ( this.optbase , "-iterations" ); 118 funevals = optimbase_get ( this.optbase , "-funevals" ); 119 ssize = optimsimplex_size ( simplex ) 120 this = neldermead_log (this,sprintf("=================================================================")); 121 this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter)); 122 this = neldermead_log (this,sprintf("Function Eval #%d",funevals)); 123 this = neldermead_log (this,sprintf("Xopt : [%s]",_strvec(xlow))); 124 this = neldermead_log (this,sprintf("Fopt : %e",flow)); 125 this = neldermead_log (this,sprintf("DeltaFv : %e",deltafv)); 126 this = neldermead_log (this,sprintf("Center : [%s]",_strvec(currentcenter))); 127 this = neldermead_log (this,sprintf("Size : %e",ssize)); 128 str = optimsimplex_tostring ( simplex ) 129 for i = 1:nbve 130 this = neldermead_log (this,str(i)); 131 end 132 end 133 neldermead_outputcmd ( this, "iter" , simplex , step ) 134 // 135 // Update termination flag 136 // 137 if ( iter > 1 ) then 138 [ this , terminate , status] = neldermead_termination (this , ... 139 fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex ); 140 if ( terminate ) then 141 if ( verbose == 1 ) then 142 this = neldermead_log (this,sprintf("Terminate with status : %s",status)); 143 end 144 break 145 end 146 end 147 // 148 // Compute xbar, center of better vertices 149 // 150 if ( verbose == 1 ) then 151 this = neldermead_log (this,sprintf("Reflect")); 152 end 153 xbar = optimsimplex_xbar ( simplex ); 154 if ( verbose == 1 ) then 155 this = neldermead_log (this,sprintf("xbar=[%s]" , _strvec(xbar))); 156 end 157 // 158 // Reflect the worst point with respect to center 159 // 160 xr = neldermead_interpolate ( xbar , xhigh , this.rho ); 161 if ( verbose == 1 ) then 162 this = neldermead_log (this,sprintf("xr=[%s]" , _strvec(xr))); 163 end 164 // Adjust point to satisfy bounds and nonlinear inequality constraints 165 if ( hasbounds | nbnlc > 0 ) then 166 [ this , status , xr ] = _scaleinboundsandcons ( this , xr , xbar ); 167 if ( ~status ) then 168 status = "impossibleconstr" 169 break 170 end 171 end 172 [ this.optbase , fr , cr , index ] = optimbase_function ( this.optbase , xr , 2 ); 173 if ( verbose == 1 ) then 174 this = neldermead_log (this,sprintf("xr=[%s], f(xr)=%f", _strvec(xr) , fr)); 175 end 176 if ( fr >= flow & fr < fn ) then 177 if ( verbose == 1 ) then 178 this = neldermead_log (this,sprintf(" > Perform reflection")); 179 end 180 simplex = optimsimplex_setve ( simplex , ihigh , fr , xr ) 181 step = "reflection"; 182 elseif ( fr < flow ) then 183 // Expand 184 if ( verbose == 1 ) then 185 this = neldermead_log (this,sprintf("Expand")); 186 end 187 xe = neldermead_interpolate ( xbar , xhigh , this.rho*this.chi ); 188 // Adjust point to satisfy bounds and nonlinear inequality constraints 189 if ( hasbounds | nbnlc > 0 ) then 190 [ this , status , xe ] = _scaleinboundsandcons ( this , xe , xbar ); 191 if ( ~status ) then 192 status = "impossibleconstr" 193 break 194 end 195 end 196 [ this.optbase , fe , ce , index ] = optimbase_function ( this.optbase , xe , 2 ); 197 if ( verbose == 1 ) then 198 this = neldermead_log (this,sprintf("xe=[%s], f(xe)=%f", strcat(string(xe)," ") , fe )); 199 end 200 if (fe < fr) then 201 if ( verbose == 1 ) then 202 this = neldermead_log (this,sprintf(" > Perform Expansion")); 203 end 204 simplex = optimsimplex_setve ( simplex , ihigh , fe , xe ) 205 step = "expansion"; 206 else 207 if ( verbose == 1 ) then 208 this = neldermead_log (this,sprintf(" > Perform reflection")); 209 end 210 simplex = optimsimplex_setve ( simplex , ihigh , fr , xr ) 211 step = "reflection"; 212 end 213 elseif ( fr >= fn & fr < fhigh ) then 214 // Outside contraction 215 if ( verbose == 1 ) then 216 this = neldermead_log (this,sprintf("Contract - outside")); 217 end 218 xc = neldermead_interpolate ( xbar , xhigh , this.rho*this.gamma ); 219 // Adjust point to satisfy bounds and nonlinear inequality constraints 220 if ( hasbounds | nbnlc > 0 ) then 221 [ this , status , xc ] = _scaleinboundsandcons ( this , xc , xbar ); 222 if ( ~status ) then 223 status = "impossibleconstr" 224 break 225 end 226 end 227 [ this.optbase , fc , cc , index ] = optimbase_function ( this.optbase , xc , 2 ); 228 if ( verbose == 1 ) then 229 this = neldermead_log (this,sprintf("xc=[%s], f(xc)=%f", strcat(string(xc)," ") , fc)); 230 end 231 if ( fc <= fr ) then 232 if ( verbose == 1 ) then 233 this = neldermead_log (this,sprintf(" > Perform Outside Contraction")); 234 end 235 simplex = optimsimplex_setve ( simplex , ihigh , fc , xc ) 236 step = "outsidecontraction"; 237 else 238 // Shrink 239 if ( verbose == 1 ) then 240 this = neldermead_log (this,sprintf(" > Perform Shrink")); 241 end 242 [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this ); 243 step = "shrink"; 244 end 245 else 246 // ( fr >= fn & fr >= fhigh ) 247 // Inside contraction 248 if ( verbose == 1 ) then 249 this = neldermead_log (this,sprintf("Contract - inside")); 250 end 251 xc = neldermead_interpolate ( xbar , xhigh , -this.gamma ); 252 // Adjust point to satisfy bounds and nonlinear inequality constraints 253 if ( hasbounds | nbnlc > 0 ) then 254 [ this , status , xc ] = _scaleinboundsandcons ( this , xc , xbar ); 255 if ( ~status ) then 256 status = "impossibleconstr" 257 break 258 end 259 end 260 [ this.optbase , fc , cc , index ] = optimbase_function ( this.optbase , xc , 2 ); 261 if ( verbose == 1 ) then 262 this = neldermead_log (this,sprintf("xc=[%s], f(xc)=%f", strcat(string(xc)," ") , fc)); 263 end 264 if ( fc < fhigh ) then 265 if ( verbose == 1 ) then 266 this = neldermead_log (this,sprintf(" > Perform Inside Contraction")); 267 end 268 simplex = optimsimplex_setve ( simplex , ihigh , fc , xc ) 269 step = "insidecontraction"; 270 else 271 // Shrink 272 if ( verbose == 1 ) then 273 this = neldermead_log (this,sprintf(" > Perform Shrink")); 274 end 275 [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this ) 276 step = "shrink"; 277 end 278 end 279 // 280 // Sort simplex 281 // 282 if ( verbose == 1 ) then 283 this = neldermead_log (this,sprintf("Sort")); 284 end 285 simplex = optimsimplex_sort ( simplex ); 286 end 287 this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow.' ); 288 this.optbase = optimbase_set ( this.optbase , "-fopt" , flow ); 289 this.optbase = optimbase_set ( this.optbase , "-status" , status ); 290 this.simplexopt = simplex; 291endfunction 292 // 293 // _scaleinboundsandcons -- 294 // Given a point to scale and a reference point which satisfies the constraints, 295 // scale the point towards the reference point until it satisfies all the constraints, 296 // including boun constraints. 297 // Returns isscaled = %T if the procedure has succeded before -boxnbnlloops 298 // Returns isscaled = %F if the procedure has failed after -boxnbnlloops 299 // iterations. 300 // Arguments 301 // x : the point to scale 302 // xref : the reference point 303 // isscaled : %T or %F 304 // p : scaled point 305 // 306function [ this , isscaled , p ] = _scaleinboundsandcons ( this , x , xref ) 307 p = x 308 [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase ); 309 nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" ) 310 // 311 // 1. No bounds, no nonlinear inequality constraints 312 // => no problem 313 // 314 if ( ( hasbounds == %f ) & ( nbnlc == 0 ) ) then 315 isscaled = %T 316 return; 317 end 318 isscaled = %F 319 // 320 // 2. Scale into bounds 321 // 322 if ( hasbounds ) then 323 [ this.optbase , p ] = optimbase_proj2bnds ( this.optbase , p ); 324 this = neldermead_log (this,sprintf(" > After projection into bounds p = [%s]" , ... 325 _strvec(p))); 326 end 327 // 328 // 2. Scale into nonlinear constraints 329 // Try the current point and see if the constraints are satisfied. 330 // If not, move the point "halfway" to the centroid, 331 // which should satisfy the constraints, if 332 // the constraints are convex. 333 // Perform this loop until the constraints are satisfied. 334 // If all loops have been performed without success, the scaling 335 // has failed. 336 // 337 alpha = 1.0 338 p0 = p 339 while ( alpha > this.guinalphamin ) 340 [ this.optbase , feasible ] = optimbase_isinnonlincons ( this.optbase , p ); 341 if ( feasible ) then 342 isscaled = %T; 343 break; 344 end 345 alpha = alpha / 2.0 346 this = neldermead_log (this,sprintf("Scaling inequality constraint with alpha = %e", ... 347 alpha )); 348 p = ( 1.0 - alpha ) * xref + alpha * p0; 349 end 350 this = neldermead_log (this,sprintf(" > After scaling into inequality constraints p = [%s]" , ... 351 _strvec(p) ) ); 352 if ( ~isscaled ) then 353 this = neldermead_log (this,sprintf(" > Impossible to scale into constraints." )); 354 end 355endfunction 356// 357// Test with my own algorithm, 358// the "Mega Super Ultra Modified Simplex Method" !!! 359// 360nm = neldermead_new (); 361nm = neldermead_configure(nm,"-numberofvariables",4); 362nm = neldermead_configure(nm,"-function",optimtestcase); 363nm = neldermead_configure(nm,"-x0",[0.0 0.0 0.0 0.0]'); 364nm = neldermead_configure(nm,"-maxiter",200); 365nm = neldermead_configure(nm,"-maxfunevals",400); 366nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-3); 367nm = neldermead_configure(nm,"-simplex0method","axes"); 368nm = neldermead_configure(nm,"-nbineqconst",3); 369nm = neldermead_configure(nm,"-method","mine"); 370nm = neldermead_configure(nm,"-mymethod",neldermead_constraints); 371nm = neldermead_search(nm); 372// Check optimum point 373xopt = neldermead_get(nm,"-xopt"); 374assert_checkalmostequal ( xopt , [0.0 1.0 2.0 -1.0]', 1e-3, 1e-3 ); 375// Check optimum point value 376fopt = neldermead_get(nm,"-fopt"); 377assert_checkalmostequal ( fopt , -44.0 , 1e-5 ); 378// Check status 379status = neldermead_get(nm,"-status"); 380assert_checkequal ( status , "tolsize" ); 381nm = neldermead_destroy(nm); 382