Matlab 蚁群算法

clear;
clc;

G=textread（‘IEEE33.txt‘）;%支路号—节点号矩阵
alpha=1;%蚂蚁在运动过程中所积累信息（即信息素）在蚂蚁选择路径时的相对重要程度，alpha过大时，算法迭代到一定代数后将出现停滞现象
beta=5;%启发式因子在蚂蚁选择路径时的相对重要程度
rho=0.5;%0Q=100;%蚂蚁释放的信息素量，对本算法的性能影响不大
Tau=ones（137）*Q;%路径上的初始信息素
R=G（:4）‘;
eta=1./R;%对应LN中的支路能见度

Nc=0;
Nc_max=50;
m=25;
trip=zeros（m32）;%m只蚂蚁走过的m条路径
trip（:1）=1;
P=zeros（1m）;%同一次迭代中每只蚂蚁所走路径的网损值
Ploss=zeros（1Nc_max）;%每次迭代的网损值


PP=zeros（1m）;


while Nc    Nc=Nc+1    
    for i=1:m
        t=1;
        S=[1];%蚂蚁k连入树的节点集合
        W=[234567891011121314151617181920212223242526272829303132];%蚂蚁k连入树的节点集合
        E=[218];%两节点集合间所有可选路径
        while t<32
            t=t+1;
            long=length（E）;
            tao=zeros（1long）;
            eta1=zeros（1long）;
            for j=1:long 
                tao（j）=Tau（E（j））;%返回可选路径上的初始信息素
                eta1（j）=eta（E（j））;%返回可选路径的能见度
            end
            tao;
            eta1;
            he=sum（（tao.^alpha）.*（eta1.^beta））;%求和
            
%******************************按概率选择路径******************************&            
            q0=0.5; 
            probability=zeros（1long）;
            a=rand;
            if a>q0
                for k=1:long
                    probability（k）=tao（k）^alpha*eta1（k）^beta;
                end
                probability;
                for o=1:long  
                    chance=sort（probability‘ascend‘）;
                    Locate=find（probability==chance（end））;
                    l=E（Locate（1））;%若rand大于0.5，在E中选取概率最大的路径
                    N=G（l2:3）;%所选路径的首、末端节点编号
                    P=intersect（NS）;%N与已连入树的节点的公共节点
                    if P==N
                        Locate=find（E==l）;%E是存储数据的数组名，find是找到等于l的数的位置
                        E（Locate）=[];%删除数组E中等于l的元素
                        E=unique（E）;
                        probability（probability==chance（end））=[];
                    else
                        break
                    end
                end
                %修正3个集合
                if intersect（G（l2）S）==G（l2）
                    S=[SG（l3）];
                    Locate=find（W==G（l3））;%W是存储数据的数组名，find是找到等于G（l3）的数的位置
                    W（Locate）=[];%删除数组W中等于G（l3）的元素
                    if isempty（W）==1
                        trip（it）=l;%将已走过的路径存放在trip中
                        break
                    end
                    [rowcol]=find（G（:2:3）==G（l3））;%与节点j相关联的新增支路所在位置
                    E=[Erow‘];
                    Locate=find（E==l）;%E是存储数据的数组名，find是找到等于l的数的位置
                    E（Locate）=[];%删除数组a中等于l的元素
                    E=unique（E）;
                    trip（it）=l;%将已走过的路径存放在trip中                    
                else            
                    if intersect（G（l3）S）==G（l3）
                        S=[SG（l2）];
                        Locate=find（W==G（l2））;%W是存储数据的数组名，find是找到等于G（l2）的数的位