Moselを使って、複数のモデルを管理し、パラレルに解く∗
Dash Optimization Whitepaper
Last update 18 October, 2006

Multiple models and parallel solving
with Mosel
Y. Colombani and S. Heipcke
Dash Optimization, Blisworth House, Blisworth, Northants NN7 3BX, UK
Yves.Colombani@dashoptimization.com,
Susanne.Heipcke@dashoptimization.com
18 October, 2006

アブストラクト
この論文では、Moselを使い、複数のモデルを、順次、そして、並行的に解く方法を、いくつかの例で説明します。可能なコンフィギュレーションのすべてを網羅的にリストすることができないので、この小論では、例を使って、Moselモジュールmmjobsの使い方を紹介します。例としては、あるモデルのいくつかのインスタンスをコンカレントに実行する、マスターに、順次、サブモデルを埋め込む、ディコンポジション・
アルゴリズム(Dantzig-Wolfe and Benders decomposition)を実行するなどを使います。また、より専門的な観点から、小論では、モデル管理、コンカレント・モデルのシンクロナイゼーション、共有メモリIOドライ
バーの使い方なども議論します。

∗

この翻訳は、あくまでも、皆様のご参考に供するためのものですので、内容的の厳密性は、保障できません。したがって、英文のWhitepaperもお読みなさることをお勧めいたします。英文は、下記より入手できます。
http://jp.dashoptimization.com/home/services/publications/support̲papers.html

1目次
1 イントロダクション.......................................................................................................... 42 基本的なタスク.................................................................................................................. 6
2.1 サブモデルを実行する ................................................................................................ 6
2.2 サブモデルの実行を停止する...................................................................................... 7
2.3 サブモデルのアウトプット ......................................................................................... 8
2.4 メモリへのコンパイル ................................................................................................ 9
2.5 ランタイム・パラメータ............................................................................................. 9
2.6 複数のサブモデルをランする.................................................................................... 10
2.6.1 Sequential サブモデル ...................................................................................... 10
2.6.2 Parallelサブモデル............................................................................................. 11
2.7 異なるモデル間でのデータのやり取り ..................................................................... 122.7.1 共有メモリドライバーを使う ............................................................................. 122.7.2 メモリ･パイプドライバーを使う........................................................................ 13
3 カラムの生成：異なるモデルを、順番に解く ................................................................. 16
3.1 例題: cutting stock.................................................................................................. 16
3.2 インプリメンテーション........................................................................................... 17
3.2.1 Master model..................................................................................................... 18
3.2.2 Knapsack model ................................................................................................ 19
3.3 結果 .......................................................................................................................... 20
4 いくつかのモデル・インスタンスをパラレルに解く....................................................... 22
4.1 例題：economic lot sizing ........................................................................................ 22
4.1.1 Cutting plane algorithm ................................................................................... 23
4.2 インプリメンテーション .......................................................................................... 24
4.2.1 Master model..................................................................................................... 24
4.2.2 ELS model ......................................................................................................... 26
4.3 Results...................................................................................................................... 29
5 Dantzig-Wolfe decomposition: sequentialとparallel解法を組み合わせる..................... 31
5.1 例題: multi-item, multi-period production planning ............................................. 32
5.1.1 もともとのモデル .............................................................................................. 32
5.1.2 問題の分解.......................................................................................................... 34
5.2 インプリメンテーション .......................................................................................... 36
5.2.1 Master model..................................................................................................... 37
5.2.2 単一工場問題 ..................................................................................................... 40
5.2.3 Master problemのsubroutines .......................................................................... 422

5.3 結果 .......................................................................................................................... 45
6 Benders decomposition: working with several different submodels ............................ 46
6.1 小さい例題 ................................................................................................................ 48
6.2 インプリメンテーション........................................................................................... 48
6.2.1 Master model..................................................................................................... 48
6.2.2 Submodel 1: fixed continuous variables ........................................................... 51
6.2.3 Submodel 2: fixed integer variables ................................................................. 53
6.2.4 Submodel 0: start solution ................................................................................ 55
6.3 Results...................................................................................................................... 57
7 Summary ........................................................................................................................ 58
Bibliography ...................................................................................................................... 59

3

1 イントロダクション
モーゼルのリリース1.6では、モーゼル言語を使い、直接、複数のモデルを扱うことができるようになりました。モジュールmmjobsによって、新たに提供されるようになった機能には、モデル管理機能、イベン
トの待ち行列に基づくコンカレント・モデルのシンクロナイゼーション機能、および、共有メモリIOドライバーが含まれています。以下のリストは、利用可能な機能の概要を記したものです。このモジュールについての詳細は、Mosel Reference Manualのmmjobsの部分をご参照ください。
•

モデル管理:

モデルのソースファイルのコンパイル、bim ファイルのロード、モデルの実行と中断、モデル情報(status, exit code, ID)のリトリーバル、IO ストリームの出力先変更。
•

シンクロナイゼーション・メカニズム:

イベントの送信、リトリーブ、イベント、もしくは、イベントクラスの待ち、イベント情報(class, value, sender model)のリト
リーブ•

共有メモリIOドライバー:

コンカレント・モデル間でデータをやりとりするための、

memドライバーの共有メモリ・バージョン(write access by a single model, read accessby several models simultaneously)、これは、Moselがファイルネームを予期している場所ならどこでも、とくに、initializations blockで使用できる。
•

メモリ・パイプIOドライバー:

コンカレント・モデル間でデータをやりとりするためのメモリ・パイプ(write access by several models, read access by a single model)
で、これは、Moselがファイルネームを予期している場所ならどこでも、とくに、
initializations blockで使用できる。

mmjobsにより、新たに二つの「タイプ」が導入されました。すなわち、ModelとEventです。Model
は、Moselモデルの参照に使われます。モデルの参照を行う前に、bimファイルをロードしてイニ
シアライズしておかなければなりません。
Eventは、Mosel言語の中で、一つのeventを意味します。Eventは、classとvalueで特徴づけられ、モデルとparentモデルの間でやりとりできます。各モデルには、event queueが取り付けられ、このモデルに送られてくるすべてのeventを受け取ります。その管理は、FIFO (First In ‒ First Out)で行なわれます。
この論文の最初の章では、読者に、Moselで、下記のようなモデルを扱うときに、通常行なわれる基本的なタスクを説明します。
•

マスターモデルからサブモデルを実行する：その順序は、コンパイル、ロード、ラン、
ウエイトの順です。

•

サブモデルの実行を止める•

サブモデルからアウトプットを得る： (特に、IVEの下では)サブモデルを静かにする

4

ため、ファイルに出力する•メモリへコンパイルする•

ランタイムパラメータを渡す•

複数のサブモデルのランを行なう‒ 順繰りに‒ 同時並行的に•モデル間でデータをやり取りする：共有メモリIOドライバー、メモリパイプIOドライ
バーを使う

この論文の残りの部分で、mmjobsを使って行なう、もう少し進んだ例題について議論します。
そこでは、それらのインプリメンテーションの仕方の詳細も説明します。これらの例題は、すべて、Dash websiteからダウンロードできます。
•

カラムの生成： Moselユーザガイドに掲載されている二つの別個のカラム生成のモデ
ルを、共有メモリ経由でデータを渡し順番に解く。

•

パラレルで解く：

同一のモデルのいくつかのインスタンスを、別個の解法アルゴリ

ズムパラメータで、コンカレントにランする。改善された解の値は、ランが行なわれているすべてのモデルに送られ、バウンドを更新する。そして、最も早く最適解を得たモデルは、他のモデルをすべて停止させる。
•

Dantzig-Wolfeディコンポジション： コンカレントな方法でサブ問題の一組のインス
タンス（subproblem instances）を解くときの反復手順、更新マスター問題を解く。
データのやりとりは、共有メモリ経由で行なう。

•

Benders ディコンポジション： それぞれに異なる一連のサブ問題を繰り返して解く。
データのやりとりは、共有メモリ経由で行なう。

ここで強調しておきたいのは、multiple modelsとmultiple problemsの違いです。 — Moselは、複数のモデルを扱えますが、これらの個々のモデルは、皆、単一の問題に対応したものです。これは、例えば､ある一つのモデルが、Xpress-Optimizerのようなソルバーを、何回かcallするような場合、 ソルバーは、単一の問題を表現したものを扱い、常に、最後に行なった最適化ランのソリューションのみが得られることを意味します。

5

2基本的なタスク
この章では、Mosel でいくつかのモデルを扱うときに、mmjobsモジュールの機能を使って行な
う基本的なタスクについて説明します。

2.1 サブモデルを実行する
ここで、下記の簡単なtestsub.mosというサブモデルを、そのmaster modelから実行したいとします。
model "Test submodel"
forall(i in 10..20) write(i^2, " ")writelnend-model

ここで、読者の皆さんは、標準的なcompile-load-runという手順に慣れていると前提します。こ
の手順は、Mosel command line、ホストのアプリケーション、別のMoselモデルなど、どこから実行するに関りなく、Moselモデルの実行に、常に、必要な手順です。Moselモデルをmaster modelから実行する場合、compile-load-run-wait という標準的な手順を補強する必要があります。これの背後にある理論的説明は、サブモデルを別の一つのスレッドとして開始し、こうすることで、そのサブモデルが、master model が終わる前に終了する（もしくは、master modelが、サブモデ
ルの結果を利用して、実行を継続する−後述の2.7節を参照）ようにさせるためです。

下のモデルruntestsub.mosは、上のtestsub.mosというモデルを実行するため、compile-load-run-wait という手順を、その基本的な形式で使います。それに対応するMosel のサブルーチンは、mmjobsで定義しますが、これは、usesステートメントの中に含まれている必要があります。サブモデルはそのままです（すなわち、サブモデル自体がこのモジュールの持つ機能を何も使わない限り、サブにmmjobsを入れる必要はありません）。master modelの最後のステートメントdropnexteventについては、もう少し説明が必要でしょう。サブモデルの終了は、
master modelに送られるevent of class EVENT̲ENDで知らされます。waitステートメントは、
master modelの実行を、eventを受け取るまで、一時、休止させます。ここでの例の場合、サブモデルが送る唯一のevent がこの終了メッセージなので、その種類を確認することなく、mastermodelのevent queue から、メッセージを削除します。
model "Run model testsub"
uses "mmjobs"
declarations
modSub: Model
end-declarations! Compile the model file6

if compile("testsub.mos")<>0 then exit(1); end-if
load(modSub, "testsub.bim")run(modSub)
wait
dropnextevent! Load the bim file
! Start model execution! Wait for model termination! Ignore termination event message

end-model

この二つのモデルは、Mosel modelの標準中的な実行方法(Mosel command line, within
Xpress-IVE, or from a host application)のいずれからも、master model を実行することでランされます。Mosel command line から実行する場合、例えば､下のようなコマンドを使います。
mosel -c "exec runtestsub.mos"

その結果として、サブモデルにより、下のようなプリントが出てきます（Xpress-IVEを使っている読者は、2.3節を参照してください）。

100 121 144

169 196 225

256 289 324

361 400

サブモデルの終了ステータスについて、より正確な情報を得たいときは、dropnexteventというステートメントを、下記のステートメントに置き換えてください。これにより、サブモデルにより送られるevent がリトリーブされ、そのクラス(the termination event has the predefinedclass EVENT̲END) 、および、そのバリュー(default value 0)がプリントされます。それに加え、モデルのexit code が表示されます(value sent by an exit statement terminating the model
execution or the default value 0)。
declarations
ev: Eventend-declarations
ev:=getnextevent
writeln("Event class: ", getclass(ev))
writeln("Event value: ", getvalue(ev))
writeln("Exit code : ", getexitcode(modSub))

2.2 サブモデルの実行を停止する
仮に、サブモデルの実行に長い時間が掛かるような場合、master modelを停止することなく、サブモデルを停止したほうが望ましいことがあります。master modelを停止することなく、サブモ
デルを停止するには、下記のようにmaster model(file runsubwait.mos)を修正し、waitステー
トメントにduration (in seconds)を加えてください。 サブモデルが、1秒、実行した後、termination event message を送っていない場合、サブモデルは、stop のcall により停止されます。このとき、model reference が表示されます。

7

model "Run model testsub"
uses "mmjobs"
declarations
modSub: Model
end-declarations! Compile the model file
if compile("testsub.mos")<>0 then exit(1); end-if
load(modSub, "testsub.bim")! Load the bim filerun(modSub)! Start model execution
wait(1)! Wait 1 second for an event
if isqueueempty thenwriteln("Stopping the submodel")stop(modSub)
wait
end-ifdropnextevent! No event has been sent: model still runs
! Stop the model! Wait for model termination! Ignore termination event message

end-model

2.3 サブモデルのアウトプット
これまでの説明に基づき、Xpress-IVEでmaster modelを実行した読者の皆さんは、サブモデルからのアウトプットを何もご覧になれなかったと思います。これは、IVE がmaster model のアウ
トプットしか捉えないからです。この場合、サブモデルのアウトプットにアクセスする一番よい
方法は、このアウトプットの出力先をファイルに変更することです。サブモデルのアウトプットの出力先を変更することは、いくつかの（サブ）モデルを同時並行的にランしている場合にも必要です。なぜなら、スクリーン上に表示されるアウトプットは、これらのモデルからのアウトプットが混ざって表示されるからです。
下記のようなステートメントを使い、サブモデルの中で、アウトプットがプリントされる前に、
直接、アウトプットの出力先を変更するのがよいでしょう。
fopen("testout.txt", F_OUTPUT+F_APPEND)
fopen("tee:testout.txt&", F_OUTPUT)
fopen("null:", F_OUTPUT)! Output to file (in append mode)! Output to file and on screen! Disable all output

ここで、最初のステートメントは、アウトプットをtestout.txtというファイルに出力させ、2行目のステートメントは、testout.txtというファイルに書き込み中、アウトプットをスクリーンに保持し、そして3番目のステートメントは、モデルを完全にsilentにします。アウトプット・
ファイルは、モデルで、プリント・ステートメントの後に、fclose(F̲-OUTPUT)というステート
メントを追加することによりクローズされます。
同様のことが、master modelからも行なえます。その方法は、対応するサブモデルのrun ステー
トメント前に、下記のように、アウトプットの出力先変更を加えます。
setdefstream(modSub, F_OUTPUT, "testout.txt") ! Output to filesetdefstream(modSub, F_OUTPUT, "tee:testout.txt&")! Output to file and on screen

8setdefstream(modSub, F_OUTPUT, "null:")! Disable all output

サブモデルでのアウトプットの出力先変更は、アウトプット・ストリームをリセットすることで終了できます。
setdefstream(modSub, F_OUTPUT, "")2.4 メモリへのコンパイル
Moselファイルfilename.mosのディフォルトのコンパイルにより、filename.bimというバイナリーのモデル・ファイルが生成されます。サブモデルの物理的なBIM fileの生成を避けるため、下の例runsubmem.mosに示されているように、サブモデルをメモリにコンパイルすることもできます。物理的なBIM fileにアクセスするよりも、メモリを使うと、通常、効率性がよくなります。
さらに、master modelが実行されている場所にwrite accessができない場合、この機能は便利です。
model "Run model testsub"
uses "mmjobs", "mmsystem"
declarations
modSub: Model
end-declarations! Compile the model file
if compile("", "testsub.mos", "shmem:testsubbim")<>0 thenexit(1)
end-ifload(modSub, "shmem:testsubbim")! Load the bim file from memory
fdelete("shmem:testsubbim")! ... and release the memory blockrun(modSub)! Start model execution
wait! Wait for model termination
dropnextevent! Ignore termination event message
end-model

compileは3つのアーギュメントが必要です。すなわち、compilation flags (例えば､デバッグには"g")、model file name、および、output file name (ここで、共有メモリドライバーの名前
により、レイベルがプリフィックスされている。モデルをレイベルしたので、コンパイルされたモデルで使われているメモリを解放できますが、それは、fdeleteをコールすることで行ないます。このサブルーチンは、mmjobsに加え、mmsystemというモジュールをロードし、そこで得られます。

2.5 ランタイム・パラメータモデルをランしているときに、（すなわち、モデル自体を修正して、リコンパイルすることなく）
Moselモデルのデータを修正する便利な方法は、runtime parametersを使う方法です。このようなパラメータは、モデルの最初に部分のparametersブロックでdeclareすることです。parameters9ブロックでは、すべてのパラメータにたいしてディフォルト値が与えられており、ディフォルト
値以外の値が指定されていない場合は、モデルの実行時には、これらのディフォルト値が使われます。
下のモデルrtparams.mosについて考えてみてください。このモデルは、4種類のパラメータ（integer、real、string、Boolean）を受け取り、それらの値をプリントします。
model "Runtime parameters"
parameters
PARAM1 = 0
PARAM2 = 0.5
PARAM3 = ’’
PARAM4 = false
end-parameters
writeln(PARAM1, " ", PARAM2, " ", PARAM3, " ", PARAM4)
end-model
Master model“runrtparam.mos”が、この（sub）modelを起動させるには、以下のモデルで行えます。−全てのランタイム パラメータで新しい値が与えられます。
model "Run model rtparams"
uses "mmjobs"
declarations
modPar: Model
end-declarations! Compile the model file
if compile("rtparams.mos")<>0 then exit(1); end-if
load(modPar, "rtparams.bim")! Load the bim file
! Start model execution
run(modPar, "PARAM1=" + 2 + ",PARAM2=" + 3.4 +",PARAM3=’a string’" + ",PARAM4=" + true)
wait! Wait for model termination
dropnextevent! Ignore termination event message
end-model

2.6 複数のサブモデルをランする
master model から、どのようにしてパラメータ化されたモデルをランするかを理解できれば、
master model から、異なるサブモデルのインスタンスを実行するのは簡単です。次の二つの節
で、サブモデルをsequential に、そして、parallel に実行することについて説明します。 説
明を判りやすくするため、ここでの例では、サブモデルは、すべて、単一モデルのパラメータ化されたバージョンを使います。もちろん、単一のmaster model から、異なるサブモデルをコンパイルし、ロードし、ランすることもできます。

2.6.1 Sequential サブモデル
サブモデルのいくつかのインスタンスをシークエンシャルにランするには、下記の例で見るように(file runrtparamseq.mos)、単一モデルのインスタンスに使ったmaster modelに、小さな修正を加えるだけです。ものごとを単純にするため、各々の実行ごとに、単一のパラメータだけをリ

10

セットします。
model "Run model rtparams in sequence"
uses "mmjobs"
declarations
A = 1..10
modPar: Model
end-declarations! Compile the model file
if compile("rtparams.mos")<>0 then exit(1); end-if
load(modPar, "rtparams.bim")! Load the bim file
forall(i in A) dorun(modPar, "PARAM1=" + i)
wait
dropnextevent
end-doend-model! Start model execution! Wait for model termination! Ignore termination event message

サブモデルがコンパイルされ、一度、ロードされ、サブモデル・インスタンスの実行が始まると、現在、実行されているインスタンスが終了するまで、次のインスタンスの開始を待ちます。

2.6.2 Parallel サブモデル
サブモデルのインスタンスをパラレルに実行するには、もう少し、修正を行なうことが必要です。
この場合も、サブモデルのコンパイルは、一回、行なえばよいのですが、しかし、インスタンス
をパラレルでランする数だけ、ロードしなければなりません。今回は、サブモデルをすべて開始
し、それらの終了を待つことになりますから、wait ステートメントを別個のループに移します。
model "Run model rtparams in parallel"
uses "mmjobs"
declarations
A = 1..10
modPar: array(A) of Model
end-declarations! Compile the model file
if compile("rtparams.mos")<>0 then exit(1); end-if
forall(i in A) doload(modPar(i), "rtparams.bim")run(modPar(i), "PARAM1=" + i)
end-do! Load the bim file
! Start model execution

forall(i in A) do
wait
dropnextevent
end-do! Wait for model termination! Ignore termination event message

end-model

スクリーンにサブモデルのアウトプットが現れる順番は、決まっていません。なぜなら、モデル

11

がパラレルにランされているからです。しかし、サブモデルの実行は非常に早いので、よくわからないでしょう。サブモデルrtparams.mos のwritelnステートメントの直前に、wait(1)を加えてみてください。また、モデルの最初に、uses "mmjobs"ステートメントを加えることも必要です。これがすんだら、master modelのいくつかのランのアウトプットを比較してみてください。
おそらく、ランごとに、アウトプットのシークエンスが異なっていることが判るでしょう。

2.7 異なるモデル間でのデータのやり取り
ランタイム・パラメータは、個々のデータの値をサブモデルに渡す方式の一つです。しかし、この方式は、サブモデルに、例えば､データ・テーブルやデータ・セットを渡す方式ではありません。また、これは、サブモデルから情報をリトリーブしたり、実行の最中に、モデル間でデータをやり取りする方式としては使えません。これらのタスクは、すべて、モジュールmmjobsで定義される二つのIO driver、すなわち、shmem driver、および、mempipe driverで対応します。前に述べたように、shmem driver は、one-to-many communication (one model writing, manyreading) のためのものであり、mempipe driver は、many-to-one communicationのためのものです。one model writing、one model readingを行なう場合、どちらのドライバーも使えますが、
shmemのほうが、概念的に使いやすいでしょう。

2.7.1 共有メモリドライバーを使う
shmem共有メモリドライバーを使って、データ・ブロックをメモリに書いたり、メモリから読んだりします。このドライバーの使い方は、物理的なファイルを使うときの方法とよく似ています。
2.4節で、既に、その使い方の例を見ましたが、そこでは、filename はlabel で置き換えられ、"mmjobs.shmem:aLabel"のようなドライバーの名前がprefixされています。モデルにmmjobsモジ
ュールがロードされていたり、メモリに、master modelのような、別のモデルがあるような場合
は、短縮形"shmem:aLabel"を使うこともできます。
異なるモデル間でのデータをやりとりは、initializationsブロックを通して行なわれます。通常、shmem driver は、データをバイナリー形式でセーブするためにraw driver と結合されます。
ここで、2.1節で見た最初のテスト・サブモデルの修正版であるtestsubshm.mosについて考えてみましょう。このモデルは、メモリからindex range を読み、resulting array をメモリに書き出します。
model "Test submodel"
declarations
A: range
B: array(A) of real
end-declarations
initializations from "raw:"
A as "shmem:A"
end-initializations
forall(i in A) B(i):= i^212initializations to "raw:"
B as "shmem:B"
end-initializations
end-model

このサブモデルをランするmasterモデルrunsubshm.mosは、下記のようなものになるでしょう。
model "Run model testsubshm"
uses "mmjobs"
declarations
modSub: Model
A = 30..40
B: array(A) of real
end-declarations! Compile the model file
if compile("testsubshm.mos")<>0 then exit(1); end-if
load(modSub, "testsubshm.bim")! Load the bim fileinitializations to "raw:"
A as "shmem:A"
end-initializations
run(modSub)
wait
dropnextevent! Start model execution! Wait for model termination! Ignore termination event message

initializations from "raw:"
B as "shmem:B"
end-initializations
writeln(B)end-model

サブモデルのランが開始される前に、メモリにindex range が書かれ、ランが終わるとresult
array をリトリーブして、master modelからプリントします。
メモリ･ブロックが使用されなくなったら、とくに、データ・ブロックが大きい場合、fdelete(subroutine provided by module mmsystem)をコールして、これらのブロックを解放することをお勧めします。なぜなら、これらのブロックは、それらのブロックを生成したモデルが終了しても、そして、明確にモデルがアンロードされたとしても、mmjobsモジュールが、Moselセッショ
ンの終了により、アンロードされるまで、残っているからです。このような理由から、mastermodelの終わりに、下記のステートメントを加えます。
fdelete("shmem:A")
fdelete("shmem:B")

2.7.2 メモリ･パイプドライバーを使う
メモリ･パイプIO ドライバーmempipe は、共有メモリドライバーについて見てきたこととは逆に

13

動きます。パイプは、そこに書くことができるようにするためには、まず、openされなければなりません。これを行なうには、initializations to の前に、initializations fromをコールします。サブモデル(file testsubpip.mos)は、次のようになるでしょう。
model "Test submodel"
declarations
A: range
B: array(A) of real
end-declarations
initializations from "mempipe:indata"A
end-initializations
forall(i in A) B(i):= i^2
initializations to "mempipe:resdata"B
end-initializations
end-model

これは、確かに、前のサブモデルと大きく変わるものではありません。しかし、master model(filerunsubpip.mos)の変更は、もっと大きいものです。すなわち、ここでは、インプット･データは、サブモデルが開始されてから、master modelによって書かれます。master modelは、次いで、サブモデルがその実行を終了する前に、結果のデータを読むための新しいパイプをopenします。
model "Run model testsubpip"
uses "mmjobs"
declarations
modSub: Model
A = 30..40
B: array(A) of real
end-declarations! Compile the model file
if compile("testsubpip.mos")<>0 then exit(1); end-if
load(modSub, "testsubpip.bim")! Load the bim filerun(modSub)! Start model execution

initializations to "mempipe:indata"A
end-initializations
initializations from "mempipe:resdata"B
end-initializations
wait
dropnextevent! Wait for model termination! Ignore termination event message

writeln(B)]end-model

14

ファイルが、読むためにパイプをオープンすると、このパイプを通して、要求したデータを受け取るまで、ブロックされたままの状態のままです。したがって、このケースの場合、プログラム・
コントロールのフローは、下記のようになります。
1.

master modelは、サブモデルをスタートする。

2.

サブモデルは、インプット・データ・パイプをオープンして、master modelが、それに、書くのを待つ。

3.

インプット・データが渡されるとサブモデルは実行を続け、master modelは、アウトプット・データ・パイプをオープンして、アウトプットが渡されるのを待つ。

4.

アウトプット・データ・パイプがオープンされると、サブモデルは、アウトプット・データ・パイプに書き、終了する。

5.

master modelは、結果をプリントする。

15

3 カラムの生成：異なるモデルを、順番に解く
この章で使うcutting stock example は、ʻMosel User Guideʼ記載のものです。カラム生成のアルゴリズムと、それをMoselでどのように行なったらよいかについては、読者の皆さんは、このマニュアルを参照してください。
カラム生成のアルゴリズムは、たくさんの変数を持つ線形計画モデルでよく使われます。たくさんの変数があると、問題マトリックスのすべてのカラムをはっきりと示すのが困難だからです。
最初、カラムの数を制限したマトリックスから始め、問題のソリューションを得てから、カラム
生成アルゴリズムを使って、手持ちのソリューションを改善するカラムを追加して行きます。
cutting stock問題のカラム生成アルゴリズムでは、新しいカラム(= cutting pattern)を決定す
るのに、現在のソリューションのdual valueに基づき、knapsack問題を解く必要があります。User
Guideのインプリメンテーションと下記の相違は、このknapsack（サブ）問題を扱う方法が異なる、ということにあります。User Guideのインプリメンテーションでは、Moselのconstraint

hiding functionalityが使い、制約式の部分集合を滑らかにします（blend out）が、下記では、モデルの中で、サブ問題が、そのままで実行されます。この両方が実行することは、まったく同一のアルゴリズムで、また、パーフォーマンスも同程度です。しかし、インスタンスの数が多いと、two-model 方式の方が、少し、効率性がよいようです。なぜなら、各モデルは、(un)hidden
constraintsを選択することなく、正確に、解くべき問題を定義するからです。
ここの例では、問題に加える変化は、すべての最適化のランごとに、問題を完全に、再度、ロ−ドすることが必要になるような変化です。このマルチモデル方式の利点は、main (cutting stock)
problemへ加えられる変化(bound updates)が小さく、すべての最適化のランごとに、この問題を、再度、ロ−ドすることが必要でない場合に、明白に現れてきます。

3.1 例題: cutting stock
ある製紙工場では、固定幅MAXWIDTH のロールを生産し、その後、顧客オーダーに基づき、小さなロールにカットしています。ロールは、NWIDTHS だけの種類の異なるサイズにカットできます。
各サイズi (DEMANDi)ごとのオーダーは、需要として与えられます。この製紙工場の目的関数は、
需要を、できるだけ少ない数の紙ロールで満たし、ロスを最小化することです。
ロールの数を最小にするという目的関数は、現在の手持ちの需要を満たすのに最も良いカット・パターンの部分集合を選択する、として表現できます。手動で、どのようにしたら、すべての可能なカット・パターンを計算したらよいかは定かではないので、パターンの基本セット(PATTERNS1,...,PATTERNSNWIDTH) から始めます。これは、大きい元のロールから、同じ幅の、各小さいロールを、できるだけ多くカットするパターンとなっています。ここで変数usejを定義して、カット・パターンj (j ∈ WIDTHS = {1, … ,NWIDTH})が使われる回数を示すとすると、目的関数は、これらの変数の合計となり、制約式は、各サイズへの需要が満たされることとして定式化されます。

16∑ useminimizej∈WIDTHS

Subject toj∑ PATTERNS ·use ≥ DEMAND
ij

j∈WIDTHSji∀j ∈ WIDTHS : use j ≤ [ DEMAND j / PATTERN ij ], use j ∈ IN

この製紙工場は、カット・パターンの基本セットで、需要を満たせますが、やり方によっては、
大きなロールを必要以上に使うこと、および、オーダー数量以上の小さなロールを生産してしまうことで大きなムダを生み出してしまうこともありえます。したがって、カラム生成ヒューリス
ティックを使い、ムダを小さくするカット・パターンを見つけることが必要です。
その行ない方は、MIP サーチを行なう前に、トップ・ノードで、ヒューリスティックにカラム生
成ループを行なうことです。カラム生成ループを繰り返すごとに、下記のステップを行ないます。

1.

LP問題を解き、ベイシスをセーブする

2.

ソリューションの値を得る

3.

手許にある現在のソリューションに基づき、より有利なカット・パターンを計算する

4.

新しいカラム(= cutting pattern)を生成する: 目的関数、および、対応する需要制約
式に項を追加する

5.

修正した問題とセーブしておいたベイシスをロードする

このループのステップ3 では、下記のインテジャーknapsack問題を解くことが必要です。

∑ C ·xmaximize z =j∈WIDTHSij∑ A ·x ≤ B

j∈WIDTHSjj∀ j ∈ WIDTHS : x j integer
この二番目の最適化問題は、主問題であるcutting stock 問題とは変数を共有していないので、主問題とは独立しています。

3.2 インプリメンテーション
インプリメンテーションは、二つの部分かなっています。すなわち、「cutting stock problem の
定義、および、カラム生成アルゴリズムを持つmaster model (file paperp.mos)」、そして、masterからランされる｢knapsack model (file knapsack.mos)｣です。

17

3.2.1 Master model
cutting stockモデルの主要部分は、下記のようなものです。
model "Papermill (multi-prob)"uses "mmxprs", "mmjobs"
forward procedure column_gen
forward function knapsack(C:array(range) of real,
A:array(range) of real,
B:real,xbest:array(range) of integer): realdeclarations
NWIDTHS = 5
WIDTHS = 1..NWIDTHS
RP: range
MAXWIDTH = 94
EPS = 1e-6

! Number of different widths! Range of widths! Range of cutting patterns! Maximum roll width! Zero tolerance

WIDTH: array(WIDTHS) of real
! Possible widths
DEMAND: array(WIDTHS) of integer! Demand per width
PATTERNS: array(WIDTHS, WIDTHS) of integer ! (Basic) cutting patterns
use: array(RP) of mpvar
soluse: array(RP) of real
Dem: array(WIDTHS) of linctr
MinRolls: linctr ! Objective function! Rolls per pattern! Solution values for variables ‘use’! Demand constraints

Knapsack: Model
end-declarations! Reference to the knapsack model

WIDTH:: [ 17, 21, 22.5, 24, 29.5]
DEMAND:: [150, 96, 48, 108, 227]
! Make basic patterns
forall(j in WIDTHS) PATTERNS(j,j) := floor(MAXWIDTH/WIDTH(j))
forall(j in WIDTHS) do
create(use(j))! Create NWIDTHS variables ‘use’
use(j) is_integer! Variables are integer and boundeduse(j) <= integer(ceil(DEMAND(j)/PATTERNS(j,j)))end-doMinRolls:= sum(j in WIDTHS) use(j)! Objective: minimize no. of rolls! Satisfy all demands

forall(i in WIDTHS)
Dem(i):= sum(j in WIDTHS) PATTERNS(i,j) * use(j) >= DEMAND(i)res:= compile("knapsack.mos")load(Knapsack, "knapsack.bim")
column_gen! Compile the knapsack model! Load the knapsack model! Column generation at top node

minimize(MinRolls)! Compute the best integer solution! for the current problem (including! the new columns)
writeln("Best integer solution: ", getobjval, " rolls")
write(" Rolls per pattern: ")forall(i in RP) write(getsol(use(i)),", ")writeln18

end-model

カラム生成ループごとに、新しいデータでランだけすればよいように、カラム生成ヒューリステ
ィックを開始する前に、knapsack モデルがコンパイルされ、ロードされます。（手順column̲gen
は、User Guide で示されている例と同じですので、ここでは、省略します。）knapsack モデル
は、function knapsackからランされますが、function knapsackは、そのパラメータとして、
knapsack 問題とそのソリューションの値のデータを使います。function knapsackは、すべてのデータを共有メモリにセーブし、次いで、knapsackモデルのランを行なって、共有メモリからソリューションをリトリーブします。返す値は、knapsack 問題の目的関数の値(zbest)です。
function knapsack(C:array(range) of real,
A:array(range) of real,
B:real,xbest:array(range) of integer):real
initializations to "raw:noindex"
A as "shmem:A" B as "shmem:B" C as "shmem:C"
end-initializations
run(Knapsack, "NWIDTHS="+NWIDTHS)
wait
dropnextevent! Start solving knapsack subproblem! Wait until subproblem finishes! Ignore termination message

initializations from "raw:"
xbest as "shmem:xbest" returned as "shmem:zbest"
end-initializations
end-function

二つのモデルを、順繰りに行なうには、（master modelを継続するには、knapsack問題の結果をリトリーブする必要があるので）、runステートメントのすぐ後ろに、手順waitへのコールを追加しなければなりません。そうしておかないと、master modelの実行は、child modelとコンカ
レントリーに継続してしまいます。child modelは、終了すると、
ʻterminationʼ event (an eventof class EVENT̲END)を送ります。しかし、アルゴリズムは、このevent を必要としないので、dropnexteventをコールして、モデルのevent queue から、このeventを取り除きます。

3.2.2 Knapsack model
knapsackモデルの実行は、簡単でステートメント。問題データを、すべて、共有メモリから取り、問題を解いてから、ソリューションを共有メモリにセーブします。
model "Knapsack"uses "mmxprs"
parameters
NWIDTHS=5
end-parameters! Number of different widths

declarations

19

WIDTHS = 1..NWIDTHS
A,C: array(WIDTHS) of real
B: real
KnapCtr, KnapObj: linctr
x: array(WIDTHS) of mpvar
xbest: array(WIDTHS) of integer
end-declarations! Range of widths! Constraint + obj. coefficients! RHS value of knapsack constraint! Knapsack constraint+objective! Knapsack variables! Solution values

initializations from "raw:noindex"
A as "mmjobs.shmem:A" B as "mmjobs.shmem:B" C as "mmjobs.shmem:C"
end-initializations! Define the knapsack problem
KnapCtr:= sum(j in WIDTHS) A(j)*x(j) <= B
KnapObj:= sum(j in WIDTHS) C(j)*x(j)forall(j in WIDTHS) x(j) is_integer! Solve the problem and retrieve the solution
maximize(KnapObj)
z:=getobjval
forall(j in WIDTHS) xbest(j):=round(getsol(x(j)))
initializations to "raw:"
xbest as "mmjobs.shmem:xbest" z as "mmjobs.shmem:zbest"
end-initializations
end-model

このモデルでは、共有メモリドライバー shmem を、モジュールの名前mmjobsでprefixしてあります。このようにするのは、knapsack問題を切り離してランしたい、すなわち、メモリにmmjobs
モジュールをロードするcutting stock master modelなしでランしたい場合のみです。
ドライバーにモジュールの名前を明示的に加えなければならないケースは、それぞれに異なるモジュールで定義された、いくつかのshmemドライバーを区別しなければならないときです。

3.3 結果上のモデルのデータを使い、カラム生成アルゴリズムは、cutting stock問題のLP緩和の値をもともとの177.67から160.95に降ろし、6つの新しいパターンを生成します。MIPは、下記のパターンを使って、161ロールというソリューションを見つけます。

WidthsPattern

17

21

22.5

24

29.5

Usage

3

0

0

4

0

0

1

5

0

0

0

0

3

15

6

0

1

0

3

0

32

8

2

0

0

0

2

75

10

0

2

1

0

1

32

11

0

0

2

2

0

6

20

21

4 いくつかのモデル・インスタンスをパラレルに解く
この章では、どのようにして、いくつかのモデルをパラレルに実行し、これらのモデル間で、ソリューション情報をやり取りするかについて説明します。この仕組みは、Moselをマルチプロセ
ッサー・マシンで使う、つまり、利用可能なプロセッサーの数に対応する数のモデルを解こうとするような場合、特に、興味深いものです。
ここでのアイディアは、（ソリューション・アルゴリズムのパラメータ設定のみが異なる）同一のMIPモデルのいくつかのインスタンスをコンカレントリーにランし、最初のモデルが終了した
ときに、すべてのランを停止するというアイディアです。異なるソリューション・アルゴリズム
が、あるアルゴリズムは（良い）ソリューションを素早く見つける、別のソリューション・アル
ゴリズムは、最も良いソリューションが見つかったとき、オプティマリティの判定に優れているというような意味で、相互補完的であるとするならば、MIPサーチの間に、ソリューション上を交換することによるシナジー効果を期待できると考えてもよいでしょう。
この仕組みを実行するには、まず、master model とparameterizable child modelを定義します。
master modelは、モデルのランを開始し、ソリューションの更新を調整、統合し、また、child model
は、ロードされ、好きなだけの変形（version）でランされます。child modelは、すべて、同一のソルバー(Xpress-Optimizer)を使いますが、別のソルバーも、そのソルバーがサーチとインタラクトするのに必要な機能を定義できれば、使うことができます。

4.1 例題：economic lot sizing
Economic lot sizing (ELS)では、所与の計画期間を対象に生産計画を作成します。この例では、計画期間はTIMES = 1, ･･･ , Tです。期間tの需要はDEMANDpt、製品の種類はp (p ∈ PRODUCTS)
で、需要DEMANDptは、この期間の生産と、前の期間からの繰越在庫で満たされなければなりません。
段取りコストSETUPCOSTtはその期間の生産量と関連し、各期間の生産キャパシティCAPtには制限があります。在庫保有のコストはありません。
期間tでの、製品pの生産量を示すデシジョン変数producept、および、期間tで製品pの段取りが発生するかどうかを示すバイナリー変数setupptを導入します。このバイナリー変数は、段取りがあるときは(setuppt = 1) 、ないときは(setuppt = 0)となります。
この問題の数学的なモデルは下記のようになります。

minimize∑ ( SETUPCOST · ∑ setupt∈TIMEStp∈PRODUCTS

pt +∑ PRODCOST ·produce ）ptp∈PRODUCTS

t

t

s =1

s =1

∀p ∈ PRODUCTS , t ∈ TIMES : ∑ produce ps ≥ ∑ DEMAND ps∀p ∈ PRODUCTS , t ∈ TIMES : produce pt ≤ D ptT ·setup pt
22pt

∀t ∈ TIMES :∑ produce ≤ CAP

p∈PRODUCTS

ptt∀p ∈ PRODUCTS , t ∈ TIMES : setup pt ∈ {0,1}, produce pt ≥ 0

目的関数は、総コストを最小にすることです。2行目の制約式で、期間0かたら期間tまでの製品pの生産が、この期間の製品pの需要を満たさなければならないことを定式化します。それに続く
制約式は、「期間tに生産が行なわれれば、期間tには段取りもある」こと、ここで、Dptlは、期間tかたら期間lまでの製品pの需要を意味します。また、期間tの生産キャパシティには制限があること、そして、変数setupptはバイナリー変数であることを意味しています。

4.1.1 Cutting plane algorithm
ELS でよく知られた不等式は、いわゆる(l, S)-inequalities [PW94]です。Dptlが期間tから期間lまでの製品pの需要を意味するき、各期間l、および、期間1かたら期間lの各部分集合期間Sに対して、(l, S)-inequalityは、下記のようになります。
l

l

t =1
t∈s

t =1
t∉s

∑ produce pt + ∑ D ptl·setup pt ≥ D p1l
この不等式は、Sに含まれる期間の実際の生産producept＋残りの期間（Sに含まれていない期間）の潜在的な最大生産Dptl · setupptは、少なくとも、期間1かたら期間lまでの総需要に等しくなけ
ればならないということを言っています。
これらの(l, S)-inequalities に基づき、下記のcutting plane algorithm を作成することがで
きます。

1.

LP問題を解く

2.

下記の不等式を調べて、(l, S)-inequalitiesの違反を識別するl∑ min( produce , D ·setup ) ≥ D
t =1

pt

pti

pt

p1l

3.

違反している（満たされていない）不等式を、新たなcutとして問題に加える

4.

LP問題を、再度、解く

このアルゴリズムをどのように作るかについては、無数のオプションがあります。例えば､•Cutをルート・ノードの中のみに生成するか、もしくは、サーチの中でも生成するか(Cut-and-Branch vs Branch-and-Cut)•

ノードで生成するcut generation passの数（例えば､passを一つ、もしくは、上の2.
から4.のステップをループさせ、cutが生成されないようになるまで繰り返す）

23•

cut生成を、サーチ･ツリーの中で、どれくらいまで深く行なうか (ある程度の深さまでか、それとも、すべてのノード？)•(l, S)-cutsだけを使うか、それとも、他のもの（例えば､ソルバーにより生成される
default cuts）と組み合わせるか

下記の(l, S)-cut生成アルゴリズムの実行では、トップ・ノードのみでcut生成を(TOPONLY =true)、一つ、もしくは、いくつかのrounds of cuts (SEVERALROUNDS = true)を行なうように設定しましょう。

4.2 インプリメンテーション

mmjobsでは、eventは、あるchildから別のchildに、直接、送られることはなく、必ず、「parentとchildの対」の間でやりとりします。したがって、「solution found」というメッセージは、
parent modelに送られる必要があり、次いで、parent modelは、このメッセージをその他のすべてのchild modelに送ります。
ここで強調しておきたいもう一つの点は、ELS model fileのコンパイルは1回のみ行われるということ、しかし、メモリにロードされるインスタンスの数は、ランしたいchild modelの数と対応していなければならないということです。

4.2.1 Master model
master modelは、child modelをコンパイルし、ロードし、ランをします。そして、ソリューシ
ョンの更新を調整します。ソリューションの更新には、注意が必要です。なぜなら、新しいソリューションは、必ずしも、child modelにより、これまでに見つけられ、報告されているソリューションよりも良いとは限らないからです。
例えば､2つのモデルが、ほとんど同時に、ソリューションを見つけたとすると、master modelに早く着いたソリューションがよりよいソリューションかも知れません。そのような場合、早く着いたソリューションが、続いてきたソリューションにより、書き換えられてはいけません。
ランの終わりにソリューション・ディスプレーへの表示を良いものにするために、master model
は、problem data from fileの一部を読み込みます。
model "Els master"uses "mmjobs"
parameters
DATAFILE = "els5.dat"
T = 45P=4
end-parameters
declarations
RM = 0..5
TIMES = 1..TPRODUCTS = 1..P

! Range of models! Time periods! Set of products

24solprod: array(PRODUCTS,TIMES) of real ! Sol. values for var.s produce
solsetup: array(TIMES) of real
! Sol. values for var.s setup
DEMAND: array(PRODUCTS,TIMES) of integer ! Demand per period
modELS: array(RM) of Modelmodid: array(set of integer) of integer
NEWSOL = 2
Msg: Event
end-declarations! Models! Model indices! Identifier for "sol. found" event! Messages sent by models! Compile, load, and run models M1 and M2
M1:= 1; M2:=2
res:= compile("elsp.mos")load(modELS(M1), "elsp.bim")load(modELS(M2), "elsp.bim")
forall(m in RM) modid(getid(modELS(m))):= mrun(modELS(M1), "ALG="+M1+",DATAFILE="+DATAFILE+",T="+T+",P="+P)
run(modELS(M2), "ALG="+M2+",DATAFILE="+DATAFILE+",T="+T+",P="+P)
objval:= MAX_REAL
algsol:= -1; algopt:= -1
repeatwait! Wait for the next event
Msg:= getnextevent! Get the event
if getclass(Msg)=NEWSOL then! Get the event class
if getvalue(Msg) < objval then
! Value of the event (= obj. value)
algsol:= modid(getfromid(Msg))! ID of model sending the event
objval:= getvalue(Msg)
writeln("Improved solution ", objval, " found by model ", algsol)
forall(m in RM | m <> algsol) send(modELS(m), NEWSOL, objval)
else
writeln("Solution ", getvalue(Msg), " found by model ",modid(getfromid(Msg)))
end-ifend-ifuntil getclass(Msg)=EVENT_END ! A model has finishedalgopt:= modid(getfromid(Msg))
forall(m in RM) stop(modELS(m))! Retrieve ID of terminated model! Stop all running models! Retrieve the best solution from shared memory
initializations from "raw:noindex"
solprod as "shmem:solprod"+algsol
solsetup as "shmem:solsetup"+algsol
end-initializations
initializations from DATAFILE
DEMAND
end-initializations! Solution printing
writeln("Best solution found by model ", algsol)
writeln("Optimality proven by model ", algopt)
writeln("Objective value: ", objval)
write("Period setup ")forall(p in PRODUCTS) write(strfmt(p,-7))
forall(t in TIMES) do
write("¥n ", strfmt(t,2), strfmt(solsetup(t),8), " ")forall(p in PRODUCTS) write(strfmt(solprod(p,t),3), " (",DEMAND(p,t),")")
end-dowriteln25

end-model

ここでのインプリメンテーションでは、array modidを定義します。これは、モデルで使われているmodel index と、Moselの中のモデル内部IDの間の通信のやり取りを確立するためのものです。child modelがmaster modelにeventを送ると、function getfromidにより、かならず、そのIDがリトリーブされ、対応するmodel indexがストアされ、これを使って、後で、ソリューショ
ンのプリントを行ないます。

4.2.2 ELS model
ELSのchild modelは、モデルが分離して実行できるように書かれます。特に、モデルは、それぞれ、ファイルからのデータの完全なinitializationを行ないますが、これは、master modelの大きな効率性を保持するためで、データは共有メモリ経由で、child modelに渡されます。（しかし、ここでの例では、データ・ハンドリングの時間は、ソリューション・アルゴリズムのランに必要な時間と較べ、ごく僅かです。）
ELS modelの主要部分には、モデル自体の定義とソリューション・アルゴリズムの選択が入っています。
model Els
uses "mmxprs","mmjobs"
parameters
ALG = 0
DATAFILE = "els4.dat"
T = 60P=4
end-parameters! Default algorithm: no user cuts

forward procedure tree_cut_gen
forward public function cb_node: boolean
forward public function cb_updatebnd(node:integer): integer
forward public procedure cb_intsol
declarations
NEWSOL = 2
EPS = 1e-6
TIMES = 1..TPRODUCTS = 1..P! "New solution" event class! Zero tolerance! Time periods! Set of products

DEMAND: array(PRODUCTS,TIMES) of integer ! Demand per period
SETUPCOST: array(TIMES) of integer! Setup cost per period
PRODCOST: array(PRODUCTS,TIMES) of real ! Production cost per period
CAP: array(TIMES) of integer! Production capacity per period
D: array(PRODUCTS,TIMES,TIMES) of integer ! Total demand in periods t1 - t2
produce: array(PRODUCTS,TIMES) of mpvar ! Production in period t
setup: array(TIMES) of mpvar! Setup in period t
solprod: array(PRODUCTS,TIMES) of real
solsetup: array(TIMES) of real

! Sol. values for var.s produce! Sol. values for var.s setup

Msg: Event! An event26

end-declarations
initializations from DATAFILE
DEMAND SETUPCOST PRODCOST CAP
end-initializations
forall(p in PRODUCTS,s,t in TIMES) D(p,s,t):= sum(k in s..t) DEMAND(p,k)! Objective: minimize total cost
MinCost:= sum(t in TIMES) (SETUPCOST(t) * setup(t) +sum(p in PRODUCTS) PRODCOST(p,t) * produce(p,t) )! Production in period t must not exceed the total demand for the! remaining periods; if there is production during t then there! is a setup in t
forall(t in TIMES)sum(p in PRODUCTS) produce(p,t) <= sum(p in PRODUCTS) D(p,t,T) * setup(t)! Production in periods 0 to t must satisfy the total demand! during this period of time
forall(p in PRODUCTS,t in TIMES)sum(s in 1..t) produce(p,s) >= sum (s in 1..t) DEMAND(p,s)! Capacity limits
forall(t in TIMES) sum(p in PRODUCTS) produce(p,t) <= CAP(t)forall(t in TIMES) setup(t) is_binary! Variables setup are 0/1

setparam("zerotol", EPS/100)! Set Mosel comparison tolerance
SEVERALROUNDS:=false; TOPONLY:=false
case ALG of
1: do
setparam("XPRS_CUTSTRATEGY", 0)
setparam("XPRS_HEURSTRATEGY", 0)
end-do2: do
setparam("XPRS_CUTSTRATEGY", 0)
setparam("XPRS_HEURSTRATEGY", 0)
setparam("XPRS_PRESOLVE", 0)
end-do3: tree_cut_gen
4: do
tree_cut_gen
SEVERALROUNDS:=true
end-do5: dotree_cut_gen
SEVERALROUNDS:=true
TOPONLY:=true
end-doend-case! No cuts! No heuristics! No cuts! No heuristics! No presolve! User branch&cut (single round)! User branch&cut (several rounds)! User cut&branch (several rounds)! Parallel setup
setcallback(XPRS_CB_PRENODE, "cb_updatebnd") ! Node pre-treatment callback
setcallback(XPRS_CB_INTSOL, "cb_intsol") ! Integer solution callback
setparam("XPRS_SOLUTIONFILE",0)! Do not save solutions to file! Solve the problem
minimize(MinCost)
end-model

27

手順tree̲cut̲genにより、ユーザのcut生成ルーチンがセットアップできます。これは、branch-and-bound searchのトップノードにおいてのみcutを生成するか(TOPONLY) 、ノードあたり、いくつかのcutを生成するか(SEVERALROUNDS)のどちらかです。ここでは、cut生成ルーチンcb̲node自体の定義は省かれています。
procedure tree_cut_gen
setparam("XPRS_HEURSTRATEGY", 0)
setparam("XPRS_CUTSTRATEGY", 0)
setparam("XPRS_PRESOLVE", 0)
setparam("XPRS_EXTRAROWS", 5000)
setcallback(XPRS_CB_CUTMGR, "cb_node")
end-procedure! Switch heuristics off! Switch automatic cuts off! Switch presolve off! Reserve extra rows in matrix! Define the cut manager callback

コンカレントにランされているchild model間の情報のやり取りは、2つの部分を持っています。
すなわち、(a)見つけられたインテジャー・ソリューションは、セーブされ、master modelに伝
えられなければならない。(b)master modelから送られるbound updates は、サーチに組み込まれなければならない、です。Xpress-Optimizer は、ユーザ構造の中にソリューションをセーブするための固有なinteger solution callbackを提供しています。An obvious place for boundupdates in nodes is the cut-manager callback function. ノードでのbound updatesの場所は、

cut-manager callback functionです。しかし、このfunctionが、もし、アルゴリズムの何らかの設定で、既に、使われてしまっているときは、すべてのノードで呼び出せる別のcallback
function、すなわち、node pre-treatment callbackを使います。
! Update cutoff value
public function cb_updatebnd(node:integer): integer
if not isqueueempty thenrepeatMsg:= getnextevent
until isqueueempty
newcutoff:= getvalue(Msg)
setparam("XPRS_MIPABSCUTOFF", newcutoff)
if (newcutoff < getparam("XPRS_LPOBJVAL")) thenreturned:= 1! Node becomes infeasible
end-ifend-ifend-function! Store and communicate new solution
public procedure cb_intsol
objval:= getparam("XPRS_LPOBJVAL")! Retrieve current objective value
cutoff:= getparam("XPRS_MIPABSCUTOFF")if(cutoff > objval) thensetparam("XPRS_MIPABSCUTOFF", objval)
end-if! Get the solution values
forall(t in TIMES) do
forall(p in PRODUCTS) solprod(p,t):=getsol(produce(p,t))solsetup(t):=getsol(setup(t))
end-do28

! Store the solution in shared memory
initializations to "raw:noindex"
solprod as "shmem:solprod"+ALG
solsetup as "shmem:solsetup"+ALG
end-initializations! Send "solution found" signal
send(NEWSOL, objval)
end-procedure

Figure 1:主要部分が、アンギュラー構造を持つ係数マトリックス」

bound update callback function は、event queueにeventが入っているかどうかをチェックします。そして、event が入っていると、event queueのeventを、すべて、取り出して、最後のeventの値を新しいcutoff valueとします。event queue を空にするループの背後にある論拠は、mastermodelは、最後にチェックした時点以後に、いくつかの改善されたソリューション値を送信してきているかもしれないが、最も良い値は最後に送信されてきたものである、すなわち、event
queue の最後の値だ、ということです。integer solution callbackは、共有メモリにソリュー
ションの値を書き、identifier of the model (= value of ALG)を加えます。これにより、同時に、ソリューションの値を書き込んだかも知れない2つのchild modelが、同意ツリーのメモリ領域を使わないようにします。

4.3 Results
2つのモデルを使って行なったランは、下記のようなログを生成するでしょう。ここで、サーチ

29

を終わらせたモデルは、最適解を見つけたモデルとは異なることに注目してください。
Improved solution 1283 found by model 2
Improved solution 1250 found by model 2
Improved solution 1242 found by model 1
Improved solution 1236 found by model 2
Improved solution 1234 found by model 2
Best solution found by model 2
Optimality proven by model 1
Objective value: 123430

5 Dantzig-Wolfe decomposition: sequential と parallel 解法を組み合わせる
Dantzig-Wolfe decomposition (詳しくは [Teb01] をお読みください)は、一つの解法です。こ
の解法は、比較的、少数の制約式を取り除いてやると、問題が、それぞれ独立した問題に分解できるような問題で有効な解法です。このことは、Figure 1のように、制約式マトリックスの行と列を並べ替えることができることを意味し、そこでは、ノンゼロ係数は、灰色の領域のみに位置します。このようなprimal block angular structure は、Xpress-IVEで問題マトリックスを視覚化すると、即座に、はっきり見えるようになります。しかし、多くの場合、期間、製品、工場の位置などのようなcommon index (sub)setでグループ化して、制約式の定義を組み立てなおすことが必要になるでしょう。
いくつかの、もしくは、すべてのsubproblemの変数を含む制約式（目的関数を含む）は、global

constraints (also: common, linking, or master constraints)と呼ばれます。これらの制約式
は、master problemを定式化するのに使われます。係数マトリックスの対角線上にある個々のブ
ロックは、master problemの統制の下に、pricing subproblemsとして解かれます。そして、masterproblem を解くことにより、もともとの問題の解を得ます。master problem は、(pricing
problemsの基本的な実行可能解とun-bounded directionsの集合で定義される)たくさんの変数を持つので、私たちは、これらの変数の小さな部分集合を対象に、restricted master problemを使います。restricted master problemの変数の「active set of variablesに入れる変数」は、
pricing subproblemsを解くことで決められます。restricted master problemのdual valueに基づくpricing problemの目的関数により、各extreme pointでの目的関数の値は、そのextreme
pointに対応するmaster problemの変数のreduced cost (or price)になります。
最大化問題では、修正pricing problemを解くことにより、maximum reduced costを持つ基本実
行可能解が生成されます。あるextreme pointでの目的関数の値が正ならば、master problemに、対応するmaster problem変数が追加されます。そして、すべてのextreme pointでのminimum
objective valueが負であるならば、master problemの変数で、その時点でのmaster problemのソリューションを改善できる変数は存在しません。
Dantzig-Wolfe decomposition の計算上の利点は、pricing problemの大量の計算量をこなすときに生れます。pricing problemは、もともとの問題よりも大幅に小さく、したがって、解くのが容易になります。この小論のテーマとの関連で興味深い一つの側面は、subproblemが相互に独立であることで、したがって、subproblemを同時並行的に解くことができる、ということです。
このdecompositionによるアプローチの潜在的な欠点は、master problemのサイズがとても大きいことです。master problemは、もともとの問題よりも制約式の数は少ないのですが、変数の数が大幅に多くなります。通常、すべての変数を明示的に生成することは求められませんが、masterproblem の実行可能領域は、もともとの問題よりも複雑なので、ソリューション・パスが長くなります。さらに、master problemの変数をダイナミックに生成することから、数値計算上の問題（numerical problems）が発生することがあります。
decompositionアプローチのパーフォーマンスには、多数の要因が影響します。したがって、特31

定のアプリケーションごとに、このソリューションが適しているかどうかを実験的に計算してみることが必要でしょう。このようなテストには、与えられた問題を異なった方法で分解してみることも含まれています。問題の分解の定義ための一般的なルールの一つは、global constraints
を少なくすることを狙うべきである、ということです。なぜなら、master problemを素早く解けることが重要だからです。さらに、pricing problemは、degeneracyによる計算上の問題を回避するため、pricing problem それ自体で、良く定式化された問題であるように構築されるべきです。

5.1 例題: multi-item, multi-period production planning
Coco社は、二つの種類の粉末のココアの素を生産できる工場を二つ持っています。最初の工場は、1ヶ月あたり、合計400トンのキャパシティを持ち、二番目の工場は、合計500トンのキャパシテ
ィを持っています。マーケッティング部門は、次の4ヶ月で売れる製品別の最大販売量の推定値、および、予想販売価格を生産部門に寄越します。生産にはいくつかの原材料が必要ですが、最終製品1トンあたりに必要な原材料所要量は判っています。期末に、最終製品、および、原材料を工場で貯蔵しておき、次の期に繰り越せますが、その保有には、1期」、1トンあたり、一定のコストが発生します。工場での原材料の貯蔵能力は、最大300トンです。利益を最大にするには、
この2つの工場を、どのように操業したらよいでしょうか。

5.1.1 もともとのモデル
それぞれ、PRODSで最終製品の、FACTで工場の、RAWで原材料の、そして、TIMES = {1, . . . , NT}
で対象とする期間の集合を表すものとします。
また、makepftというデシジョン変数を定義して、製品pの、工場fにての、期間tでの生産量を表すものとします。さらに、ある期から次の期への移行を定式化し、異なるタイプのコストを考慮するため、いくつかの変数を追加します。sellpft で、製品pの、工場fでの、期間tでの販売量を表し、buyrftで、原材料rの、工場fで購入される、期間tでの購入量を表し、また、それぞれ、期間t = 1, . . . , NT +1を対象に定義されるpstockpft 、rstockrft で、期間tの期初に、工場fで保有されている最終製品、および、原材料を表すものとします。
さらに、MXSELLptで、期間tでの、製品pの最大販売可能量を表し、MXMAKEfで、工場fのキャパシ
ティの制限を、MXRSTOCKで原材料貯蔵キャパシティを表すものとします。

IPSTOCKpfで、計画初期に製品pの、工場fで、貯蔵されている数量を、また、IRSTOCKrfで、計画
初期に原材料rの、工場fで、貯蔵されている数量を表すものとします。さらに、CPSTOCK 、CRSTOCK
で、それぞれ、最終製品、および、原材料の単位あたりの貯蔵コストを表します。
利益最大の目的関数は、売上(REVp)から生産コスト(CMAKEpf )、原材料購入コスト(CBUYrt)、最終
製品、および、原材料の貯蔵コストを差し引いたものです。

32∑ ∑ ∑ REV ·sell

maximizep∈PRODS f ∈FACT t∈TIME−∑ ∑ ∑ CMAKE ·makeptp∈PRODS f ∈FACT t∈TIME

−

∑ ∑

pft

pt

− ∑

∑ ∑ CBUY ·buyrtr∈RAW f ∈FACT t∈TIME

NT +1

rft

NT +1

∑ CPSTOCK·pstock pft − ∑p∈PRODS f ∈FACT t = 2pft∑ ∑ CRSTOCK·rstockr∈RAW f ∈FACT t = 2

rft

期間を跨いで製品を貯蔵することを可能にするには、次の 3 種類の制約式が必要となります。す
なわち、最終製品、原材料の在庫バランスの制約式(PBalpft)、および、(RBalrft)、原材料貯蔵キャパシティの制限の制約式(MxRStockft)です。
製品pのt + 1 期の期初在庫pstockp,f,t+1は、製品pのt 期の期初在庫＋t 期の生産数量−t 期の販売数量です。原材料rのt + 1 期の期初在庫rstockr,f,t+1は、原材料rのt 期の期初在庫＋t 期の購入数量−t 期の使用数量です。

∀p ∈ PRODS , ∀f ∈ FACT , ∀t ∈ TIME : PBal pft := pstock p , f ,t +1 = pstock pft + make pft − sell pft∀r ∈ RAW , ∀f ∈ FACT , ∀t ∈ Time :RBal rft := rstock r , f ,t +1 = rstock rft + buy rft −∑ REQ ·makep∈PRODSprpft∀f ∈ FACT , ∀t ∈ {2,..., NT + 1 } : MxRStock ft := ∑ rstock rft ≤ MXRSTOCKr∈RAW

さらに、2種類のキャパシティ制約式が必要です。すなわち、工場の生産キャパシティは制限されており(制約式 MxMakeft)、各期間で販売できる最終製品の最大数量は制限されています(制約式MxSellft)。
下記は、この問題を数学的に表現したものです。最終製品、および、原材料のt = 1期の期初在
庫レベルは、与えられた値に固定されています。

maximize∑ ∑ ∑ REV ·sellp∈PRODS f ∈FACT t∈TIME−∑ ∑ ∑ CMAKE ·makep∈PRODS f ∈FACT t∈TIME

−

∑ ∑

pt

NT +1

pft

pt

− ∑

∑ ∑ CBUY ·buyr∈RAW f ∈FACT t∈TIME∑ CPSTOCK·pstock pft − ∑p∈PRODS f ∈FACT t = 2

pft

rft

NT +1

∑ ∑ CRSTOCK·rstockr∈RAW f ∈FACT t = 2

33

rt

rft

∀p ∈ PRODS , ∀f ∈ FACT , ∀t ∈ TIME : PBal pft := pstock p , f ,t +1 = pstock pft + make pft − sell pft∀r ∈ RAW , ∀f ∈ FACT , ∀t ∈ Time :RBal rft := rstock r , f ,t +1 = rstock rft + buy rft −∀p ∈ PRODS , ∀t ∈ TIME : MxSell pt :=∀f ∈ FACT , ∀t ∈ TIME : MxMake ft :=∑ REQ ·makep∈PRODSprpft∑ sellpft≤ MXSELL p∑ makepft≤ MXMAKE f

f ∈FACTp∈PRODS∀f ∈ FACT , ∀t ∈ {2,..., NT + 1 } : MxRStock ft := ∑ rstock rft ≤ MXRSTOCKr∈RAW∀p ∈ PRODS , ∀f ∈ FACT : pstock pf 1 = IPSTOCK pf∀r ∈ RAW , ∀f ∈ FACT : rstock rf 1 = IRSTOCK rf∀p ∈ PRODS , ∀f ∈ FACT , ∀t ∈ TIME : make pft ≥ 0, sell pft ≥ 0∀r ∈ RAW , ∀f ∈ FACT , ∀t ∈ TIME : buy rft ≥ 0∀p ∈ PRODS , ∀f ∈ FACT , ∀t ∈ {1,..., NT + 1 } : pstock pft ≥ 0∀r ∈ RAW , ∀f ∈ FACT , ∀t ∈ {1,..., NT + 1 } : rstock rft ≥ 0

5.1.2 問題の分解
これから、上で説明したモデルを、生産拠点で分解します。この問題を分解する方法は、ここに示す方法だけではないことに注意してください。この問題は、製品ごとに、または、期間で、分解することもできます。しかし、これらの二つの分解の方法は、工場による分解に比較して、
global constraints の数が多くなります。このことは、master problem を解くのが、もっと、難しくなることを意味しています。
工場fごとに、下記のサブモデルを得ます（boundsの形式での販売数量の制約式MxSellをサブモデルにいれることは必須ではありませんが、もっと正確な、もしくは、もっと速いソリューションが得られるかもしれません。）

maximize∑ ∑ REV ·sellp∈PRODS t∈TIME

34

pt

pt

∑ ∑ CMAKE ·make − ∑ ∑ CBUY ·buy−pp∈PRODS t∈TIME

−

pt

NT +1

NT +1

p∈PRODS t = 2

r∈RAW t = 2∑

rtr∈RAW t∈TIME∑ CPSTOCK·pstock pt − ∑rt∑ CRSTOCK·rstockrt∀p ∈ PRODS , ∀t ∈ TIME : PBal pt := pstock p ,,t +1 = pstock pt + make pt − sell pt∀r ∈ RAW , ∀t ∈ Time : RBal rt := rstock r ,t +1 = rstock rt + buy rt −∀t ∈ TIME : MxMaket :=∑ REQ ·makep∈PRODSprpt∑ make ≤ MXMAKE
pt

p∈PRODS∀t ∈ {2,..., NT + 1 } : MxRStock t := ∑ rstock rt ≤ MXRSTOCKr∈RAW∀p ∈ PRODS , ∀t ∈ TIME : MxSell pt := sell pt ≤ MXSELL p∀p ∈ PRODS , : pstock p1 = IPSTOCK p∀r ∈ RAW : rstock r1 = IRSTOCK r∀p ∈ PRODS , ∀t ∈ TIME : make pt ≥ 0, sell pt ≥ 0∀r ∈ RAW , ∀t ∈ TIME : buy rt ≥ 0∀p ∈ PRODS , ∀t ∈ {1,..., NT + 1 } : pstock pt ≥ 0∀r ∈ RAW , ∀t ∈ {1,..., NT + 1 } : rstock rt ≥ 0

master problemは、global constraints、すなわち、販売数量を制限する制約式MxSellのみを持っています。（ここでは、変数が非負であることを示す制約式は省いてあります。）

maximize∑ ∑ ∑ REV ·sellp∈PRODS f ∈FACT t∈TIME−∑ ∑ ∑ CMAKE ·makep∈PRODS f ∈FACT t∈TIME

−

∑ ∑

pt

pft

pt

− ∑

∑ ∑ CBUY ·buyr∈RAW f ∈FACT t∈TIME

NT +1∑ CPSTOCK·pstock pft − ∑p∈PRODS f ∈FACT t = 2pft∑ ∑ CRSTOCK·rstock35∑ sellf ∈FACT

rft

NT +1

r∈RAW f ∈FACT t = 2∀p ∈ PRODS , ∀t ∈ TIME : MxSell pt :=

rt

pft

≤ MXSELL prft

分解アルゴリズムでは、master problemのデシジョン変数は、下記のようなサブ問題で生成されるソリューション（proposal)で表現されます。

∀p ∈ PRODS , ∀f ∈ FACT , ∀t ∈ TIME : sell pft =nPROPf∑ Pr op _ sell
k =1

pftk

weight fk

ここで、Prop̲sellpftkは、サブ問題fにより生成されたk番目のproposalの中の変数sellptの解の値を意味しています。そして、Prop̲costfkは、このproposalのobjective valueです。すべてのサブ問題fにたいして、convexity制約式Convexfを変数weightfkに追加する必要があります。

max imize ∑nPRop ff ∈FACT∀p ∈ PRODS , ∀t ∈ TIME : MxSell pt :=∑ Pr op _ cos t ·weight
fk

k =1

fk

nPROPf∑ ∑ Pr op _ sellf ∈FACT

k =1∀f ∈ FACT : Convex f :=pftk

·weight fk ≤ MXSELL pnPROPf∑ weight = 1
k =1

fk

∀f ∈ FACT , ∀k ∈ 1,..., nPROPf : weight fk ≥ 0

以下、この問題をmodified master problem と呼びます。これ以上、技術的な詳細に入ることはしませんが、ここでは、もともとの問題とmodified master problem との間には、情報のやりとりがある、ということだけを述べておきます。

5.2 インプリメンテーション
分解アルゴリズムには、下記のフェーズがあります。
•

Phase 1: modified master problemの実行可能解を得るために、一組のproposalを生成する。

•

Phase 2: modified master problemを最適化する。

•

Phase 3:もともとの問題のソリューションを計算する。

フェーズ1と2でかれるサブ問題解は、modified master problem のdual valueの合計を目的関数としてとります。空のmaster problemから出発したり、ミスリーディングなrandom dual
informationに影響されないようにしたりするために、フェーズ0として、crashフェーズを追加します。これにより、各サブ問題に、もともとの目的関数を使って、一つのproposalを生成します。

36

5.2.1 Master model
Below follows the body of the master model file. The definitions of the function bodies will be shown
later in Section 5.2.3.
model "Coco3 Master"uses "mmxprs","mmjobs","mmsystem"
parameters
DATAFILE = "coco2.dat"
ALG = 0 ! 0: stop phase with 1st failed subpb.
! 1: stop when all subprob.s failend-parameters
forward procedure process_sub_result
forward procedure solve_master(phase:integer)
forward procedure process_master_result
forward function calc_solution:real
forward procedure print_solution
declarations
PHASE_0=2
PHASE_1=3
PHASE_2=4
PHASE_3=5
EVENT_SOLVED=6
EVENT_FAILED=7
EVENT_READY=8
NPROD, NFACT, NRAW, NT: integer
end-declarations! Event codes sent to submodels! Event codes sent by submodels

initializations from DATAFILE
NPROD NFACT NRAW NT
end-initializations
declarations
PRODS = 1..NPRODFACT = 1..NFACTRAW = 1..NRAWTIME = 1..NT! Range of products (p)! factories (f)! raw materials (r)! time periods (t)nIter: integer
nPROP: array(FACT) of integer
end-declarations! Iteration counter! Counters of proposals from subprob.s!**** Master problem ****

declarations
MXSELL: array(PRODS,TIME) of realexcessS: mpvar
weight: array(FACT,range) of mpvar
MxSell: array(PRODS,TIME) of linctr
Convex: array(FACT) of linctr
Price_convex: array(FACT) of realPrice_sell: array(PRODS,TIME) of real
end-declarations! Max. amount of p that can be sold! Violation of sales/buying limits! weights for proposals! Sales limit constraints! Convexity constraints! Dual price on convexity constraints! Dual price on sales limits

initializations from DATAFILEMXSELL
end-initializations37

!**** Submodels ****

declarations
submod: array(FACT) of Model
Stopped: set of integer
modid: array(set of integer) of integer
end-declarations! One subproblem per factory! Model indicesres:= compile("g","cocoSubF.mos")! Compile the submodel file
forall(f in FACT) do! Load & run one submodel per product
Price_convex(f):= 1
load(submod(f), "cocoSubF.bim")
modid(getid(submod(f))):= frun(submod(f), "Factory=" + f + ",DATAFILE=" + DATAFILE)wait! Wait for child model to be readydropnextevent
end-do!**** Phase 0: Crash ****

nIter:=1; finished:=false
writeln("¥nPHASE 0 -- Iteration ", nIter); fflush
forall(f in FACT)
send(submod(f), PHASE_0, 0)! Start solving all submodels (Phase 1)

forall(f in FACT) do
wait
ev:= getnextevent
if getclass(ev)=EVENT_SOLVED thenprocess_sub_result
elif getclass(ev)=EVENT_FAILED thenfinished:= true
end-ifend-do! Wait for child (termination) events! Add new proposal to master problem

if finished thenwriteln("Problem is infeasible")exit(1)
end-ifsolve_master(1)process_master_result! Solve the updated Ph. 1 master problem! Store initial pricing data for submodels!**** Phase 1: proposal generation (feasibility) ****

repeatnoimprove:= 0
nIter+=1
writeln("¥nPHASE 1 -- Iteration ", nIter); fflush
forall(f in FACT)
send(submod(f), PHASE_1, Price_convex(f))
forall(f in FACT) do
wait
ev:= getnextevent
if getclass(ev)=EVENT_SOLVED thenprocess_sub_result
elif getclass(ev)=EVENT_FAILED thennoimprove += 1
end-ifend-do! Start solving all submodels (Phase 1)! Wait for child (termination) events! Add new proposal to master problem

if noimprove = NFACT then38

writeln("Problem is infeasible")exit(2)
end-ifif ALG=0 and noimprove > 0 thenwriteln("No improvement by some subproblem(s)")
break
end-ifsolve_master(1)if getobjval>0.00001 thenprocess_master_result
end-ifuntil getobjval <= 0.00001! Solve the updated Ph. 1 master problem! Store new pricing data for submodels!**** Phase 2: proposal generation (optimization) ****
writeln("¥n**** PHASE 2 ****")
finished:=false
repeatsolve_master(2)process_master_result! Solve Phase 2 master problem! Store new pricing data for submodels

nIter+=1
writeln("¥nPHASE 2 -- Iteration ", nIter); fflush
forall(f in FACT)
send(submod(f), PHASE_2, Price_convex(f))! Start solving all submodels (Phase 2)

forall(f in FACT) do
wait! Wait for child (termination) events
ev:= getnextevent
if getclass(ev)=EVENT_SOLVED thenprocess_sub_result ! Add new proposal to master problem
elif getclass(ev)=EVENT_FAILED thenif ALG=0 thenfinished:=true ! 1st submodel w/o prop. stops phase 2
else
Stopped += {modid(getfromid(ev))} ! Stop phase 2 only when no submodel! generates a new proposal
end-ifend-ifend-doif getsize(Stopped) = NFACT then finished:= true; end-if
until finishedsolve_master(2)! Re-solve master to integrate! proposal(s) from last ph. 2 iteration!**** Phase 3: solution to the original problem ****

writeln("¥n**** PHASE 3 ****")
forall(f in FACT) dosend(submod(f), PHASE_3, 0)
wait
dropnextevent
end-do! Stop all submodels

writeln("Total Profit=", calc_solution)
print_solution
end-model

このモデルの、最初のdeclarations blockにより、このmaster modelとchild (sub)modelとの間

39

でやりとりされるメッセージを識別するのに使うevent codeが定義されます。child modelでも、同一のdeclarationを行なう必要があります。

5.2.2 単一工場問題
Subproblemを持つモデル・ファイルcocoSubF.mos は、下記の内容を持っています。
model "Coco Subproblem (factory based decomp.)"uses "mmxprs", "mmjobs"
parameters
Factory = 0
TOL = 0.00001
DATAFILE = "coco3.dat"
end-parameters
forward procedure process_solution
declarations
PHASE_0=2
PHASE_1=3
PHASE_2=4
PHASE_3=5
EVENT_SOLVED=6
EVENT_FAILED=7
EVENT_READY=8
NPROD, NFACT, NRAW, NT: integer
end-declarations! Event codes sent to submodels! Event codes sent by submodelssend(EVENT_READY,0)! Model is ready (= running)

initializations from DATAFILE
NPROD NFACT NRAW NT
end-initializations
declarations
PRODS = 1..NPRODFACT = 1..NFACTRAW = 1..NRAWTIME = 1..NT! Range of products (p)!
factories (f)!
raw materials (r)!
time periods (t)REV: array(PRODS,TIME) of real
CMAKE: array(PRODS,FACT) of real

! Unit selling price of products! Unit cost to make product p! at factory f! Unit cost to buy raw materials! Requirement by unit of product p! for raw material r! Max. amount of p that can be sold! Max. amount factory f can make! over all products! Initial product stock levels! Initial raw material stock levels! Unit cost to store any product p! Unit cost to store any raw mat. r! Raw material storage capacity! Amount of products made at factory! Amount of product sold from factory! Amount of raw material bought! Product stock levels at start! of period t

CBUY: array(RAW,TIME) of real
REQ: array(PRODS,RAW) of real
MXSELL: array(PRODS,TIME) of real
MXMAKE: array(FACT) of real
IPSTOCK: array(PRODS,FACT) of real
IRSTOCK: array(RAW,FACT) of real
CPSTOCK: real
CRSTOCK: real
MXRSTOCK: real
make: array(PRODS,TIME) of mpvar
sell: array(PRODS,TIME) of mpvar
buy: array(RAW,TIME) of mpvar
pstock: array(PRODS,1..NT+1) of mpvar

40

rstock: array(RAW,1..NT+1) of mpvar! Raw material stock levels! at start of period t! Amount of products made! Amount of product sold! Amount of raw mat. bought! Product stock levels! Raw mat. stock levels! Profit of proposal! Dual price on sales limits

sol_make: array(PRODS,TIME) of real
sol_sell: array(PRODS,TIME) of real
sol_buy: array(RAW,TIME) of real
sol_pstock: array(PRODS,1..NT+1) of real
sol_rstock: array(RAW,1..NT+1) of real
Profit: linctr
Price_sell: array(PRODS,TIME) of real
end-declarations
initializations from DATAFILECMAKE REV CBUY REQ MXSELL MXMAKEIPSTOCK IRSTOCK MXRSTOCK CPSTOCK CRSTOCK
end-initializations! Product stock balance
forall(p in PRODS,t in TIME)
PBal(p,t):= pstock(p,t+1) = pstock(p,t) + make(p,t) - sell(p,t)! Raw material stock balance
forall(r in RAW,t in TIME)
RBal(r,t):= rstock(r,t+1) =rstock(r,t) + buy(r,t) - sum(p in PRODS) REQ(p,r)*make(p,t)! Capacity limit
forall(t in TIME)
MxMake(t):= sum(p in PRODS) make(p,t) <= MXMAKE(Factory)! Limit on the amount of prod. p to be soldforall(p in PRODS,t in TIME) sell(p,t) <= MXSELL(p,t)! Raw material stock limit
forall(t in 2..NT+1)
MxRStock(t):= sum(r in RAW) rstock(r,t) <= MXRSTOCK! Initial product and raw material stock levels
forall(p in PRODS) pstock(p,1) = IPSTOCK(p,Factory)
forall(r in RAW) rstock(r,1) = IRSTOCK(r,Factory)! Total profit
Profit:=sum(p in PRODS,t in TIME) REV(p,t) * sell(p,t) ! revenuesum(p in PRODS,t in TIME) CMAKE(p,Factory) * make(p,t) - ! prod. costsum(r in RAW,t in TIME) CBUY(r,t) * buy(r,t) ! raw mat.
sum(p in PRODS,t in 2..NT+1) CPSTOCK * pstock(p,t) ! p storagesum(r in RAW,t in 2..NT+1) CRSTOCK * rstock(r,t)! r storage! (Re)solve this model until it is stopped by event "PHASE_3"
repeatwait
ev:= getnextevent
Phase:= getclass(ev)if Phase=PHASE_3 thenbreak
end-ifPrice_convex:= getvalue(ev)if Phase<>PHASE_0 theninitializations from "raw:noindex"
Price_sell as "shmem:Price_sell"
end-initializations
end-if41! Stop the execution of this model! Get new pricing data

! (Re)solve this model
if Phase=PHASE_0 thenmaximize(Profit)
elif Phase=PHASE_1 thenmaximize(sum(p in PRODS,t in TIME) Price_sell(p,t)*sell(p,t) + Price_convex)
else! PHASE 2
maximize(
Profit - sum(p in PRODS,t in TIME) Price_sell(p,t)*sell(p,t) Price_convex)
end-ifwriteln("Factory ", Factory, " - Obj: ", getobjval," Profit: ", getsol(Profit), " Price_sell: ",getsol(sum(p in PRODS,t in TIME) Price_sell(p,t)*sell(p,t) )," Price_convex: ", Price_convex)
Fflushif getobjval > TOL thenprocess_solution
elif getobjval <= TOL thensend(EVENT_FAILED,0)
else
send(EVENT_READY,0)
end-ifuntil false! Solution found: send values to master! Problem is infeasible (Phase 0/1) or! no improved solution found (Phase 2)!----------------------------------------------------------! Process solution dataprocedure process_solution
forall(p in PRODS,t in TIME) do
sol_make(p,t):= getsol(make(p,t))sol_sell(p,t):= getsol(sell(p,t))
end-doforall(r in RAW,t in TIME) sol_buy(r,t):= getsol(buy(r,t))
forall(p in PRODS,t in 1..NT+1) sol_pstock(p,t):= getsol(pstock(p,t))
forall(r in RAW,t in 1..NT+1) sol_rstock(r,t):= getsol(rstock(r,t))
Prop_cost:= getsol(Profit)
send(EVENT_SOLVED,0)
initializations to "mempipe:noindex,sol"
Factorysol_make sol_sell sol_buy sol_pstock sol_rstock
Prop_cost
end-initializations
end-procedure
end-modelchild modelは、PHASE̲3 codeを得るまで、何度も、繰り返し解かれます。各々の繰り返しで、
child modelはメモリにソリューションを書き、master modelがソリューションを処理できるよ
うにします。

5.2.3 Master problem の subroutines
下記のmaster modelの三つのサブルーチンにより、master modelは、サブルーチン(process̲sub̲result)が得たソリューションを得て、master問題(solve̲master),を繰り返し解き、master modelのソリューションをchild model(process̲master̲result)に渡します。

42

declarations
Prop_make: array(PRODS,FACT,TIME,range) of real ! Amount of products made
Prop_sell: array(PRODS,FACT,TIME,range) of real ! Amount of product sold
Prop_buy: array(RAW,FACT,TIME,range) of real ! Amount of raw mat. bought
Prop_pstock: array(PRODS,FACT,1..NT+1,range) of real ! Product stock levels
Prop_rstock: array(RAW,FACT,1..NT+1,range) of real ! Raw mat. stock levels
Prop_cost: array(FACT,range) of real ! Cost/profit of each proposal
end-declarations
procedure process_sub_result
declarations
f: integer! Factory index! Solution values of the proposal:! Amount of products made! Amount of product sold! Amount of raw mat. bought! Product stock levels! Raw mat. stock levels! Cost of the proposal

sol_make: array(PRODS,TIME) of real
sol_sell: array(PRODS,TIME) of real
sol_buy: array(RAW,TIME) of real
sol_pstock: array(PRODS,1..NT+1) of real
sol_rstock: array(RAW,1..NT+1) of real
pc: real
end-declarations! Read proposal data from memory
initializations from "mempipe:noindex,sol"fsol_make sol_sell sol_buy sol_pstock sol_rstockpc
end-initializations! Add the new proposal to the master problem
nPROP(f)+=1
create(weight(f,nPROP(f)))
forall(p in PRODS,t in TIME) do
Prop_make(p,f,t,nPROP(f)):= sol_make(p,t)
Prop_sell(p,f,t,nPROP(f)):= sol_sell(p,t)
end-doforall(r in RAW,t in TIME) Prop_buy(r,f,t,nPROP(f)):= sol_buy(r,t)
forall(p in PRODS,t in 1..NT+1) Prop_pstock(p,f,t,nPROP(f)):= sol_pstock(p,t)
forall(r in RAW,t in 1..NT+1) Prop_rstock(r,f,t,nPROP(f)):= sol_rstock(r,t)
Prop_cost(f,nPROP(f)):= pc
writeln("Sol. for factory ", f, ":¥n make: ", sol_make, "¥n sell: ",sol_sell, "¥n buy: ", sol_buy, "¥n pstock: ", sol_pstock,"¥n rstock: ", sol_rstock)end-procedure!----------------------------------------------------------procedure solve_master(phase: integer)
forall(f in FACT)
Convex(f):= sum (k in 1..nPROP(f)) weight(f,k) = 1
if phase=1 thenforall(p in PRODS,t in TIME)
MxSell(p,t):=sum(f in FACT,k in 1..nPROP(f)) Prop_sell(p,f,t,k)*weight(f,k) excessS <= MXSELL(p,t)
minimize(excessS)
else
forall(p in PRODS,t in TIME)
MxSell(p,t):=sum(f in FACT,k in 1..nPROP(f)) Prop_sell(p,f,t,k)*weight(f,k) <=MXSELL(p,t)
maximize(sum(f in FACT, k in 1..nPROP(f)) Prop_cost(f,k) * weight(f,k))
end-if43

writeln("Master problem objective: ", getobjval)
write(" Weights:")
forall(f in FACT,k in 1..nPROP(f)) write(" ", getsol(weight(f,k)))writeln
end-procedure!----------------------------------------------------------procedure process_master_result
forall(p in PRODS,t in TIME) Price_sell(p,t):=getdual(MxSell(p,t))
forall(f in FACT) Price_convex(f):=getdual(Convex(f))
initializations to "raw:noindex"
Price_sell as "shmem:Price_sell"
end-initializations
end-procedure

最後に、master modelは、もともとの問題(calc̲solution)へのソリューションを計算する二つのサブルーチンにより完了し（decomposition algorithmのフェーズ3）、ソリューションがプリントされます(print̲solution)。もともとの問題のソリューションは、modified master problemの値とsubproblemで生成されるproposalとから得ます。
declarations
REV: array(PRODS,TIME) of real
CMAKE: array(PRODS,FACT) of real

! Unit selling price of products! Unit cost to make product p! at factory f! Unit cost to buy raw materials! Fixed cost of factory f being! open for one period! Unit cost to store any product p! Unit cost to store any raw mat. R

CBUY: array(RAW,TIME) of real
COPEN: array(FACT) of real
CPSTOCK: real
CRSTOCK: real

Sol_make: array(PRODS,FACT,TIME) of real ! Solution value (products made)Sol_sell: array(PRODS,FACT,TIME) of real ! Solution value (product sold)Sol_buy: array(RAW,FACT,TIME) of real ! Solution value (raw mat. bought)
Sol_pstock: array(PRODS,FACT,1..NT+1) of real ! Sol. value (prod. stock)
Sol_rstock: array(RAW,FACT,1..NT+1) of real ! Sol. value (raw mat. stock)end-declarations
initializations from DATAFILECMAKE REV CBUY CPSTOCK CRSTOCK COPEN
end-initializations
function calc_solution: real
forall(p in PRODS,f in FACT,t in TIME) do
Sol_sell(p,f,t):=sum(k in 1..nPROP(f)) Prop_sell(p,f,t,k) * getsol(weight(f,k))
Sol_make(p,f,t):=sum(k in 1..nPROP(f)) Prop_make(p,f,t,k) * getsol(weight(f,k))
end-doforall(r in RAW,f in FACT,t in TIME) Sol_buy(r,f,t):=sum(k in 1..nPROP(f)) Prop_buy(r,f,t,k) * getsol(weight(f,k))
forall(p in PRODS,f in FACT,t in 1..NT+1) Sol_pstock(p,f,t):=sum(k in 1..nPROP(f)) Prop_pstock(p,f,t,k) * getsol(weight(f,k))
forall(r in RAW,f in FACT,t in 1..NT+1) Sol_rstock(r,f,t):=sum(k in 1..nPROP(f)) Prop_rstock(r,f,t,k) * getsol(weight(f,k))
returned:=sum(p in PRODS,f in FACT,t in TIME) REV(p,t) * Sol_sell(p,f,t) -44

sum(p in PRODS,f in FACT,t in TIME) CMAKE(p,f) * Sol_make(p,f,t) sum(r in RAW,f in FACT,t in TIME) CBUY(r,t) * Sol_buy(r,f,t) sum(p in PRODS,f in FACT,t in 2..NT+1) CPSTOCK * Sol_pstock(p,f,t) sum(r in RAW,f in FACT,t in 2..NT+1) CRSTOCK * Sol_rstock(r,f,t)
end-function!----------------------------------------------------------procedure print_solution
writeln("Finished products:")
forall(f in FACT) do
writeln("Factory ", f, ":")
forall(p in PRODS) do
write(" ", p, ": ")forall(t in TIME) write(strfmt(Sol_make(p,f,t),6,1), "(",
strfmt(Sol_pstock(p,f,t+1),5,1), ")")
writeln
end-doend-dowriteln("Raw material:")forall(f in FACT) do
writeln("Factory ", f, ":")
forall(r in RAW) do
write(" ", r, ": ")forall(t in TIME) write(strfmt(Sol_buy(r,f,t),6,1), "(",
strfmt(Sol_rstock(r,f,t+1),5,1), ")")
writeln
end-doend-dowriteln("Sales:")
forall(f in FACT) do
writeln("Factory ", f, ":")
forall(p in PRODS) do
write(" ", p, ": ")forall(t in TIME) write(strfmt(Sol_sell(p,f,t),4))
writeln
end-doend-dowriteln("¥nComputation time: ", gettime)
end-procedure

5.3 結果ここでのテスト・ケースに関しては、decomposition algorithm を使って得られたソリューショ
ンは、最適解に近いものですが、最適解は、必ずしも得られるとは限りません。その理由は、
decomposition algorithm が繰り返しの連続で、これにより、いろいろなレベルでエラーが蓄積し、submodelでの停止判定基準（stopping criterion ）として使われるトレランスを低め、多くの場合、解を改善しなくなるからです。しかし、decomposition algorithm のコンフィギュレーションそれ自体により、解に影響することもできます。例えば､フェーズ1、フェーズ2で、最初のsubmodel が改善された解を返して寄越さなくなったときに停止する、つまり、新しい
proposalが生成されなくなるまで続けるという方針です。たくさんのproposalを生成することで、より良いソリューションが見つかるかもしれませんが、これにより、(sub)problemsを解く回数も増加し、したがって、解く時間が長くなります。

45

6 Benders decomposition: working with several different submodels
Benders decompositionは、大きなMIP問題を解くためのdecomposition methodの一つです。 標
準的な解法を使って解くには大きすぎるMIP問題を、そのまま解く代わりに、一連の線形、および、pure integer subproblemsを使いますが、このとき、後者の制約式の数は、もともとの問題のそれよりも少なくなります。下記のdecomposition algorithmについての説明は、[Hu69]から取ったものです。ここには、この方法の有限性（finiteness）、および、この方法が、常に、最適解をもたらすことの証明がありますので、興味のある読者はご覧ください。
下記のMIP問題(problem I)の標準形式について考えましょう。

minimize z = Cx + Dy
Ax + By ≥ bx, y ≥ 0, y integer上の定式化、および、本章を通じて、ベクトル、および、マトリックスには太字を使います。例
NCTVAR

えば､Cx +Dy は、

NINTVAR∑ C ·x + ∑ D ·y を意味し、ここで、NCTVARは連続変数の数を表し、ま
i =1

i

i

i =1

i

i

た、NINTVARはインテジャー変数の数を表します。yが特定の値yʼを取るとき、上の問題は線形問題になり、下記のように表現できます(problem II)。
minimize Cx
Ax ≥ b − Byʼx≥ 0problem IIの双対問題は、下記のproblem IIdのようになります。
maximize u(b − Byʼ)uA ≤ Cu ≥ 0
双対問題IId の興味深い特徴は、uの実行可能領域が、yに独立であるということです。さらに、双対定理から、もし、問題IId がインフィージブルであったり、有限最適解を持たなかったりす
るならば、もともとの問題problem Iも有限最適解を持たない、ということがいえます。また、双対定理から、問題IId が有限最適解u*を持つならば、このソリューションは、主問題(すなわち、u* (b − Byʼ) = Cx*)の最適解と同じ値を持ち、実行可能領域の頂点pでの解upは、up(b − Byʼ)≤ Cx*という値を持つということがいえます。もともとのMIP問題の目的関数に、実行可能領域の頂点pでの解up を代入すると、下記のような制約式を得ます。
z ≥ up (b − By) + Dy

46

もともとのMIP問題(problem I)の最適解z*を得るには、Benders decomposition algorithmとして知られている、下記のpartitioning algorithm を使うことができます。

Step 0Find a uʼ that satisfies uʼA ≤ C
if no such uʼ exists
then the original problem (I) has no feasible solution
else continue with Step 1
end-ifStep 1

Solve the pure integer programminimize zz ≥ uʼ(b − By) + Dyy ≥ 0, y integer
If z is unbounded from below, take a z to be any small value zʼ.

Step 2

With the solution yʼ of Step 1, solve the linear programmaximize u (b − Byʼ)uA ≤ Cu≥0
If u goes to inifinity with u(b − Byʼ)then add the constraintPi∑ u ≤ M, where M is a large positive constant, and resolveiithe problem
end-if
Let the solution of this program be uʼʼ.
if zʼ − Dyʼ ≤ uʼʼ(b − Byʼ)then continue with Step 3
else return to Step 1 and add the constraint zʼ ≥ Dy + uʼʼ (b − By)
end-ifStep 3

With the solution yʼ of Step 1, solve the linear programminimize Cx
Ax ≥ b − Byʼx ≥0xʼ and yʼ are the optimum solution z* = Cxʼ + Dyʼ

47

このアルゴリズムは、有限の解をもたらすことが証明できます（provably finite）。そして、最適解をもたらし、実行の中のどの時点でも、最適解z*の上限、下限がわかります。

6.1 小さい例題
ここで、Benders decompositionを使って、下記の小さい問題を解いてみましょう。この問題の連続変数xiの数はNCTVAR = 3、インテジャー変数yiの数はNINTVAR = 3です。また、不等式の制約式はNC = 4です。

maximize 5 · x1 + 4 · x2 + 3 · x3 + 2 · y1 + 2 · y2 + 3 · y3
x1 −+ 2 · y3 ≤ 5

x2 + 5 · x3 + 3 · y1

≤ 7−2·x + y −y+ y ≤ 4+ 5 · y + 5 · y + 5 · y ≤ −2

4 · x2 + 3 · x3 − y1
2 · x13 · x1

3

+

y2

1

2

1

2

3

3

x, y ≥ 0, y integer

6.2 インプリメンテーション
6.2.1 Master model
master modelはデータを読み、ソリューション・アルゴリズムを定義し、サブモデルとのやり取りを調整し、最後にソリューションをプリントします。 (fixed integer variables を使って双対問題を解く)アルゴリズムのステップ2では、主問題を解き、双対問題の解の値をOptimizer から得るか、もしくは、自分の手で双対問題を定義して、それを解くか、という選択ができます。
Model parameter ALGにより、ユーザは、どちらを採るかを選択します。
master modelの主要部分は、下記のようなものになります。ソリューション・アルゴリズムをス
タートさせる前に、サブモデルを、すべて、コンパイルし、ロードし、スタートします。それは、アルゴリズムのそれぞれのステップで、対応するサブモデルの(re)solvingをトリガーするだけで済むようにするためです。
model "Benders (master model)"uses "mmxprs", "mmjobs"
parameters
NCTVAR = 3
NINTVAR = 3
NC = 4
BIGM = 1000
ALG = 1

! 1: Use Benders decomposition (dual)! 2: Use Benders decomposition (primal)

DATAFILE = "bprob33.dat"
end-parameters
forward procedure start_solution
forward procedure solve_primal_int(ct: integer)
forward procedure solve_cont48

forward function eval_solution: boolean
forward procedure print_solution
declarations
STEP_0=2
STEP_1=3
STEP_2=4
EVENT_SOLVED=6
EVENT_INFEAS=7
EVENT_READY=8! Event codes sent to submodels! Event codes sent by submodels

CtVars = 1..NCTVAR
IntVars = 1..NINTVAR
Ctrs = 1..NC! Continuous variables! Discrete variables! Set of constraints (orig. problem)

A: array(Ctrs,CtVars) of integer
B: array(Ctrs,IntVars) of integerb: array(Ctrs) of integer
C: array(CtVars) of integer
D: array(IntVars) of integer
Ctr: array(Ctrs) of linctr
CtrD: array(CtVars) of linctr
MC: array(range) of linctr! Coeff.s of continuous variables! Coeff.s of discrete variables! RHS values! Obj. coeff.s of continuous variables! Obj. coeff.s of discrete variables! Constraints of orig. problem! Constraints of dual problem! Constraints generated by alg.

sol_u: array(Ctrs) of real
sol_x: array(CtVars) of real
sol_y: array(IntVars) of real
sol_obj: real
RM: range
stepmod: array(RM) of Model
end-declarations! Solution of dual problem! Solution of primal prob. (cont.)! Solution of primal prob. (integers)! Objective function value (primal)! Model indices! Submodels

initializations from DATAFILE
ABbCD
end-initializations! **** Submodels ****
initializations to "raw:noindex"! Save data for submodels
A as "shmem:A" B as "shmem:B"b as "shmem:b" C as "shmem:C" D as "shmem:D"
end-initializations! Compile + load all submodels
if compile("benders_int.mos")<>0 then exit(1); end-if
create(stepmod(1)); load(stepmod(1), "benders_int.bim")
if compile("benders_dual.mos")<>0 then exit(2); end-if
if ALG=1 thencreate(stepmod(2)); load(stepmod(2), "benders_dual.bim")elsecreate(stepmod(0)); load(stepmod(0), "benders_dual.bim")
if compile("benders_cont.mos")<>0 then exit(3); end-if
create(stepmod(2)); load(stepmod(2), "benders_cont.bim")run(stepmod(0), "NCTVAR=" + NCTVAR + ",NINTVAR=" + NINTVAR + ",NC=" + NC)
end-if! Start the execution of the submodels
run(stepmod(1), "NINTVAR=" + NINTVAR + ",NC=" + NC)run(stepmod(2), "NCTVAR=" + NCTVAR + ",NINTVAR=" + NINTVAR + ",NC=" + NC)
forall(m in RM) dowait! Wait for "Ready" messages

49ev:= getnextevent
if getclass(ev) <> EVENT_READY thenwriteln("Error occurred in a subproblem")exit(4)
end-ifend-do! **** Solution algorithm ****
start_solution
ct:= 1
repeatwriteln("¥n**** Iteration: ", ct)
solve_primal_int(ct)solve_cont
ct+=1
until eval_solution
print_solution! 0. Initial solution for cont. var.s! 1. Solve problem with fixed cont.
! 2. Solve problem with fixed int.
!
Test for optimality! 3. Retrieve and display the solution

異なるサブモデルを開始したサブルーチンは、'start solvingʼeventを送信し、サブモデルが解けたら、ソリューションをリトリーブします。start solution の生成には、パラメータALG のセッティングに従い、正しいサブモデルを選択する必要があります。もしこの問題がインフィー
ジブルであることが判明したら、問題全体がインフィージブルであるので、モデルの実行を停止します。
! Produce an initial solution for the dual problem using a dummy objective
procedure start_solution
if ALG=1 then
! Start the problem solving
send(stepmod(2), STEP_0, 0)
else
send(stepmod(0), STEP_0, 0)
end-ifwait! Wait for the solution
ev:=getnextevent
if getclass(ev)=EVENT_INFEAS thenwriteln("Problem is infeasible")exit(6)
end-ifend-procedure!----------------------------------------------------------! Solve a modified version of the primal problem, replacing continuous! variables by the solution of the dual
procedure solve_primal_int(ct: integer)
send(stepmod(1), STEP_1, ct)! Start the problem solving
wait! Wait for the solution
ev:=getnextevent
sol_obj:= getvalue(ev)! Store objective function value
initializations from "raw:noindex"
sol_y as "shmem:y"end-initializations
end-procedure! Retrieve the solution!----------------------------------------------------------! Solve the Step 2 problem (dual or primal depending on value of ALG)
! for given solution values of yprocedure solve_contsend(stepmod(2), STEP_2, 0)! Start the problem solving

50

wait
dropnextevent! Wait for the solution

initializations from "raw:noindex"
sol_u as "shmem:u"
end-initializations
end-procedure! Retrieve the solution

master modelは、termination criterion (function eval̲- solution)が満たされたかどうかをテストし、最終ソリューションをプリントします(procedure print̲solution)。この手順は、また、すべてのサブモデルを停止します。
function eval_solution: boolean
write("Test optimality: ", sol_obj - sum(i in IntVars) D(i)*sol_y(i)," = ", sum(j in Ctrs) sol_u(j)* (b(j) sum(i in IntVars) B(j,i)*sol_y(i)) )returned:= ( sol_obj - sum(i in IntVars) D(i)*sol_y(i) =sum(j in Ctrs) sol_u(j)* (b(j) - sum(i in IntVars) B(j,i)*sol_y(i)) )writeln(if(returned, " : true", " : false"))end-function!----------------------------------------------------------procedure print_solution! Retrieve results
initializations from "raw:noindex"
sol_x as "shmem:x"
end-initializations
forall(m in RM) stop(stepmod(m))! Stop all submodels

write("¥n**** Solution (Benders): ", sol_obj, "¥n x: ")forall(i in CtVars) write(sol_x(i), " ")write(" y: ")forall(i in IntVars) write(sol_y(i), " ")writeln
end-procedure

6.2.2 Submodel 1: fixed continuous variables
decomposition algorithmの最初のステップで、pure integer problemを解く必要があります。
このモデルの実行が開始されるとき、モデルはinvariant data を読み込み、変数をセットアッ
プします。次いで、wait statement (first line of the repeat-until loop) により、master model
が(solving) eventを送ってくるまで休止します。この問題を解く度に、メモリから読んでくる連続問題の、前回のランの解の値を使って、integer problemの新しい制約式MC(k)が定義されます。モデル全体、および､解法が埋め込まれているendless loop は、master modelからのʻstop
modelʼcommandによってのみ、終了されます。このサブモデル(file benders̲- int.mos)のソース全体を下記にリストします。
model "Benders (integer problem)"uses "mmxprs", "mmjobs"
parameters

51

NINTVAR = 3
NC = 4
BIGM = 1000
end-parameters
declarations
STEP_0=2
STEP_1=3
EVENT_SOLVED=6
EVENT_READY=8! Event codes sent to submodels! Event codes sent by submodels

IntVars = 1..NINTVAR
Ctrs = 1..NC! Discrete variables! Set of constraints (orig. problem)

B: array(Ctrs,IntVars) of integerb: array(Ctrs) of integer
D: array(IntVars) of integer
MC: array(range) of linctr! Coeff.s of discrete variables! RHS values! Obj. coeff.s of discrete variables! Constraints generated by alg.

sol_u: array(Ctrs) of real
sol_y: array(IntVars) of real

! Solution of dual problem! Solution of primal prob.

y: array(IntVars) of mpvar
z: mpvar
end-declarations! Discrete variables! Objective function variable

initializations from "raw:noindex"
B as "shmem:B" b as "shmem:b" D as "shmem:D"
end-initializations
z is_free
forall(i in IntVars) y(i) is_integer
forall(i in IntVars) y(i) <= BIGM

! Objective is a free variable! Integrality condition! Set (large) upper boundsend(EVENT_READY,0)! Model is ready (= running)

repeatwait
ev:= getnextevent
ct:= integer(getvalue(ev))
initializations from "raw:noindex"
sol_u as "shmem:u"
end-initializations! Add a new constraint
MC(ct):= z >= sum(i in IntVars) D(i)*y(i) +sum(j in Ctrs) sol_u(j)*(b(j) - sum(i in IntVars) B(j,i)*y(i))
minimize(z)! Store solution values of yforall(i in IntVars) sol_y(i):= getsol(y(i))
initializations to "raw:noindex"
sol_y as "shmem:y"end-initializations
send(EVENT_SOLVED, getobjval)
write("Step 1: ", getobjval, "¥n y: ")forall(i in IntVars) write(sol_y(i), " ")write("¥n Slack: ")52

forall(j in 1..ct) write(getslack(MC(j)), " ")writelnfflushuntil false
end-model

我われが解いているこの問題は、イタレーションの次の問題とたった一つの制約式が違うだけですから、root LP-緩和のソリューションのベイシス(MIP solutionのベイシスではない)をセーブして置き、次の最適化のランのときにリロードすることは行う価値があります。しかし、ここでの小さなテスト用のモデルでは、実行スピードの改善にあまり役立ちません。ベイシスをセーブし、読み込むには、call to minimizeを、下記の一連のステートメントで置き換えます。
declarations
bas: basis
end-declarations
loadprob(z)loadbasis(bas)minimize(XPRS_TOP, z)
savebasis(bas)minimize(XPRS_GLB, z)

6.2.3 Submodel 2: fixed integer variables
decomposition algorithm の二番目のステップは、すべてのインテジャー変数を、一番目のステップで見つかったソリューションの値に固定されたサブモデルを解くことです。このステップを実行するモデルの構造は、前のサブモデルと非常によく似ています。モデルがランされると、モデルは、メモリからinvariant dataを読み、目的関数をセットアップします。次いで、モデルは、
ループの最初で、line waitで休止し、master modelから、step 2 solving event が送られてくるのを待ちます。解の反復のたびごとに、メモリから読み込まれる係数を使って、constraints
CTR は再定義され、ソリューションは、メモリ返されます。下記は、このモデルのソースです(filebenders̲cont.mos)。
model "Benders (continuous problem)"uses "mmxprs", "mmjobs"
parameters
NCTVAR = 3
NINTVAR = 3
NC = 4
BIGM = 1000
end-parameters
declarations
STEP_0=2
STEP_2=4
STEP_3=5
EVENT_SOLVED=6
EVENT_READY=8CtVars = 1..NCTVAR! Event codes sent to submodels! Event codes sent by submodels! Continuous variables

53

IntVars = 1..NINTVAR
Ctrs = 1..NC! Discrete variables! Set of constraints (orig. problem)

A: array(Ctrs,CtVars) of integer
B: array(Ctrs,IntVars) of integerb: array(Ctrs) of integer
C: array(CtVars) of integer
Ctr: array(Ctrs) of linctr! Coeff.s of continuous variables! Coeff.s of discrete variables! RHS values! Obj. coeff.s of continuous variables! Constraints of orig. problem

sol_u: array(Ctrs) of real
sol_x: array(CtVars) of real
sol_y: array(IntVars) of real

! Solution of dual problem! Solution of primal prob. (cont.)! Solution of primal prob. (integers)

x: array(CtVars) of mpvar
end-declarations! Continuous variables

initializations from "raw:noindex"
A as "shmem:A" B as "shmem:B"b as "shmem:b" C as "shmem:C"
end-initializations
Obj:= sum(i in CtVars) C(i)*x(i)send(EVENT_READY,0)! Model is ready (= running)! (Re)solve this model until it is stopped by event "STEP_3"
repeatwait
dropnextevent
initializations from "raw:noindex"
sol_y as "shmem:y"end-initializations
forall(j in Ctrs)
Ctr(j):= sum(i in CtVars) A(j,i)*x(i) +sum(i in IntVars) B(j,i)*sol_y(i) >= b(j)minimize(Obj)! Solve the problem! Store values of u and x
forall（j in Ctrs) sol_u(j):= getdual(Ctr(j))
forall(i in CtVars) sol_x(i):= getsol(x(i))
initializations to "raw:noindex"
sol_u as "shmem:u" sol_x as "shmem:x"
end-initializations
send(EVENT_SOLVED, getobjval)
write("Step 2: ", getobjval, "¥n u: ")forall(j in Ctrs) write(sol_u(j), " ")write("¥n x: ")forall(i in CtVars) write(getsol(x(i)), " ")writelnfflushuntil false
end-model

546.2.4 Submodel 0: start solution
decomposition algorithm を開始するには、連続変数に、初期値を生成する必要があります。これは、ダミー目的関数を使い、双対問題を解くことで行ないます。このアルゴリズムのステップ
2でも、双対問題を使いますが、これは、前の章で見た主問題を置き換えるときです。このサブモデルのインプリメンテーションでは、二つのケースを考えます。solving loop の中で、mastermodelにより送られたtype (class) of eventをテストし、それによって解くべき問題を選択します。
このモデルの主要部分は、下記のMosel code (file benders̲- dual.mos)で実行されます。
model "Benders (dual problem)"uses "mmxprs", "mmjobs"
parameters
NCTVAR = 3
NINTVAR = 3
NC = 4
BIGM = 1000
end-parameters
forward procedure save_solution
declarations
STEP_0=2
STEP_2=4
EVENT_SOLVED=6
EVENT_INFEAS=7
EVENT_READY=8! Event codes sent to submodels! Event codes sent by submodels

CtVars = 1..NCTVAR
IntVars = 1..NINTVAR
Ctrs = 1..NC! Continuous variables! Discrete variables! Set of constraints (orig. problem)

A: array(Ctrs,CtVars) of integer
B: array(Ctrs,IntVars) of integerb: array(Ctrs) of integer
C: array(CtVars) of integer! Coeff.s of continuous variables! Coeff.s of discrete variables! RHS values! Obj. coeff.s of continuous variables

sol_u: array(Ctrs) of real
sol_x: array(CtVars) of real
sol_y: array(IntVars) of real

! Solution of dual problem! Solution of primal prob. (cont.)! Solution of primal prob. (integers)

u: array(Ctrs) of mpvar
end-declarations! Dual variables

initializations from "raw:noindex"
A as "shmem:A" B as "shmem:B"b as "shmem:b" C as "shmem:C"
end-initializations
forall(i in CtVars) CtrD(i):= sum(j in Ctrs) u(j)*A(j,i) <= C(i)send(EVENT_READY,0)! Model is ready (= running)! (Re)solve this model until it is stopped by event "STEP_3"
repeatwait

55

ev:= getnextevent
Alg:= getclass(ev)if Alg=STEP_0 then

! Produce an initial solution for the! dual problem using a dummy objective

maximize(sum(j in Ctrs) u(j))
if(getprobstat = XPRS_INF) thenwriteln("Problem is infeasible")send(EVENT_INFEAS,0)
else
write("**** Start solution: ")save_solution
end-if! Problem is infeasible

else! STEP 2: Solve the dual problem for! given solution values of yinitializations from "raw:noindex"
sol_y as "shmem:y"end-initializations
Obj:= sum(j in Ctrs) u(j)* (b(j) - sum(i in IntVars) B(j,i)*sol_y(i))
maximize(XPRS_PRI, Obj)if(getprobstat=XPRS_UNB) then
BigM:= sum(j in Ctrs) u(j) <= BIGM
maximize(XPRS_PRI, Obj)
end-ifwrite("Step 2: ")save_solution
BigM:= 0
end-if! Write solution to memory! Reset the ’BigM’ constraint

until false

このモデルは、subroutine save̲solutionで終わります。subroutine save̲solutionは、メモリにソリューションを書き、EVENT̲- SOLVED message を送ることで、master modelにソリューシ
ョンが利用可能であることを知らせます。
! Process solution dataprocedure save_solution! Store values of u and x
forall(j in Ctrs) sol_u(j):= getsol(u(j))
forall(i in CtVars) sol_x(i):= getdual(CtrD(i))
initializations to "raw:noindex"
sol_u as "shmem:u" sol_x as "shmem:x"
end-initializations
send(EVENT_SOLVED, getobjval)
write(getobjval, "¥n u: ")forall(j in Ctrs) write(sol_u(j), " ")write("¥n x: ")forall(i in CtVars) write(getdual(CtrD(i)), " ")writelnfflush
end-procedure
end-model

566.3 Results
この小さなテストモデルの最適解では、目的関数の値は18.1852です。ここで使ったプログラム
により、下記のアウトプットが得られ、そこから、この問題は、decomposition algorithmを3回繰り返して得られたことが判ります(looping around steps 1 and 2)。
**** Start solution: 4.05556
u: 0.740741 1.18519
x: 0.611111 0.166667

2.12963 00.111111**** Iteration: 1
Step 1: -1146.15
y: 1000 0 0
Slack: 0
Step 2: 1007
u: 0 1 0 0
x: 0 251.75 0
Test optimality: -3146.15 = 1007 : false**** Iteration: 2
Step 1: 17.0185
y: 3 0 0
Slack: 0 -1.01852
Step 2: 12.5
u: 0 1 2.5 0
x: 0.5 2.5 0
Test optimality: 11.0185 = 12.5 : false**** Iteration: 3
Step 1: 18.1852
y: 2 0 0
Slack: 0 -5.18519 -0.185185
Step 2: 14.1852
u: 0.740741 1.18519 2.12963 0
x: 1.03704 2.22222 0.037037
Test optimality: 14.1852 = 14.1852 : true**** Solution (Benders): 18.1852
x: 1.03704
2.222220.037037

y: 2 0 0

57

７要約
このホワイトペーパーの中で示されている例は、mmjobs モジュールで用意した機能をどのように使うかを示しております。これらの例は、読者の皆さんが、技術的な可能性を網羅的に概観することなしに、問題をパラレルに解いたり、その他の複数問題解法アプローチ
を実際に使ったり、試してみたりするためのスターティングポイントとするためのものです。
Mosel で使えるどの solver モジュールも mmjobs.と連動して使えます。ただし、同一のsolver で複数問題をパラレルに解くことは、使おうとする solver が multi-threaded 環境
で動く場合においてのみ可能です。Xpress-Optimizer、それから派生した QP、SLP、SP、および、Xpress-Kalis は、この条件を満たしています。

58

Bibliography[Hu69] T. C. Hu. Integer programming and network flows. Addison-Wesley, Reading, MA, 1969.
[PW94] Y. Pochet and L. A. Wolsey. Algorithms and Reformulations for Lot Sizing Problems. Technical Report Discussion
Paper 9427, CORE, Université Catholique de Louvain, 1994.
[Teb01] J. R. Tebboth. A computational study of Dantzig-Wolfe decomposition. PhD thesis, University of Buckingham,
2001.

59