群知能〜蟻コロニー最適化法の実装

蟻コロニー最適化法をCで実装します。
参考図書は「機械学習と深層学習(小高知宏,オーム社)」です。

群知能(Swarm Intelligence)とは

生物の集団が見せる知的な行動を模した機械学習手法です。
蟻コロニー最適化法の他に
・粒子群最適化法：群れ全体が効率的に餌を見つける挙動をシミュレート
・AFSA：魚の群れが見せる捕食や追尾などの行動特性をシミュレート
などがあります。

蟻コロニー最適化法とは

蟻が巣穴と餌場の間の最短経路を見つける挙動をシミュレートしたものです。

蟻は自分が歩いた道筋にフェロモンを残します。
他の蟻はそのフェロモンの跡を辿って移動しようとしますが、フェロモンは時間の経過と共に蒸発して消えてしまいます。フェロモンが蒸発する前に他の蟻が通れば追加上書きできます。
巣穴と餌場の距離が短いとフェロモンがすぐに上書きされるので、その道筋に多くの蟻が群がることになります。
この結果、最短経路に群れが誘導されるのです。

学習手続き

1.設定したすべての経路のフェロモン濃度を０に初期化する。
2.以下を適切な回数繰り返す
　2.1.すべての蟻を巣穴から餌場まで歩かせる
　 　 経路は、確率εでランダムに、確率(1-ε)でフェロモン濃度に応じて選択する。
　2.2.すべての経路のフェロモン濃度に定数ρ(0~1)を乗じて蒸発させる。
　2.3.各蟻の移動距離Lを求め、個体の選択した経路に対応するフェロモン濃度に、Q*(1/L)を加える。

コードの概要

巣穴と餌場の間に図のような９つの中間分岐点があるような道を考えます。
それぞれの経路の距離は、簡単のため5と1としています。
巣穴から餌場までの最短経路をたどる行動知識を獲得します。

コード解説１〜変数について

cost[][]:道筋の長さを格納
pheromone[][]:それぞれの道のフェロモン濃度を格納
ー１次元目は分岐の方向、２次元目はステップ数です

mstep[][]:実際に歩いた経路の記録
ー１次元目は蟻の個体番号、２次元目はステップ数です。０か１でそれぞれの経路の選択を表します。

コード解説２〜蟻を歩かせる

学習手続き2.1を実現するために、ε-greedy法により経路を選択します。
確率εでランダムに、確率(1-ε)でフェロモン濃度が高い道を選択する

if((rand1()<EPSILON)||(abs(pheromone[0][s]-pheromone[1][s])<1e-9)){
    mstep[m][s]=rand01(); //ランダムに行動
}
else{ //濃度が高い方を選択
    if(pheromone[0][s]>pheromone[1][s])mstep[m][s]=0;
    else mstep[m][s]=1; 
}

コード解説３〜フェロモンの更新

###1.フェロモンの蒸発　
一律にすべての道のフェロモンをRHO(ここでは0.8)倍に薄めます

for(int i=0;i<2;++i){
    for(int j=0;j<STEP;++j){
        pheromone[i][j]*=RHO; //すべてをRHO倍に薄めている
    }
}

###2.フェロモンの上書き
各個体の移動距離Lを求め、個体の選択した経路に対応するフェロモン濃度に、Q*(1/L)を加えます。

for(int m=0;m<NOA;++m){
    /*個体mの移動距離lmの計算*/
    lm=0;
    for(int i=0;i<STEP;++i){
        lm+=cost[mstep[m][i]][i];
    }
    /*フェロモンの上塗り*/
    for(int i=0;i<STEP;++i){
        pheromone[mstep[m][i]][i]+=Q*(1.0/lm);
    }
    sum_lm+=lm;
}
printf("Average Distance=%lf\n\n",sum_lm/NOA); //蟻の歩いた平均距離を出力

結果

最初フェロモン濃度を０に初期化し、経路もランダムに選んでいたため、最初は平均移動距離が
25.8(1試行目), 28.0(2試行目)と大きくなっています。

学習が進むにつれ、どんどん小さくなり、(8試行目：22.4, 9試行目:22.8)

19試行目で16.8となり、かなり効率よく移動できるようになったことがわかります。

全コード

aco.c

#include<stdio.h>
#include<stdlib.h>

#define _CRT_SECURE_NO_WARNINGS

#define NOA 10       //蟻の個体数
#define ILIMIT 20    //繰り返しの回数
#define Q 3          //フェロモン更新の定数
#define RHO 0.8      //蒸発の定数
#define STEP 10      //道のりのステップ数
#define EPSILON 0.15 //行動選択のランダム性を決定
#define SEED 31      //乱数のシード

void printp(double pheromone[2][STEP]); //表示
void printmstep(int mstep[NOA][STEP]);  //蟻の挙動の表示
void walk(int cost[2][STEP],double phenomenone[2][STEP],
        int mstep[NOA][STEP]);  //蟻を歩かせる
void update(int cost[2][STEP],double pheromone[2][STEP],
        int mstep[NOA][STEP]);  //フェロモンの更新
double rand1();
int rand01(); 

/**************************************************************/
int main(){
    int cost[2][STEP]={ //各ステップのコスト(距離が短い方を１、長い方を５とする)
        {1,1,1,1,1,1,1,1,1,1},
        {5,5,5,5,5,5,5,5,5,5}
    };
    double pheromone[2][STEP]={0}; //各ステップのフェロモン量
    int mstep[NOA][STEP]={0};      //蟻が歩いた過程
    int i;                         //繰り返し回数の制御

    srand(SEED);

    for(i=0;i<ILIMIT;++i){
        /*フェロモンの状態出力*/
        printf("%d:\n",i); 
        printp(pheromone);
        /*蟻を歩かせる*/
        walk(cost,pheromone,mstep);
        /*フェロモンの更新*/
        update(cost,pheromone,mstep);
    }
    /*フェロモンの状態出力*/
    printf("%d:\n",i); 
    printp(pheromone);

    return 0;
}

/****************************************************************/
/*フェロモンの更新*/
void update(int cost[2][STEP],double pheromone[2][STEP],int mstep[NOA][STEP]){
    int m;           //蟻の個体番号
    int lm;          //蟻の歩いた距離
    double sum_lm=0; //蟻の歩いた距離の合計

    /*フェロモンの蒸発*/
    for(int i=0;i<2;++i){
        for(int j=0;j<STEP;++j){
            pheromone[i][j]*=RHO; //すべてを0.8倍に薄めている
        }
    }

    /*蟻による上塗り*/
    for(int m=0;m<NOA;++m){
        /*個体mの移動距離lmの計算*/
        lm=0;
        for(int i=0;i<STEP;++i){
            lm+=cost[mstep[m][i]][i];
        }
        /*フェロモンの上塗り*/
        for(int i=0;i<STEP;++i){
            pheromone[mstep[m][i]][i]+=Q*(1.0/lm);
        }
        sum_lm+=lm;
    }
    printf("Average Distance=%lf\n\n",sum_lm/NOA); //蟻の歩いた平均距離を出力
}

/*蟻を歩かせる*/
void walk(int cost[2][STEP],double pheromone[2][STEP],int mstep[NOA][STEP]){
    int m; //蟻の個体番号
    int s; //蟻の現在ステップ位置

    for(m=0;m<NOA;++m){
        for(s=0;s<STEP;++s){
            /*ε-greedyによる行動選択*/
            if((rand1()<EPSILON)||(abs(pheromone[0][s]-pheromone[1][s])<1e-9)){
                mstep[m][s]=rand01(); //ランダムに行動
            }
            else{ //濃度が高い方を選択
                if(pheromone[0][s]>pheromone[1][s])mstep[m][s]=0;
                else mstep[m][s]=1; 
            }
        }
    } 
    printmstep(mstep);//蟻の挙動の出力
}

/*蟻の挙動の表示*/
void printmstep(int mstep[NOA][STEP]){
    printf("*mstep\n");
    for(int i=0;i<NOA;++i){
        for(int j=0;j<STEP;++j){
            printf("%d",mstep[i][j]);
        }
        printf("\n");
    }
}

/*フェロモンの表示*/
void printp(double pheromone[2][STEP]){
    for(int i=0;i<2;++i){
        for(int j=0;j<STEP;++j){
            printf("%4.2lf ",pheromone[i][j]);
        }
        printf("\n");
    }
}

/*乱数生成*/
double rand1(){
    return (double)rand()/RAND_MAX;
}

int rand01(){
    int rnd;
    while((rnd=rand())==RAND_MAX);
    return (int)((double)rnd/RAND_MAX*2);
}