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