home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 4.ddi / OPTIM.DI$ / FMINU.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  7.2 KB  |  251 lines
  1. function [x,OPTIONS] = fminu(FUN,x,OPTIONS,GRADFUN,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10)
  2. %FMINU    Finds the minimum of a function of several variables.
  3. %    X=FMINU('FUN',X0) starts at the matrix X0 and finds a minimum to the
  4. %    function which is described in FUN (usually an M-file: FUN.M).
  5. %    The function 'FUN' should return a scalar function value: F=FUN(X).
  6. %
  7. %    X=FMINU('FUN',X0,OPTIONS) allows a vector of optional parameters to
  8. %    be defined. OPTIONS(1) controls how much display output is given; set 
  9. %    to 1 for a tabular display of results, (default is no display: 0). 
  10. %    OPTIONS(2) is a measure of the precision required for the values of 
  11. %    X at the solution. OPTIONS(3) is a measure of the precision
  12. %    required of the objective function at the solution.
  13. %    For more information type HELP FOPTIONS. 
  14. %
  15. %    X=FMINU('FUN',X0,OPTIONS,'GRADFUN') enables a function'GRADFUN'
  16. %    to be entered which returns the partial derivatives of the function,
  17. %    df/dX, at the point X: gf = GRADFUN(X).
  18. %
  19. %    The default algorithm is the BFGS Quasi-Newton method with a 
  20. %    mixed quadratic and cubic line search procedure. 
  21.  
  22. %    Copyright (c) 1990 by the MathWorks, Inc.
  23. %    Andy Grace 7-9-90.
  24.  
  25.  
  26. % ------------Initialization----------------
  27. XOUT=x(:);
  28. nvars=length(XOUT);
  29.  
  30. evalstr = [FUN];
  31. if ~any(FUN<48)
  32.     evalstr=[evalstr, '(x'];
  33.     for i=1:nargin - 4
  34.         evalstr = [evalstr,',P',int2str(i)];
  35.     end
  36.     evalstr = [evalstr, ')'];
  37. end
  38.  
  39. if nargin < 3, OPTIONS=[]; end
  40. if nargin < 4, GRADFUN=[]; end
  41.  
  42.  
  43. if length(GRADFUN)
  44.     evalstr2 = [GRADFUN];
  45.     if ~any(GRADFUN<48) 
  46.         evalstr2 = [evalstr2, '(x'];
  47.         for i=1:nargin - 4
  48.             evalstr2 = [evalstr2,',P',int2str(i)];
  49.         end
  50.         evalstr2 = [evalstr2, ')'];
  51.     end
  52. end
  53.  
  54. f = eval(evalstr); 
  55. n = length(XOUT);
  56. GRAD=zeros(nvars,1);
  57. OLDX=XOUT;
  58. MATX=zeros(3,1);
  59. MATL=[f;0;0];
  60. OLDF=f;
  61. FIRSTF=f;
  62. [OLDX,OLDF,HESS,OPTIONS]=optint(XOUT,f,OPTIONS);
  63. CHG = 1e-7*abs(XOUT)+1e-7*ones(nvars,1);
  64. SD = zeros(nvars,1);
  65. diff = zeros(nvars,1);
  66. PCNT = 0;
  67.  
  68.  
  69. OPTIONS(10)=2; % Iteration count (add 1 for last evaluation)
  70. status =-1;
  71.  
  72. while status ~= 1
  73. % Work Out Gradients
  74.     if ~length(GRADFUN) | OPTIONS(9) 
  75.         OLDF=f;
  76. % Finite difference perturbation levels
  77. % First check perturbation level is not less than search direction.
  78.         f = find(10*abs(CHG)>abs(SD)); 
  79.         CHG(f) = -0.1*SD(f);
  80. % Ensure within user-defined limits
  81.         CHG = sign(CHG+eps).*min(max(abs(CHG),OPTIONS(16)),OPTIONS(17));
  82.         for gcnt=1:nvars
  83.             XOUT(gcnt,1)=XOUT(gcnt)+CHG(gcnt);
  84.             x(:) = XOUT; f = eval(evalstr); 
  85.             GRAD(gcnt)=(f-OLDF)/(CHG(gcnt));
  86.             if f < OLDF
  87.                 OLDF=f;
  88.             else
  89.                  XOUT(gcnt)=XOUT(gcnt)-CHG(gcnt);
  90.             end
  91.         end
  92. % Try to set difference to 1e-8 for next iteration
  93. % Add eps for machines that can't handle divide by zero.
  94.         CHG = 1e-8./(GRAD + eps); 
  95.         f = OLDF;
  96.         OPTIONS(10)=OPTIONS(10)+nvars;
  97. % Gradient check 
  98.         if OPTIONS(9) == 1 
  99.             GRADFD = GRAD; 
  100.             x(:)=XOUT;  GRAD(:) = eval(evalstr2); 
  101.             graderr(GRADFD, GRAD, evalstr2);
  102.             OPTIONS(9) = 0; 
  103.         end
  104.                 
  105.     else
  106.         OPTIONS(11)=OPTIONS(11)+1;
  107.         x(:)=XOUT; GRAD(:) = eval(evalstr2);
  108.     end
  109. %---------------Initialization of Search Direction------------------
  110. if status == -1
  111.     SD=-GRAD;
  112.     FIRSTF=f;
  113.     OLDG=GRAD;
  114.     GDOLD=GRAD'*SD;
  115. % For initial step-size guess assume the minimum is at zero. 
  116.     OPTIONS(18) = max(0.001, min([1,2*abs(f/GDOLD)]));
  117.     if OPTIONS(1)>0
  118.         disp([sprintf('%5.0f %12.6g %12.6g ',OPTIONS(10),f,OPTIONS(18)),sprintf('%12.3g  ',GDOLD)]);
  119.     end
  120.     XOUT=XOUT+OPTIONS(18)*SD;
  121.     status=4; 
  122.     if OPTIONS(7)==0; PCNT=1; end
  123.          
  124. else
  125. %-------------Direction Update------------------
  126.     gdnew=GRAD'*SD;
  127.     if OPTIONS(1)>0, 
  128.         num=[sprintf('%5.0f %12.6g %12.6g ',OPTIONS(10),f,OPTIONS(18)),sprintf('%12.3g  ',gdnew)];
  129.     end
  130.     if (gdnew>0 & f>FIRSTF)|~finite(f)
  131. % Case 1: New function is bigger than last and gradient w.r.t. SD -ve
  132. %    ...interpolate.
  133.         how='inter';
  134.         [stepsize]=cubici1(f,FIRSTF,gdnew,GDOLD,OPTIONS(18));
  135.         if stepsize<0|isnan(stepsize), stepsize=OPTIONS(18)/2; how='C1f '; end
  136.         if OPTIONS(18)<0.1&OPTIONS(6)==0 
  137.             if stepsize*norm(SD)<eps
  138.                 stepsize=exp(rand(1,1)-1)-0.1;
  139.                 how='RANDOM STEPLENGTH';
  140.                 status=0;
  141.             else        
  142.                 stepsize=stepsize/2;
  143.             end   
  144.         end      
  145.         OPTIONS(18)=stepsize;
  146.         XOUT=OLDX;
  147.     elseif f<FIRSTF
  148.         [newstep,fbest] =cubici3(f,FIRSTF,gdnew,GDOLD,OPTIONS(18));
  149.         sk=(XOUT-OLDX)'*(GRAD-OLDG);
  150.         if sk>1e-20
  151. % Case 2: New function less than old fun. and OK for updating HESS
  152. %         .... update and calculate new direction.
  153.         how='';   
  154.             if gdnew<0
  155.                 how='incstep';
  156.                 if newstep<OPTIONS(18),  newstep=2*OPTIONS(18)+1e-5; how=[how,' IF']; end
  157.                 OPTIONS(18)=min([max([2,1.5*OPTIONS(18)]),1+sk+abs(gdnew)+max([0,OPTIONS(18)-1]), (1.2+0.3*(~OPTIONS(7)))*abs(newstep)]);
  158.             else % gdnew>0
  159.                 if OPTIONS(18)>0.9
  160.                     how='int_st';
  161.                     OPTIONS(18)=min([1,abs(newstep)]);
  162.                 end
  163.             end %if gdnew
  164.             [HESS,SD]=updhess(XOUT,OLDX,GRAD,OLDG,HESS,OPTIONS);
  165.             gdnew=GRAD'*SD;
  166.             OLDX=XOUT;
  167.             status=4;
  168. % Save Variables for next update
  169.             FIRSTF=f;
  170.             OLDG=GRAD;
  171.             GDOLD=gdnew;
  172. % If mixed interpolation set PCNT
  173.             if OPTIONS(7)==0, PCNT=1; MATX=zeros(3,1);  MATL(1)=f; end
  174.     elseif gdnew>0 %sk<=0 
  175. % Case 3: No good for updating HESSIAN .. interpolate or halve step length.
  176.             how='inter_st'; 
  177.             if OPTIONS(18)>0.01
  178.                 OPTIONS(18)=0.9*newstep;
  179.                 XOUT=OLDX;
  180.             end
  181.             if OPTIONS(18)>1, OPTIONS(18)=1; end
  182.         else  
  183. % Increase step, replace starting point
  184.             OPTIONS(18)=max([min([newstep-OPTIONS(18),3]),0.5*OPTIONS(18)]);
  185.             how='incst2';
  186.             OLDX=XOUT;
  187.             FIRSTF=f;
  188.             OLDG=GRAD;
  189.             GDOLD=GRAD'*SD;
  190.                  OLDX=XOUT;
  191.         end % if sk>
  192. % Case 4: New function bigger than old but gradient in on
  193. %         ...reduce step length.
  194.     else %gdnew<0 & F>FIRSTF
  195.         if gdnew<0&f>FIRSTF
  196.             how='red_step';  
  197.             if norm(GRAD-OLDG)<1e-10; HESS=eye(nvars); end
  198.             if abs(OPTIONS(18))<eps
  199.                 SD=norm(nvars,1)*(rand(nvars,1)-0.5)
  200.                 OPTIONS(18)=abs(rand(1,1)-0.5)*1e-6;
  201.                           how='RANDOM SD';
  202.                      else
  203.                           OPTIONS(18)=-OPTIONS(18)/2;
  204.             end
  205.             XOUT=OLDX;
  206.         end %gdnew>0    
  207.     end % if (gdnew>0 & F>FIRSTF)|~finite(F)
  208.     XOUT=XOUT+OPTIONS(18)*SD;
  209.     if OPTIONS(1)>0, disp([num,how]),end
  210. end %----------End of Direction Update-------------------
  211.  
  212. % Check Termination 
  213.     if max(abs(SD))<2*OPTIONS(2) & (-GRAD'*SD) < 2*OPTIONS(3)
  214.         if OPTIONS(1) > 0
  215.             disp('Optimization Terminated Successfully')
  216.             disp('Gradient less than options(2)')     
  217.             disp([' NO OF ITERATIONS=', int2str(OPTIONS(10))]);
  218.         end
  219.         status=1; 
  220.     elseif OPTIONS(10)>OPTIONS(14) 
  221.             if OPTIONS(1)>=0
  222.                 disp('Warning: Maximum number of iterations has been exceeded');
  223.                 disp('       - increase options(14) for more iterations.')
  224.             end
  225.             status=1;
  226.     else
  227.  
  228. % Line search using mixed polynomial interpolation and extrapolation.
  229.         if PCNT~=0 
  230.             while PCNT > 0
  231.                 x(:) = XOUT; f = eval(evalstr); OPTIONS(10)=OPTIONS(10)+1;
  232.                 [PCNT,MATL,MATX,steplen,f, how]=searchq(PCNT,f,OLDX,MATL,MATX,SD,GDOLD,OPTIONS(18), how);
  233.                 OPTIONS(18)=steplen;
  234.                 XOUT=OLDX+steplen*SD;
  235.             end
  236.  
  237.         else
  238.             x(:)=XOUT; f = eval(evalstr); OPTIONS(10)=OPTIONS(10)+1;
  239.         end
  240.     end
  241. end
  242.  
  243. x(:)=XOUT; 
  244. f = eval(evalstr);
  245. if f > FIRSTF
  246.     OPTIONS(8) = FIRSTF; 
  247.     x(:)=OLDX;
  248. else
  249.     OPTIONS(8) = f; 
  250. end
  251.