ソルバ

ZIMPL モデリング言語の使い方

集合の例(ZIMPL User's guideより)

set A := {1,2,3};
set B := {"hi","ha","ho"};
set C := {<1,2,"x">,<6,5,"y">,<787,12.6,"oh">};
ZIMPL User's guide

集合の演算も可能.

条件付きの集合も使える

set F:= { in Q with i>j and i<5};
ZIMPL User's guide

添字付き集合

set I := {1..3};
set A[I] := <1>{"a","b"},<2>{"c","e"},<3>{"f"};
set B[ in I] := {3*i};
set P[] := powerset(I); //Iの全ての部分集合
set J := indexset(P); //Pの添字集合、すなわち、{1,2,\ldots,|A|}
set S[] := subsets(I,2); //Iの部分集合で、要素数が2のものの集合
set T[] := subsets(I,1,2); //Iの部分集合で、要素数が1以上2以下のもの2のものの集合
ZIMPL User's guide

パラメータ。添字集合をつかってもつかわなくても宣言可能。デフォルト値が設定可能.

set A:= {12 .. 30};
set C:= {<1,2,"x">,<6,5,"y">,<3,7,"z">}; //集合Cを定義。それぞれの要素は、3つ組
param q:=5;
param u[A]:=<13> 17, <17> 29, <23> 14 default 99;
param amin := min A; #=12. 集合A の要素の最小値? だろう
param umin := min in A: u[a]; #=14. uの最小値。 u[a]を最小にするような, Aの要素a
param mmax := max in {1 .. 10} : i mod 5; //
param w[C] := <1,2,"x"> 1/2, <6,5,"y"> 2/3; //
param x[ in {1..8} with i mod 2==0] := 3*i; // 1..8の間の偶数、すなわち,i=2,4,6,8に、3*iを設定する.
x[2]=6, x[4]=12, x[6]=18, x[8]=24.
ZIMPL User's guide

下から2番目w[C]は要注意。Cの要素<3,7,"z">に対する値が指定されていないけれど、 w[<3,7,"z">が参照されない限り問題がない、ということ。いわゆる、疎なデータを保持することができるということだろう。しかし、メモリ管理はどう実装しているのか、という疑問がわく.

多次元の添字付パラメータは表から初期化可能. 列(カラム)の添字は1次元でなければならない。

set I:={1..10};
set J:={"a","b","c","x","y","z"};
param h[I*J] := | "a", "c","x", "z" |
|1| 12, 17, 99, 23 |
|3| 4, 3, -17, 66*5.5|
|5| 2/3, -.4, 3, abs(-4)|
|9| 1, 2, 0, 3| default -99;
ZIMPL User's guide

各行の添字は、多次元でもよい。

param g[I*I*I] := | 1, 2, 3|
|1,3|0,0,1|
|2,1|1,0,1|;
ZIMPL User's guide

表形式と、エントリのリストの併用も可

param k[I*I] := | 7,8,9|
|4|89,67,55|
|5|12,13,14|,<1,2> 17, <3,4> 99;
ZIMPL User's guide

変数

変数のタイプは、real, binary, integerのどれか. デフォルトのタイプは、real. 上限と下限も設定可能で、デフォルト値はそれぞれ0と無限大. Binaryまたはintegerの変数は、implicitと宣言できるらしい。このimplicitと宣言することによってその情報を活かして効率のよい計算をするソルバがあるようだ,といっても、Zimple2.0の時点では、 SCIP のみが対応しているらしい。 integerの変数には, startvalで初期値を設定可能キーワードpriorityは、分枝の優先度を指定する。

implicitとは?

Binary and integer variables can be declared implicit,i.e. the variable can be assumed to have an integral value in any optimal solution to the integer program, even if it is declared continous.
ZIMPL User's guide

var x1;
var x2 binary;
var y[A] real>=2<=18;
var z[ in C] integer
var z[ in C] integer
>=a*10 <= if b <= 3 then p[b] else infinity end; // z[]の下限はa*10, 上限は、b<=3のときp[b]で、それ以外のとき無限大.
var w implicit binary;
var t[k in K] integer >= 1 <= 3 * k startval 2 * k priority 50;
ZIMPL User's guide

目的

ひとつのモデルに"objective statement"はひとつ. 目的は最大化(maximize)か最小化(minimize). 目的関数は線形でなければならない.

minimize cost: 12 * x1 - 4.4 * x2 + 5
+ sum <a> in A: u[a] * y[a] + sum <a,b,c> in C with a in E and b > 3: -a/2 * z[a,b,c];
ZIMPL User's guide

Texで書いてみると、

と書けるだろう.

制約

一般的な形式は,

subto 名前: 項 意味 項

subtoはキーワードなのでこのとおり書く. 名前は制約につける名前. 意味は、<=, >= ==のうちのいずれか。

subto time: 3 * x1 + 4 * x2 <= 7;
ZIMPL User's guide

timeが名前、3*x1+4*x2が項, <=が意味, 7が項. "範囲指定"タイプの制約も書ける. この場合は、下の例のように意味(sense)が2つ出てくる. 範囲指定タイプの制約の場合は、2つの意味(sense) はおなじでなければならない。下の例では、両方とも">=".

subto space: 50 >= sum in A: 2 * u[a] * y[a] >= 5;
ZIMPL User's guide

これもtexで書いてみると,

複数の制約式をループで作ることができる.

subto weird: forall <a> in A: sum <a,b,c> in C: z[a,b,c]==55;
ZIMPL User's guide

たとえば A:={1,2,3}であれば, 上の1つの記述で、

sum <1,b,c> in C : z[1,b,c] == 55;
sum <2,b,c> in C : z[2,b,c] == 55;
sum <3,b,c> in C : z[3,b,c] == 55;

と3本の制約式を生成することができる.

subto c21: 6*(sum in A : x[i] + sum <j> in B: y[j]) >= 2;
ZIMPL User's guide

これをTeXで書くと、

forallとsumの進んだ使い方

forall <i,j> in X cross { 1 to 5 } without { <2,3> }
with i > 5 and j < 2 do sum in X cross { 1 to 3 } cross Z do p[i] * q[j] * w[j,k] >= if i == 2 then 17 else 53;
ZIMPL User's guide

集合A とBについて、A cross Bは直積を与える. forallで制約式を生成するときに、if文を使うことができる.

subto c1: forlal in I do
if (i mod 2 == 0) then 3 * x[i] >= 4
else -2 * y[i] <= 3 end;

subto c2: sum in I :
if (i mod 2 == 0) then 3 * x[i] else -2 * y[i] end <= 3;

ZIMPL User's guide

2番めの制約式では、左辺だけにif文が効いている. iが偶数なら 3*x_i <=3, 奇数なら-2*y_i<=3;

Special Ordered Set

ZIMPLはSpecial Ordered Setもサポートするようです.

集合とパラメータをファイルからよみこむ

ファイルから集合とパラメータを読み込む形式は、

read filename as template [skip n] [use n] [fs s] [match s] [comment s]
ZIMPL User's guide

filenameは、よみこむファイル名, templateは生成するtupleの文字列よみこむファイル内の各行はフィールドに分けられる。 templateから後ろにある、skip nは、ファイルの最初のn行をスキップさせる. comment sは、ファイル内でコメント行の先頭に使われる文字のリストを指定する。 use nは、ファイル内のn行だけを使うように指定する. 各行のフィールドへの分割のルールは、

スペース、タブ、コンマ、セミコロン、コロンで区切る
double quotesで囲まれたテキストは分割されず、quoteは削除される
フィールドが分割されたはには、その分割した点のまわりのスペースとタブは全て削除される
If the split is due to comma, semicolon or colon, each occurrence of these characters starts a new field.

入力ファイルの例

Hallo;12;3
Moin 7 2
"Hallo, Peter:"; Nice to meet you" 77
ZIMPL User's guide

上の3行は、いずれも3つのフィールドとして読まれる.

set P := { read "nodes.txt" as "<1s>" };
nodes.txt:
Hamburg -> <"Hamburg">
München -> <"München">
Berlin -> <"Berlin">
ZIMPL User's guide

set Q := { read "blabla.txt" as "<1s,5n,2n>" skip 1 use 2 };
blabla.txt:
Name;Nr;X;Y;No -> skip
Hamburg;12;x;y;7 -> <"Hamburg",7,12>
Bremen;4;x;y;5 -> <"Bremen",5,4>
Berlin;2;x;y;8 -> skip
ZIMPL User's guide

<1s,5n,2n>というのは、各行の 1番目のフィールドを文字列として,5番目のフィールドを数字として、 2番目のフィールドを数字として読みこんで、集合Q の1つの要素としなさい、という意味.

set A := { read "test.txt" as "<2n>",<5>,<6> };
param winniepoh[X] :=
read "values.txt" as "<1n,2n> 3n", <1,2> 17, <3,4> 29;
ZIMPL User's guide

ファイル中の全ての値を集合の要素として読み込む場合は、文字列に対しては"<s+>"、数値に対しては"<n+>"を用いる.

set X := { read "stream.txt" as "" };
ZIMPL User's guide

ファイル stream.txtから数値の集合を読み込む.stream.txtの中身は以下のとおり.

1 2 3 7 9 5 6 23 63 37 88 4 87 27
ZIMPL User's guide

この読み込みを実行すると、集合Xは、
X:={1,2,3,7,9,5,6,23,63,37,88,4,87,27};
となる

関数定義

関数定義も可能。関数が返す値は、数値、文字列、boolean, 集合、のいずれかでなくてはならない。関数の引数は、数値か文字列のみ. 数値を返す関数の定義の冒頭は、defnumb, 文字列を返すならdefstrg, booleanを返すならdefbool, 集合を返すならdefset

defnumb dist(a,b) := sqrt(a*a+b*b);
defstrg huehott(a) := if a<0 then "hue" else "hott" end;
defbool wirklich(a,b) := a=10 or b<5;
defset bigger(i) := { <j> in K with j>i };

ZIMPL User's guide

モデリングの例

The diet problem:
食料の集合F, 栄養素の集合N, 食料fに含まれる栄養素nの量を\pi_{in}とする。栄養素nの必要量を\PI_nとする. \delta_fは、食料fがサーブ可能な量とする。食料fの購入代金をc_fとする.この条件で、必要な栄養をとる条件の中で、最も購入代金が安くなる食料の購入量をなにか? 購入する食料の量は整数とする. この問題は、以下の整数計画問題に定式化できる.

ZIMPL User's guide

最初の制約式は、栄養素nを必要量\PI_n以上摂取することをあらわしている. この問題は、ZIMPLを用いて下のように書ける(example/chvatal_diet.zpl)。

set Food := { "Oatmeal","Chicken","Eggs",
		"Milk", "Pie", "Pork"};
set Nutrients := {"Energy", "Protein", "Calcium"};
set Attr := Nutrients + {"Servings", "Price"};
param needed[Nutrients] :=
<"Energy"> 2000, <"Protein"> 55, <"Calcium"> 800;

param data[Food * Attr] :=
	|"Servings","Energy","Protein","Calcium","Price"|
|"Oatmeal"|	4,	110,	4,	2,	3 |
|"Chicken"|	3,	205,	32,	12,	24|
|"Eggs"	  |	2,	160,	13,	54,	13|
|"Milk"	  |	8,	160,	8,	284,	9 |
|"Pie"	  |	2,	420,	4,	22,	20|
|"Pork"   |     2,	260,	14,	80,	19|;

var x[<f> in Food] integer >= 0 <= data[f, "Servings"]; -- (1)
minimize cost: sum < f > in Food : data[f, "Price"] * x[f]; --(2)
subto need: forall <n> in Nutrients do
   sum < f > in Food : data[f, n] * x[f] >= needed[n]; --(3)

ZIMPL User's guide, 番号は追加

このモデルでは、各食料の栄養素を定めるテーブル\pi_{in}の列に、その食料の価格とサーブできる量を持たせている. (1)は、上の数式の2番目の制約に対応している。for all fのところを、xの添え字として書いてしまうのが定石だろう。数式では、変数の上下限も"subject to"の中に書いてしまう場合が多いけれど、 ZIMPLでは、"var"の行に変数の定義と同時に書いてしまう. (2)はsummationの書き方の定石,すなわち、\sum_{f \in F}t[f] は、sum <f> in Food: t[f]と書く.(3)は制約式の書き方の定石. 各栄養素nについて、1本制約式を定義するとき、数式では、 forall n \in Nとなるが、ZIMPLでは、

forall < n > in N do 

   (nに関する制約式)

となる.

これを、

$ zimpl chvatal_diet.zpl

と実行すると、chvatal_diet.lpファイルが生成される。

$ zimpl chvatal_diet.zpl 
****************************************************
* Zuse Institute Mathematical Programming Language *
* Release 3.0.0 Copyright (C)2009 by Thorsten Koch *
****************************************************
*   This is free software and you are welcome to   *
*     redistribute it under certain conditions     *
*      ZIMPL comes with ABSOLUTELY NO WARRANTY     *
****************************************************

Reading chvatal_diet.zpl
Instructions evaluated: 670
Name: chvatal_diet.zpl   Variables: 6   Constraints: 3   Non Zeros: 18
Writing [chvatal_diet.tbl]
Writing [chvatal_diet.lp]

これをsoplexに与えると, diet problemが解かれる.

$soplex chvatal_diet.lp
************************************************************************
* SoPlex --- the Sequential object-oriented simPlex. Release 1.4.2     *
* Copyright (C)  1997-2009 Zuse Institute Berlin                       *
************************************************************************

SoPlex parameters: 
Delta          =     1.000000e-06
Epsilon Zero   =     1.000000e-16
Epsilon Factor =     1.000000e-20
Epsilon Update =     1.000000e-16

algorithm      = Leaving
representation = Column
update         = Forest-Tomlin
pricing        = Steep
starter        = no
simplifier     = MainSM
ratiotest      = Fast
scaling        = bi-Equilibrium / no

Loading LP file chvatal_diet.lp
LP has 3 rows 6 columns 18 nonzeros
LP reading time: 0.00

Solving LP ...
IEQUSC01 Equilibrium scaling LP
IMAISM69 Main simplifier removed 0 rows, 0 columns, 0 nonzeros, 1 col bounds, 0 row bounds
IMAISM74 Reduced LP has 3 rows 6 columns 18 nonzeros
ISOLVE01 iteration =        0	lastUpdate =    0	value = 0.000000e+00
ISOLVE01 iteration =        2	lastUpdate =    2	value = 2.520837e+01
ISOLVE01 iteration =        3	lastUpdate =    1	value = 9.250000e+01
ISOLVE02 Finished solving (status=OPTIMAL, iters=3, leave=2, enter=1, objValue=9.250000e+01)

SoPlex statistics:
  Factorizations : 3
      Time spent : 0.00
  Solves         : 14
      Time spent : 0.00
  solution time  : 0.00
  iterations     : 3

The symmetric travelling salesman problem
G=(V,E)を完全グラフとする. Vは都市の集合で、Eは都市間の枝の集合とする. 巡回路に枝(i,j) \in Eが含まれるとき1、それ以外のとき0をとる0-1変数

ZIMPL User's guide

最初の制約式で出てくる\delta_vは、都市vにつながっている枝を表す. つまりこの和は、都市vに入る枝と出る枝は1本ずつ、ということを表している. 下で示すZIMPファイルでは、出る枝と入る枝の和に分けている:

ZIMPL User's guide

データとしては、都市の数と、それぞれのx,y座標を与える。都市間の距離はユークリッド距離で計算する.データでは与えない. データファイルの中身は以下のように与える

#City 	X	Y
Berlin	5251	1340
Frankfurt	5011	864
Leipzig	5133	1237
Heidelberg	4941	867
Karlsruhe	4901	840
Hamburg	5356	998
Bayreuth	4993	1159
Trier	4974	668
Hannover	5237	972
Stuttgart	4874	909
Passau	4856	1344
Augsburg	4833	1089
Koblenz	5033	759
Dortmund	5148	741
Bochum	5145	728
Duisburg	5142	679
Wuppertal	5124	715
Essen	5145	701
Juena	5093	1158

ZIMPL User's guide

set V := { read "tsp.dat" as "<1s>" comment "#"};
set E := { < i,j > in V * V with i < j };
set P[] := powerset(V);
set K := indexset(P);
param px[V] := read "tsp.dat" as "<1s> 2n" comment "#";
param py[V] := read "tsp.dat" as "<1s> 3n" comment "#";
defnumb dist(a,b) := sqrt((px[a]-px[b])^2 + (py[a]-py[b])^2;
var x[E] binary; //変数の添字は集合でよい 
minimize cost: sum < i,j > in E: dist(i,j) * x[i,j] 
subto two_connected : forall < v > in V do 
   (sum < v,j > in E: x[v,j] ) + (sum <i,v> in E: x[i,v]==2;
subto no_subtuor:
   forall < k > in K with
	card(P[k]) > 2 and card(P[k]) < card(V) -2 do --- (1)
	sum < i, j > in E with < i > in P[k] and < j > in P[k] : x[i,j] < = card(P[k])-1; ---(2)

ZIMPL User's guide, (1),(2)は追加

集合Vは、文字列(都市名を表す)の集合. 変数xの添字は集合でよい。ここでは(i,j)の二つ組を添字にしている。参照するときは、x[i,j]と、括弧<>がないことに注意. px[V], py[V]は,文字列に対して数値を関連付ける"辞書". set P[]:=powerset(V);とset K:=indexset(P)の組に注意. 制約式no_subtourで、"P[k]"と, Kの要素を添字にしてPの要素にアクセスしている. ということは、おそらくP[k]はpowerset(V)のk番目の要素(集合)を表しているのだろう. 集合P[k]に含まれる枝についての和をとるときは、

sum <i, j> in E with <i> in P[k] and <i> in P[k];

とかく.つまり、まずEの中の和をとる、と書いておいて、その後でそのEに条件をつける.

The capacitated facility location problem

工場建設の候補地の集合Pと、商店の集合Sがある。各商店では需要\delta_sがある。候補地のうちのいくつかに工場を建てて、それらの工場から商店に商品を運ぶことで、需要\delta_sを満たしたい。 1つの商店に商品を収める工場は1つである。 pに工場を建てる費用はc_pで、工場pから商店sに荷物をはこぶ費用は一個あたりc_{ps}とする. 工場pには出荷上限があって、それを\kappa_pとする。このとき、各商店の需要を満たすなかで、工場の建設費用と輸送費用の和を最小にする問題は、以下の整数計画問題に定式化できる.

ZIMPL User's guide

最初の制約式は、各商店はちょうど一つの工場に割り当てられることを要求する。 2番目の制約式は、工場pが建設されなければ(z_p=0)、pから出荷される量はどの商店に対しても0である(x_{ps}=0),という制約.3番目の制約は、工場pから出荷される商品の合計量は、上限\kappa_p以下である、という制約. これをZIMPLモデルで書くと、下のようになる.

set PLANTS := { "A", "B", "C", "D"}; // P
set STORES := { 1 .. 9 }; //S
set PS := PLANTS * STORES; 
param building[PLANTS]:= <"A"> 500, <"B"> 600, <"C"> 700, <"D"> 800; //c_p
param capacity[PLANTS]:= <"A"> 40, <"B"> 55, <"C"> 73, <"D"> 90; //\kappa_p
param demand[STORES]:= // \delta_s
<1>10, <2>14,
<3>17, <4>8,
<5> 9, <6>12,
<7>11, <8>15,
<9>16; 
param transport[PS]:=  // c_{ps}
    | 1, 2, 3, 4, 5, 6, 7, 8, 9 |
|"A"| 55, 4, 17, 33, 47, 98, 19, 10, 6 |
|"B"| 42, 12, 4, 23, 16, 78, 47, 9, 82 |
|"C"| 17, 34, 65, 25, 7, 67, 45, 13, 54 |
|"D"| 60, 8, 79, 24, 28, 19, 62, 18, 45 |; 
var x[PS] binary;
var z[PLANTS] binary;

minimize cost: sum < p > in PLANTS : building[p] * z[p]
	     + sum < p,s > in PS : transport[p,s] * x[p,s];

subto assign:
  forall < s > in STORES do 
     sum < p > in PLANTS : x[p,s] == 1;

subto build:
  forall < p,s > in PS do
    x[p,s] <= z[p];

subto limit:
  forall < p > in PLANTS do
    sum < s > in S : demand[s] * x[p,s] <= capacity[p];
ZIMPL User's guide, ただし//以下のコメントは追加

c_{ps}が、transport[PS]に対応する. \delta_sが, demand[STORES]に対応する. c_pが、building[PLANTS]に対応する. \kappa_pが、capacity[PLANTS]に対応する.

ここからは、自分でつくったモデルの例

最小費用流問題

始点sから終点tまでの需要dをみたすようにネットワーク上にものを流したい。各アークには流せる上限cap_{ij}と、単位量あたりの流す費用c_{ij}が与えられている. このとき、もっとも費用を少なく流す流しかたを求めよ.

アーク(i,j)に流す量をx_{ij}と表すと、次の線形計画問題に表せる.

これをZIMPLで書くと、下のようになる.

set V:= {1..4};
set E:={  in V * V};

param cost[E]:= |1,2,3,4|
|1|2,5,1,3|
|2|5,3,4,10|
|3|1,4,7,6|
|4|3,10,6,9|;

param cap[E]:= |1,2,3,4|
|1|20,12,10, 3|
|2|12,9 ,8,10|
|3|10,8 ,7, 6|
|4|3,10 ,6, 9|;

param p[V]:=
<1>-10,<2>0,<3>0,<4>10;

var x[ in E] >= 0 <= cap[i,j];

minimize mincost: sum  in E: cost[i,j]*x[i,j];
subto flow:
  forall  in V do
        (sum  in E: x[i,n])-(sum  in E: x[n,j])==p[n];

SCIPでの解き方

コマンドラインで、scipにlpファイルをオプション-fで与える。すると、問題を解いて結果を表示する

$scip -f facility_location.lp

対話的にも実行可能. 引数なしでscipを実行すると、対話式のSCIPが起動して、SCIPプロンプトが表れる.

SCIP version 1.2.0 [precision: 8 byte] [memory: block] [mode: optimized] [LP solver: SoPlex 1.4.2]
Copyright (c) 2002-2009 Konrad-Zuse-Zentrum fuer Informationstechnik Berlin (ZIB)

user parameter file  not found - using default parameters

SCIP>

ここで、

SCIP>help

を実行すると、命令の一覧とその説明が表示される.対話式SCIPでは、

SCIP> read facility_location.lp
SCIP> optimize

とすると、facility_location.lpに書かれた問題が解かれる. 対話式SCIPでは、情報を見ることができる。たとえば

SCIP> read facility_location.lp

と問題を読みこんでおいて、

SCIP> display problem

とすると、読み込んだ問題の情報が表示される.

運航・物流系トップページへ
 海上技術安全研究所トップページへ

美容整形外科アンチエイジング