【机器学习笔记36】蚁群算法-聚类分析

【参考资料】
【1】《模式识别与智能计算的MATLAB技术实现》
【2】《蚁群聚类算法综述》

1 算法概述

聚类数已知的算法流程
  1. 初始化蚁群参数,如蚂蚁数量、聚类数量等;每只蚂蚁对应一个解集:
样品号12N
蚂蚁 S i S_i Si类别1类别2类别1

上述表表示蚂蚁 S i S_i Si把N个样本分归属到分类中。

  1. 构建信息素矩阵
样本\类别类别1类别2类别3类别4
样本10.10.10.10.1
样本20.10.10.10.1
样本30.10.10.10.1
0.10.10.10.1
样本N0.10.10.10.1

上述信息素 τ i j \tau_{ij} τij表示把第i个样本归类到第j个类别时的信息素

  1. 构建目标函数

假设: N个样品、M个模式分类 ∣ S j , j = 1 , 2 , . . . , M ∣ |S_j,j=1,2,..., M| Sj,j=1,2,...,M,每个样品有n个特征。
目标: 将每个样品到聚类中心的距离和达到最小

m i n J ( w , c ) = ∑ j = 1 m ∑ i = 1 N ∑ p = 1 n w i j ∣ ∣ x i p − c j p ∣ ∣ 2 minJ(w, c)=\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{N}\sum\limits_{p=1}^{n} w_{ij}||x_{ip} - c_{jp}||^2 minJ(w,c)=j=1mi=1Np=1nwijxipcjp2

c j p = ∑ i = 1 N w i j x i p ∑ i = 1 N w i j c_{jp}=\dfrac{\sum\limits_{i=1}^{N}w_{ij}x_{ip}}{\sum\limits_{i=1}^{N}w_{ij}} cjp=i=1Nwiji=1Nwijxip

w i j = { 1 , N i ∈ S j 0 , N i ∉ S j w_{ij} = \begin{cases} 1, & N_i \in S_j \\ 0, & N_i \notin S_j \end{cases} wij={1,0,NiSjNi/Sj样本属于该类则为1,否则为0

其中 x i p x_{ip} xip为第i个样本的p属性、 c j p c_{jp} cjp为第j个分类的p属性

  1. 更新蚁群

每一只蚂蚁在对自己解集中样本归属判断时采用两种策略(随机选择):

  1. 根据当前时刻信息素表,选择信息素高的;(直接基于既有知识)
  2. 按照当前信息素的概率,即信息素高的类别有更大的可能性被选择;(也是基于既有知识,但非直接,而是一定程度上的随机选择)

4.1 参考第一步每只蚂蚁所具备的解集;根据目标函数公式计算每只蚂蚁的目标值;
4.2 根据目标值将蚂蚁进行排序;
4.3 取最优的L只蚂蚁进行“局部搜索”,遍历这L只蚂蚁:

4.3.1 随机选择其中第i个样本,重新计算目标函数,确定一个新的分类;
4.3.2 当这只蚂蚁的若干个样本被重新分类后,再计算一次该蚂蚁的总目标函数,若解更优;则替换原解;
4.3.3 遍历L只蚂蚁后,前L只蚂蚁中具备最优解的作为当前全局最优解;

  1. 更新信息素表
    τ i j ( t + 1 ) = ( 1 − ρ ) τ i j ( t ) + ∑ s = 1 l Δ τ i j s \tau_{ij}(t+1)=(1-\rho)\tau_{ij}(t) + \sum\limits_{s=1}^{l}\Delta \tau_{ij}^s τij(t+1)=(1ρ)τij(t)+s=1lΔτijs
    其中若蚂蚁s中的样品i属于分类j,则 Δ τ i j s = Q J \Delta \tau_{ij}^s=\dfrac{Q}{J} Δτijs=JQ,否则为0。这里Q是一个超参数;J是蚂蚁s的目标函数值。

  2. 多次迭代后达到全局最优解

2 算法实现

实际在实现过程中感觉蚁群算法还是非常收到超参数影响。包括信息素的挥发参数,以及超参数Q等等,而且效果差距非常大。当前代码中的参数是一个具备比较好效果的设置。

# -*- coding: utf-8 -*-
import numpy  as np
import sklearn.datasets as ds
import matplotlib.pyplot as plt
import random
import math
import operator

SAMPLE_NUM  = 18    #样本数量
FEATURE_NUM = 2     #每个样本的特征数量
CLASS_NUM   = 2     #分类数量
ANT_NUM     = 200    #蚂蚁数量


"""
初始化测试样本,sample为样本,target_classify为目标分类结果用于对比算法效果
"""
sample, target_classify = ds.make_blobs(SAMPLE_NUM, n_features=FEATURE_NUM, centers=CLASS_NUM, random_state=3)

"""
信息素矩阵
"""
tao_array =  [[random.random() for col in range(FEATURE_NUM)] for row in range(SAMPLE_NUM)] 

"""
蚁群解集
"""
ant_array = [[0 for col in range(SAMPLE_NUM)] for row in range(ANT_NUM)] 

t_ant_array = [[0 for col in range(SAMPLE_NUM)] for row in range(ANT_NUM)] #存储局部搜索时的临时解

"""
聚类中心点
"""
center_array = [[0 for col in range(FEATURE_NUM)] for row in range(CLASS_NUM)] 

"""
当前轮次蚂蚁的目标函数值,前者是蚂蚁编号、后者是目标函数值
"""
ant_target = [(0, 0) for col in range(ANT_NUM)] 




change_q   = 0.3       #更新蚁群时的转换规则参数,表示何种比例直接根据信息素矩阵进行更新
L          = 2         #局部搜索的蚂蚁数量
change_jp  = 0.03      #局部搜索时该样本是否变动
change_rho = 0.02      #挥发参数
Q          = 0.1      #信息素浓度参数


def _init_test_data():

    """
    初始化蚁群解集,随机确认每只蚂蚁下每个样本的分类为1或者0
    """
    for i in range(0, ANT_NUM):
        for j in range(0, SAMPLE_NUM):

            tmp = random.randint(0, FEATURE_NUM - 1)

            ant_array[i][j] = tmp

    """
    将前两个样本作为聚类中心点的初始值
    """
    for i in range(0, CLASS_NUM):

        center_array[i][0] = sample[random.randint(0, SAMPLE_NUM-1)][0]
        center_array[i][1] = sample[random.randint(0, SAMPLE_NUM-1)][1]

def _get_best_class_by_tao_value(sampleid):

    max_value = np.max(tao_array[sampleid])

    for i in range(0, CLASS_NUM):

        if max_value == tao_array[sampleid][i]:

            return i


def random_pick(some_list, probabilities):
    x = random.uniform(0,1)
    cumulative_probability = 0.0
    for item, item_probability in zip(some_list, probabilities):
        cumulative_probability += item_probability
        if x < cumulative_probability:break
    return item


def _get_best_class_by_tao_probablity(sampleid):

    tarray = np.array(tao_array[sampleid])

    parray = tarray/np.sum(tarray)

    return random_pick([0,1], parray)


def _update_ant():

    """
    更新蚁群步骤
    """

    #产生一个随机数矩阵
    r = np.random.random((ANT_NUM, SAMPLE_NUM))

    for i in range(0, ANT_NUM):
        for j in range(0, SAMPLE_NUM):

            if r[i][j] > change_q:

                tmp_index = _get_best_class_by_tao_value(j)

                #选择该样本中信息素最高的做为分类
                ant_array[i][j] = tmp_index

            else:
                #计算概率值,根据概率的大小来确定一个选项
                tmp_index = _get_best_class_by_tao_probablity(j)

                ant_array[i][j] = tmp_index

    #print(ant_array[i])
    #1. 确定一个新的聚类中心
    f_value_feature_0 = 0
    f_value_feature_1 = 0

    for i in range(0, CLASS_NUM):

        
        f_num   = 0

        for j in range(0, ANT_NUM):

            for k in range(0, SAMPLE_NUM):

                if ant_array[j][k] == 0:

                    f_num   += 1
                    f_value_feature_0 += sample[k][0] #特征1

                else:

                    f_num   += 1
                    f_value_feature_1 += sample[k][1] #特征2

        if i == 0:
            center_array[i][0] = f_value_feature_0/f_num
        else:
            center_array[i][1] = f_value_feature_1/f_num
    
    

        #print(center_array[i], f_num)


def _judge_sample(sampleid):

    """
    计算与当前聚类点的举例,判断该sample应所属的归类
    """
    target_value_0 = 0
    target_value_1 = 0

    f1 = math.pow((sample[sampleid][0] - center_array[0][0]),2)
    f2 = math.pow((sample[sampleid][1] - center_array[0][1]),2)
    target_value_0 = math.sqrt(f1 + f2)


    f1 = math.pow((sample[sampleid][0] - center_array[1][0]),2)
    f2 = math.pow((sample[sampleid][1] - center_array[1][1]),2)
    target_value_1 = math.sqrt(f1 + f2)

    if target_value_0 > target_value_1:
        return 1
    else:
        return 0


def _local_search():

    """
    局部搜索逻辑
    """
    


    #2. 根据新的聚类中心计算每个蚂蚁的目标函数

    for i in range(0, ANT_NUM):

        target_value = 0

        for j in range(0, SAMPLE_NUM):

            if ant_array[i][j] == 0:

                #与分类0的聚类点计算距离
                f1 = math.pow((sample[j][0] - center_array[0][0]),2)
                f2 = math.pow((sample[j][1] - center_array[0][1]),2)
                target_value += math.sqrt(f1 + f2)

            else:
                #与分类1的聚类点计算距离
                f1 = math.pow((sample[j][0] - center_array[1][0]),2)
                f2 = math.pow((sample[j][1] - center_array[1][1]),2)
                target_value += math.sqrt(f1 + f2) 

        #保存蚂蚁i当前的目标函数
        ant_target[i] = (i, target_value)

    #3. 对全部蚂蚁的目标进行排序,选择最优的L只蚂蚁

    ant_target.sort(key= operator.itemgetter(1)) #对ant进行排序

    #4. 对这L只蚂蚁进行解的优化
    for i in range(0, L):

        ant_id = ant_target[i][0]

        target_value = 0

        for j in range(0, SAMPLE_NUM):

            #对于该蚂蚁解集中的每一个样本
            if random.random() < change_jp:

                #将该样本调整到与当前某个聚类点最近的位置
                t_ant_array[ant_id][j] = _judge_sample(j)


        #判断是否保留这个临时解
        for j in range(0, SAMPLE_NUM):

            if t_ant_array[ant_id][j] == 0:

                #与分类0的聚类点计算距离
                f1 = math.pow((sample[j][0] - center_array[0][0]),2)
                f2 = math.pow((sample[j][1] - center_array[0][1]),2)
                target_value += math.sqrt(f1 + f2)

            else:
                #与分类1的聚类点计算距离
                f1 = math.pow((sample[j][0] - center_array[1][0]),2)
                f2 = math.pow((sample[j][1] - center_array[1][1]),2)
                target_value += math.sqrt(f1 + f2) 


        if target_value < ant_target[i][1]:

            #更新最优解
            ant_array[ant_id] = t_ant_array[ant_id]


def _update_tau_array():

    """
    更新信息素表
    """
    for i in range(0, SAMPLE_NUM):

        for j in range(0, CLASS_NUM):

            tmp = tao_array[i][j] #当前的信息素

            tmp = (1 - change_rho) * tmp #处理信息素挥发

            J = 0

            #处理信息素浓度增加
            for k in range(0, ANT_NUM):

                if ant_array[k][i] == j:

                    f1 = math.pow((sample[i][0] - center_array[j][0]),2)
                    f2 = math.pow((sample[i][1] - center_array[j][1]),2)
                    J += math.sqrt(f1 + f2)

            if J != 0:
                tmp += Q/J 

                #print(tmp, Q/J)

            tao_array[i][j] = tmp


    #print(np.var(tao_array))






"""
说明:

简单蚁群算法解决聚类问题,参考笔记《蚁群算法-聚类算法》

作者:fredric

日期:2018-12-21

"""
if __name__ == "__main__":

    _init_test_data();

    for i in range(0, 100):

        print("iterate No. {} target {}".format(i, ant_target[0][1]))

        _update_ant()

        _local_search()

        _update_tau_array()



    #画出分类
    pre = ant_array[ant_target[0][0]]

    plt.figure(figsize=(5, 6), facecolor='w')
    plt.subplot(211)
    plt.title('origin classfication')
    plt.scatter(sample[:, 0], sample[:, 1], c=target_classify, s=20, edgecolors='none')

    plt.subplot(212)
    plt.title('ant classfication')
    plt.scatter(sample[:, 0], sample[:, 1], c=pre, s=20, edgecolors='none')

    plt.plot(center_array[0][0], center_array[0][1],'ro')
    plt.plot(center_array[1][0], center_array[1][1],'bo')

    plt.show()

在这里插入图片描述

  • 8
    点赞
  • 49
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
%-- Unknown date --% else p(:,j)=0; end; if maxp(1)<p(1,j) maxp(1)=p(1,j); end; linear_index=find(maxp(1)==p(1,:)); size1=[1,n]; [r_index,c_index]=ind2sub(size1,linear_index(1)); solution_medium(k,1)=distance(g(NC,k),c_index(1)); route(k,1)=c_index(1); tabu(k,c_index(1))=1; tao(g(NC,k),c_index(1))=(1-rou)*tao(g(NC,k),c_index(1))+rou*initao(g(NC,k),c_index(1));%local updating for the first itertion tao(c_index(1),g(NC,k))=(1-rou)*tao(c_index(1),g(NC,k)))+rou*initao(c_index(1),g(NC,k)); solution(k)=solution_medium(k,1); for s=2:(n-1) c_index(s)=0; r_index(s)=0; linear_index(s)=0; maxp(s)=0; for j=1:n if tabu(k,j)==0 psum_medium0(s,j)=(tao(route(k,s-1),j)^alpha).*(yita(route(k,s-1),j)^beta); else psum_medium0(s,j)=0; end; psum_medium=psum_medium0.'; psum(k,s)=sum(psum_medium(:,s)); for j=1:n if tabu(k,j)==0; p(s,j)=(tao(route(k,s-1),j)^alpha).*(yita(route(k,s-1),j)^beta)/psum(k,s); else p(s:(n-1),j)=0; end; if maxp(s)<p(s,j) maxp(s)=p(s,j); end; linear_index=find(maxp(s)==p); size2=[n-1,n]; [r_index(s),c_index(s)]=ind2sub(size2,linear_index(1)); solution_medium(k,s)=distance(c_index(s-1),c_index(s)); route(k,s)=c_index(s); tabu(k,c_index(s))=1; tao(c_index(s-1),c_index(s))=(1-rou)*tao(c_index(s-1),c_index(s))+rou*initao(c_index(s-1),c_index(s)); tao(c_index(s),c_index(s-1))=(1-rou)*tao(c_index(s),c_index(s-1)))+rou*initao(c_index(s),c_index(s-1)); solution(k)=solution_medium(k,s); end; tao(c_index(n-1),g(NC,k))=(1-rou)*tao(c_index(n-1),g(NC,k)))+rou*initao(c_index(n-1),g(NC,k)); tao(g(NC,k),c_index(n-1))=(1-rou)*tao(g(NC,k),c_index(n-1))+rou*initao(g(NC,k),c_index(n-1)); solution_medium(k,n)=distance(c_index(n-1),g(NC,k)); solution(k)=solution(k)+solution_medium(k,n); solution_NC(NC,k)=solution(k); besolution_NC(NC,k)=solution_NC(NC,1); end; for k=1:m if besolution_NC(NC)=solution_NC(NC,k) besolution_NC(NC)=solution_NC(NC,k); end; %% ******** In this phase globle updating occurs ********%% linear_index1=find(besolution_NC(NC)==solution_NC(NC,:)); size3=[NC,m]; [r_index1(NC),c_index1(NC)]=ind2sub(size3,linear_index(1)); bestroute_NC(NC,:)=route(c_index(NC),:); [aa,bb]=size(linaer_index1); for i=1:bb [r_index1_tao(NC,i),c_index1_t(NC,i)]=ind2sub(size3,linear_index(i)); detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))=Q/solution(c_index1_t(NC,i)); detatao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))=Q/solution(c_index1_t(NC,i)); tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))=(1-rou).* tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))+rou*detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1)); tao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))=(1-rou).* tao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))+rou*detatao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i))); detatao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))=Q/solution(c_index1_t(NC,i)); detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))=Q/solution(c_index1_t(NC,i)); tao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))=(1-rou).* tao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))+rou*detatao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i))); tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))=(1-rou).* tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))+rou*detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1)); for s=2:n-1 detatao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))=Q/solution(c_index1_t(NC,i)); detatao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))=Q/solution(c_index1_t(NC,i)); tao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))=(1-rou).*tao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))+rou.*detatao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1)); tao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))=(1-rou),*tao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))+rou.*detatao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s)); end; if bestsolution>bestsolution_NC(NC) bestsolution=bestsolution_NC(NC); end; linear_index2=find(bestsolution==bestsolution_NC); size4[1,NC]; [r_index2,c_index2]=ind2sub(size4,linear_index2(1)); bestroute(1,:)=bestroute_NC(c_index2,:); initao=tao; end; %% ******** Otuput the best rusult of TSP ******** %% for NC-1:NCmax bestroute_NC(NC,n)=g(NC,c_index1(NC)); t(NC)=NC; end; bestroute(1,n)=bestroute_NC(c_index2,n); polt(x(bestroute(n)),y(bestroute(n)),'*'); hold on; end; line([x(bestroute(n)),x(bestroute(1))],[y(bestroute(n)),y(bestroute(1))]); hold on; for j=1:(n-1) line([x(bestroute(j)),x(bestroute(j+1))],[y(bestroute(j)),y(bestroute(j+1))]); hold on; end; hold off; xlabel('X coordinate'); ylabel('Y coordinate'); title('best tour in NCmax iterations'); %% ********Funtion: Open and improt file of city coordinate in TSP %% ******** %% [fname,pname,fiterindex]=uigetfile('*.*','All Files(*.*)'); set(handles.text_filename,'string',filename); fid=fopen(filename,'r')' if fid==-1; warndlg('can't open the file','WARN'); fclose(fid); end; [A0,COUNT]=fscanf(fid,'%g'); n=COUNT/3; set(handles.edit_citysum,'string',n); fclose(fid); m=str2double(get(handles.edit_citysum,'string')); NCmax=str2double(get(handles.edit_ncmax,'string')); for NC=1:NCmax for k=1:m g(NC,k)=fix((n-1)*rand(1))+1; end; end %-- 10-9-1 下午4:54 --% %-- 10-9-1 下午8:09 --% load('E:\苏雪敬个人资料\苏雪敬-开题报告\程序代码\matlab.mat') %-- 10-9-2 上午9:01 --% load('E:\苏雪敬个人资料\苏雪敬-开题报告\程序代码\matlab.mat') %% ******Initialization phase ****** %% global A0 global n; %city number global g; m=str2double(get(handles.edit_antsum, 'String')); %set ant number by using Matlab GUI initao=str2double(get(handles.edit_tao, 'String')); alpha=str2double(get(handles.edit_alpha, 'String')); beta=str2double(get(handles.edit_beta, 'String')); Q=str2double(get(handles.edit_q, 'String')); rou=str2double(get(handles.edit_rou, 'String')); NCmax=str2double(get(handles.edit_ncmax, 'String')); A=reshape(A0,3,n); %get the city coordinates in TSP x=A(2,:); %get x=coordinate y=A(3,:); %get y=coordinate for i=1:n for j=1:n distance(i,j)=sqrt((x(i)-x(j)+(y(i)-y(j))*y(i)-y(j))); end; for i=1:n for j=1:n if j!=i tao(i,j)=initao; yita(i,j)=1/distance(i,j); end; initao=initao.*ones(n,n); detatao=zeors(n,n); bestsolution=10000000000000; %% ******This is in phase in which ants build tours ****** %% for NC=1:NCmax tabu=zeros(m.n)l for k=1:m tabu(k,g(NC,k))=1; maxp(1)=0; for j=1:n if tabu(k,j)==0 psum_medium0(1,j)=(tao(g(NC,k)^alpha).*(yita(g(NC,k),j)^beta); else psum_medium0(1,j)=0; end; psum_medium=psum_medium0.'; psum(k,1)=sum(psum_medium(:,1)); for j=1:n if tabu(k,j)==0 p(1,j)=(tao(g(NC,k)^alpha).*(yita(g(NC,k),j)^beta/psm(k,1); else p(:,j)=0; end; if maxp(1)<p(1,j) maxp(1)=p(1,j); end; linear_index=find(maxp(1)==p(1,:)); size1=[1,n]; [r_index,c_index]=ind2sub(size1,linear_index(1)); solution_medium(k,1)=distance(g(NC,k),c_index(1)); route(k,1)=c_index(1); tabu(k,c_index(1))=1; tao(g(NC,k),c_index(1))=(1-rou)*tao(g(NC,k),c_index(1))+rou*initao(g(NC,k),c_index(1));%local updating for the first itertion tao(c_index(1),g(NC,k))=(1-rou)*tao(c_index(1),g(NC,k)))+rou*initao(c_index(1),g(NC,k)); solution(k)=solution_medium(k,1); for s=2:(n-1) c_index(s)=0; r_index(s)=0; linear_index(s)=0; maxp(s)=0; for j=1:n if tabu(k,j)==0 psum_medium0(s,j)=(tao(route(k,s-1),j)^alpha).*(yita(route(k,s-1),j)^beta); else psum_medium0(s,j)=0; end; psum_medium=psum_medium0.'; psum(k,s)=sum(psum_medium(:,s)); for j=1:n if tabu(k,j)==0; p(s,j)=(tao(route(k,s-1),j)^alpha).*(yita(route(k,s-1),j)^beta)/psum(k,s); else p(s:(n-1),j)=0; end; if maxp(s)<p(s,j) maxp(s)=p(s,j); end; linear_index=find(maxp(s)==p); size2=[n-1,n]; [r_index(s),c_index(s)]=ind2sub(size2,linear_index(1)); solution_medium(k,s)=distance(c_index(s-1),c_index(s)); route(k,s)=c_index(s); tabu(k,c_index(s))=1; tao(c_index(s-1),c_index(s))=(1-rou)*tao(c_index(s-1),c_index(s))+rou*initao(c_index(s-1),c_index(s)); tao(c_index(s),c_index(s-1))=(1-rou)*tao(c_index(s),c_index(s-1)))+rou*initao(c_index(s),c_index(s-1)); solution(k)=solution_medium(k,s); end; tao(c_index(n-1),g(NC,k))=(1-rou)*tao(c_index(n-1),g(NC,k)))+rou*initao(c_index(n-1),g(NC,k)); tao(g(NC,k),c_index(n-1))=(1-rou)*tao(g(NC,k),c_index(n-1))+rou*initao(g(NC,k),c_index(n-1)); solution_medium(k,n)=distance(c_index(n-1),g(NC,k)); solution(k)=solution(k)+solution_medium(k,n); solution_NC(NC,k)=solution(k); besolution_NC(NC,k)=solution_NC(NC,1); end; for k=1:m if besolution_NC(NC)=solution_NC(NC,k) besolution_NC(NC)=solution_NC(NC,k); end; %% ******** In this phase globle updating occurs ********%% linear_index1=find(besolution_NC(NC)==solution_NC(NC,:)); size3=[NC,m]; [r_index1(NC),c_index1(NC)]=ind2sub(size3,linear_index(1)); bestroute_NC(NC,:)=route(c_index(NC),:); [aa,bb]=size(linaer_index1); for i=1:bb [r_index1_tao(NC,i),c_index1_t(NC,i)]=ind2sub(size3,linear_index(i)); detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))=Q/solution(c_index1_t(NC,i)); detatao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))=Q/solution(c_index1_t(NC,i)); tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))=(1-rou).* tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))+rou*detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1)); tao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))=(1-rou).* tao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))+rou*detatao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i))); detatao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))=Q/solution(c_index1_t(NC,i)); detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))=Q/solution(c_index1_t(NC,i)); tao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))=(1-rou).* tao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))+rou*detatao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i))); tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))=(1-rou).* tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))+rou*detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1)); for s=2:n-1 detatao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))=Q/solution(c_index1_t(NC,i)); detatao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))=Q/solution(c_index1_t(NC,i)); tao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))=(1-rou).*tao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))+rou.*detatao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1)); tao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))=(1-rou),*tao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))+rou.*detatao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s)); end; if bestsolution>bestsolution_NC(NC) bestsolution=bestsolution_NC(NC); end; linear_index2=find(bestsolution==bestsolution_NC); size4[1,NC]; [r_index2,c_index2]=ind2sub(size4,linear_index2(1)); bestroute(1,:)=bestroute_NC(c_index2,:); initao=tao; end; %% ******** Otuput the best rusult of TSP ******** %% for NC-1:NCmax bestroute_NC(NC,n)=g(NC,c_index1(NC)); t(NC)=NC; end; bestroute(1,n)=bestroute_NC(c_index2,n); polt(x(bestroute(n)),y(bestroute(n)),'*'); hold on; end; line([x(bestroute(n)),x(bestroute(1))],[y(bestroute(n)),y(bestroute(1))]); hold on; for j=1:(n-1) line([x(bestroute(j)),x(bestroute(j+1))],[y(bestroute(j)),y(bestroute(j+1))]); hold on; end; hold off; xlabel('X coordinate'); ylabel('Y coordinate'); title('best tour in NCmax iterations'); %% ********Funtion: Open and improt file of city coordinate in TSP %% ******** %% [fname,pname,fiterindex]=uigetfile('*.*','All Files(*.*)'); set(handles.text_filename,'string',filename); fid=fopen(filename,'r')' if fid==-1; warndlg('can't open the file','WARN'); fclose(fid); end; [A0,COUNT]=fscanf(fid,'%g'); n=COUNT/3; set(handles.edit_citysum,'string',n); fclose(fid); m=str2double(get(handles.edit_citysum,'string')); NCmax=str2double(get(handles.edit_ncmax,'string')); for NC=1:NCmax for k=1:m g(NC,k)=fix((n-1)*rand(1))+1; end; end %-- 10-9-2 下午4:47 --% %-- 10-9-2 下午8:02 --% %{ ? 当前蚂蚁移到邻近区域内的没有被其他蚂蚁占据的节点 a) 输入参数 ? 蚂蚁编号ant_no ? S局部查找范围 ? ant_matrix(蚂蚁的窗格矩阵) ? Z平面窗格总区域 b) 输出参数 ? 蚂蚁前往的结点坐标 %} function ant_matrix=ant_move(ant_no,S,ant_matrix,Z) %1、获得ant_no点坐标 %[x,y,z]=ant_matrix(ant_no,:); one_ant=ant_matrix(ant_no,:); x=one_ant(1); y=one_ant(2); z=one_ant(3); %2、获得上下界限 x_low=x-S/2; x_high=x+S/2; if(x_low<0) x_low=0; end if(x_high>Z) x_high=Z; end y_low=y-S/2; y_high=y+S/2; if(y_low<0) y_low=0; end if(y_high>Z) y_high=Z; end %获得所有邻接结点下标数组 %{ allneighbours=[]; [row,col]=size(ant_matrix); for i=1:row if i~=ant_no [x,y,z]=ant_matrix(i,:); if x>=x_low && x<=x_high && y>=y_low && y<=y_high allneighbours=[allneighbours i]; end end end allneighbours=[allneighbours ant_no]; %} [row,col]=size(ant_matrix); positions=[]; %遍历所有上下界之间的点 for i=x_low:x_high for j=y_low:y_high status=0; for k=1:row % [x,y,z]=ant_matrix(k,:); one_ant=ant_matrix(k,:); x=one_ant(1); y=one_ant(2); z=one_ant(3); if x==i && y==j status=1; end end if status==0 positions=[positions;i j]; end end end %随机产生这个位置 [row,col]=size(positions); i=ceil(rand*row); position=positions(i,:); ant_matrix(ant_no,1)=position(1); ant_matrix(ant_no,2)=position(2); %{ ? 求Oi点在平面窗格内S*S区域内的邻接点的下标矩阵 a) 输入参数: ? Oi点 ? S局部查找范围 ? Item_Window(结点的窗格矩阵) ? Z平面窗格总区域 b) 输出参数 ? 所有在S*S区域内的结点的下标 具体算法: 1、S区域的上界 2、S区域的下界 3、比较一下点在哪些区域内 %} function allneighbours=get_allneighbours(Oi,S,item_window,Z) %1、获得Oi点坐标 %[x,y]=item_window(Oi,:); one_item=item_window(Oi,:); x=one_item(1); y=one_item(2); %2、获得上下界限 x_low=x-S/2; x_high=x+S/2; if(x_low<0) x_low=0; end if(x_high>Z) x_high=Z; end y_low=y-S/2; y_high=y+S/2; if(y_low<0) y_low=0; end if(y_high>Z) y_high=Z; end %获得所有邻接结点下标数组 allneighbours=[]; [row,col]=size(item_window); for i=1:row if i~=Oi % [x,y]=item_window(i,:); one_item=item_window(i,:); x=one_item(1); y=one_item(2); if x>=x_low && x<=x_high && y>=y_low && y<=y_high allneighbours=[allneighbours i]; end end end %{ ? 求两点空间的距离函数 a) 输入参数: ? Oi点 ? Oj点 ? 点的空间矩阵Item_Space b) 输出参数: ? 两点之间的距离 %} function distance=get_distance(Oi,Oj,item_space) distance=0.0; X1=item_space(Oi,:); X2=item_space(Oj,:); size=length(X1); for i=1:size distance=distance+(X1(i)-X2(i))^2; end distance=sqrt(distance); %{ ? 设计函数计算群体相似度计算 a) 输入参数: ? Oi点(需要计算的点)、 ? S(S*S区域内)、 ? Alpha、 ? Item_Window(结点的窗格矩阵) ? Z(Z*Z区域) ? item_space(结点的空间矩阵) b) 输出参数: ? Oi点的相似度 程序流程: 1、计算Oi的所有邻居 2、所有邻居的距离相应的和 3、判断和值和0的关系 4、最后决定值的大小 %} function fi=get_Fi(Oi,S,Alpha,item_window,Z,item_space) %1、计算Oi的所有邻居 allneighbours=get_allneighbours(Oi,S,item_window,Z); %2、所有邻居的距离相应的和 len=length(allneighbours); fi=0; for i=1:len distance=get_distance(Oi,allneighbours(i),item_space); fi=fi+(1-distance)/Alpha; end if(fi>0) fi=fi/(S^2); else fi=0; end %{ ? 放下概率的计算 a) 输入参数: ? Oi点(需要计算的点)、 ? S(S*S区域内)、 ? Alpha、 ? Item_Window(结点的窗格矩阵) ? K2 ? Z(Z*Z区域) ? item_space(结点的空间矩阵) b) 输出参数 ? 放下概率 程序流程 1、先计算群体相似概率 2、计算放下概率 %} function Pd=get_Pd(Oi,S,Alpha,item_window,K2,Z,item_space) %1、计算群体相似概率 fi=get_Fi(Oi,S,Alpha,item_window,Z,item_space); %2、计算放下概率 if fi %{ ? 拾起概率的计算 a) 输入参数: ? Oi点(需要计算的点)、 ? S(S*S区域内)、 ? Alpha、 ? Item_Window(结点的窗格矩阵) ? K2 b) 输出参数 ? 拾起概率 程序流程: 1、首先计算群体相似概率 2、计算拾起概率 %} function Pp=get_Pp(Oi,S,Alpha,item_window,K1,Z,item_space) %1 fi=get_Fi(Oi,S,Alpha,item_window,Z,item_space); %2 Pp=(K1/(K1+fi))^2; %{ ? 判断蚂蚁所在处是否有点 a) 输入参数 ? x蚂蚁横坐标 ? y蚂蚁纵坐标 ? item_window所有点的坐标矩阵 b) 输出参数 ? 蚂蚁所在处是否有点,有点返回这个点在item_window中的行号,无点0 %} function position=has_item(x,y,item_window) %{ [row,col]=size(item_window); for i=1:row end %} %position=find(item_window(find(item_window(:,1)==x),2)==y) position=0; [row,col]=size(item_window); for i=1:row if ~isempty(find(x==item_window(i,1))) && ~isempty(find(y==item_window(i,2))) position=i; end end %{ ? 初始化函数 a) 输入参数 ? 蚂蚁的数目ant_number ? 点的数目item_number ? 空间尺寸Z,空间大小Z*Z b) 输出参数 ? 蚂蚁的平面窗格矩阵 ? 点的平面窗格矩阵 %} function [ant_matrix,item_window]=initialize(ant_number,item_number,Z) ant_matrix_x=randperm(Z); ant_matrix_y=randperm(Z); item_matrix_x=randperm(Z); item_matrix_y=randperm(Z); %生成蚂蚁矩阵 for i=1:ant_number ant_matrix(i,1)=ant_matrix_x(i); ant_matrix(i,2)=ant_matrix_y(i); ant_matrix(i,3)=0; end for i=1:item_number item_window(i,1)=item_matrix_x(i); item_window(i,2)=item_matrix_y(i); end %{ test 测试程序 整个程序流程 1、初始化: 2、迭代tmax次 3、所有的蚂蚁运动一次 4、产生一个0-1之间的随机数R 5、如果当前蚂蚁处于未负载状态,而且当前蚂蚁所在处的有点Oi 5.1、计算群体相似度f(Oi)和拾起概率Pp(Oi) 5.2、如果拾起概率Pp(Oi)》R 5.2.1、当前蚂蚁拾起点Oi(注意Oi在窗格中的位置是不断变动的) 5.3 5.2结束 6、如果条件5不成立,如果当前蚂蚁处于负载状态,持有点Oi,而且当前位置没有其他点 6.1计算群体相似度f(Oi)和放下概率Pd(Oi) 6.2如果放下概率Pd(Oi)》R 6.2.1放下节点Oi(注意点Oi在窗格中的位置是不断变动的) 6.2.2修改蚂蚁状态为卸载 6.3 6.2结束 7、5结束 8、当前蚂蚁移到邻近区域内的没有被其他蚂蚁占据的节点 9、所有的蚂蚁运动一次结束 10、迭代tmax次结束 输入参数: 1、ant_number 蚂蚁数 2、item_number 点数 3、Z 区域范围 4、tmax 最大迭代次数 5、S 邻接区域 6、Alpha 相异度参数 7、item_space 点的空间矩阵 8、K1 拾起概率参数 9、K2 放下概率参数 %} clear; clc; K1=0.1; K2=0.15; Alpha=0.15; %S=6; S=30; tmax=300; ant_number=20; Z=50; axis([0 Z 0 Z]); item_space=[ 5.1 3.5 1.4 0.2; 4.9 3.0 1.4 0.2; 4.7 3.2 1.3 0.2; 4.6 3.1 1.5 0.2; 5.0 3.6 1.4 0.2; 5.4 3.9 1.7 0.4; 4.6 3.4 1.4 0.3; 5.0 3.4 1.5 0.2; ]; [row,col]=size(item_space); item_number=row; %1、初始化: [ant_matrix,item_window]=initialize(ant_number,item_number,Z); item_window %2、迭代tmax次 for i=1:tmax %3、所有的蚂蚁运动一次 for j=1:ant_number %4、产生一个0-1之间的随机数R r=rand; %5、如果当前蚂蚁处于未负载状态,而且当前蚂蚁所在处的有点Oi %[x,y,z]=ant_matrix(j,:); one_ant=ant_matrix(j,:); x=one_ant(1); y=one_ant(2); z=one_ant(3); %rows1=find(item_window(:,1)==x); %rows2=find(item_window(:,2)==y); %Oi=find(item_window(find(item_window(:,1)==x),2)==y); Oi=has_item(x,y,item_window); if z==0 && Oi>0 %5.1、计算群体相似度f(Oi)和拾起概率Pp(Oi) %fi=get_Fi(Oi,S,Alpha,item_window,Z,item_space); Pp=get_Pp(Oi,S,Alpha,item_window,K1,Z,item_space); Pp; %5.2、如果拾起概率Pp(Oi)》R if Pp>=r %5.2.1、当前蚂蚁拾起点Oi(注意Oi在窗格中的位置是不断变动的) ant_matrix(j,3)=Oi; end elseif ant_matrix(j,3)>0 && Oi==0 %6、如果条件5不成立,如果当前蚂蚁处于负载状态,持有点Oi,而且当前位置没有其他点 %6.1计算群体相似度f(Oi)和放下概率Pd(Oi) Oi=ant_matrix(j,3);%代表当前蚂蚁所携带的点的行号 Pd=get_Pd(Oi,S,Alpha,item_window,K2,Z,item_space); Pd; %6.2如果放下概率Pd(Oi)》R if Pd>=r %6.2.1放下节点Oi(注意点Oi在窗格中的位置是不断变动的) item_window(Oi,1)=x; item_window(Oi,2)=y; ant_matrix(j,3)=0; end end %8、当前蚂蚁移到邻近区域内的没有被其他蚂蚁占据的节点 ant_matrix=ant_move(j,S,ant_matrix,Z); end end item_window axis([0 Z 0 Z]); plot(item_window(:,1),item_window(:,2),'--rs'); %-- 10-9-2 下午9:11 --% load('E:\苏雪敬个人资料\苏雪敬-开题报告\程序代码\matlab.mat') %-- 10-9-3 下午2:54 --% load('E:\苏雪敬个人资料\苏雪敬-开题报告\程序代码\lf蚁群算法.mat') matlab7

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值