nlconst.m

matpower软件下

         if stepsize < 1e-4,  
            stepsize = -stepsize; 
         
            % Semi-infinite may have changing sampling interval
            % so avoid too stringent check for improvement
            if meritFunctionType == 5, 
               stepsize = -stepsize; 
               MATL2 = MATL2 + 10; 
            end
         end
         if hasBoundOnStep && (iter <= relLineSrchBndDuration)
           % Bound total displacement:
           % |stepsize*SD(i)| <= relLineSrchBnd*max(|x(i)|, |typicalx(i)|)
           % for all i.
           indxViol = abs(stepsize*SD) > relLineSrchBnd*max(abs(MATX),abs(typicalx));
           if any(indxViol)
             stepsize = sign(stepsize)*min(  min( abs(stepsize), ...
                  relLineSrchBnd*max(abs(MATX(indxViol)),abs(typicalx(indxViol))) ...
                  ./abs(SD(indxViol)) )  );
           end   
         end
         
         XOUT = MATX + stepsize*SD;
         x(:)=XOUT; 
      
         if strcmp(funfcn{2},'fseminf')
            f = feval(funfcn{3},x,varargin{3:end});
         
            [nctmp,nceqtmp,NPOINT,NEWLAMBDA,OLDLAMBDA,LOLD,s] = ...
               semicon(x,LAMBDA,NEWLAMBDA,OLDLAMBDA,POINT,FLAG,s,varargin{:});
            nctmp = nctmp(:); nceqtmp = nceqtmp(:);
            non_ineq = length(nctmp);  % the length of nctmp can change
            ineq = non_ineq + lin_ineq;
            ncstr = ineq + eq;
            % Possibly changed constraints, even if same number,
            % so ACTIND may be invalid.
            ACTIND = 1:eq;
         else
            f = feval(funfcn{3},x,varargin{:});
            if constflag
               [nctmp,nceqtmp] = feval(confcn{3},x,varargin{:});
               nctmp = nctmp(:); nceqtmp = nceqtmp(:);
            else
               nctmp = []; nceqtmp=[];
            end
         end
            
         nc = [nceqtmp(:); nctmp(:)];
         c = [Aeq*XOUT-Beq; nceqtmp(:); A*XOUT-B; nctmp(:)];  

         numFunEvals = numFunEvals + 1;
         ga=[abs(c( (1:eq)' )) ; c( (eq+1:length(c))' )];
         if ~isempty(c)
            mg=max(ga);
         else
            mg = 0;
         end

         MERIT = f+sum(OLDLAMBDA.*(ga>0).*ga);
         if meritFunctionType == 0 || meritFunctionType == 5
            if mg > 0
               MERIT2 = mg;
            elseif f >=0 
               MERIT2 = -1/(f+1);
            else 
               MERIT2 = 0;
            end
            if ~infeasIllPosedMaxSQPIter && f < 0
               MERIT2 = MERIT2 + f - 1;
            end
         else
            MERIT2=mg+f;
         end
         if haveoutputfcn
             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

                                                                                                                                                                                                                            end  % line search loop
      %------------Finished Line Search-------------
   
      if meritFunctionType~=5
         mf=abs(stepsize);
         LAMBDA=mf*LAMBDA+(1-mf)*LOLD;
      end

      x(:) = XOUT;
      switch funfcn{1} % evaluate function gradients
      case 'fun'
         ;  % do nothing...will use finite difference.
      case 'fungrad'
         [f,gf_user] = feval(funfcn{3},x,varargin{:});
         gf_user = gf_user(:);
         numGradEvals=numGradEvals+1;
      case 'fun_then_grad'
         gf_user = feval(funfcn{4},x,varargin{:});
         gf_user = gf_user(:);
         numGradEvals=numGradEvals+1;
      otherwise
         error('optim:nlconst:UndefinedCalltypeInFMINCON', ...
               'Undefined calltype in FMINCON.');
      end
      numFunEvals=numFunEvals+1;
   
   
      % Evaluate constraint gradients
      switch confcn{1}
      case 'fun'
         gnceq=[]; gncineq=[];
      case 'fungrad'
         [nctmp,nceqtmp,gncineq,gnceq] = feval(confcn{3},x,varargin{:});
         nctmp = nctmp(:); nceqtmp = nceqtmp(:);
         numGradEvals=numGradEvals+1;
      case 'fun_then_grad'
         [gncineq,gnceq] = feval(confcn{4},x,varargin{:});
         numGradEvals=numGradEvals+1;
      case ''
         nctmp=[]; nceqtmp =[];
         gncineq = zeros(numberOfVariables,length(nctmp));
         gnceq = zeros(numberOfVariables,length(nceqtmp));
      otherwise
         error('optim:nlconst:UndefinedCalltypeInFMINCON', ...
               'Undefined calltype in FMINCON.');
      end
      gnc_user = [gnceq, gncineq];
      gc = [Aeq', gnceq, A', gncineq];
   
   end % if ~done   
end % while ~done


% Update 
numConstrEvals = numGradEvals;

% Gradient is in the variable gf
GRADIENT = gf;

% If a better solution was found earlier, use it:
if f > bestf 
   XOUT = bestx;
   f = bestf;
   HESS = bestHess;
   GRADIENT = bestgrad;
   lambda = bestlambda;
   mg = bestmg;
   gf = bestgrad;
   optimError = bestOptimError;
end

FVAL = f;
x(:) = XOUT;

if haveoutputfcn
    [xOutputfcn, optimValues] = callOutputFcn(outputfcn,XOUT,xOutputfcn,'done',iter,numFunEvals, ...
        f,mg,stepsize,gf,SD,meritFunctionType,optimError,how,howqp,outputfcnVarargin{:});
    % Do not check value of 'stop' as we are done with the optimization
    % already.
end

OUTPUT.iterations = iter;
OUTPUT.funcCount = numFunEvals;
OUTPUT.stepsize = stepsize;
OUTPUT.algorithm = 'medium-scale: SQP, Quasi-Newton, line-search';
if meritFunctionType == 1
   OUTPUT.firstorderopt = [];
else   
   OUTPUT.firstorderopt = optimError;
end
OUTPUT.cgiterations = [];
OUTPUT.message = outMessage;

[lin_ineq,Acol] = size(Ain);  % excludes upper and lower

lambda_out.lower=zeros(lenlb,1);
lambda_out.upper=zeros(lenub,1);

lambda_out.eqlin = lambdaNLP(1:lin_eq);
ii = lin_eq ;
lambda_out.eqnonlin = lambdaNLP(ii+1: ii+ non_eq);
ii = ii+non_eq;
lambda_out.lower(arglb) = lambdaNLP(ii+1 :ii+nnz(arglb));
ii = ii + nnz(arglb) ;
lambda_out.upper(argub) = lambdaNLP(ii+1 :ii+nnz(argub));
ii = ii + nnz(argub);
lambda_out.ineqlin = lambdaNLP(ii+1: ii + lin_ineq);
ii = ii + lin_ineq ;
lambda_out.ineqnonlin = lambdaNLP(ii+1 : end);

% NLCONST finished
%--------------------------------------------------------------------------
function [xOutputfcn, optimValues, stop] = callOutputFcn(outputfcn,x,xOutputfcn,state,iter,numFunEvals, ...
    f,mg,stepsize,gf,SD,meritFunctionType,optimError,how,howqp,varargin)
% CALLOUTPUTFCN assigns values to the struct OptimValues and then calls the
% outputfcn.  
%
% The input STATE can have the values 'init','iter', or 'done'. 
% We do not handle the case 'interrupt' because we do not want to update
% xOutputfcn or optimValues (since the values could be inconsistent) before calling
% the outputfcn; in that case the outputfcn is called directly rather than
% calling it inside callOutputFcn.
%
% For the 'done' state we do not check the value of 'stop' because the
% optimization is already done.

optimValues.iteration = iter;
optimValues.funccount = numFunEvals;
optimValues.fval = f;
optimValues.constrviolation = mg;
optimValues.stepsize = stepsize;
if ~isempty(SD)
    optimValues.directionalderivative = gf'*SD;  
else
    optimValues.directionalderivative = [];
end
optimValues.gradient = gf;
optimValues.searchdirection = SD;
if meritFunctionType == 1
    optimValues.firstorderopt = [];
else
    optimValues.firstorderopt = optimError;
end
optimValues.procedure = [how,'  ',howqp];
xOutputfcn(:) = x;  % Set x to have user expected size
switch state
    case {'iter','init'}
        stop = feval(outputfcn,xOutputfcn,optimValues,state,varargin{:});
    case 'done'
        stop = false;
        feval(outputfcn,xOutputfcn,optimValues,state,varargin{:});
    otherwise
        error('optim:nlconst:UnknownStateInCALLOUTPUTFCN', ...
              'Unknown state in CALLOUTPUTFCN.')
end

%--------------------------------------------------------------------------
function [x,FVAL,lambda_out,EXITFLAG,OUTPUT,GRADIENT,HESS] = cleanUpInterrupt(xOutputfcn,optimValues)
% CLEANUPINTERRUPT updates or sets all the output arguments of NLCONST when the optimization 
% is interrupted.  The HESSIAN and LAMBDA are set to [] as they may be in a
% state that is inconsistent with the other values since we are
% interrupting mid-iteration.

x = xOutputfcn;
FVAL = optimValues.fval;
EXITFLAG = -1; 
OUTPUT.iterations = optimValues.iteration;
OUTPUT.funcCount = optimValues.funccount;
OUTPUT.stepsize = optimValues.stepsize;
OUTPUT.algorithm = 'medium-scale: SQP, Quasi-Newton, line-search';
OUTPUT.firstorderopt = optimValues.firstorderopt; 
OUTPUT.cgiterations = [];
OUTPUT.message = 'Optimization terminated prematurely by user.';
GRADIENT = optimValues.gradient;
HESS = []; % May be in an inconsistent state
lambda_out = []; % May be in an inconsistent state

%--------------------------------------------------------------------------
function [activeLb,activeUb,activeIneqLin,activeIneqNonlin] = ...
    activeInequalities(c,tol,arglb,argub,linEq,nonlinEq,linIneq)
% ACTIVEINEQUALITIES returns the indices of the active inequalities
% and bounds.
% INPUT:
% c                 vector of constraints and bounds (see nlconst main code)
% tol               tolerance to determine when an inequality is active
% arglb, argub      boolean vectors indicating finite bounds (see nlconst
%                   main code)
% linEq             number of linear equalities
% nonlinEq          number of nonlinear equalities
% linIneq           number of linear inequalities
%
% OUTPUT
% activeLB          indices of active lower bounds
% activeUb          indices of active upper bounds  
% activeIneqLin     indices of active linear inequalities
% activeIneqNonlin  indices of active nonlinear inequalities
%

% We check wether a constraint is active or not using '< tol'
% instead of '<= tol' to be onsistent with nlconst main code, 
% where feasibility is checked using '<'.
finiteLb = nnz(arglb); % number of finite lower bounds
finiteUb = nnz(argub); % number of finite upper bounds

indexFiniteLb = find(arglb); % indices of variables with LB
indexFiniteUb = find(argub); % indices of variables with UB

% lower bounds
i = linEq + nonlinEq; % skip equalities

% Boolean vector that indicates which among the finite
% bounds is active
activeFiniteLb = abs(c(i + 1 : i + finiteLb)) < tol;

% indices of the finite bounds that are active
activeLb = indexFiniteLb(activeFiniteLb);

% upper bounds
i = i + finiteLb;

% Boolean vector that indicates which among the finite
% bounds is active
activeFiniteUb = abs(c(i + 1 : i + finiteUb)) < tol;

% indices of the finite bounds that are active
activeUb = indexFiniteUb(activeFiniteUb);

% linear inequalities
i = i + finiteUb;
activeIneqLin = find(abs(c(i + 1 : i + linIneq)) < tol); 
% nonlinear inequalities
i = i + linIneq;
activeIneqNonlin = find(abs(c(i + 1 : end)) < tol);   

%--------------------------------------------------------------------------
function printColumnwise(a,b,c,d)
% PRINTCOLUMNWISE prints vectors a, b, c, d (which
% in general have different lengths) column-wise.
% 
% Example: if a = [1 2], b = [4 6 7], c = [], d = [8 11 13 15]
% then this function will produce the output (without the headers):
%
% a  b  c   d
%-------------
% 1  4      8
% 2  6     11
%    7     13
%          15
%
length1 = length(a); length2 = length(b);
length3 = length(c); length4 = length(d);

for k = 1:max([length1,length2,length3,length4])
    % fprintf stops printing numbers as soon as it encounters [].
    % To avoid this, we convert all numbers to string
    % (fprintf doesn't stop when it comes across the blank
    % string ' '.)
   if k <= length1
      value1 = num2str(a(k));
   else
      value1 = ' ';
   end
   if k <= length2
      value2 = num2str(b(k));
   else
      value2 = ' ';
   end   
   if k <= length3
      value3 = num2str(c(k));
   else
      value3 = ' ';
   end      
   if k <= length4
      value4 = num2str(d(k));
   else
      value4 = ' ';
   end  
   fprintf('%5s %10s %10s %10s\n',value1,value2,value3,value4);
end