SCIP Optimization Suiteを使ってMINLPを解く

SCIPは,http://scip.zib.de/から入手可能。非線形計画問題を解くために,CppADと,lpoptを使う.

ZIBの提供するSCIP Optimization Suite(以下,SCIP)は,整数変数を含む数理計画問題を解くことができます。 非商用には無料で使用可能です。(詳しくはライセンスを確認ください。) このSCIPを使って,混合整数非線形計画問題を解いてみます。
SCIPにはZIMPLというモデリング言語が含まれています。 ZIMPLを使うと,数式に近い形で数理計画問題を記述してソルバに与えることができます。 ZIMPLは非線形計画問題を記述することはできますが, SCIPにはそれらを解く機能が含まれていません.非線形計画問題を解くには, COIN-ORによって提供されているソフトウェアIPOPTをインストールし,SCIPのコンパイル時にIPOPTを使うことを 示す必要があります。そうすることで,SCIPによって非線形計画問題を解くことができるようになります。
IPOPTとSCIPのインストールはそれなりの手間がかかるかもしれませんが,ここでは,それらのインストールについては 述べません.コマンドラインから

	$scip
によってSCIPが実行可能な状態になっている,ということを仮定して話を進めます。

準備体操として,次の二次関数の最小化問題を,ZIMPLで書いて解きます。

min f(v):=0.0036*v^2-0.1015*v+0.8848 
subject to 11 <= v <= 19
これを解くためのZIMPLファイルは,次の通りです。
param vmin:=11;
param vmax:=19;
param alpha:=0.0036;
param beta:=-0.1015;
param gamma:=0.8848;
var v;
var q;
minimize cost: alpha*q+beta*v+gamma;
subto vmin:
		vmin<=v;
subto vmax:
		v<=vmax;
subto quad:
		q>=v^2;
		
これをex1.zplとして保存して
$ scip -f ex1.zpl
とすると,v=14.0972222235149が得られます。この解は,二次関数を平方完成すれば得られることは,高校の数学の範囲内です。

ZIMPLのモデルでは,目的関数が線形でなければなりません。いまは,二次関数を最小化しようとしているので,補助変数qを導入して,

	min. alpha*v^2+beta*v+gamma
	
の代わりに
	min. alpha*q+beta*v+gamma
	s.t. q>=v^2
	
としていることに注意してください。これで,二次の要素は制約式に移し,目的関数は線形にすることができます.

次に,最短路問題に毛を生やした問題を解いてみます。

まず,基本的な最短路問題の例題として,ZIMPL User Guide version 3.3.1 Section 2 Introductionのネットワークを 使います。そのZIMPLモデルはUser Guideにあるとおり, 次のように書けます。

set V := {"a","b","s","t"};
set A := {<"s","a">,<"s","b">,<"a","b">,<"a","t">,<"b","t">};
param c[A] := <"s","a">17,<"s","b">47,<"a","b">19,<"a","t">53,<"b","t">23;
defset dminus(v):={<i,v> in A};
defset dplus(v):={<v,j>in A};
var x[A] binary;
minimize cost: sum in A: c[i,j] * x[i,j];
subto fc:
        forall  in V - {"s","t"}:
                sum in dminus(v):x[i,v] == sum in dplus(v): x[v,i];
subto uf:
        sum in dplus("s"):x[s,i]==1;
ここでは、アークの費用c[i,j]は定数として与えられています。そして, 変数は,各アークに対して定義された0-1変数x[i,j]のみです。

この最短路問題に毛を生やしてみましょう。
目的関数を定数ではなく,アーク上を移動する 際の速度vの関数とします。 速度vで移動するときの単位距離あたりの移動コストが, 上の二次関数f(v)で得られるとします。c[i,j]は,アーク上の距離と解釈することにします。 アーク(i,j)上の速度は,実数変数sp[i,j]と表しています。

set V := {"a","b","s","t"};
set A := {<"s","a">,<"s","b">,<"a","b">,<"a","t">,<"b","t">};
param c[A] := <"s","a">17,<"s","b">47,<"a","b">19,<"a","t">53,<"b","t">23;
param alpha:=0.0036;
param beta:=-0.1015;
param gamma:=0.8848;
param vmin:=11;
param vmax:=19;
defset dminus(v):={<i,v>in A};
defset dplus(v):={<v,j>in A};
var x[A] binary;
var sp[A];
var ecost[A];
minimize cost: sum in A: ecost[i,j]*c[i,j];
subto fc:
        forall  in V - {"s","t"}:
                sum in dminus(v):x[i,v] == sum in dplus(v): x[v,i];
subto uf:
        sum in dplus("s"):x[s,i]==1;
subto ecost:
        forall  in A:
                ecost[i,j]>=alpha*sp[i,j]^2+beta*sp[i,j]+gamma*x[i,j];
subto vmin:
        forall  in A:
                vmin*x[i,j]<=sp[i,j];
subto vmax:
        forall  in A:
                sp[i,j]<=vmax*x[i,j];
このZIMPLモデルをたとえばsp.zplという名前のファイルとして保存して
$scip -f sp.zpl
とすると,s->a->b->tと行くのが最適で,s->a,a->b, b->tのいずれも14.1の速度で行くのが最適である,と いう解が得られます。これは最初の最短路と同じです。当たり前です。距離が最も短い経路を,最もコストの小さい速度14.1で行くのが最適だといっているのですから。

次に,各ノードを訪問できる時刻に制限をつけます。これを時間枠制約といいます。 時間枠制約を表すために,各ノードへの訪問時刻(ノードを出発する時刻)を表す変数が必要です。 ノードiからの出発時刻をvtime[i]とします。 アーク(i,j)を使用するときは,次の関係が成り立つ必要があります。

	vtime[i]+c[i,j]/sp[i,j]<=vtime[j]
ここで,c[i,j]はアーク(i,j)の距離,sp[i,j]はアーク(i,j)上の移動速度とします。つまりこの式は,ノードiを出発する時刻に,アーク(i,j)上の移動時間を足したものより,ノードjを出発する時刻の方が大きい,ということを表しています。 この関係式は,アーク(i,j)を使用するときだけ成り立てばよいものです。 ただし,ZIMPLはc[i,j]/sp[i,j]という分数の表現が使えないので,かわりに両辺にsp[i,j]をかけた
	vtime[i]*sp[i,j]+c[i,j]<=vtime[j]*sp[i,j]
を使います。
アーク(i,j)を使用することは,x[i,j]が1であることで示されます。ZIMPLでは,変数の値を用いた分岐処理を書くことができます。この場合の処理は,具体的には次の通りに書けます。
subto tm_trans:
	forall  in A:
		vif x[i,j]>0
		then vtime[i]*sp[i,j]+c[i,j]<=vtime[j]*sp[i,j] end;
ここで,x[i,j]が1であるという条件は,x[i,j]>0と表しています。
変数の値による分岐処理が気持ち悪い,という方は,上の制約のかわりに次の書き方があります。
subto tm_trans:
	forall  in A:
		vtime[i]*sp[i,j]+c[i,j]*x[i,j]<=vtime[j]*sp[i,j];
この書き方では,c[i,j]の代わりにc[i,j]*x[i,j]を用いています。 こうすると,x[i,j]が0のときには両辺がともに0で成り立ち, x[i,j]が1のときにはvtime[i],vtime[j],sp[i,j]間の制約を意図通りに表します。
ノードiへの最早訪問時刻(earliest arrival time)をearliest[i], 最遅訪問時刻(latest arrival time)をlatest[i]と表します これらの値はパラメータとして指定するのですが,その方法は次のとおりです。
param latest[V]:=<"s"> 0,<"a"> 0.8,<"b">2,<"t"> 1.8 default 10000;
param earliest[V]:=<"s"> 0,<"a"> 0,<"b">0,<"t"> 0 default 0;
これらを用いて,時間枠制約は次のように表すことができます。
subto tm_latest:
        forall  in V:
                vtime[v]<=latest[v];
subto tm_earliest:
        forall  in V:
                earliest[v]<=vtime[v];
これらをまとめたZIMPLモデルは次のとおりです。これを例えばsp_tw.zplという名前のファイルに保存し,
$scip -f sp_tw.zpl
を実行すれば,答えが得られます。
set V:={"s","a","b","t"};
set A := {<"s","a">,<"s","b">,<"a","b">,<"a","t">,<"b","t">};
param c[A] := <"s","a">14,<"s","b">14,<"a","b">14,<"a","t">42,<"b","t">14;
param latest[V]:=<"s"> 0,<"a"> 0.8,<"b">2,<"t"> 1.8 default 10000;
param earliest[V]:=<"s"> 0,<"a"> 0,<"b">0,<"t"> 0 default 0;
param alpha:=0.0036;
param beta:=-0.1015;
param gamma:=0.8848;
param vmin:=11;
param vmax:=19;
defset dminus(v):={ in A};
defset dplus(v):={ in A};
var x[A] binary;
var sp[A];
var ecost[A];
var vtime[V] >=0 <=1000;
minimize cost: sum in A: ecost[i,j]*c[i,j];
subto fc:
        forall  in V - {"s","t"}:
                sum in dminus(v):x[i,v] == sum in dplus(v): x[v,i];
subto uf:
        sum in dplus("s"):x[s,i]==1;
subto ecost:
        forall  in A:
                ecost[i,j]>=alpha*sp[i,j]^2+beta*sp[i,j]+gamma*x[i,j];
subto vmin:
        forall  in A:
                vmin*x[i,j]<=sp[i,j];
subto vmax:
        forall  in A:
                sp[i,j]<=vmax*x[i,j];
subto tm_trans:
        forall  in A:
                vtime[i]*sp[i,j]+c[i,j]*x[i,j]<=vtime[j]*sp[i,j];
subto tm_latest:
        forall  in V:
                vtime[v]<=latest[v];
subto tm_earliest:
        forall  in V:
                earliest[v]<=vtime[v];
このモデルは,非線形のコストをもった時間枠制約付き最短路問題です。数理計画問題のクラスとしては,整数変数を含んでいて,非線形の項をもったものです。このクラスの問題を解くことのできるソルバは,ありがたいです。
(2014.1.9)