nlconst.m

matpower软件下

               if verbosity > 0
                   disp(OUTPUT.message)
               end
               return;
           end
       end
       
       %-------------TEST CONVERGENCE---------------
       % If NoStopIfFlatInfeas option is on, in addition to the objective looking
       % flat, also require that the iterate be feasible (among other things) to 
       % detect that no further progress can be made.
       if ~noStopIfFlatInfeas
         noFurtherProgress = ( max(abs(SD)) < 2*tolX || abs(gf'*SD) < 2*tolFun ) && ...
               (mg < tolCon || infeasIllPosedMaxSQPIter);
       else
         noFurtherProgress = ( max(abs(SD)) < 2*tolX || (abs(gf'*SD) < 2*tolFun && ...
               feasError < tolCon*feasScal) ) && ( mg < tolCon || infeasIllPosedMaxSQPIter );
       end
         
       if optimError < tolFun*optimScal && feasError < tolCon*feasScal
           outMessage = ...
             sprintf(['Optimization terminated: first-order optimality measure less\n' ...
                      ' than options.TolFun and maximum constraint violation is less\n' ...
                      ' than options.TolCon.']);
           if verbosity > 0
               disp(outMessage) 
           end
           EXITFLAG = 1;
           done = true;

           if inequalitiesExist
              % Report active inequalities
              [activeLb,activeUb,activeIneqLin,activeIneqNonlin] = ...
                  activeInequalities(c,tolCon,arglb,argub,lin_eq,non_eq,size(Ain));           

              if any([activeLb;activeUb;activeIneqLin;activeIneqNonlin])
                 if verbosity > 0
                    fprintf('Active inequalities (to within options.TolCon = %g):\n',tolCon)
                    disp('  lower      upper     ineqlin   ineqnonlin')
                    printColumnwise(activeLb,activeUb,activeIneqLin,activeIneqNonlin);
                 end
              else
                 if verbosity > 0
                    disp('No active inequalities')
                 end 
              end 
           end   
       elseif noFurtherProgress
           % The algorithm can make no more progress.  If feasible, compute 
           % the new up-to-date Lagrange multipliers (with new gradients) 
           % and recompute the KKT error.  Then output appropriate termination
           % message.
           if mg < tolCon
               if meritFunctionType == 1
                   optimError = inf;
               else
                   lambdaNLP(:,1) = 0;
                   [Q,R] = qr(AN(ACTIND,:)');
                   ws = warning('off');
                   lambdaNLP(ACTIND) = -R\Q'*gf;
                   warning(ws);
                   lambdaNLP(eq+1:ncstr) = max(0,lambdaNLP(eq+1:ncstr));
                   if strcmp(funfcn{2},'fseminf'), lambdaNLP = NEWLAMBDA; end
                   normgradLag = norm(gf + AN'*lambdaNLP,inf);
                   normcomp = norm(lambdaNLP(eq+1:ncstr).*c(eq+1:ncstr),inf);
                   if isfinite(normgradLag) && isfinite(normcomp)
                       optimError = max(normgradLag, normcomp);
                   else
                       optimError = inf;
                   end
               end
               optimScal = 1;
               if optimError < tolFun*optimScal
                   outMessage = ...
                     sprintf(['Optimization terminated: first-order optimality ' ...
                              'measure less than options.TolFun\n and maximum ' ...
                              'constraint violation is less than options.TolCon.']);
                   if verbosity > 0
                       disp(outMessage)
                   end
                   EXITFLAG = 1;
               elseif max(abs(SD)) < 2*tolX
                   outMessage = ...
                     sprintf(['Optimization terminated: Search direction less than 2*options.TolX\n' ...
                              ' and maximum constraint violation is less than options.TolCon.']);
                   if verbosity > 0
                       disp(outMessage)
                   end
                   EXITFLAG = 4;
               else 
                   outMessage = ...
                      sprintf(['Optimization terminated: Magnitude of directional derivative in search\n' ... 
                               ' direction less than 2*options.TolFun and maximum constraint violation\n' ...
                               '  is less than options.TolCon.']);
                   if verbosity > 0 
                       disp(outMessage) 
                   end
                   EXITFLAG = 5;
               end 
               
               if inequalitiesExist
                  % Report active inequalities
                  [activeLb,activeUb,activeIneqLin,activeIneqNonlin] = ...
                      activeInequalities(c,tolCon,arglb,argub,lin_eq,non_eq,size(Ain));  

                  if any([activeLb;activeUb;activeIneqLin;activeIneqNonlin])
                     if verbosity > 0
                        fprintf('Active inequalities (to within options.TolCon = %g):\n', tolCon)
                        disp('  lower      upper     ineqlin   ineqnonlin')
                        printColumnwise(activeLb,activeUb,activeIneqLin,activeIneqNonlin);
                     end
                  else
                     if verbosity > 0
                        disp('No active inequalities')
                     end
                  end
               end
           else                         % if mg >= tolCon
               if max(abs(SD)) < 2*tolX
                   outMessage = sprintf(['Optimization terminated: No feasible solution found.\n', ...
                      ' Search direction less than 2*options.TolX but constraints are not satisfied.']);
               else
                   outMessage = sprintf(['Optimization terminated: No feasible solution found.\n' ...
                                        '  Magnitude of directional derivative in search direction\n', ...
                                        '  less than 2*options.TolFun but constraints are not satisfied.']);
               end
               if strcmp(howqp,'MaxSQPIter')
                   outMessage = sprintf(['Optimization terminated: No feasible solution found.\n' ...
                           ' During the solution to the last quadratic programming subproblem, the\n' ...
                           ' maximum number of iterations was reached. Increase options.MaxSQPIter.']);
               end 
               EXITFLAG = -2;
               if verbosity > 0
                 disp(outMessage)
               end
           end                          % of "if mg < tolCon"
           done = true;
       else % continue
           % NEED=[LAMBDA>0] | G>0
           if numFunEvals > maxFunEvals
               XOUT = MATX;
               f = OLDF;
               gf = OLDgf;
               outMessage = sprintf(['Maximum number of function evaluations exceeded;\n' ...
                                         ' increase OPTIONS.MaxFunEvals.']);
               if verbosity > 0
                   disp(outMessage)
               end
               EXITFLAG = 0;
               done = true;
           end
           if iter > maxIter
               XOUT = MATX;
               f = OLDF;
               gf = OLDgf;
               outMessage = sprintf(['Maximum number of iterations exceeded;\n' ...
                                         ' increase OPTIONS.MaxIter.']);
               if verbosity > 0
                   disp(outMessage)
               end
               EXITFLAG = 0;
               done = true;
           end
       end 
   else % ------------------------0th Iteration----------------------------------
       if verbosity > 1
           disp(header)
           % Print 0th iteration information (some columns left blank)
           if meritFunctionType == 1,
               gamma = mg+f;
               CurrOutput = sprintf(formatstrFirstIter,iter,numFunEvals,gamma,how); 
               disp(CurrOutput)
           else
               if mg > 0
                   how = 'Infeasible start point';
               else
                   how = '';
               end
               CurrOutput = sprintf(formatstrFirstIter,iter,numFunEvals,f,mg,how); 
               disp(CurrOutput)
           end
       end
       
       % Initialize the output function.
       if haveoutputfcn
           [xOutputfcn, optimValues, stop] = callOutputFcn(outputfcn,XOUT,xOutputfcn,'init', iter,numFunEvals, ...
               f,mg,[],gf,[],meritFunctionType,[],[],[],outputfcnVarargin{:});
           if stop
               [x,FVAL,lambda_out,EXITFLAG,OUTPUT,GRADIENT,HESS] = cleanUpInterrupt(xOutputfcn,optimValues);
               if verbosity > 0
                   disp(OUTPUT.message)
               end
               return;
           end
           
           % OutputFcn call for 0th iteration
           [xOutputfcn, optimValues, stop] = callOutputFcn(outputfcn,XOUT,xOutputfcn,'iter',iter,numFunEvals, ...
               f,mg,[],gf,[],meritFunctionType,[],how,'',outputfcnVarargin{:});
           if stop  % Stop per user request.
               [x,FVAL,lambda_out,EXITFLAG,OUTPUT,GRADIENT,HESS] = ...
                   cleanUpInterrupt(xOutputfcn,optimValues);
               if verbosity > 0
                   disp(OUTPUT.message)
               end
               return;
           end
           
       end % if haveoutputfcn
   end % if iter > 0
   
   % Continue if termination criteria do not hold or it is the 0th iteration-------------------------------------------
   if ~done 
      how=''; 
      iter = iter + 1;

      %-------------SEARCH DIRECTION---------------
      % For equality constraints make gradient face in 
      % opposite direction to function gradient.
      for i=1:eq 
         schg=AN(i,:)*gf;
         if schg>0
            AN(i,:)=-AN(i,:);
            c(i)=-c(i);
         end
      end
   
      if numGradEvals>1  % Check for first call    
         if meritFunctionType~=5,   
            NEWLAMBDA=LAMBDA; 
         end
         [ma,na] = size(AN);
         GNEW=gf+AN'*NEWLAMBDA;
         GOLD=OLDgf+OLDAN'*LAMBDA;
         YL=GNEW-GOLD;
         sdiff=XOUT-OLDX;

         % Make sure Hessian is positive definite in update.
         if YL'*sdiff<stepsize^2*1e-3
            while YL'*sdiff<-1e-5
               [YMAX,YIND]=min(YL.*sdiff);
               YL(YIND)=YL(YIND)/2;
            end
            if YL'*sdiff < (eps*norm(HESS,'fro'));
               how=' Hessian modified twice';
               FACTOR=AN'*c - OLDAN'*OLDC;
               FACTOR=FACTOR.*(sdiff.*FACTOR>0).*(YL.*sdiff<=eps);
               WT=1e-2;
               if max(abs(FACTOR))==0; FACTOR=1e-5*sign(sdiff); end
               while YL'*sdiff < (eps*norm(HESS,'fro')) && WT < 1/eps
                  YL=YL+WT*FACTOR;
                  WT=WT*2;
               end
            else
               how=' Hessian modified';
            end
         end
         
         if haveoutputfcn
             % Use the xOutputfcn and optimValues from last call to outputfcn (do not call
             % callOutputFcn)
             stop = feval(outputfcn,xOutputfcn,optimValues,'interrupt',outputfcnVarargin{:});
             if stop
                 [x,FVAL,lambda_out,EXITFLAG,OUTPUT,GRADIENT,HESS] = ...
                     cleanUpInterrupt(xOutputfcn,optimValues);
                 if verbosity > 0
                     disp(OUTPUT.message)
                 end
                 return;
             end
         end
         
         %----------Perform BFGS Update If YL'S Is Positive---------
         if YL'*sdiff>eps
             HESS=HESS ...
                 +(YL*YL')/(YL'*sdiff)-((HESS*sdiff)*(sdiff'*HESS'))/(sdiff'*HESS*sdiff);
             % BFGS Update using Cholesky factorization  of Gill, Murray and Wright.
             % In practice this was less robust than above method and slower. 
             %   R=chol(HESS); 
             %   s2=R*S; y=R'\YL; 
             %   W=eye(numberOfVariables,numberOfVariables)-(s2'*s2)\(s2*s2') + (y'*s2)\(y*y');
             %   HESS=R'*W*R;
         else
            how=' Hessian not updated';
         end
      else % First call
         OLDLAMBDA=repmat(eps+gf'*gf,ncstr,1)./(sum(AN'.*AN')'+eps);
         ACTIND = 1:eq;     
      end % if numGradEvals>1
      numGradEvals=numGradEvals+1;
   
      LOLD=LAMBDA;
      OLDAN=AN;
      OLDgf=gf;
      OLDC=c;
      OLDF=f;
      OLDX=XOUT;
      XN=zeros(numberOfVariables,1);
      if (meritFunctionType>0 && meritFunctionType<5)
         % Minimax and attgoal problems have special Hessian:
         HESS(numberOfVariables,1:numberOfVariables)=zeros(1,numberOfVariables);
         HESS(1:numberOfVariables,numberOfVariables)=zeros(numberOfVariables,1);
         HESS(numberOfVariables,numberOfVariables)=1e-8*norm(HESS,'inf');
         XN(numberOfVariables)=max(c); % Make a feasible solution for qp
      end
   
      GT =c;
   
      HESS = (HESS + HESS')*0.5;
   
      [SD,lambda,exitflagqp,outputqp,howqp,ACTIND] ...
         = qpsub(HESS,gf,AN,-GT,[],[],XN,eq,-1, ...
         Nlconst,size(AN,1),numberOfVariables,OPTIONS,defaultopt,ACTIND,phaseOneTotalScaling);
    
      lambdaNLP(:,1) = 0;
      lambdaNLP(ACTIND) = lambda(ACTIND);
      lambda((1:eq)') = abs(lambda( (1:eq)' ));
      ga=[abs(c( (1:eq)' )) ; c( (eq+1:ncstr)' ) ];
      if ~isempty(c)
          mg = max(ga);
      else
          mg = 0;
      end

      if strncmp(howqp,'ok',2); 
          howqp =''; 
      end
      if ~isempty(how) && ~isempty(howqp) 
          how = [how,'; '];
      end

      LAMBDA=lambda((1:ncstr)');
      OLDLAMBDA=max([LAMBDA';0.5*(LAMBDA+OLDLAMBDA)'])' ;

      %---------------LINESEARCH--------------------
      MATX=XOUT;
      MATL = f+sum(OLDLAMBDA.*(ga>0).*ga) + 1e-30;

      infeasIllPosedMaxSQPIter = strcmp(howqp,'infeasible') || ...
          strcmp(howqp,'ill posed') || strcmp(howqp,'MaxSQPIter');
      if meritFunctionType==0 || meritFunctionType == 5
         % This merit function looks for improvement in either the constraint
         % or the objective function unless the sub-problem is infeasible in which
         % case only a reduction in the maximum constraint is tolerated.
         % This less "stringent" merit function has produced faster convergence in
         % a large number of problems.
         if mg > 0
            MATL2 = mg;
         elseif f >=0 
            MATL2 = -1/(f+1);
         else 
            MATL2 = 0;
         end
         if ~infeasIllPosedMaxSQPIter && f < 0
            MATL2 = MATL2 + f - 1;
         end
      else
         % Merit function used for MINIMAX or ATTGOAL problems.
         MATL2=mg+f;
      end
      if mg < eps && f < bestf
         bestf = f;
         bestx = XOUT;
         bestHess = HESS;
         bestgrad = gf;
         bestlambda = lambda;
         bestmg = mg;
         bestOptimError = optimError;
      end
      MERIT = MATL + 1;
      MERIT2 = MATL2 + 1; 
      stepsize=2;
      while  (MERIT2 > MATL2) && (MERIT > MATL) ...
            && numFunEvals < maxFunEvals
         stepsize=stepsize/2;