基于蚁群算法的路由问题研究

在组播路由中，用matlab画图，组播总费用提前收敛，到达各个节点的费用后收敛，甚至有些节点费用一直没有收敛。不应该是每一个节点收敛后，组播树总费用才收敛吗？ 是不是我的代码有问题

以下是我的代码：

function [MRT,EDGES,cost]=ACA_QoS_MR(C,D,S,E,Dmax,K,M,Alpha,Beta,Gamma,Tau,Rho,Q)

%% Ant Colony Algorithm for QoS Multicast Routing

% QoS组播路由蚁群算法

%% 输入参数列表

% C 费用邻接矩阵（N×N）

% D 延时邻接矩阵（N×N）

% S 源节点

% E 组播目的节点（行向量）

% Dmax 延时约束

% K 迭代次数（指蚂蚁出动多少波）

% M 蚂蚁个数（每一波蚂蚁有多少个）

% Alpha 表征信息素重要程度的参数

% Beta 表征启发式因子（费用）重要程度的参数

% Gamma 表征启发式因子（延时）重要程度的参数

% Tau 初始信息素矩阵

% Rho 信息素蒸发系数

% Q 信息素增加强度系数

%% 输出参数列表

% MRT 最优组播树（01邻接矩阵表示）

% EDGES 最优组播树所有的边

% cost 最优组播树的费用

%%

C=[inf,2,inf,inf,inf,3,inf,2,inf,inf,inf,inf,inf;

2,inf,inf,inf,1,inf,inf,inf,inf,inf,inf,inf,inf;

inf,inf,inf,inf,inf,inf,2,inf,2,2,inf,inf,inf;

inf,inf,inf,inf,1,3,1,inf,inf,inf,inf,inf,inf;

inf,1,inf,1,inf,4,inf,inf,inf,inf,inf,inf,inf;

3,inf,inf,3,4,inf,inf,inf,inf,2,inf,inf,3;

inf,inf,2,1,inf,inf,inf,inf,inf,2,inf,inf,inf;

2,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,2,inf;

inf,inf,2,inf,inf,inf,inf,inf,inf,inf,inf,2,inf;

inf,inf,2,inf,inf,2,2,inf,inf,inf,3,inf,inf;

inf,inf,inf,inf,inf,inf,inf,inf,inf,3,inf,inf,2;

inf,inf,inf,inf,inf,inf,inf,2,2,inf,inf,inf,2;

inf,inf,inf,inf,inf,3,inf,inf,inf,inf,2,2,inf;];

D=[inf,3,inf,inf,inf,1,inf,1,inf,inf,inf,inf,inf;

3,inf,inf,inf,4,inf,inf,inf,inf,inf,inf,inf,inf;

inf,inf,inf,inf,inf,inf,4,inf,5,4,inf,inf,inf;

inf,inf,inf,inf,3,1,3,inf,inf,inf,inf,inf,inf;

inf,4,inf,3,inf,1,inf,inf,inf,inf,inf,inf,inf;

1,inf,inf,1,1,inf,inf,inf,inf,4,inf,inf,3;

inf,inf,4,3,inf,inf,inf,inf,inf,3,inf,inf,inf;

1,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,1,inf;

inf,inf,5,inf,inf,inf,inf,inf,inf,inf,inf,1,inf;

inf,inf,4,inf,inf,4,3,inf,inf,inf,3,inf,inf;

inf,inf,inf,inf,inf,inf,inf,inf,inf,3,inf,inf,1;

inf,inf,inf,inf,inf,inf,inf,1,1,inf,inf,inf,3;

inf,inf,inf,inf,inf,3,inf,inf,inf,inf,1,3,inf;];

B=[inf,90,inf,inf,inf,80,inf,90,inf,inf,inf,inf,inf;

90,inf,inf,inf,100,inf,inf,inf,inf,inf,inf,inf,inf;

inf,inf,inf,inf,inf,inf,90,inf,80,90,inf,inf,inf;

inf,inf,inf,inf,110,80,100,inf,inf,inf,inf,inf,inf;

inf,100,inf,110,inf,80,inf,inf,inf,inf,inf,inf,inf;

80,inf,inf,80,80,inf,inf,inf,inf,90,inf,inf,45;

inf,inf,90,100,inf,inf,inf,inf,inf,50,inf,inf,inf;

90,inf,inf,inf,inf,inf,inf,inf,inf,inf,inf,90,inf;

inf,inf,80,inf,inf,inf,inf,inf,inf,inf,inf,90,inf;

inf,inf,90,inf,inf,90,50,inf,inf,inf,80,inf,inf;

inf,inf,inf,inf,inf,inf,inf,inf,inf,80,inf,inf,90;

inf,inf,inf,inf,inf,inf,inf,90,90,inf,inf,inf,90;

inf,inf,inf,inf,inf,45,inf,inf,inf,inf,90,90,inf;];

S=4;

E=[9,13,1];

Dmax=900;

K=200;

M=8;

Alpha=0.61;

Beta=0.41;

Gamma=0.8;

Tau=ones(13,13);

Rho=0.5*ones(13,13);

Q=2;

%% 第一步：变量初始化

N=size(C,1);%网络节点个数为N

P=length(E);%目的节点个数为M

MRT=zeros(N,N);

cost=inf;

ROUTES=cell(P,K,M);%用细胞结构存储到每一个目的节点的每一代的每一只蚂蚁的爬行路线

DELAYS=inf*ones(P,K,M);%用三维数组存储每代每个蚂蚁爬行到各个目的节点的延时

COSTS=inf*ones(P,K,M);%用三维数组存储每代每个蚂蚁爬行到各个目的节点的费用

%% 第二步：启动到P个目的节点的K轮蚂蚁觅食活动，每轮派出M只蚂蚁

for p=1:P

Tau=ones(N,N);

for k=1:K

for m=1:M

%% 第三步：状态初始化

W=S;%当前节点初始化为起始点

Path=S;%爬行路线初始化

PD=0;%爬行路线延时初始化

PC=0;%爬行路线费用初始化

TABU=ones(1,N);%禁忌表初始化

TABU(S)=0;%S已经在初始点了，因此要排除

CC=C;%费用邻接矩阵备份

DD=D;%延时邻接矩阵备份

%% 第四步：下一步可以前往的节点

DW=DD(W,:);

DW1=find(DW<inf);

for j=1:length(DW1)

if TABU(DW1(j))==0

DW(j)=inf;

end

end

LJD=find(DW<inf);%可选节点集

Len_LJD=length(LJD);%可选节点的个数

%% 觅食停止条件：蚂蚁未遇到食物或者陷入死胡同

while (W~=E(p))&&(Len_LJD>=1)

%% 第五步：转轮赌法选择下一步怎么走

PP=zeros(1,Len_LJD);

for i=1:Len_LJD

PP(i)=(Tau(W,LJD(i))^Alpha)*(C(W,LJD(i))^Beta)*(D(W,LJD(i))^Gamma);

end

PP=PP/(sum(PP));%建立概率分布

Pcum=cumsum(PP);

Select=find(Pcum>=rand);

to_visit=LJD(Select(1));%下一步将要前往的节点

%% 第六步：状态更新和记录

Path=[Path,to_visit];%路径增加

PD=PD+DD(W,to_visit);%路径延时累计

PC=PC+CC(W,to_visit);%路径费用累计

W=to_visit;%蚂蚁移到下一个节点

for kk=1:N

if TABU(kk)==0

CC(W,kk)=inf;

CC(kk,W)=inf;

DD(W,kk)=inf;

DD(kk,W)=inf;

end

end

TABU(W)=0;%已访问过的节点从禁忌表中删除

DW=DD(W,:);

DW1=find(DW<inf);

for j=1:length(DW1)

if TABU(DW1(j))==0

DW(j)=inf;

end

end

LJD=find(DW<inf);%可选节点集

Len_LJD=length(LJD);%可选节点的个数

%%

end

%% 第七步：记下每一代每一只蚂蚁的觅食路线和路线长度

ROUTES{p,k,m}=Path;

if Path(end)==E(p)&&PD<=Dmax

DELAYS(p,k,m)=PD;

COSTS(p,k,m)=PC;

else

DELAYS(p,k,m)=inf;

COSTS(p,k,m)=inf;

end

end

%% 第八步：更新信息素

Delta_Tau=zeros(N,N);%更新量初始化

for m=1:M

if COSTS(p,k,m)<inf&&DELAYS(p,k,m)<Dmax

ROUT=ROUTES{p,k,m};

TS=length(ROUT)-1;%跳数

Cpkm=COSTS(p,k,m);

for s=1:TS

x=ROUT(s);

y=ROUT(s+1);

Delta_Tau(x,y)=Delta_Tau(x,y)+Q/Cpkm;

Delta_Tau(y,x)=Delta_Tau(y,x)+Q/Cpkm;

end

end

end

Tau=(1-Rho).*Tau+(1-Rho).*Delta_Tau;%信息素挥发一部分，新增加一部分

end

end

%% 第九步：整理输出结果

MINCOSTS=NaN*ones(1,K);

allcost=zeros(1,0);

for k=1:K

for m=1:M

COSTkm=COSTS(:,k,m);

DELAYkm=DELAYS(:,k,m);

if sum(COSTkm)<inf&&sum(DELAYkm)<inf

Tree=zeros(N,N);

for p=1:P

path=ROUTES{p,k,m};

RLen=length(path);

for i=1:(RLen-1)

Tree(path(i),path(i+1))=1;

Tree(path(i+1),path(i))=1;

end

end

TC=Tree.*C;

for ii=1:N

for jj=1:N

if isnan(TC(ii,jj))

TC(ii,jj)=0;

end

end

end

mincost=0.5*sum(sum(TC));

if mincost<cost

MINCOSTS(1,k)=mincost;

MRT=Tree;

cost=mincost;

end

allcost=[allcost,cost];

end

end

end

MM=triu(MRT);

T1=find(MM==1);

T2=ceil(T1/N);

T3=mod(T1,N);

EDGES=[T3,T2];

%% 绘收敛曲线

figure(1)

COSTS2=zeros(M,K,P);

DELAYS2=zeros(M,K,P);

for p=1:P

for k=1:K

for m=1:M

if COSTS(p,k,m)<inf

COSTS2(m,k,p)=COSTS(p,k,m);

DELAYS2(m,k,p)=DELAYS(p,k,m);

end

end

end

end

LC1=zeros(1,K);

LC2=zeros(1,K);

for k=1:K

costs=COSTS2(:,k,1);

delays=DELAYS2(:,k,1);

pos1=find(costs>0);

pos2=find(delays>0);

len1=length(pos1);

len2=length(pos2);

LC1(k)=sum(costs)/len1;

LC2(k)=sum(delays)/len2;

end

plot(LC1,'ko-');

hold on

plot(LC2,'bs-');

legend('费用','延时')

title('到达路径9的费用延时变化情况')

figure(2)

COSTS2=zeros(M,K,P);

DELAYS2=zeros(M,K,P);

for p=1:P

for k=1:K

for m=1:M

if COSTS(p,k,m)<inf

COSTS2(m,k,p)=COSTS(p,k,m);

DELAYS2(m,k,p)=DELAYS(p,k,m);

end

end

end

end

LC1=zeros(1,K);

LC2=zeros(1,K);

for k=1:K

costs=COSTS2(:,k,2);

delays=DELAYS2(:,k,2);

pos1=find(costs>0);

pos2=find(delays>0);

len1=length(pos1);

len2=length(pos2);

LC1(k)=sum(costs)/len1;

LC2(k)=sum(delays)/len2;

end

plot(LC1,'ko-');

hold on

plot(LC2,'bs-');

legend('费用','延时')

title('到达路径13的费用延时变化情况')

figure(3)

COSTS2=zeros(M,K,P);

DELAYS2=zeros(M,K,P);

for p=1:P

for k=1:K

for m=1:M

if COSTS(p,k,m)<inf

COSTS2(m,k,p)=COSTS(p,k,m);

DELAYS2(m,k,p)=DELAYS(p,k,m);

end

end

end

end

LC1=zeros(1,K);

LC2=zeros(1,K);

for k=1:K

costs=COSTS2(:,k,3);

delays=DELAYS2(:,k,3);

pos1=find(costs>0);

pos2=find(delays>0);

len1=length(pos1);

len2=length(pos2);

LC1(k)=sum(costs)/len1;

LC2(k)=sum(delays)/len2;

end

plot(LC1,'ko-');

hold on

plot(LC2,'bs-');

legend('费用','延时')

title('到达路径1的费用延时变化情况')

figure(4)

plot(allcost,'b-')

title('组播树费用收敛曲线')