QoS组播路由蚁群算法

matlab的旋律

具体算法程序如下：
function [MRT,EDGES,cost,allcost,LC1,LC2]=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         最优组播树的费用
%%
%% 第一步：变量初始化
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) + 1e-10;
                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+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];
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

对寻优参数进行加权计算：
function [MRT,EDGES,cost,allcost,LC1,LC2]= MACA_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         最优组播树的费用
%%
%% 第一步：变量初始化
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) + 1e-10;
                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+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];

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

调用函数的脚本程序如下：
clear
close all
clc

N = 20;
C=rand(N);
D=rand(N);
S=3;
E=[6,1];
Dmax=10;
K = 6;
M=3;
Alpha=0.23;
Beta=0.4;
Gamma=0.8;
Tau=ones(N);
Rho=ones(N);
Q=6;
[MRT,EDGES,cost,allcost,LC1,LC2] = ACA_QoS_MR(C,D,S,E,Dmax,K,M,Alpha,Beta,Gamma,Tau,Rho,Q);

figure(1)
subplot(2,1,1)
plot(LC1,'ko-');
hold on
plot(LC2,'bs-');
legend('费用','延时')
title('路径的费用延时变化情况')

figure(2)
subplot(2,1,1)
plot(allcost,'b-')
title('组播树费用收敛曲线')
[MMRT,MEDGES,Mcost,Mallcost,MLC1,MLC2]= MACA_QoS_MR(C,D,S,E,Dmax,K,M,Alpha,Beta,Gamma,Tau,Rho,Q);


figure(1)
subplot(2,1,2)
plot(MLC1,'ro-');
hold on
plot(MLC2,'gs-');
legend('费用','延时')
title('改进算法路径的费用延时变化情况')

figure(2)
subplot(2,1,2)
plot(Mallcost,'r-')
title('改进算法组播树费用收敛曲线')

程序运行生成的图形由于参数的随机性，每次都有变化，其中一次结果如下：