home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / FUNFUN.DI$ / FMINS.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  4.3 KB  |  168 lines
  1. function [x, options] = fmins(funfcn,x,options,grad,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10)
  2. %FMINS    Minimize a function of several variables.
  3. %    FMINS('F',X0) attempts to return a vector x which is a local minimizer 
  4. %    of F(x) near the starting vector X0.  'F' is a string containing the
  5. %    name of the objective function to be minimized.  F(x) should be a
  6. %    scalar valued function of a vector variable.
  7. %
  8. %    FMINS('F',X0,OPTIONS) uses a vector of control parameters.
  9. %    If OPTIONS(1) is nonzero, intermediate steps in the solution are
  10. %    displayed; the default is OPTIONS(1) = 0.  OPTIONS(2) is the termination
  11. %    tolerance for x; the default is 1.e-4.  OPTIONS(3) is the termination
  12. %    tolerance for F(x); the default is 1.e-4.  OPTIONS(14) is the maximum
  13. %    number of steps; the default is OPTIONS(14) = 500.  The other components
  14. %    of OPTIONS are not used as input control parameters by FMIN.  For more
  15. %    information, see FOPTIONS.
  16. %
  17. %    FMINS('F',X0,OPTIONS,[],P1,P2,...) provides for up to 10 additional
  18. %    arguments which are passed to the objective function, F(X,P1,P2,...)
  19. %
  20. %    FMINS uses a Simplex search method.
  21. %
  22. %    See also FMIN. 
  23.  
  24. %    Reference: J. E. Dennis, Jr. and D. J. Woods, New Computing
  25. %    Environments: Microcomputers in Large-Scale Computing,
  26. %    edited by A. Wouk, SIAM, 1987, pp. 116-122.
  27.  
  28. %    C. Moler, 8-19-86
  29. %    Revised Andy Grace, 6-22-90, 1-17-92 CBM.
  30. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  31.  
  32. if nargin<3, options = []; end
  33. options = foptions(options);
  34. prnt = options(1);
  35. tol = options(2);
  36. tol2 = options(3);
  37. % The input argument grad is there for compatability with FMINU in
  38. % the Optimization Toolbox, but is not used by this function.
  39.  
  40. evalstr = [funfcn];
  41. if ~any(funfcn<48)
  42.     evalstr=[evalstr, '(x'];
  43.     for i=1:nargin - 4
  44.         evalstr = [evalstr,',P',int2str(i)];
  45.     end
  46.     evalstr = [evalstr, ')'];
  47. end
  48.  
  49. n = prod(size(x));
  50. if (~options(14)) 
  51.     options(14) = 200*n; 
  52. end
  53.  
  54. % Set up a simplex near the initial guess.
  55. xin = x(:);
  56. v = 0.9*xin;
  57. x(:) = v; fv = eval(evalstr); 
  58. for j = 1:n
  59.    y = xin;
  60.    if y(j) ~= 0
  61.       y(j) = 1.1*y(j);
  62.    else
  63.       y(j) = 0.1;
  64.    end
  65.    v = [v y];
  66.    x(:) = y; f = eval(evalstr);
  67.    fv = [fv  f];
  68. end
  69. [fv,j] = sort(fv);
  70. v = v(:,j);
  71.  
  72. cnt = n+1;
  73. if prnt
  74.    clc
  75.    format compact
  76.    format short e
  77.    home
  78.    cnt
  79.    disp('initial ')
  80.    disp(' ')
  81.    v
  82.    f
  83. end
  84.  
  85. alpha = 1;  beta = 1/2;  gamma = 2;
  86. [n,np1] = size(v);
  87. onesn = ones(1,n); 
  88. ot = 2:n+1;
  89. on = 1:n;
  90.  
  91. % Iterate until the diameter of the simplex is less than tol.
  92. while cnt < options(14)
  93.     if max(max(abs(v(:,ot)-v(:,onesn)))) <= tol & ...
  94.            max(abs(fv(1)-fv(ot))) <= tol2
  95.         break
  96.     end
  97.  
  98.     % One step of the Nelder-Mead simplex algorithm
  99.  
  100.     vbar = (sum(v(:,on)')/n)';
  101.     vr = (1 + alpha)*vbar - alpha*v(:,n+1);
  102.     x(:) = vr;
  103.     fr = eval(evalstr); 
  104.     cnt = cnt + 1; 
  105.     vk = vr;  fk = fr; how = 'reflect ';
  106.     if fr < fv(n)
  107.         if fr < fv(1)
  108.             ve = gamma*vr + (1-gamma)*vbar;
  109.             x(:) = ve;
  110.             fe = eval(evalstr);
  111.             cnt = cnt + 1;
  112.             if fe < fv(1)
  113.                 vk = ve; fk = fe;
  114.                 how = 'expand  ';
  115.             end
  116.         end
  117.     else
  118.         vt = v(:,n+1); ft = fv(n+1);
  119.         if fr < ft
  120.             vt = vr; ft = fr;
  121.         end
  122.         vc = beta*vt + (1-beta)*vbar;
  123.         x(:) = vc;
  124.         fc = eval(evalstr); 
  125.         cnt = cnt + 1;
  126.         if fc < fv(n)
  127.             vk = vc; fk = fc;
  128.             how = 'contract';
  129.         else
  130.             for j = 2:n
  131.                 v(:,j) = (v(:,1) + v(:,j))/2;
  132.                 x(:) = v(:,j);
  133.                 fv(j) = eval(evalstr); 
  134.             end
  135.         cnt = cnt + n-1;
  136.         vk = (v(:,1) + v(:,n+1))/2;
  137.         x(:) = vk;
  138.         fk = eval(evalstr); 
  139.         cnt = cnt + 1;
  140.         how = 'shrink  ';
  141.         end
  142.     end
  143.     v(:,n+1) = vk;
  144.     fv(n+1) = fk;
  145.     [fv,j] = sort(fv);
  146.     v = v(:,j);
  147.  
  148.     if prnt
  149.         home
  150.         cnt
  151.         disp(how)
  152.         disp(' ')
  153.         v
  154.         fv
  155.     end
  156. end
  157. x(:) = v(:,1);
  158. if prnt, format, end
  159. options(10)=cnt;
  160. options(8)=min(fv); 
  161. if cnt==options(14) 
  162.     if options(1) >= 0
  163.         disp(['Warning: Maximum number of iterations (', ...
  164.                int2str(options(14)),') has been exceeded']);
  165.         disp( '         (increase OPTIONS(14)).')
  166.     end
  167. end
  168.