進化的手法〜遺伝的アルゴリズムでナップサック問題を解く

遺伝的アルゴリズム(GA)により、ナップサック問題を解くプログラムをCで実装します。
データセットはテキストファイルにして、それを読み込みます。

進化的手法とは
--
生物の進化を参考に考えられた手法で、その代表が遺伝的アルゴリズム(GA)です。
まず、1または0の並び（＝染色体）でパターンの表現をします。
そして、その染色体を評価・淘汰・生成を繰り返して、より良い遺伝子を持つ染色体が増えていきます。

（ちなみに、記号列による染色体のデータ表現を「遺伝子型」、染色体を解釈して具体的な情報に戻したデータ表現を「表現型」といいます。）

染色体に対して遺伝的操作を加えることで、より良い染色体を生み出していきます。
遺伝的操作には以下のようなものがあります。
・交叉：親の遺伝子を組み合わせて子供の染色体をつくる。
　（例）0000と1111が親なら子は0011と1100
・突然変異：遺伝子をランダムに反転(0と1を逆にする)して、染色体のある部分(=遺伝子座)を書き換える。
　（例）0000が突如0100になる
・選択：染色体の価値の評価に基づいて、選り分ける。評価値の高い染色体を選んで後世代に残す。

###遺伝的アルゴリズムの手順
1.染色体の集合を乱数で初期化する
2.以下を繰り返す
　2.1.以下を適切な回数繰り返し、次の世代の染色体を適当な個数生み出す
　　・親を評価値を利用して選択
　　・交叉により子を生み出す
　2.2.以下を適切な回数繰り返し、突然変異を起こす
　　・ある確率で遺伝子座を選ぶ
　　・突然変異
　2.3.親世代と同じ数だけ、子世代の染色体を選択

ナップサック問題
--
ある容量のナップサックと，数種類の品物(それぞれ価格と容量を持つ)が与えられ，ナップサックの容量の範囲内でいくつかの品物を詰めようとする時，ナップザックに入れる品物の価格の和が最大となる組み合わせを求める、という問題です。

ポイントは、
最大容量の範囲に収められる中で、最も評価値、つまり価格の総和が大きくなるような選び方をすること。
そして、最大容量を超えてしまうと、評価値は0になってしまうことです。

コード解説
--
###変数について
POOLSIZEを染色体の個数,Nを品物の個数として定義します。
そして、染色体の集団は以下のように表します。

```
int pool[POOLSIZE][N];     //染色体プール
int ngpool[POOLSIZE*2][N]; //次世代染色体プール(次世代の候補)
```
ngpoolの１次元目が２倍なのは、子世代を多めに作っておき、選択で良いものを選ぶためです。


###交叉
以下のようにして、ランダムに交叉点を決めてそこを境に反転させることで、一点交叉を実現します。
（c1,c2は子、mとpは親です）

```
  cp=rndn(N);  //交叉点の決定
  for(j=0;j<cp;++j){
    c1[j]=m[j];
    c2[j]=p[j];
  }
  for(;j<N;++j){
    c2[j]=m[j];
    c1[j]=p[j];
  }
```

###選択（ルーレット選択）
染色体の評価値に比例した面積が割り当てられたルーレットを用いて、選択をします。
これにより、評価値の高い染色体が選択される可能性を高めつつ、評価の低い染色体が選ばれる可能性も残すことができます。
ルーレットを以下のように作成します。関数evalfitで評価値を計算し、rouletteに代入します。評価値の総和はtotalfitnessに代入され、各染色体が選択される確率は(roulette/totalfitness)になります。
そして染色体を選択後、それを子世代に引き継ぎます。

```c
    /*ルーレットの作成*/
    totalfitness=0;
    for(c=0;c<POOLSIZE*2;++c){
      roulette[c]=evalfit(ngpool[c]);
      totalfitness+=roulette[c]; //適応度の合計値を計算
    }
    /*染色体を一つ選ぶ*/
    ball=rndn(totalfitness);
    acc=0;
    for(c=0;c<POOLSIZE*2;++c){
      acc+=roulette[c];  //適応度を積算
      if(acc>ball)break; //対応する遺伝子
    }
    /*染色体のコピー*/
    for(j=0;j<N;++j){
      pool[i][j]=ngpool[c][j];
 
```

###突然変異
以下のようにして、確率MRATEで突然変異をおこし、遺伝子を反転させます。

```
 if((double)rndn(100)/100.0<=MRATE){   //反転させる
    if(ngpool[i][j]==1)ngpool[i][j]=0;    
    else if(ngpool[i][j]==0)ngpool[i][j]=1;
 }
```

結果
--
最初はランダムに選ぶため、染色体集合の評価値の平均は低く第０世代は619.7でした
<img src=https://qiita-image-store.s3.amazonaws.com/0/184212/302af1fd-e1e1-c8a1-a585-841f86d2a303.png width=70%>
第20世代になると959.1まで上がります。
<img src=https://qiita-image-store.s3.amazonaws.com/0/184212/73ee63b1-0d4f-edda-187d-19b1dc43c56c.png width=38%>
最後第49世代になると990.1となりました。
<img src=https://qiita-image-store.s3.amazonaws.com/0/184212/4cc18afa-4c5c-5183-678d-a9399082c0d3.png width=38%>

コード全文
--
```c:ga.c
/*ナップサック問題のGAによる求解プログラム*/
#define _CRT_SECURE_NO_WARNINGS
#include<stdio.h>
#include<stdlib.h>
#include<limits.h>

#define N 30            //品物の個数
#define LIMIT 750 //ナップサックの容量制限
#define POOLSIZE 30     //プールサイズ
#define LASTG 50        //繰り返しを打ち切る世代
#define MRATE 0.01      //突然変異が起こる確率
#define SEED 33

void initparcel();       //荷物の初期化
int evalfit(int gene[]); //適応度の計算
void mating(int pool[POOLSIZE][N],int ngpool[POOLSIZE*2][N]);  //交叉
void mutation(int ngpool[POOLSIZE*2][N]);                      //突然変異
void printp(int pool[POOLSIZE][N]);                            //結果出力
void initpool(int pool[POOLSIZE][N]);                          //集団の初期化
int rndn();              //n未満の乱数の生成
int selectp(int roulette[POOLSIZE],int totalfitness);          //親の選択
void crossing(int m[],int p[],int c1[],int c2[]);              //特定の２染色体の交叉
void selectng(int ngpool[POOLSIZE*2][N],int pool[POOLSIZE][N]);//次世代の選択

int parcel[N][2];        //荷物データ

/**********************************************************************/
int main(int argc,char *argv[]){
 int pool[POOLSIZE][N];     //染色体プール
 int ngpool[POOLSIZE*2][N]; //次世代染色体プール
 int generation;            //現在の世代数

 srand(SEED);

 initparcel();   //荷物データの読み込み
 initpool(pool); //集団の初期化

 for(generation=0;generation<LASTG;++generation){
    printf("%d世代\n",generation);
    mating(pool,ngpool);   //交叉
    mutation(ngpool);      //突然変異
    selectng(ngpool,pool); //次世代の選択
    printp(pool);          //結果出力
 }
 return 0;
}
/************************************************************************/
/*荷物の初期化 データの読み取り*/
void initparcel(){
  int i=0;
  while((i<N)&&(scanf("%d %d",&parcel[i][0],&parcel[i][1])!=EOF))++i;
}

/*集団の初期化*/
void initpool(int pool[POOLSIZE][N]){
  for(int i=0;i<POOLSIZE;++i){
    for(int j=0;j<N;++j){
      pool[i][j]=rndn(2);
    }
  }
}

/*交叉*/
void mating(int pool[POOLSIZE][N],int ngpool[POOLSIZE*2][N]){
  int totalfitness=0;     //適応度の合計値
  int roulette[POOLSIZE]; //適応度を格納
  int mama,papa;          //親の遺伝子

  /*ルーレットの作成*/
  for(int i=0;i<POOLSIZE;++i){
    roulette[i]=evalfit(pool[i]);
    /*適応度の合計値を計算*/
    totalfitness+=roulette[i];
  }
  /*選択と交叉の繰り返し*/
  for(int i=0;i<POOLSIZE;++i){
    do{
      /*親の選択*/
      mama=selectp(roulette,totalfitness);
      papa=selectp(roulette,totalfitness);
    }while(mama==papa); //重複の削除

    /*特定の２染色体の交叉*/
    crossing(pool[mama],pool[papa],ngpool[i*2],ngpool[i*2+1]);
  }
}

/*突然変異*/
void mutation(int ngpool[POOLSIZE*2][N]){
  for(int i=0;i<POOLSIZE*2;++i){
    for(int j=0;j<N;++j){
      if((double)rndn(100)/100.0<=MRATE){   //反転させる
	if(ngpool[i][j]==1)ngpool[i][j]=0;    
	else if(ngpool[i][j]==0)ngpool[i][j]=1;
      }
    }
  }
}


/*次世代の選択*/
void selectng(int ngpool[POOLSIZE*2][N],int pool[POOLSIZE][N]){
  int i,j,c;
  int totalfitness=0;       //適応度の合計値
  int roulette[POOLSIZE*2]; //適応度を格納
  int ball;                 //玉（選択位置の数値）
  int acc=0;                //適応度の積算値

  /*選択を繰り返す*/
  for(i=0;i<POOLSIZE;++i){
    /*ルーレットの作成*/
    totalfitness=0;
    for(c=0;c<POOLSIZE*2;++c){
      roulette[c]=evalfit(ngpool[c]);
      totalfitness+=roulette[c]; //適応度の合計値を計算
    }
    /*染色体を一つ選ぶ*/
    ball=rndn(totalfitness);
    acc=0;
    for(c=0;c<POOLSIZE*2;++c){
      acc+=roulette[c];  //適応度を積算
      if(acc>ball)break; //対応する遺伝子
    }
    /*染色体のコピー*/
    for(j=0;j<N;++j){
      pool[i][j]=ngpool[c][j];
    }
  }
}

/*結果の出力*/
void printp(int pool[POOLSIZE][N]){
  int fitness;           //適応度
  double totalfitness=0; //適応度の合計値
  int elite,bestfit=0;   //エリート遺伝子の処理用変数

  for(int i=0;i<POOLSIZE;++i){
    for(int j=0;j<N;++j){
      printf("%1d",pool[i][j]);
    }
    fitness=evalfit(pool[i]);
    printf("\t%d\n",fitness);
    if(fitness>bestfit){ //エリート解
      bestfit=fitness;
      elite=i;
    }
    totalfitness+=fitness;
  }
  /*エリート解の適応度を出力*/
  printf("%d\t%d \t",elite,bestfit);
  /*平均の適応度を出力*/
  printf("%lf\n",totalfitness/POOLSIZE);
}

/*親の選択*/
int selectp(int roulette[POOLSIZE],int totalfitness){
  int ball;  //玉（選択位置）の制御
  int acc=0; //適応度の積算値
  int i;

  ball=rndn(totalfitness);
  for(i=0;i<POOLSIZE;++i){
    acc+=roulette[i];
    if(acc>ball)break; //対応する遺伝子
  }
  return i;
}

/*特定の２染色体の交叉*/
void crossing(int m[],int p[],int c1[], int c2[]){
  int j;
  int cp; //交叉する点

  /*交叉点の決定*/
  cp=rndn(N);

  /*前半部のコピー*/
  for(j=0;j<cp;++j){
    c1[j]=m[j];
    c2[j]=p[j];
  }
  for(;j<N;++j){
    c2[j]=m[j];
    c1[j]=p[j];
  }
}

/*適応度の計算*/
int evalfit(int g[]){
  int pos;      //遺伝子座の指定
  int value=0;  //評価値
  int capa=0;   //容量

  /*各遺伝子座を調べて容量と評価値を計算*/
  for(pos=0;pos<N;++pos){
    capa+=parcel[pos][0]*g[pos];
    value+=parcel[pos][1]*g[pos];
  }
  /*致死遺伝子の処理(容量制限を超えてしまった)*/
  if(capa>=LIMIT)value=0;
  return value;
}

/*乱数*/
int rndn(int l){
  int rndno; //生成した乱数
  while((rndno=((double)rand()/RAND_MAX)*l)==l);
  return rndno;
}
```

使ったデータセット
--
品物の容量と価格のリストです
ga <data.txt で読み込んで使います。

```text:data.txt
65 27
39 82
9 85
72 71
87 91
91 28
34 92
58 79
3 27
12 82
92 15
39 49
83 54
76 43
6 26
77 2
68 6
24 60
60 47
6 40
91 58
44 68
50 33
91 92
57 62
97 49
96 68
39 77
89 6
24 97
```