Lin-Kernighan法で巡回セールスマン問題を解く

[さざなみ]

以下のサイトを見て巡回セールスマン問題をLin-Kernighan法で解くプログラムを実装ようとしています。

巡回セールスマン問題におけるLin-Kernighan法の実装
http://www.th.cs.meiji.ac.jp/researches/2007/takigawa/

ファイルから座標を読み込んだあと、nn法で巡回路を作り
Lin-Kernighan法で解を近似していくプログラムを作ったのですが
実行した結果はLin-Kernighan法を実行するとnn法で作った巡回路よりも長い巡回路長になってしまいます。

そこで、参考にしたサイトのアルゴリズムと自分が作成したlk関数とで
違うところなどがありましたら教えてほしいです。

どうかお力を貸してください。
よろしくお願いします。

コード:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAX_N 10000   // 点の数の最大値
#define INF 100000000 // 無限大の定義
#define EPSILON 0.00000001 //ε 小さい正の値
#define SWAP(a,b){int temp; temp=(a); (a)=(b); (b)=temp; }   


struct point {
  int x;
  int y;
};



double dist(struct point p, struct point q) { // pとq の間の距離を計算 
  return sqrt((p.x-q.x)*(p.x-q.x)+(p.y-q.y)*(p.y-q.y));
}

double tour_length(struct point p[MAX_N], int n, int tour[MAX_N]) {
  int i;
  double sum=0.0;
  for(i=0;i<n;i++) sum+=dist(p[tour[i]],p[tour[(i+1)%n]]);
  return sum;// 総距離が関数の戻り値
}

void read_tsp_data(char *filename, struct point p[MAX_N],int *np) {
  FILE *fp;
  char buff[100];
  int i;

  if ((fp=fopen(filename,"rt")) == NULL) {// 指定ファイルを読み込み用に開く
    fprintf(stderr,"Error: File %s open failed.\n",filename);
    exit(0);
  }   

  while((fgets(buff,sizeof(buff),fp)!=NULL)   // DIMENSION で始まる行に出会う
	&&(strncmp("DIMENSION",buff,9)!=0)) ; // まで読み飛ばす. 
  sscanf(buff,"DIMENSION : %d",np);           // 点の数 *np を読み込む

  while((fgets(buff,sizeof(buff),fp)!=NULL)   // NODE_COORD_SECTIONで始まる
	&&(strncmp("NODE_COORD_SECTION",buff,18)!=0)) ; // 行に出会うまで, 
                                                        // 読み飛ばす. 
  for(i=0;i<*np;i++) {                       // i=0 から i=(*np)-1まで
    if(fgets(buff,sizeof(buff),fp)!=NULL) 
      sscanf(buff,"%*d %d %d",&(p[i].x),&(p[i].y)); // i番目の点の座標を読み込む
  }                                 

  fclose(fp);
}

void nn(struct point p[MAX_N], int n, int tour[MAX_N]){
  FILE *fp;
  int i,j,nearest;
  int visited[MAX_N]; // 都市iが訪問済みならば1そうでなければ0
  int min;

  for(i=0;i<n;i++) visited[i]=0; // 最初は全都市は未訪問
  tour[0]=0;         // 最初に訪問する都市は 0 
  visited[0]=1;      // 都市0は訪問済み

  for(i=0;i<n-1;i++) {
    //最後に訪問した都市tour[i]から最短距離にある未訪問都市nearestを
    //見つける
    min = INF;          
    for(j=1;j<n;j++) { 
      if(!visited[j] && dist(p[tour[i]],p[j])<min){ 
	      nearest=j;      //都市tour[i]から暫定的に最短距離にある未訪問都市をnearestとする
	      min = dist(p[tour[i]],p[j]); // その暫定最短距離をminとする
      }
    }
    tour[i+1]=nearest; // i+1 番目に訪問する都市を nearest にして, 
    visited[nearest]=1;// nearest を訪問済みとする. 
  }

}

void lk(struct point p[MAX_N], int n, int tour[MAX_N]){

  int i, j, min = INF, base1, base2, retsu, step = 0, g, h, sucess;
  struct point base, nbase, probe, nprobe; 
  double t_dist, omega;

  //初めの順回路はnnで作ったやつを使う。
  
  do{
  //最も短い経路長を持つ3点を総当りで求める
    for(i=0; i<n-3; i++){
      t_dist = dist(p[tour[i]],p[tour[i+1]])+dist(p[tour[i+1]],p[tour[i+2]]);
      if(t_dist < min){
        min = t_dist;
        base1 = i;
      }
    } 
  
    //その3点をツアー配列の先頭に移動させた配列をつくる
    retsu = 0;
    for(base1=i; base1<n; base1++){     
      base2 = base1;
      for(; base2>retsu; base2=base2-1){
        SWAP(tour[base2-1],tour[base2]);
      }
      retsu++;
    }
     
    //1回目の動作はループ外 
    base = p[tour[2]];
    nbase = p[tour[3]];
   
    for(i=4; i<n-1; i++){
      probe = p[tour[i]];
      nprobe = p[tour[i+1]];
      if(dist(nbase, nprobe) < dist(base, nbase)){
        omega = dist(base, nbase) - dist(nbase, nprobe);
        g = 3; h = i;
        while(g < h){
          SWAP(tour[g], tour[h]);
          g++; h--;
        }
        break;
      }
    }
  
   //2回目以降の動作
  
   do{
      sucess = 0;
      base = p[tour[2]];
      nbase = p[tour[3]];
      
      for(i=4; i<n-1; i++){
        probe = p[tour[i]];
        nprobe = p[tour[i+1]];
        if(dist(nbase, nprobe) < dist(base, nbase) + omega){
          omega = dist(base, nbase) + dist(probe, nprobe) - dist(base, probe) - dist(nbase, nprobe) + omega;
          g = 3; h = i;
          while(g < h){
            SWAP(tour[g], tour[h]);
            g++; h--;
          }
          break;
          sucess = 1;
        }
      }
    }while(sucess == 1);
    step++;
      
  }while(step != 1000);  
}

void write_tour_data(char *filename, int n, int tour[MAX_N]){
  FILE *fp; 
  int i;
 
 // 構築した巡回路をfilenameという名前のファイルに書き出すためにopen
  if((fp=fopen(filename,"wt"))==NULL){ 
    fprintf(stderr,"Error: File %s open failed.\n",filename);
    exit(EXIT_FAILURE);
  }
  fprintf(fp,"%d\n",n);
  for(i=0;i<n; i++){
   fprintf(fp,"%d ",tour[i]);
  }
  fprintf(fp,"\n");
  fclose(fp);
}

main(int argc, char *argv[]) {
  int  n;                   // 点の数 
  struct point  p[MAX_N];   // 各点の座標を表す配列 
  int tour[MAX_N];   // 巡回路を表現する配列

  if(argc != 2) {
    fprintf(stderr,"Usage: %s <tsp_filename>\n",argv[0]);
    exit(EXIT_FAILURE);
  }

  // 点の数と各点の座標を1番目のコマンドライン引数で指定されたファイルから読み込む
  read_tsp_data(argv[1],p, &n);
  // 最近近傍法による巡回路構築
  nn(p,n,tour);
  // 巡回路をテキストファイルとして出力
  write_tour_data("tour1.dat",n,tour);
  // 巡回路長を画面に出力
  printf("%lf\n",tour_length(p,n,tour));
  // Lin-Kernighan法
  lk(p,n,tour);
  // 巡回路長を画面に出力
  printf("%lf\n",tour_length(p,n,tour));
  // 巡回路をテキストファイルとして出力
  write_tour_data("tour2.dat",n,tour);
  
  exit(EXIT_SUCCESS);
}

[box]

提示のコードの142～143行目

コード:

          break;
          sucess = 1;

ものすごく違和感があります。ふつうは逆じゃないかと思うのですが、どうなんでしょうか。

[さざなみ]

boxさん

ご指摘ありがとうございます。
ですが、その部分を入れ替えたら今度はプログラムを実行しても終わらなくなってしまいました。

もう少し、自分で考えてみます。

ありがとうございました。

[NakaTaka]

すでに気づいているかもしれませんが、
>ですが、その部分を入れ替えたら今度はプログラムを実行しても終わらなくなってしまいました。
これに関して、breakする前にsucessの値を１にするので、ifの条件を満たす限り

コード:

} while (sucess == 1);

を抜けられません。
回答でなくて、本当に申し訳ございません。

[f8]

そこで、参考にしたサイトのアルゴリズムと自分が作成したlk関数とで
違うところなどがありましたら教えてほしいです。

アルゴリズムのページを見ましたがあの情報だけで実装するのは私には無理っぽいです。(^_^;
したがって、正しいのか間違っているのか判断はできないことをご理解ください。

コード中の私のコメントには f8 を付けています。

コード:

void lk(struct point p[MAX_N], int n, int tour[MAX_N]){

  int i, j, min = INF, base1, base2, retsu, step = 0, g, h, sucess;
  struct point base, nbase, probe, nprobe; 
  double t_dist, omega;

  //初めの順回路はnnで作ったやつを使う。
  
  do{  // f8: ループ開始はここでいいの？
  //最も短い経路長を持つ3点を総当りで求める
    for(i=0; i<n-3; i++){
      t_dist = dist(p[tour[i]],p[tour[i+1]])+dist(p[tour[i+1]],p[tour[i+2]]);
      if(t_dist < min){
        min = t_dist;
        base1 = i;
      }
    } 
  
    //その3点をツアー配列の先頭に移動させた配列をつくる
    retsu = 0;
    for(base1=i; base1<n; base1++){  // f8: 15行目でbase1を求めているので再初期化は不要では？
      base2 = base1;
      for(; base2>retsu; base2=base2-1){
        SWAP(tour[base2-1],tour[base2]);
      }
      retsu++;
    }
     
    //1回目の動作はループ外 
    base = p[tour[2]];
    nbase = p[tour[3]];
   
    for(i=4; i<n-1; i++){
      probe = p[tour[i]];
      nprobe = p[tour[i+1]];
      if(dist(nbase, nprobe) < dist(base, nbase)){
        omega = dist(base, nbase) - dist(nbase, nprobe);
        g = 3; h = i;
        while(g < h){
          SWAP(tour[g], tour[h]);
          g++; h--;
        }
        break;
      }
    }
  
   //2回目以降の動作
  
   do{
      sucess = 0;
      base = p[tour[2]];
      nbase = p[tour[3]];
      
      for(i=4; i<n-1; i++){
        probe = p[tour[i]];
        nprobe = p[tour[i+1]];
        if(dist(nbase, nprobe) < dist(base, nbase) + omega){
          omega = dist(base, nbase) + dist(probe, nprobe) - dist(base, probe) - dist(nbase, nprobe) + omega;
          g = 3; h = i;
          while(g < h){
            SWAP(tour[g], tour[h]);
            g++; h--;
          }
          break;       // f8: sucess = 1;
          sucess = 1;  // f8: break;
        }
      }
    }while(sucess == 1);
    step++;
      
  }while(step != 1000);  
}

[さざなみ]

f8さん，NakaTakaさん

ご指摘ありがとうございます．

お二方にアドバイスいただいた事もふまえて自分でも考えてみたのですが
なかなか上手くいかないため、現在は焼きなまし法を実装しようと考えています。

今回はありがとうございました。

焼きなまし法でもうまくいかないところがあり、また別のトピックでお世話になるかもしれません。
そのときはよろしくお願いします。

[f8]

なかなか上手くいかないため、現在は焼きなまし法を実装しようと考えています。

あそこまで実装してやめるのはもったいないですねー。(^_^;

クローズになっていますが、私が試したコードを貼っておきます。
追加・修正したところには f8 付きのコメントを残しています。

別トピックで提示されたテストデータをいくつか試しました。
全部は試してないのであれですが、この実装では "Drilling problem" のデータは解けないようです。(^_^;

サンプル数が大きめの fnl4461.tsp は解けました(正しいかは不明)。

コード:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAX_N 10000   // 点の数の最大値
#define INF 100000000 // 無限大の定義
#define EPSILON 0.00000001 //ε 小さい正の値
#define SWAP(a,b){int temp; temp=(a); (a)=(b); (b)=temp; }


struct point {
    int x;
    int y;
};



double dist(struct point p, struct point q) {   // pとq の間の距離を計算
    return sqrt((p.x-q.x)*(p.x-q.x)+(p.y-q.y)*(p.y-q.y));
}

double tour_length(struct point p[MAX_N], int n, int tour[MAX_N]) {
    int i;
    double sum=0.0;
    for(i=0;i<n;i++) sum+=dist(p[tour[i]],p[tour[(i+1)%n]]);
    return sum;// 総距離が関数の戻り値
}

void read_tsp_data(char *filename, struct point p[MAX_N],int *np) {
    FILE *fp;
    char buff[100];
    int i;

    if ((fp=fopen(filename,"rt")) == NULL) {    // 指定ファイルを読み込み用に開く
        fprintf(stderr,"Error: File %s open failed.\n",filename);
        exit(0);
    }

    while((fgets(buff,sizeof(buff),fp)!=NULL)   // DIMENSION で始まる行に出会う
        &&(strncmp("DIMENSION",buff,9)!=0)) ;   // まで読み飛ばす. 
    sscanf(buff,"DIMENSION : %d",np);           // 点の数 *np を読み込む

    while((fgets(buff,sizeof(buff),fp)!=NULL)           // NODE_COORD_SECTIONで始まる
        &&(strncmp("NODE_COORD_SECTION",buff,18)!=0)) ; // 行に出会うまで, 
                                                        // 読み飛ばす. 
    for(i=0;i<*np;i++) {                       // i=0 から i=(*np)-1まで
        if(fgets(buff,sizeof(buff),fp)!=NULL) 
            sscanf(buff,"%*d %d %d",&(p[i].x),&(p[i].y));   // i番目の点の座標を読み込む
    }                                 

    fclose(fp);
}

// f8: tourの中身を標準出力へ
void print_tour(const int* tour, int n)
{
    int i;
    printf("--tour:");
    for(i = 0; i < n; i++) {
        printf(" %d", *(tour+i));
    }
    printf("\n");
}

// f8: pの中身を標準出力へ
void print_point(const struct point* p, int n)
{
    int i;
    printf("--point\n");
    for(i = 0; i < n; i++) {
        printf("%d: (%d,%d)\n", i, (p+i)->x, (p+i)->y);
    }
}

void nn(struct point p[MAX_N], int n, int tour[MAX_N]){
    FILE *fp;
    int i,j,nearest;
    int visited[MAX_N]; // 都市iが訪問済みならば1そうでなければ0
    int min;

    for(i=0;i<n;i++) visited[i]=0; // 最初は全都市は未訪問
    tour[0]=0;         // 最初に訪問する都市は 0 
    visited[0]=1;      // 都市0は訪問済み

    for(i=0;i<n-1;i++) {
        //最後に訪問した都市tour[i]から最短距離にある未訪問都市nearestを
        //見つける
        min = INF;
        for(j=1;j<n;j++) { 
            if(!visited[j] && dist(p[tour[i]],p[j])<min){ 
                nearest=j;      //都市tour[i]から暫定的に最短距離にある未訪問都市をnearestとする
                min = dist(p[tour[i]],p[j]); // その暫定最短距離をminとする
            }
        }
        tour[i+1]=nearest; // i+1 番目に訪問する都市を nearest にして,
        visited[nearest]=1;// nearest を訪問済みとする.
    }

}

void lk(struct point p[MAX_N], int n, int tour[MAX_N])
{
    int i, j, base1, base2, retsu, step = 0, g, h, sucess;
    struct point base, nbase, probe, nprobe;
    double min = INF, t_dist, omega;    // f8: minをdouble型へ変更

    //初めの順回路はnnで作ったやつを使う。

    //do{   // f8: コメントアウト
        //最も短い経路長を持つ3点を総当りで求める
        for(i=0; i<n-3; i++){
            t_dist = dist(p[tour[i]],p[tour[i+1]])+dist(p[tour[i+1]],p[tour[i+2]]);
            printf("dist((%d,%d),(%d,%d))+dist((%d,%d),(%d,%d))=%lf\n", p[tour[i]], p[tour[i+1]], p[tour[i+1]], p[tour[i+2]], t_dist);  // f8:
            if(t_dist < min){
                min = t_dist;
                base1 = i;
            }
        }

        printf("min=%lf, base1=%d\n", min, base1);  // f8:
        print_tour(tour, n);                        // f8:

        //その3点をツアー配列の先頭に移動させた配列をつくる
        retsu = 0;
        for(/*base1=i*/; base1<n; base1++){     // f8: base1の再初期化しない
            base2 = base1;
            for(; base2>retsu; base2=base2-1){
                SWAP(tour[base2-1],tour[base2]);
            }
            retsu++;
        }

    do{     // f8: ループ開始はここでいいんじゃね？
        print_tour(tour, n);    // f8:

        //1回目の動作はループ外
        base = p[tour[2]];
        nbase = p[tour[3]];

        for(i=4; i<n-1; i++){
            probe = p[tour[i]];
            nprobe = p[tour[i+1]];
            if(dist(nbase, nprobe) < dist(base, nbase)){
                omega = dist(base, nbase) - dist(nbase, nprobe);
                g = 3; h = i;
                while(g < h){
                    SWAP(tour[g], tour[h]);
                    g++; h--;
                }
                break;
            }
        }

        //2回目以降の動作
        do{
            sucess = 0;
            base = p[tour[2]];
            nbase = p[tour[3]];

            for(i=4; i<n-1; i++){
                probe = p[tour[i]];
                nprobe = p[tour[i+1]];
                if(dist(nbase, nprobe) < dist(base, nbase) + omega){
                    omega = dist(base, nbase) + dist(probe, nprobe) - dist(base, probe) - dist(nbase, nprobe) + omega;
                    g = 3; h = i;
                    while(g < h){
                        SWAP(tour[g], tour[h]);
                        g++; h--;
                    }
                    sucess = 1;     // f8:
                    break;          // f8:
                }
            }
            printf("loop2: sucess=%d\n", sucess);   // f8:
        }while(sucess == 1);
        step++;

        printf("loop1: step=%d\n", step);   // f8:
    }while(step != 1000);
}
 
void write_tour_data(char *filename, int n, int tour[MAX_N]){
    FILE *fp; 
    int i;

    // 構築した巡回路をfilenameという名前のファイルに書き出すためにopen
    if((fp=fopen(filename,"wt"))==NULL){ 
        fprintf(stderr,"Error: File %s open failed.\n",filename);
        exit(EXIT_FAILURE);
    }
    fprintf(fp,"%d\n",n);
    for(i=0;i<n; i++){
        fprintf(fp,"%d ",tour[i]);
    }
    fprintf(fp,"\n");
    fclose(fp);
}
 
main(int argc, char *argv[]) {
    int  n;                   // 点の数
    struct point  p[MAX_N];   // 各点の座標を表す配列
    int tour[MAX_N];   // 巡回路を表現する配列

    if(argc != 2) {
        fprintf(stderr,"Usage: %s <tsp_filename>\n",argv[0]);
        exit(EXIT_FAILURE);
    }

    // 点の数と各点の座標を1番目のコマンドライン引数で指定されたファイルから読み込む
    read_tsp_data(argv[1],p, &n);
    print_point(p, n);      // f8:

    // 最近近傍法による巡回路構築
    nn(p,n,tour);
    // 巡回路をテキストファイルとして出力
    write_tour_data("tour1.dat",n,tour);
    // 巡回路長を画面に出力
    printf("%lf\n",tour_length(p,n,tour));
    // Lin-Kernighan法
    lk(p,n,tour);
    // 巡回路長を画面に出力
    printf("%lf\n",tour_length(p,n,tour));
    // 巡回路をテキストファイルとして出力
    write_tour_data("tour2.dat",n,tour);

    exit(EXIT_SUCCESS);
}

[さざなみ]

f8さん

わざわざコードの添付までありがとうございます。

nnの関数を実行したあとの巡回路長の長さと，lkの関数を実行したあとの巡回路長の長さが変わらないことが多いので
やはり参考にしたサイトだけでは情報が足りなかったかなと思います。

自分でも他に改善できる個所はないか、探してみることにします。

本当にありがとうございました。