有效集 matlab代码

%有效集
function activeset 
H=[2 -1; -1 4]; 
c=[-1 -10]'; 
Ae=[ ]; be=[ ]; 
Ai=[-3 -2; 1 0; 0 1]; 
bi=[-6 0 0]'; 
x0=[0 0]'; 
[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0)


function [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) 

% 初始化 
epsilon=1.0e-9; err=1.0e-6; 
k=0; x=x0; 
n=length(x); kmax=1.0e3; 
ne=length(be); ni=length(bi); 
lamk=zeros(ne+ni,1); 
index=ones(ni,1); 
for i=1:ni 
    if (Ai(i,:)*x>bi(i)+epsilon)
        index(i)=0;
    end
end 
while (k<=kmax) 
    %求解子问题 
    Aee=[]; if(ne>0), Aee=Ae; end 
    for j=1:ni
        if(index(j)>0), Aee=[Aee; Ai(j,:)]; end
    end
    gk=H*x+c; 
    [m1,n1] = size(Aee); 
    [dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1)); 
    if(norm(dk)<=err) 
        y=0.0;
        if(length(lamk)>ne) 
            [y,jk]=min(lamk(ne+1:length(lamk))); end 
        if(y>=0) 
            exitflag=0; 
        else
            exitflag=1;
            for i=1:ni 
                if(index(i)&&(ne+sum(index(1:i)))==jk) 
                    index(i)=0; 
                    break; 
                end
            end
        end
        k=k+1; 
    else
        exitflag=1; 
        %求步长 
        alpha=1.0; tm=1.0; 
        for i=1:ni 
            if((index(i)==0)&&(Ai(i,:)*dk<0)) 
                tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk); 
                if(tm1<tm) 
                    tm=tm1; ti=i; 
                end
            end
        end
        alpha=min(alpha,tm); 
        x=x+alpha*dk; 
        %修正有效集 
        if(tm<1), index(ti)=1; 
        end
    end
    if(exitflag==0), break; end
    k=k+1;
end
output.fval=0.5*x'*H*x+c'*x; 
output.iter=k;

% 求解子问题 
function [x,lambda]=qsubp(H,c,Ae,be) 
ginvH=pinv(H); 
[m,n]=size(Ae); 
if(m>0) 
    rb=Ae*ginvH*c + be; 
    lambda=pinv(Ae*ginvH*Ae')*rb; 
    x=ginvH*(Ae'*lambda-c); 
else
    x=-ginvH*c;
    lambda=0; 
end