simple.m

线性规划 包括单纯形法

function simple = simple(N,a,x)
Fun = @fun;
eps = 1/(10.^12);
Sig = 10000;
alpha = 1;
beta = 0.5;
gama = 2;
count =0;
p = a*(sqrt(N+1)+N-1)/(sqrt(2)*N);
q = a*(sqrt(N+1)-1)/(sqrt(2)*N);
simp = zeros(N,N+1);
reflect = ones(0,N);
f_worst = 0;
i_worst = 0;
f_best = 0;
i_best = 0;

%--------------------建立单纯行矩阵------------------%
for i = 1:N
    for j = 1:N+1;
        simp(i,j) = x(i);
    end
end

for i =1:N
    for j =2:N+1
        if( ( j - i ) ~= 1)
        simp(i,j) = x(i) + q;
        end
    end
end

for i = 2:N+1
   simp(i-1,i) = x(i-1) +p;
end

%---------------获取x()时函数f()的值------------------%

for i = 1:N+1
        f(i)= feval(Fun,(simp(:,i)));
end

%----------------开始迭代----------------%
while( Sig > eps)
    count = count + 1;
    f_worst = f(1);
    i_worst = 1;
    f_best = f(1);
    i_best = 1;
%----------------求得最大值---------------%    
    for i = 2:N+1
        if(f(i) > f_worst)
            i_worst = i;
            f_worst = f(i);
        end
        if(f(i)< f_best)
            i_best = i;
            f_best = f(i);
        end
    end
  
  %------------------求得次最大值----------------%  
    i_subworst = 1;
    f_subworst = f(1);
    if(i_worst == 1)
        i_subworst = 2;
        f_subworst = f(2);
    end
    
    for i = (i_subworst+1):N+1
        if((i~=i_worst) & (f(i)>f_subworst))
            i_subworst = i;
            f_subworst = f(i);
        end
    end
    
    %---------------------求得形心-----------------------%
   for j = 1:N
     x_center(j) = 0;
     for i = 1:N+1
         if (i~=i_worst)
             x_center(j) = x_center(j) + simp(j,i);
         end
     end
     x_center(j) = x_center(j)/N;
    end

    %--------------------得到投影点-----------------------%
    for(i = 1:N)
        x_reflect(i) = (1+alpha)*x_center(i)-alpha*simp(i,i_worst);
    end
     
    f_reflect = fun(x_reflect);
    %---------------------放大-----------------------%
    if(f_reflect < f(i_best))
        for (i = 1:N)
            x_expand(i) = gama*x_reflect(i)+(1-gama)*x_center(i);
        end
        f_expand = fun(x_expand);
        if(f_expand <f(i_best))
            for (i = 1:N)
                simp(i,i_worst) = x_expand(i);
            end
            f(i_worst) = f_expand;
        else
            for(i = 1:N)
                simp(i,i_worst)= x_reflect(i);
            end
            fun(i_worst) = f_reflect;
        end
        
    else
            if(f_reflect <= f(i_subworst))
                for (i = 1:N)
                    simp(i,i_worst) = x_reflect(i);
                end
                f(i_worst) = f_reflect;
                continue;
            else
                if(f_reflect <= f(i_worst))
                    for j = 1:N
                        simp(j,i_worst) = x_reflect(j)
                    end
                    f(i_worst) = f_reflect;
              
                end
            end
            %------------------收缩-------------------%
            for i = 1:N;
                x_compact(i) = beta*simp(i,i_worst)+(1-beta)*x_center(i);
            end
            f_compact = feval(Fun,x_compact);
            
            %----------------收缩不行就重置------------%
            if(f_compact > f(i_worst)) 
                for(i = 1:N)
                    for( j = 1:N+1)
                        simp(i,j) = (simp(i,j) + simp(i,i_best))/2;
                    end
                end
            else
                for (i = 1:N)
                    simp(i,i_worst) = x_compact(i);
                end
                f(i_worst) = f_compact;
            end
            
            %-------------------判断结束条件------------------%
            Sig = 0;
            f_center = feval(Fun,x_center);

            for (j = 1:N+1)
                Sig = Sig + (feval(Fun,(simp(:,j))) -f_center).^2;
            end
            Sig = sqrt(Sig/(N+1));
        end
    end             
end

%-------------------显示结果-------------------%
x = (simp(:,1))'
fmin = feval(Fun,(simp(:,1)))