matlab_无线传感器网络路由优化_蚁群算法_有消耗能量_生命周期

clear tic; M=100; A=[0;0]; B=[10;10]; %sink rand('state',1) G=rand(2,100)*10; %在10*10空间里面随机生成100个节点 U=[A G B]; C=U'; C0=0.2; f=5; K=100; %迭代次数（指蚂蚁出动多少波） Alpha0=1; Alpha1=1; Tau_max=88;%信息数 Tau_min=0.1; w1=0.2; w2=0.85; c2=2; d0=2; dead=0; Alpha=1; Gamma=1; Beta=5; Rh=0.5; Beta0=5;Rho=0.5;Q=500; Beta1=5; flagmax=0; flagmin=0; N=size(C,1); D=ones(N,N); Length_best=ones(K,1); Routh_best=cell(K,1); Best_path=cell(K,1); Best_length=ones(K,1); Length_average=ones(K,1); for i=1:N for j=1:N if i < j D(i,j)=sqrt((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2); end D(j,i)=D(i,j); end end Eta=ones(N,N); Eta2=ones(N,N); Eta3=ones(N,N); Eta4=ones(N,N); D_average=mean(D); R_max=D_average; for i=1:N for j=1:N Eta(i,j)=1/(w1*D(i,j)+w2*D(j,N)); %表示节点i与节点j间路径能量消耗启发因子 没有考虑到不同的能耗模型 %Eta(i,j)=1/D(j,N); %Eta(i,j)=1/D(i,j); %Eta(i,j)=1/(D(i,j)+D(j,N)); Eta2(i,j)=D(j,N); Eta3(i,j)=D(i,j); Eta4(i,j)=D(i,N); end end rand('state',5) %rand('state',0)作用在于如果指定状态，产生随机结果相同 E=ones(N,1); E1=zeros(N,1); E2=zeros(N,1); E3=ones(N,1); E10=zeros(N,1); E20=zeros(N,1); E30=ones(N,1); kd=2000; Eelec=50*10.^-9; %接收能量损耗 Eamp=100*10.^-12; %传输损耗系数 Tau = ones(N,N); Routes=cell(K,M);%%用细胞结构存储每一代的每一只蚂蚁的爬行路线 PL=zeros(K,M); %用矩阵存储每一代的每一只蚂蚁的爬行路线长度 dead_first_flag=0; %%?-----------启动K轮蚂蚁觅食活动，每轮派出M只蚂蚁---- for k=1:K disp(k); % Rh =1+(Rho*k/K); % Alpha = Alpha1*(1+exp(-C0*k)); % Beta = Beta0*(1+exp(-C0*k)); for m=1:M %%?第一步：状态初始化? W=1; %当前节点初始化为起始点 Path=1;%爬行路线初始化 PLkm=0;%爬行路线长度初始化 Tabukm=ones(1,N); %禁忌表初始化 Tabukm(1)=0;%已经在初始点了，因此要排除 DD=D;%邻接矩阵初始化 %%?第二步：下一步可以前往的节点 while W~=N for j=1:N if Tabukm(j)==0 %判断是否节点被访问过 DD(W,j)=inf; %如果访问过则另其当前位置到任意一点的距离无穷大 if(E3(i)<0) dead=dead+1; if(dead==1) dead_first_flag=k disp(k) end end end end LJD=find(DD(W,:)<c2*R_max);%簇头选取 -----------------------可改进 Len_LJD=length(LJD);%可选节点的个数 %% 觅食停止条件：蚂蚁陷入死胡同 if Len_LJD>=1 %%?第三步：转轮赌法选择下一步怎么走? PP=zeros(1,Len_LJD); for i=1:Len_LJD if 0.5*d0<Eta2(W,LJD(i))<d0 % if 1.98*Eta3(W,LJD(i))*Eta4(W,LJD(i))<Eta4(W,LJD(i))^2+Eta3(W,LJD(i))^2-Eta2(W,LJD(i))^2<2*Eta3(W,LJD(i))*Eta4(W,LJD(i)) if sqrt(3)*Eta3(W,LJD(i))*Eta4(W,LJD(i))<Eta4(W,LJD(i))^2+Eta3(W,LJD(i))^2-Eta2(W,LJD(i))^2<2*Eta3(W,LJD(i))*Eta4(W,LJD(i)) %if sqrt(2)*Eta3(W,LJD(i))*Eta4(W,LJD(i))<Eta4(W,LJD(i))^2+Eta3(W,LJD(i))^2-Eta2(W,LJD(i))^2<2*Eta3(W,LJD(i))*Eta4(W,LJD(i)) % %2016-6-12加入蚂蚁搜索方向选择-45~45 %if Eta3(W,LJD(i))*Eta4(W,LJD(i))<Eta4(W,LJD(i))^2+Eta3(W,LJD(i))^2-Eta2(W,LJD(i))^2<2*Eta3(W,LJD(i))*Eta4(W,LJD(i)) %if Eta4(W,LJD(i))^2+Eta3(W,LJD(i))^2-Eta2(W,LJD(i))^2<2*Eta3(W,LJD(i))*Eta4(W,LJD(i)) %2016-6-12加入蚂蚁搜索方向选择 -90~90 PP(i)=(Tau(W,LJD(i))^Alpha)*(Eta(W,LJD(i))^Beta)*(E3(LJD(i))^Gamma);%% -----------------------可改进 %PP(i)=(Tau(W,LJD(i))^Alpha)*(E3(LJD(i))); end end end PP=PP/(sum(PP)); %建立概率分布 Pcum=cumsum(PP); %B = cumsum(A)返回数组不同维数的累加和。 Select=find(Pcum>=rand); to_visit=LJD(Select(1)); %%?第四步：状态更新和记录? Path =[Path ,to_visit]; %路径增加 PLkm=PLkm+DD(W,to_visit);%路径长度增加 % d0=0.3; % if DD(W,to_visit)>d0 % E10(W)=E10(W)+Eelec*kd+Eamp*kd*(DD(W,to_visit)-d0)^4+d0^2; %节点总能耗 Eelec*kd为发送能耗 Eamp为传输系数 % E20(to_visit)=E20(to_visit)+Eelec*kd; %节点接受能耗 % E30=E-E10-E20; %节点剩余能量 % else %%?-------------------------------------------加入能耗只考虑-自由能量模型 没有考虑多路径衰减模型------------------------------------------------------------------------------- E1(W)=E1(W)+Eelec*kd+Eamp*kd*(DD(W,to_visit))^2; %节点总能耗 Eelec*kd为发送能耗 Eamp为传输系数 E2(to_visit)=E2(to_visit)+Eelec*kd; %节点接受能耗 E3=E-E1-E2; %节点剩余能量 % end %%?----------------------------------------------end--------------------------------------------------------------------------------------- Tau(W,to_visit)=(1-Rh)*Tau(W,to_visit)+Q/D(W,to_visit);%局部信息素更新 ————————可改进 % if Tau(W,to_visit)<Tau_min;%极大极小 2016-6-29加入信息素限制，防止过早收敛 % Tau(W,to_visit)=Tau_min; % flagmin=1+flagmin; % end if Tau(W,to_visit)>Tau_max%极大极小 2016-6-29加入信息素限制，防止过早收敛 Tau(W,to_visit)=Tau_max; flagmax=1+flagmax; end W=to_visit;%蚂蚁移到下一个节点? Tabukm(W)=0; %已访问过的节点从禁忌表中删除? else break ; end end %%?第五步：记下每一代每一只蚂蚁的觅食路线和路线长度? Routes{k,m}=Path; if Path(end)==N PL(k,m)=PLkm; else PL(k,m)=0; end end %%?第六步：更新信息素 Delta_Tau=zeros(N,N); %更新量初始化? for m=1:M if PL(k,m) Rout=Routes{k,m}; TS=length(Rout)-1;%跳数? PL_km=PL(k,m); end end if k>1 PL(k,1)=Length_best(k-1); end Length_best(k)=min(PL(k,:)); e=zeros(m,1); e_average=zeros(m,1); e_min=zeros(m,1); for m=1:M e(m)=0; g=Routes{k,m}; f=length(g); for j=1:f e(m)=e(m)+E3(g(j)); %剩余能量 end e_average(m)=e(m)/f; e_min(m)=min(e(m)); L=zeros(M,1); f=zeros(M,1); L(m)=PL(k,m); %f(m)=e_average(m)/L(m);%2016-6-18路径评价 f(m)=e_average(m)*e_min(m)/L(m); %2016-6-29加入路径评价 end B=find(f==max(f)); Best_path{k}=Routes{k,B(1)}; Best_length(k)=PL(k,B(1)); for i=1:length(B) q=B(i); Route=Routes{k,q}; for j=2:length(Route) Delta_Tau(Route(j),Route(j-1))=Q/Best_length(k); end end Position=find(PL(k,:)==Length_best(k)); Routh_best{k}=Routes{k,Position(1)}; Length_average(k)=mean(PL(k,:)); for i=1:length(Position) p=Position(i); Route=Routes{k,p}; for j=2:length(Route) Delta_Tau(Route(j),Route(j-1))=Q/Length_best(k); end end Tau=Tau+Delta_Tau; figure(1) plot(k,sum(E1+E2)/100,'+r'); xlabel('发送轮数'); ylabel('节点能耗'); legend('Leach-Ant','EEABR','LSAC','OARA'); hold on; figure(2) plot(k,sum(E3)/100-0.02,'r'); xlabel('迭代次数'); ylabel('节点平均剩余能量'); hold on; end deadn=0; %第一个节点死亡代数绘图 for i=1:K for j=1:N if (E3(j)<0&&dead_first_flag==i) deadn=deadn+