function [MRT,EDGES,mincost]=ACA_QoS_MR(C,Cc,D,Dd,Djj,Dj,Pl,B,S,E,K,M,Alpha,Beta,Rho,Q,PL,BC,DJ,DC)
%% Ant Colony Algorithm for QoS Multicast Routing
% QoS组播路由蚁群算法
% editor:zhuhuazheng
%% 输入参数列表
% C 费用邻接矩阵（N×N）
% D 延时邻接矩阵（N×N）
% S 源节点
% E 组播目的节点（行向量）
% DC 延时约束
% K 迭代次数（指蚂蚁出动多少波）
% M 蚂蚁个数（每一波蚂蚁有多少个）
% Alpha 表征信息素重要程度的参数
% Beta 表征启发式因子重要程度的参数
% Tau 初始信息素矩阵
% Rho 信息素蒸发系数
% Q 信息素增加强度系数
% C=[ 4,2,3,inf,9,inf,inf,inf;%两点之间的费用
% 2,1,6,4,inf,inf,inf,4;
% 3,6,2,8,3,8,inf,inf;
% inf,4,8,9,inf,1,inf,7;
% 9,inf,3,inf,11,21,inf,inf;
% inf,inf,8,1,21,5,10,inf;
% inf,inf,inf,inf,inf,10,7,12;
% inf,4,inf,7,inf,inf,12,3];
%
% Cc=diag(C);%每个节点的费用
% D=[1,9,3,inf,18,inf,inf,inf;%两点之间的延迟
% 9,5,7,6,inf,inf,inf,12;
% 3,7,3,10,3,6,inf,inf;
% inf,6,10,2,inf,9,inf,9;
% 18,inf,3,inf,1,8,inf,inf;
% inf,inf,6,9,8,11,4,inf;
% inf,inf,inf,inf,inf,4,2,12;
% inf,12,inf,9,inf,inf,12,7];
%
% Dd=diag(D);%每个节点的延迟
% Dj=[1,0,0,inf,3,inf,inf,inf;%两点之间的延迟抖动
% 0,1,2,0,inf,inf,inf,2;
% 0,2,1,2,1,1,inf,inf;
% inf,0,2,0,inf,1,inf,3;
% 3,inf,1,inf,1,1,inf,inf;
% inf,inf,1,1,1,1,1,inf;
% inf,inf,inf,inf,inf,1,0,2;
% inf,2,inf,3,inf,inf,2,3];
%
% Djj=diag(Dj);%每个节点的延迟抖动
% Pl=[1e-5,1e-4,1e-6,1e-7,1e-6,1e-4,1e-6,1e-3];%每个节点的丢包率
% B=[inf,80,80,inf,100,inf,inf,inf;%两点之间的带宽
% 80,inf,90,70,inf,inf,inf,120;
% 80,90,inf,75,90,30,inf,inf;
% inf,70,75,inf,inf,120,inf,40;
% 100,inf,90,inf,inf,110,inf,inf;
% inf,inf,30,120,110,inf,130,inf;
% inf,inf,inf,inf,inf,130,inf,80;
% inf,120,inf,40,inf,inf,80,inf];
%
% S=1;
% E=[2,4,5,7];
% DC=50;
% BC=70;
% DJ=6;
% PL=0.001;
% K=50;
% M=20;
% Alpha=1;
% Beta=1;
% Rho=0.5;
% Q=100;
%% 输出参数列表
% MRT 最优组播树（01邻接矩阵表示）
% EDGES 最优组播树所有的边
% cost 最优组播树的费用
%%
%% 第一步：变量初始化
N=size(C,1);%网络节点个数为N
%Net=ones(N,N);%网络节点初始化
P=length(E);%目的节点个数为M
MRT=zeros(N,N);
mincost=inf;
ROUTES=cell(P,K,M);%用细胞结构存储到每一个目的节点的每一代的每一只蚂蚁的爬行路线
DELAYS=inf*ones(P,K,M);%用三维数组存储每代每个蚂蚁爬行到各个目的节点的延时
COSTS=inf*ones(P,K,M);%用三维数组存储每代每个蚂蚁爬行到各个目的节点的费用
DELAYJ=inf*ones(P,K,M);
PLOSS=inf*ones(P,K,M);
%% 第二步：启动到P个目的节点的K轮蚂蚁觅食活动，每轮派出M只蚂蚁
for p=1:P
Tau=ones(N,N);%初始化信息素
for k=1:K %出动K波
k
%% 检查所有节点的丢包率,删除与它连接的边
NMeet=find(Pl>=PL);
for i=1:length(NMeet)
D(NMeet(i),:)=inf;
D(:,NMeet(i))=inf;
Dj(NMeet(i),:)=inf;
Dj(:,NMeet(i))=inf;
end
%% 检查所有边的带宽，删除不满足带宽要求的边
NMetB=find(B<BC);
H1=rem(NMetB,N);
H2=ceil(NMetB/N);
H3=find(H1==0);
if ~isempty(H3)
H1(H3)=N;
end
for i=1:length(H1)
D(H1(i),H2(i))=inf;
D(H2(i),H1(i))=inf;
Dj(H1(i),H2(i))=inf;
Dj(H2(i),H1(i))=inf;
end
%%
for m=1:M% M只蚂蚁搜索
%% 第三步：状态初始化
m
W=S;%当前节点初始化为起始点
Path=S;%爬行路线初始化
PD=Dd(S);%爬行路线延时初始化
PC=Cc(S);%爬行路线费用初始化
PDj=Djj(S);%爬行路线延迟抖动初始化
PPl=Pl(S);%丢包率初始化
TABU=ones(1,N);%禁忌表初始化
TABU(S)=0;%S已经在初始点了，因此要排除
CC=C;%费用邻接矩阵备份
CCc=Cc;%各节点费用备份
DD=D;%延时邻接矩阵备份
DDd=Dd;%各节点延迟备份
DDjj=Djj;%各节点延迟抖动备份
DDj=Dj;%延迟抖动矩阵备份
%% 第四步：下一步可以前往的节点
DW=DD(W,:);%取第W个点与其他节点的延迟
DW1=find(DW<inf);
for j=1:length(DW1)
if TABU(DW1(j))==0
DW(j)=inf;
end
end
LJD=find(DW<inf);%可选节点集
LJD(find(LJD==W))=[];
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)*((1/(D(W,LJD(i))+Dj(W,LJD(i))))^Beta);%转移概率
end
PP=PP/(sum(PP));%建立概率分布
Pcum=cumsum(PP);
Select=find(Pcum>=rand);
to_visit=LJD(Select(1));%下一步将要前往的节点
% [Select,to_visit]=max(PP);
%% 第六步：状态更新和记录
%判断是否满足限制条件总延迟DC，延迟抖动DJ
PPD=PD;
PPC=PC;
PPDj=PDj;
PPPl=PPl;
PD=PD+DD(W,to_visit)+DDd(to_visit);%路径延时累计
PC=PC+CC(W,to_visit)+CCc(to_visit);%路径费用累计
PDj=PDj+DDj(W,to_visit)+DDjj(to_visit);%路径延迟抖动累计
PPl=1-(1-PPl)*(1-Pl(to_visit));%路径丢包率抖动累计
count=1;
while PD>DC||PDj>DJ
%重新选择一个节点
Select=find(Pcum>=rand);
to_visit=LJD(Select(1));%下一步将要前往的节点
PD=PPD;
PC=PPC;
PDj=PPDj;
PPl=PPPl;
PD=PD+DD(W,to_visit)+DDd(to_visit);%路径延时累计
PC=PC+CC(W,to_visit)+CCc(to_visit);%路径费用累计
PDj=PDj+DDj(W,to_visit)+DDjj(to_visit);%路径延迟抖动累计
PPl=1-(1-PPl)*(1-Pl(to_visit));%路径丢包率抖动累计
count=count+1;
if count>2
break;
end
end
W=to_visit;%蚂蚁移到下一个节点
Path=[Path,to_visit];%路径增加
for kk=1:N
if TABU(kk)==0
DD(W,kk)=inf;
DD(kk,W)=inf;
DDd(kk)=inf;
DDj(W,kk)=inf;
DDj(kk,W)=inf;
DDjj(kk)=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);%可选节点集
LJD(find(LJD==W))=[];
Len_LJD=length(LJD);%可选节点的个数
%%
end
%% 第七步：记下每一代每一只蚂蚁的觅食路线和路线长度
ROUTES{p,k,m}=Path;
if Path(end)==E(p)&&PD<=DC&&PDj<=DJ
DELAYS(p,k,m)=PD;
DELAYJ(p,k,m)=PDj;
COSTS(p,k,m)=PC;
PLOSS(p,k,m)=PPl;
else
DELAYS(p,k,m)=inf;
COSTS(p,k,m)=inf;
DELAYJ(p,k,m)=inf;
PLOSS(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)<=DC&&DELAYJ(p,k,m)<=DJ&&PLOSS(p,k,m)<=PL
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/(DELAYS(p,k,m)+DELAYJ(p,k,m));
