最長片道きっぷの経路を求める
[付録3-2] 2005年4月版(普及編)

目次


ご家庭で LOP

 これまで延々と「整数計画法による最長経路の導出」について説明してきましたが、 いざ、これを追試してみようと思ったとき、 最大の問題は「ソルバーが手元にない」ということでした。 私たちが最初に最長経路を求めたとき、 利用したのは「CPLEX」という商用のソルバーでしたが、 これは性能が高いかわりにお値段も高く、 「最長経路を求めたい」という理由だけで一個人が買うのは難しいでしょう。
 本編では「LP solve」というフリーの整数(線形)計画ソルバーがあることを紹介しましたが、 性能面ではやはり、商用のソルバーに太刀打ちできません。 具体的には、解ける問題の大きさ(規模)の上限が小さいのです。 その当時、JR全線を対象とした最長経路を LP solve で求めることは、さらに問題を分割しない限りは不可能のようでした。

 ところが、本編の発表から約5年、事情は大きく変わりました。LP solve は着実にバージョンを上げ、そのたびに性能も改善されてきましたが、これとは別に GNU の整数(線形)計画処理系「GLPK」(GNU Linear Programming Kit)が登場しました。
 あくまで「Kit」なので、 本来は「自作プログラムで線形計画問題を解きたい人のためのライブラリ」 (部品のようなもの)であって、 ソルバーという「完成品」の提供が主目的ではないようですが、 いちおう簡易なソルバー(glpsol.exe)が付属してきます。 簡易といっても最長経路の算出には十分な機能を持っていますし、 キモである「整数計画問題を解く部分」の性能はなかなかのもので、 少なくとも最長経路の算出に関しては LP solve をしのぐ高性能ぶりを発揮します。
 この GLPK を利用すると、「本編」と同様の方針で(本州を分割することなく)、 フリーのソフトウェアだけでJR全線を対象とした最長経路を求めることが可能です。 とはいっても商用のソルバーにはまだ遠く及ばず、 難しい問題には歯が立たないこともあるのですが、 それでも自分の手元でそこそこの規模の問題が解けるというのは非常に便利で、 最近は私も、主に GLPK を用いて、自分のPCで最長経路を求めています。「ご家庭で LOP」という夢の(?)時代がやってきたのです。

 念のため書いておきますが、 ふだん計算に利用しているのはふつうの小型ノートPC(Pentium M 900MHz・メモリ768MB)で、特に高性能というわけではありません。 この程度の性能のPCでも、JR線を対象とするなら、 多くの場合には現実的な時間内に解答を得られます。


さっそく導入してみる

 というわけで、興味のある方はぜひ、「ホームシアター」ならぬ「ホーム LOP 環境」を整えてみましょう。
 必要になるものは以下のとおりです。ほんとうに必須なのは GLPK だけですが、「解き直し」などの作業を効率的に行ったり、GLPK の計算結果を路線図上に描画したりするための道具もあったほうが便利です。

  1. GLPK
  2. LOP toolkit(SWA 謹製)
  3. Perl
  4. Ghostscript 関係

GLPK

 まずは GLPK ですが、Windows 版のバイナリは gnuwin32.sourceforge.net の GLPK のページで入手できます。よく分からない人は「Complete package, except sources」をダウンロードしてきてセットアップしてください。 動かすだけなら「Binaries」の Zip ファイルを展開するだけでOKです。
 その他のOSの人は自力でがんばってください。ただし、LOP toolkit は Windows 環境を前提としたものなので、以下でさらにがんばる必要があります。
 2005年3月現在、最新のバージョンは4.8です。 バージョンが上がるごとに性能も上がっているようですので、 できるだけ新しいバージョンを利用するといいでしょう。

LOP toolkit(SWA 謹製)

 次に、この LOP のページから LOP toolkit (v0.01、約30KB) をダウンロードしてください。LHA で圧縮してありますので、適当なフォルダに解凍してください。 (これができない人は、以下で必ず挫折すると思いますので、 この時点であきらめてください。;-)
 LOP toolkit の用法については後述します。

 この「LOP toolkit」は「本来は自分専用だけど、 せっかくだから公開するか」という程度のスタンスで公開しているものですので、 無保証ですし、サポートも全くありません。 「この制約式は変ではないか」「スクリプトの○行目に不具合がある」 といった具体的な指摘はありがたくお受けしますが、 単に「動かないんですけど」といった問い合わせは勘弁してください。 動かなければバッチファイルを修正して、それでも動かなければ Perl スクリプトの中身を見てみるか、といったレベルのユーザを想定しています。

Perl

 LOP toolkit では Perl を利用するので、Perl もインストールしてください。 バージョンは4以降なら何でもかまいません。jperl でも特に問題はないはずですが、MARS のバイナリデータを CSV に変換するスクリプトはちょっと怪しげです。
 よく分からない人は ActivePerl の最新版を入れておけば無難でしょう。

Ghostscript 関係

 LOP toolkit では、計算結果を Postscript 形式の地図に変換するので、Postscript(日本語入り)を解釈し、 画面表示する環境が必要です。現実的には Ghostscript(+GSView)ということになるでしょう。
 Ghostscript で日本語を扱う方法はいくつかあり、TeX 関係のサイトにはたいてい詳細な説明があると思いますので、 ここでは詳しくは述べません。
 ただ、「どうすればいいか見当も付かない」という人のために、「インストール(Windows) - TeX Wiki」という Web ページ(三重大の奥村先生による)を紹介しておきます。このページは TeX のインストール方法を紹介したものですが、下のほうに「Ghostscript のインストール」「GSView のインストール」という節がありますので、 それに従って Ghostscript と GSView をインストールしてください。 うまくインストールできれば、Postscript ファイル(拡張子「.ps」)に GSView が関連づけられ、計算結果を即座に地図に変換して見ることができます。


LOP toolkit の使いかた

 以上のインストール作業が終わったところで、いよいよ「LOP toolkit」の使いかたを紹介します。
 かんたんに言えば、LOP toolkit は、以下の2つの要素から成っています。

  1. 路線データを入力すると制約式を出力する
  2. 1. で出力した制約式(=整数計画問題)を GLPK 付属のソルバーに解かせ、 「孤立ループ」が出現した場合、それを除去するための制約式を自動的に追加する

システム構成図

 たいしたシステムではありませんが、スクリプトやファイルの数は多いので、図A3-4 にそれぞれの関係をまとめておきます。 赤い枠のものがプログラム、黒い枠のものがデータで、 矢印は入出力の関係を示します。 なお、念のため書いておきますが、ソルバー(glpsol.exe)は GLPK に付属するものであって、筆者が制作したものではありません。

図A3-4

図A3-4 システム構成図

データの準備

 図A3-4 をよく見ると、「路線データ」と「色なし地図」を用意すれば、 あとは自動的に事が進むということが分かります。
 路線データ、地図データについての解説は LOP toolkit に含まれるドキュメントで行いますが、地図データ、 路線データとも、私がふだん使っているデータをサンプルとして LOP toolkit に含めてありますので、本州のJR線を対象にして計算を行う場合、 一から作らなくてもだいじょうぶです。
 また、路線データについては、JR運賃計算システム「MARS」のデータを変換して作ることもできます。MARS のデータは全国を対象としたものなので、 必要な地域のデータだけを抜き出して使うことになります。 ただ、この方法で路線データを作ると、 サンプルで添付した地図データとの対応がとれていないので、 手作業で対応をとることになります。

具体的な作業例

 具体的な作業はコマンドプロンプト(DOS 窓)で行います。路線データ edges.csv と(色なし)地図データ map.txt を準備すれば、あとは以下のようにコマンド3個で終わりです。 (以下ではタイプ Lee の最適解を求めようとしています。 「Lee」となっているところを「Pe」「B」「B8」に変えると、 それぞれのタイプの最適解が得られます。)

C:\LOP>perl edge2cs.pl edges.csv > css.csv
C:\LOP>perl makedata.pl Lee > east_Lee.mod
C:\LOP>auto east_Lee Lee map.txt

 「auto」の実体はバッチファイル「auto.bat」で、 このバッチファイルからソルバー(glpsol.exe)やループチェック用スクリプト (loopchk.pl)、計算結果の「色つき地図化」に必要な諸スクリプトを呼び出し、 計算経過を地図に描き込みながら、ループ除去が済むまで計算を繰り返します。 1回目の計算結果は「res001.txt」、 その結果にループ除去制約式を加えた2回目の計算結果は「res002.txt」…というテキストファイルに残され、 孤立ループが完全になくなると「res.txt」に計算結果が残ります。 これらのテキストファイルを保存しておくことで、 計算の過程をトレースすることができます。 また、計算経過は「res.ps」という Postscript 文書に上書きしていくので、GSView で res.ps を開き、自動更新をオンにしておくと、解き直しの過程を視覚的に追うことができます。


GLPK の出力の見かた

 最後に、GLPK の出力の見かたについて説明しておきます。 私も出力のすべてを理解しているわけではないので、 必要と思われるところだけを説明します。

線形緩和で上界をつかむ

 まず、ソルバー glpsol.exe を起動すると、 定義した変数などの情報がすさまじい勢いで流れていき、それから計算が始まります。 開始直後は以下のような表示になっているはずです。

Model has been successfully generated
lpx_simplex: original LP has 602 rows, 322 columns, 2555 non-zeros
lpx_simplex: presolved LP has 517 rows, 165 columns, 1942 non-zeros
lpx_adv_basis: size of triangular part = 517
*     1:   objval =  7.760000000e+002   infeas =  0.000000000e+000 (0)
*   171:   objval =  3.237266667e+004   infeas =  0.000000000e+000 (0)
OPTIMAL SOLUTION FOUND

 いきなり「OPTIMAL SOLUTION FOUND」と出てくるので「もう計算が終わったのか!」と思いますが、 これは、ソルバーに与えた数々の制約のうち、 「変数の値が整数である」という制約を全部無視して解いた場合の最適解が見つかったということです。 このように、整数計画問題から「この変数は整数でなくてはならない」という制約 (=整数制約)を全部取り除くことを「線形緩和」と呼びます。
 詳しくは線形計画法の教科書などを見てほしいのですが、 ある整数計画問題と、それを線形緩和した問題とを比較すると、 線形緩和した問題のほうが圧倒的に易しいのです。しかし LOP を解くうえでは、枝に対応する変数に「通る=1、 通らない=0」という意味を与えているので、 計算した結果、この変数の値が0.5とか0.221とかになってしまうと困ります。 したがって、最長経路を得るにはもとの整数計画問題をがんばって解くしかありません。

 ただ、この「線形緩和した場合の最適解」にも意味がないわけではありません。 線形緩和した場合の最適値(ここでは32372.66667)は、 本来解くべき整数計画問題の最適値の上界を示しているからです。
 今、解いている問題は「経路全体の距離を最大化したい」というものですから、 新たな制約を加えれば加えるほど、 最適値は小さくなっていきます(この部分が理解できない人は、「[2-4] 整数計画法で(戦略編)」の 「フィルターが云々…」という部分をよく読み直してください)。 今すんなり解けてしまった「線形緩和問題」と本来解くべき「整数計画問題」とでは、 整数制約が加わる分だけ、整数計画問題のほうが制約が多いことになります。 ですから、最大化問題においては、 常に「(整数計画問題の最適値)≦(それを線形緩和した問題の最適値)」ということがいえるわけです。
 このことを LOP に即して言いかえると、「これから先、 ループ除去だ何だと苦労して問題を解き進めても、 この制約式における最長経路の長さがここで出てきた値(上記の例では32372.66667)を超えることはない」ということになります。
 通常、ある範囲の路線データを用意したら、Lee、Pe、B、B8 と4つのタイプについて計算を行い、 4タイプのうち最も距離の長いものが晴れて最長経路となるわけですが、 どれか1タイプを最後まで解いて「暫定最長距離」を求められたら、 2タイプ目以降は、その「暫定最長距離」を超えられるかどうか、 ということをまずチェックするといいでしょう。たとえば最初にタイプ Lee を解いたら最長距離(この時点では「暫定最長距離」)が33000だったとします。 次にタイプ Pe を試して、線形緩和しても32372.66667にしかならないということが分かれば、 その時点でタイプ Pe の計算は(真の最長経路になる見込みがないので)やめていいわけです。

計算時間の見込みは?

 線形緩和した問題を解き終わると、いよいよ本来の整数計画問題を解き始めます。 すると、以下のように途中経過が出力されます。 最適値のとり得る範囲が更新されれば随時出力されますが、 しばらく更新されない場合でも(freeze したのでは?と心配しないために)およそ5秒おきに出力されます。

Integer optimization begins...
Objective function is integral
+   171: mip =     not found yet <=              +inf        (1; 0)
+   301: mip =  2.809800000e+004 <=  3.237200000e+004  15.2% (18; 0)
+   499: mip =  2.867900000e+004 <=  3.237200000e+004  12.9% (44; 2)
+   642: mip =  2.887500000e+004 <=  3.172400000e+004   9.9% (53; 12)
+   738: mip =  2.892500000e+004 <=  3.172400000e+004   9.7% (59; 20)
+  1591: mip =  2.893200000e+004 <=  3.076000000e+004   6.3% (126; 52)
+  2438: mip =  2.902500000e+004 <=  3.076000000e+004   6.0% (174; 98)
+  3698: mip =  2.903900000e+004 <=  3.076000000e+004   5.9% (226; 217)
+  8628: mip =  2.904200000e+004 <=  3.076000000e+004   5.9% (450; 562)
+ 13215: mip =  2.904200000e+004 <=  3.076000000e+004   5.9% (592; 1106)
+ 18257: mip =  2.904200000e+004 <=  3.039900000e+004   4.7% (647; 1793)
+ 22837: mip =  2.904200000e+004 <=  3.039900000e+004   4.7% (551; 2741)
+ 27119: mip =  2.904200000e+004 <=  3.022400000e+004   4.1% (298; 4049)

 たとえば2行目に注目してみます。

 まず、左端の「+ 301」という数字は GLPK のマニュアルによると「simplex iteration number」で、累計の計算量だと考えればいいでしょう。 増分は解く問題によりますが、おおよそ時間に比例して増えていきます。

 次の「mip = 2.809800000e+004 <= 3.237200000e+004」という部分は、見てのとおり、 「現段階では最適値が28098以上32372以下ということが分かっている」という意味です。
 より正確に言うと、「28098」は「現時点での暫定最長経路の距離」、 「32372」は上界で、「総延長2809.8kmとなる経路は実際に作れる」 「逆に、3237.2kmを超える経路は絶対に作れない」ということが判明しています。 その間の長さの経路は「実際に作れるかどうか分からない状態」です。 右側の数字(32372)が、 先に出てきた線形緩和問題の最適値(32372.66667)を上回っていないことに注意してください。
 このように最適値となり得る範囲を随時確認できるので、この欄の値によっては、 計算の途中で「このタイプには望みがない」と判断することができます。
 その右の「15.2%」は、最適値のとり得る幅(28098〜32372)、 言いかえれば「実際にその長さの経路を作れるかどうか分からない範囲」 を割合で表示したものです(計算方法は「32372÷28098−1」)。 「分からない範囲」は計算を進めていくと小さくなり、最終的には0になります。

 最後に「(18; 0)」ですが、マニュアルによると「18」が「number of subproblems in the active list」、「0」が「number of subproblems which have been solved (considered)」です。それぞれ「これから解かなければいけない子問題の数」 「これまでに解き終わった子問題の数」といった意味です。 (ソルバーが整数計画問題を解く際には、 もとの整数計画問題をいくつかの子問題に分割し、 必要ならその子問題をさらに分割し…ということを行っています)。
 正直言って私も正確に理解しているわけではありませんが、 整数計画法を道具として使うぶんには、前者が「未決裁の箱に残っている書類の数」、 後者が「決裁済の箱に入れた書類の数」といった感覚的な理解でいいと思います。
 前者は「未決裁の数」なので、 最初のうちは(子問題の分割が続いて)増えていきますが、 そのうち減少に転じて、最終的には0になります。 したがって前者の数は、計算が峠を越えたか、 それともまだ折り返し地点に到達していないか、 というのを判断する目安になります。 一方、後者は「決裁済の数」なので単純に増えていくだけです。

計算結果

 そしてお待ちかね、計算結果です。

+ 31431: mip =  2.904200000e+004 <=  2.991500000e+004   3.0% (155; 5204)
+ 33295: mip =  2.904200000e+004 <=     tree is empty   0.0% (0; 6255)
INTEGER OPTIMAL SOLUTION FOUND
Time used:   35.3 secs
Memory used: 2.6M (2744937 bytes)
lpx_print_mip: writing MIP problem solution to `res.txt'...

 この場合、得られた最適値は29042ということになります。 計算時間、使ったメモリ量などは見てのとおりです。 結果は「res.txt」に書き込まれたということなので、この中身を見てみます。

Problem:    test
Rows:       602
Columns:    322 (210 integer, 140 binary)
Non-zeros:  2555
Status:     INTEGER OPTIMAL
Objective:  distance = 29042 (MAXimum) 32372.66667 (LP)

 6行目の「Objective」(目的関数の意)を見ると、 整数計画問題の最適値(MAXimum)と、 線形緩和した場合の最適値(LP;Linear Programming =線形計画の意)が並んで書かれています。もちろん最長経路の長さは前者です。
 以下、 最適値をとるときに各変数の値がどうなっているかというのが延々書いてありますが、 重要なのは枝の使用状況なので、「e001」などという変数名を探しましょう。

   No. Column name       Activity     Lower bound   Upper bound
------ ------------    ------------- ------------- -------------
     1 isolation                   2             2             =
     2 maxv                        2             2             =
     3 maxv3                       2             2             =
     4 forksum                     0             0             =
     5 e112         *              0             0             1
     6 e243         *              1             0             1
     7 e251         *              1             0             1
(中略)
   143 e402         *              1             0             1
   144 e999         *              0             0             1

 ここでは「Activity」の欄に注目します。 たとえば112番の枝の使用状況を示す「e112」は「0」、 つまり112番の枝は通りません。以下、243番は通る、251番は通る、 …、402番は通る、999番は通らない、ということが分かります。 「色つき地図」を作るスクリプトはこの部分から計算結果を得ています。



(c) SWA / KASAI TakayaSWA へメールを送る
最終更新: 2005年 4月 9日
簡易包装にご協力ください。