function [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30]=suanfa()
%参数设定
G_max=250000;
M_0=8;
V_0=5;
T=10;
N=10;
np=5;
t=50;
n=30;
syms x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26 x27 x28 x29 x30;
for g=1:G_max
%决定分流数目
M(1:n)=[M_0];
V(1:n)=[V_0];
n(1:N)=[0];
for i=1:N
n(i)=min(max(1,floor(M(i)*V(i)/T)),np);
end
%水流演化
A(1:N,1:n)=[0];
TU(1:N,1:3)=[0];
X(1:N,1:n)=[0];
U(1:N,1:n(i))=[0];
r=rand(1)*((-2)*1)+1;
for i=1:N
    X(i,1:n)=eval('[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30]');
    for I=1:n
       q=q();%  调用子函数
       X(i,I)=X(i,I)+q;
       A(i,I)=X(i,1:n);
    end
    TU(i,3)=TU(i,2);
    TU(i,2)=TU(i,1);
    TU(i,1)=X(i);
    for k=i:n(i)
        U(i,k)=f1(X(i,1:n));%调用子函数   
    end
end
%计算分流重量与分流速度
for i=1:N
    for k=1:n(i)
        for v=1:n(i)
            a=0;
            a=a+v;
        end
        w(i,k)=((n(i)+1-k)/a)*M(i);
        S(i)=M(i)*f1(X(i,1:n));
        B(i,k)=w(i,k)*f1(U(i,k));
        M_g=M(i)/w(i,k);
        if squrt(2*9.8*(S(i)-B(i,k))/w(i,k)+M_g*V(i)^2)>0
            b(i,k)=squrt(2*9.8*(S(i)-B(i,k))/w(i,k)+M_g*V(i)^2);
        else b(i,k)=0;
        end
    end
end        
%汇流
for i=1:N
    for k=1:n(i)
        for K=k:n(i)
        if U(i,k)==U(i,K)
            m(i,k)=w(i,K)+w(i,k);
            v(i,k)=(w(i,k)*b(i,k)+w(i,K)*b(i,k))/(w(i,k)+w(i,K));
        else m(i,k)=w(i,k);
            v(i,k)=b(i,k);
        end
        end
    end
end
%判断
v=0;
for i=1:N
    for k=1:n(i)
        if v(i,k)>0
            v=1;
        end
    end
end
if v==1
    for i=1:N
        for k=1:n(i)
            m(i,k)=(1-1/t)*m(i,k);
        end
    end
    %降雨
    for i=1:N
        Xp(1:i,1:n)=[0];
        for c=1:n
            X(i,c)=X(i,c)+q;
            Xp(i,1:n)=X(i,1:n);
        end
        for k=1:i
            for c=1:n(i)
                M_k=0;
                M_k=m(k,c)+M_k;
            end
        end
        Mp(i)=(M_0-M_k)*M(i)/M_k;
        Vp(i)=V_0;
    end
    if Xp(i,1:n)<TU(i,1)
        TU(i,1)=Xp(i,1:n);
    else if Xp(i,1:n)<TU(i,2)
            TU(i,2)=Xp(i,1:n);
        else if Xp(i,1:n)<TU(i,3)
                TU(i,3)=Xp(i,1:n);
            else X(i)=U(i,1:n)
            end
        end
    end
    else for i=1:N
        for k=1:n(i)
            m(i,k)=0;
        end
    end
    %降雨
    for i=1:N
        Xp(1:i,1:n)=[0];
        for c=1:n
            X(i,c)=X(i,c)+q;
            Xp(i,1:n)=X(i,1:n);
        end
        for k=1:i
            M_k=0;
            M_k=M(k)+M_k;
        end
        Mp(i)=M(i)*M_0/M_k;
        Vp(i)=V_0;
    end
    if Xp(i,1:n)<TU(i,1)
        TU(i,1)=Xp(i,1:n);
    else if Xp(i,1:n)<TU(i,2)
            TU(i,2)=Xp(i,1:n);
        else if Xp(i,1:n)<TU(i,3)
                TU(i,3)=Xp(i,1:n);
            else X(i)=U(i,1:n)
            end
        end
    end
end
%水流演化
A(1:N,1:n)=[0];
TU(1:N,1:3)=[0];
X(1:N,1:n)=[0];
U(1:N,1:n(i))=[0];
r=rand(1)*((-2)*1)+1;
for i=1:N
X(i,1:n)=Xp(i,1:n);
    for I=1:n
       q=q( input_args );%  调用子函数
       X(i,I)=X(i,I)+q;
       A(i,I)=X(i,1:n);
    end
    TU(i,3)=TU(i,2);
    TU(i,2)=TU(i,1);
    TU(i,1)=X(i);
    for k=i:n(i)
        U(i,k)=f1(X(i,1:n));%调用子函数   
    end
end
end
for i=1:N
    [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30]=A(i,1:n);
end
end
%子函数
function [f] = f1(X)%子函数 sphere测试函数
f=0;
    for i=1:n
        f=X(i)^2+f;
    end
end
function [f] = f2(X)%子函数 rosenbrock测试函数
d=n-1;
f=0;
for i=1:d
  f=(X(i)^2-X(i+1))^2+(X(i)-1)^2+f;
end
end
function [f] = f3(X)%子函数 rastrigrin测试函数
d=0;
for i=1:n
  d=X(i)^2-10*cos(2*pi*X(i))+d;
end
f=10*n+d;
end
function [q] =q()% 子函数
q=rand(1)*((-2)*r)+r;
end