home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 4.ddi / OPTIM.DI$ / CONSTR.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  10.9 KB  |  422 lines
  1. function [x, OPTIONS,lambda, HESS]=constr(FUN,x,OPTIONS,VLB,VUB,GRADFUN,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15)
  2. %CONSTR    Finds the constrained minimum of a function of several variables.
  3. %
  4. %    X=CONSTR('FUN',X0) starts at X0 and finds a constrained minimum to 
  5. %    the function which is described in FUN (usually an M-file: FUN.M).
  6. %    The function 'FUN' should return two arguments: a scalar value of the 
  7. %    function to be minimized, F, and a matrix of constraints, G: 
  8. %    [F,G]=FUN(X). F is minimized such that G < zeros(G).
  9. %
  10. %    X=CONSTR('FUN',X,OPTIONS) allows a vector of optional parameters to 
  11. %    be defined. For more information type HELP FOPTIONS.
  12. %    
  13. %    X=CONSTR('FUN',X,OPTIONS,VLB,VUB) defines a set of lower and upper
  14. %    bounds on the design variables, X, so that the solution is always in 
  15. %    the range VLB < X < VUB. 
  16. %    
  17. %    X=CONSTR('FUN',X,OPTIONS,VLB,VUB,'GRADFUN') allows a function 
  18. %    'GRADFUN' to be entered which returns the partial derivatives of the 
  19. %    function and the  constraints at X:  [gf,GC] = GRADFUN(X).
  20.  
  21. %    Copyright (c) 1990 by the MathWorks, Inc.
  22. %    Andy Grace 7-9-90.
  23.  
  24. %
  25. %    X=CONSTR('FUN',X,OPTIONS,VLB,VUB,GRADFUN,P1,P2,..) allows
  26. %    coefficients, P1, P2, ... to be passed directly to FUN:
  27. %    [F,G]=FUN(X,P1,P2,...). Empty arguments ([]) are ignored.
  28.  
  29.  
  30. global OPT_STOP OPT_STEP;
  31. OPT_STEP = 0; 
  32. OPT_STOP = 0; 
  33.  
  34. % Set up parameters.
  35. XOUT(:)=x;
  36.  
  37. if ~any(FUN<48) % Check alphanumeric
  38.     etype = 1;
  39.     evalstr = [FUN,];
  40.         evalstr=[evalstr, '(x'];
  41.         for i=1:nargin - 6
  42.         etype = 2;
  43.                    evalstr = [evalstr,',P',int2str(i)];
  44.         end
  45.         evalstr = [evalstr, ')'];
  46. else
  47.     etype = 3;
  48.         evalstr=[FUN,'; g=g(:);'];
  49. end
  50.  
  51. if nargin < 3, OPTIONS=[]; end
  52. if nargin < 4, VLB=[]; end
  53. if nargin < 5, VUB=[]; end
  54. if nargin < 6, GRADFUN=[]; end
  55.  
  56. VLB=VLB(:); lenvlb=length(VLB);
  57. VUB=VUB(:); lenvub=length(VUB);
  58. bestf = Inf; 
  59.  
  60. nvars = length(XOUT);
  61.  
  62. CHG = 1e-7*abs(XOUT)+1e-7*ones(nvars,1);
  63. if lenvlb*lenvlb>0
  64.       if any(VLB(1:lenvub)>VUB), error('Bounds Infeasible'), end
  65. end
  66. for i=1:lenvlb
  67.        if lenvlb>0,if XOUT(i)<VLB(i),XOUT(i)=VLB(i)+1e-4; end,end
  68. end
  69. for i=1:lenvub
  70.        if lenvub>0,if XOUT(i)>VUB(i),XOUT(i)=VUB(i);CHG(i)=-CHG(i);end,end
  71. end
  72.  
  73. % Used for semi-infinite optimization:
  74. s = nan; POINT =[]; NEWLAMBDA =[]; LAMBDA = []; NPOINT =[]; FLAG = 2;
  75.  
  76. x(:) = XOUT; 
  77.  
  78.  
  79. if etype == 1, 
  80.     [f, g(:)] = feval(FUN,x);
  81. elseif etype == 2
  82.     [f, g(:)] = eval(evalstr);
  83. else 
  84.     eval(evalstr);
  85. end
  86.  
  87. ncstr = length(g);
  88. if ncstr == 0   
  89.     g = -1; 
  90.     ncstr = 1; 
  91.     if etype ~= 3 
  92.         evalstr = ['[f,g] =', evalstr, ';'];
  93.         etype = 3;
  94.     end
  95.     evalstr = [evalstr,'g=-1; ']; 
  96. end
  97.  
  98.  
  99. if length(GRADFUN)
  100.     if ~any(GRADFUN<48) % Check alphanumeric
  101.         gtype = 1;
  102.         evalstr2 = [GRADFUN,'(x'];
  103.         for i=1:nargin - 6
  104.             gtype = 2;
  105.             evalstr2 = [evalstr2,',P',int2str(i)];
  106.         end
  107.         evalstr2 = [evalstr2, ')'];
  108.     else
  109.         gtype = 3;
  110.         evalstr2=[GRADFUN,';'];
  111.     end
  112. end
  113.  
  114. OLDX=XOUT;
  115. OLDG=g;
  116. OLDgf=zeros(nvars,1);
  117. gf=zeros(nvars,1);
  118. OLDAN=zeros(ncstr,nvars);
  119. LAMBDA=zeros(ncstr,1);
  120. sizep = length(OPTIONS);
  121. OPTIONS = foptions(OPTIONS);
  122. if lenvlb*lenvlb>0
  123.       if any(VLB(1:lenvub)>VUB), error('Bounds Infeasible'), end
  124. end
  125. for i=1:lenvlb
  126.        if lenvlb>0,if XOUT(i)<VLB(i),XOUT(i)=VLB(i)+eps; end,end
  127. end
  128. OPTIONS(18)=1;
  129. if OPTIONS(1)>0
  130.     disp('')
  131.           disp('f-COUNT   FUNCTION       MAX{g}         STEP  Procedures');
  132. end
  133. HESS=eye(nvars,nvars);
  134. if sizep<1 |OPTIONS(14)==0, OPTIONS(14)=nvars*100;end
  135. OPTIONS(10)=1;
  136. OPTIONS(11)=1;
  137. GNEW=1e8*CHG;
  138.  
  139.  
  140.  
  141. %---------------------------------Main Loop-----------------------------
  142. status = 0; 
  143. while status ~= 1
  144.  
  145. %----------------GRADIENTS----------------
  146.  
  147.     if ~length(GRADFUN) | OPTIONS(9)
  148. % Finite Difference gradients
  149.         POINT = NPOINT; 
  150.         oldf = f;
  151.         oldg = g;
  152.         ncstr = length(g);
  153.         FLAG = 0; % For semi-infinite
  154.         gg = zeros(nvars, ncstr);  % For semi-infinite
  155. % Try to make the finite differences equal to 1e-8.
  156.                 CHG = -1e-8./(GNEW+eps);
  157.         CHG = sign(CHG+eps).*min(max(abs(CHG),OPTIONS(16)),OPTIONS(17));
  158.         OPT_STEP = 1;
  159.         for gcnt=1:nvars
  160.             if gcnt == nvars, FLAG = -1; end
  161.             temp = XOUT(gcnt);
  162.             XOUT(gcnt)= temp + CHG(gcnt);
  163.             x(:) =XOUT; 
  164.             if etype == 1, 
  165.             [f, g(:)] = feval(FUN,x);
  166.             elseif etype == 2
  167.                 [f, g(:)] = eval(evalstr);
  168.             else 
  169.                 eval(evalstr);
  170.             end
  171.             OPT_STEP = 0;
  172. % Next line used for problems with varying number of constraints
  173.             if ncstr~=length(g), diff=length(g); g=v2sort(oldg,g); end
  174.             gf(gcnt,1) = (f-oldf)/CHG(gcnt);
  175.             gg(gcnt,:) = (g - oldg)'/CHG(gcnt); 
  176.             XOUT(gcnt) = temp;
  177.         end
  178. % Gradient check
  179.                 if OPTIONS(9) == 1
  180.                         gfFD = gf;
  181.             ggFD = gg; 
  182.                         x(:)=XOUT; 
  183.             if gtype == 1
  184.                 [gf(:), gg] = feval(GRADFUN, x);
  185.             elseif gtype == 2
  186.                 [gf(:), gg] = eval(evalstr2);
  187.             else
  188.                 eval(evalstr2);
  189.             end
  190.             disp('Function derivative')
  191.                         graderr(gfFD, gf, evalstr2);
  192.             disp('Constraint derivative')
  193.                         graderr(ggFD, gg, evalstr2);
  194.                         OPTIONS(9) = 0;
  195.                 end
  196.         FLAG = 1; % For semi-infinite
  197.         OPTIONS(10) = OPTIONS(10) + nvars;
  198.         f=oldf;
  199.         g=oldg;
  200.     else
  201. % User-supplied gradients
  202.         if gtype == 1
  203.             [gf(:), gg] = feval(GRADFUN, x);
  204.         elseif gtype == 2
  205.             [gf(:), gg] = eval(evalstr2);
  206.         else
  207.             eval(evalstr2);
  208.         end
  209.     end
  210.     AN=gg';
  211.     how='';
  212.  
  213. %-------------SEARCH DIRECTION---------------
  214.  
  215.     for i=1:OPTIONS(13) 
  216.         schg=AN(i,:)*gf;
  217.         if schg>0
  218.             AN(i,:)=-AN(i,:);
  219.             g(i)=-g(i);
  220.         end
  221.     end
  222.  
  223.     if OPTIONS(11)>1  % Check for first call    
  224. % For equality constraints make gradient face in 
  225. % opposite direction to function gradient.
  226.         if OPTIONS(7)~=5,     
  227.             NEWLAMBDA=LAMBDA; 
  228.         end
  229.         [ma,na] = size(AN);
  230.         GNEW=gf+AN'*NEWLAMBDA;
  231.         GOLD=OLDgf+OLDAN'*LAMBDA;
  232.         YL=GNEW-GOLD;
  233.         sdiff=XOUT-OLDX;
  234. % Make sure Hessian is positive definite in update.
  235.         if YL'*sdiff<OPTIONS(18)^2*1e-3
  236.             while YL'*sdiff<-1e-5
  237.                 [YMAX,YIND]=min(YL.*sdiff);
  238.                 YL(YIND)=YL(YIND)/2;
  239.             end
  240.             if YL'*sdiff < (eps*norm(HESS,'fro'));
  241.                 how=' mod Hess(2)';
  242.                 FACTOR=AN'*g - OLDAN'*OLDG;
  243.                 FACTOR=FACTOR.*(sdiff.*FACTOR>0).*(YL.*sdiff<=eps);
  244.                 WT=1e-2;
  245.                 if max(abs(FACTOR))==0; FACTOR=1e-5*sign(sdiff); end
  246.                 while YL'*sdiff < (eps*norm(HESS,'fro')) & WT < 1/eps
  247.                     YL=YL+WT*FACTOR;
  248.                     WT=WT*2;
  249.                 end
  250.                 else
  251.                       how=' mod Hess';
  252.             end
  253.          end
  254.  
  255. %----------Perform BFGS Update If YL'S Is Positive---------
  256.         if YL'*sdiff>eps
  257.             HESS=HESS+(YL*YL')/(YL'*sdiff)-(HESS*sdiff*sdiff'*HESS')/(sdiff'*HESS*sdiff);
  258. % BFGS Update using Cholesky factorization  of Gill, Murray and Wright.
  259. % In practice this was less robust than above method and slower.
  260. %    R=chol(HESS); 
  261. %    s2=R*S; y=R'\YL; 
  262. %    W=eye(nvars,nvars)-(s2'*s2)\(s2*s2') + (y'*s2)\(y*y');
  263. %    HESS=R'*W*R;
  264.         else
  265.               how=[how,' (no update)'];
  266.         end
  267.  
  268.     else % First call
  269.             OLDLAMBDA=(eps+gf'*gf)*ones(ncstr,1)./(sum(AN'.*AN')'+eps) ;
  270.     end % if OPTIONS(11)>1
  271.     OPTIONS(11)=OPTIONS(11)+1;
  272.  
  273.     LOLD=LAMBDA;
  274.     OLDAN=AN;
  275.     OLDgf=gf;
  276.     OLDG=g;
  277.     OLDF=f;
  278.     OLDX=XOUT;
  279.     XN=zeros(nvars,1);
  280.     if (OPTIONS(7)>0&OPTIONS(7)<5)
  281.     % Minimax and attgoal problems have special Hessian:
  282.         HESS(nvars,1:nvars)=zeros(1,nvars);
  283.         HESS(1:nvars,nvars)=zeros(nvars,1);
  284.         HESS(nvars,nvars)=1e-8*norm(HESS,'inf');
  285.         XN(nvars)=max(g); % Make a feasible solution for qp
  286.     end
  287.         if lenvlb>0,
  288.                 AN=[AN;-eye(lenvlb,nvars)];
  289.                 GT=[g;-XOUT(1:lenvlb)+VLB];
  290.         else
  291.                 GT=g;
  292.         end
  293.         if lenvub>0
  294.                 AN=[AN;eye(lenvub,nvars)];
  295.                 GT=[GT;XOUT(1:lenvub)-VUB];
  296.         end
  297.     [SD,lambda,howqp]=qp(HESS,gf,AN,-GT, [], [], XN,OPTIONS(13),-1); 
  298.     lambda(1:OPTIONS(13)) = abs(lambda(1:OPTIONS(13)));
  299.     ga=[abs(g(1:OPTIONS(13)));g(OPTIONS(13)+1:ncstr)];
  300.     mg=max(ga);
  301.     if OPTIONS(1)>0
  302.         if howqp(1) == 'o'; howqp = ' '; end
  303.         disp([sprintf('%5.0f %12.6g %12.6g ',OPTIONS(10),f,mg), sprintf('%12.3g  ',OPTIONS(18)),how, ' ',howqp]);
  304.     end
  305.     LAMBDA=lambda(1:ncstr);
  306.     OLDLAMBDA=max([LAMBDA';0.5*(LAMBDA+OLDLAMBDA)'])' ;
  307.  
  308. %---------------LINESEARCH--------------------
  309.     MATX=XOUT;
  310.     MATL = f+sum(OLDLAMBDA.*(ga>0).*ga) + 1e-30;
  311.     infeas = (howqp(1) == 'i');
  312.     if OPTIONS(7)==0 | OPTIONS(7) == 5
  313. % This merit function looks for improvement in either the constraint
  314. % or the objective function unless the sub-problem is infeasible in which
  315. % case only a reduction in the maximum constraint is tolerated.
  316. % This less "stringent" merit function has produced faster convergence in
  317. % a large number of problems.
  318.         if mg > 0
  319.             MATL2 = mg;
  320.         elseif f >=0 
  321.             MATL2 = -1/(f+1);
  322.         else 
  323.             MATL2 = 0;
  324.         end
  325.         if ~infeas & f < 0
  326.             MATL2 = MATL2 + f - 1;
  327.         end
  328.     else
  329. % Merit function used for MINIMAX or ATTGOAL problems.
  330.         MATL2=mg+f;
  331.     end
  332.     if mg < eps & f < bestf
  333.         bestf = f;
  334.         bestx = XOUT;
  335.     end
  336.     MERIT = MATL + 1;
  337.     MERIT2 = MATL2 + 1; 
  338.     OPTIONS(18)=2;
  339.     while (MERIT2 > MATL2) & (MERIT > MATL) & OPTIONS(10) < OPTIONS(14)
  340.         OPTIONS(18)=OPTIONS(18)/2;
  341.         if OPTIONS(18) < 1e-4,  
  342.             OPTIONS(18) = -OPTIONS(18); 
  343.  
  344.         % Semi-infinite may have changing sampling interval
  345.         % so avoid too stringent check for improvement
  346.             if OPTIONS(7) == 5, 
  347.                 OPTIONS(18) = -OPTIONS(18); 
  348.                 MATL2 = MATL2 + 10; 
  349.             end
  350.         end
  351.         XOUT = MATX + OPTIONS(18)*SD;
  352.         x(:)=XOUT; 
  353.         if etype == 1,
  354.             [f, g(:)] = feval(FUN,x);
  355.         elseif etype == 2
  356.              [f, g(:)] = eval(evalstr);
  357.         else
  358.             eval(evalstr);
  359.         end
  360.         OPTIONS(10) = OPTIONS(10) + 1;
  361.         ga=[abs(g(1:OPTIONS(13)));g(OPTIONS(13)+1:length(g))];
  362.         mg=max(ga);
  363.         MERIT = f+sum(OLDLAMBDA.*(ga>0).*ga);
  364.         if OPTIONS(7)==0 | OPTIONS(7) == 5
  365.             if mg > 0
  366.                 MERIT2 = mg;
  367.             elseif f >=0 
  368.                 MERIT2 = -1/(f+1);
  369.             else 
  370.                 MERIT2 = 0;
  371.             end
  372.             if ~infeas & f < 0
  373.                 MERIT2 = MERIT2 + f - 1;
  374.             end
  375.         else
  376.             MERIT2=mg+f;
  377.         end
  378.     end
  379. %------------Finished Line Search-------------
  380.  
  381.     if OPTIONS(7)~=5
  382.         mf=abs(OPTIONS(18));
  383.         LAMBDA=mf*LAMBDA+(1-mf)*LOLD;
  384.     end
  385.     if max(abs(SD))<2*OPTIONS(2) & abs(gf'*SD)<2*OPTIONS(3) & (mg<OPTIONS(4) | (howqp(1) == 'i' & mg > 0 ) )
  386.         if OPTIONS(1)>0
  387.             disp([sprintf('%5.0f %12.6g %12.6g ',OPTIONS(10),f,mg),sprintf('%12.3g  ',OPTIONS(18)),how, ' ',howqp]);
  388.             if howqp(1) ~= 'i'
  389.                  disp('Optimization Terminated Successfully')
  390.                 disp('Active Constraints:'), 
  391.                 find(LAMBDA>0) 
  392.             end
  393.         end
  394.         if (howqp(1) == 'i' & mg > 0)
  395.             disp('Warning: No feasible solution found.')
  396.         end
  397.         status=1;
  398.  
  399.     else
  400.     %    NEED=[LAMBDA>0]|G>0
  401.         if OPTIONS(10) >= OPTIONS(14) | OPT_STOP
  402.             XOUT = MATX;
  403.             f = OLDF;
  404.             if ~OPT_STOP
  405.                 disp('Maximum number of iterations exceeded')
  406.                 disp('increase OPTIONS(14)')
  407.             else
  408.                 disp('Optimization terminated prematurely by user')
  409.             end
  410.             status=1;
  411.         end
  412.     end  
  413. end
  414.  
  415. % If a better unconstrained solution was found earlier, use it:
  416. if f > bestf 
  417.     XOUT = bestx;
  418.     f = bestf;
  419. end
  420. OPTIONS(8)=f;
  421. x(:) = XOUT;
  422.