direction.m

基于matlab的约束非线性规划算法库

function [Aset, Abounds, T, GNlambda, p, Jc, b, betak, lenJp, cTc, JpTf, Jp] = Direction(phase,m,n,X,f,c,k,...
                  Teq,Tineq,LB,UB,Aset,Abounds,ffun,cfun,params)
%Call: [Aset, Abounds, T, GNlambda, p, Jc, b, betak, lenJp, cTc, JpTf] = Direction(phase,m,n,X,f,c,k,...
%                 Teq,Tineq,LB,UB,Aset,Abounds,ffun,cfun,params)
%Form the current active set and compute the searchdirection p
%and the GN-Lagrange multipliers (GNlambda).
%Aset holds the indeces of c(X) that make up the current active set
%PredReduc is the predicted reduction in the minimization
%On entry:
%X is the current point
%f is F(X) - the residual vector (m*1) at the current point
%c is c(X) - the constraints value (Teq+Tineq)*1 at the current point
%LB holds the lower bounds
%UB holds the upper bounds
%Abounds indicates what variables that attain their lower or upper bound respectively
%
epsrank = sqrt(eps);
sqrteps = epsrank;
delta = 10;
if ~isempty(Aset)
    b = c(Aset); % form a vector containing the value of constraints in current working set
    T = max(size(b));
else
    b=[];
    T=0;
end
if T > n
    disp('Direction: more active constraints than n')
    pause
end
%Find indeces for free and fixed variables.
free = find(Abounds == 0);
fixed = find(Abounds ~= 0);
nfixed = max(size(fixed)); %Number of fixed variables
%Compute the Jacobians for c(X) and f(X) respectively
%Use differens approximated Jacobian if user don't provide it
[Dummy, Jc] = feval(cfun,X,1,params{:});
if isempty(Jc)
    if ~(Teq==0 & Tineq==0)% check if any constraint at all
       Jc = jacdif(cfun,X,params);
    end
end
[Dummy J] = feval(ffun,X,1,params{:});
if isempty(J)
    J = jacdif(ffun,X,params);
end
Tnfixed = T+nfixed;
%%%%%%%%%%%%%
%[GNlambda, r, p, betak, lenJp, JpTf] = SysEqu(Aset,Jc,Teq,T,J,Tnfixed,nfixed,fixed,Abounds,b,f,m,n);
%[GNlambda, lambda, r, p, betak, lenJp, JpTf] = PseudoRank(Aset,Jc,Teq,T,J,Tnfixed,nfixed,fixed,Abounds,b,f,m,n);
[GNlambda, lambda, r, p, betak, lenJp, JpTf,Jp,DeltaP1,Q1,P2,J1,J2] = ...
                    Pseudo2(Aset,Jc,Teq,T,J,Tnfixed,nfixed,fixed,Abounds,b,f,m,n);
if phase == 2 %Newton should be used
    if T == 0 %The unconstrained case
        [p, flag] = NewtonUC(Abounds,n-nfixed,J,f,X,ffun,params);
    else
        [p, flag] = NewtonC(Aset,Abounds,n-nfixed,J1,f,X,DeltaP1,Q1,P2,J2,lambda,ffun,cfun,params);
    end
    if flag ~= 0
       disp('Not pos. def. from chol()')
    end
    cTc = dot(b,b);
else %Gauss-Newton should be used
if Tnfixed == 0
    cTc = 0;
else
    cTc = dot(b,b);
    %Check Lagrange multipliers
    NegInd = find(lambda(Teq+1:Tnfixed) < -sqrteps);
    Minlambda = min(lambda(Teq+NegInd)); %value of most negative multiplier
    Ind = min(find(lambda == Minlambda)); % index of most negative multiplier
    %if MinGNl is empty there is no constraint to delete
    if ~isempty(Minlambda) & norm(r)<(-Minlambda*delta)
   %Check GN-Lagrange multipliers
    %NegInd = find(GNlambda(Teq+1:Tnfixed) < -sqrteps)
    %MinGNl = min(GNlambda(Teq+NegInd)) %value of most negative multiplier
    %Ind = min(find(GNlambda == MinGNl)) % index of most negative multiplier
    %if MinGNl is empty there is no constraint to delete
    %if ~isempty(MinGNl) & norm(r)<(-MinGNl*delta)
    %if ~isempty(MinGNl) & norm(r)<(-MinGNl*delta) & lambda(Ind)<0 & ...
            %2*abs(MinGNl-lambda(Ind))<min(abs(MinGNl),abs(lambda(Ind))) %Check if we are close enough
        disp('One constraint is deleted')
        if Ind > T
            % free a fixed variable
            Abounds(fixed(Ind-T)) = 0;
            %Find indeces for free variables.
            free = find(Abounds == 0);
            fixed = find(Abounds ~= 0);
            nfixed = max(size(fixed)); %Number of fixed variables
        else
            % delete one inequality from the working set
            Aset = setdiff(Aset,Aset(Ind));
            b = c(Aset); % form a vector containing the value of constraints in current working set
            T = max(size(b));
        end
        Tnfixed = T+nfixed;
        [GNlambda, lambda, r, p, betak, lenJp, JpTf,Jp,DeltaP1,Q1,P2,J1,J2] = Pseudo2(Aset,Jc,Teq,T,J,Tnfixed,nfixed,fixed,Abounds,b,f,m,n);
        cTc = dot(b,b);
   end
end
end