Decision Optimizationでスケジューリング(その2)

はじめに

　IBMは、CPLEXという最適化ソフトを製品として持っていますが、現在、そのクラウド対応版としてDecision Optimizationというツールを提供しています。具体的にはクラウド上のWatson Studioでも動きますし、OpenShiftベースの製品であるCP4Dでも動きます。
　下記リンク先のチュートリアルを参考に、Python APIを使ったスケジューリング問題の実装を解説していきます。

　元のチュートリアルは、かなり長いものなので、その中の問題毎に記事を分けて解説することにします。
当記事は、その中の第2弾です。
なお、一つ一つの関数の細かい解説は、下記APIリファレンスに詳しいので、こちらも適宜参照して下さい。

事前準備

　当記事の実習を動かすための事前準備に関しては、別記事Decision Optimizationでスケジューリング(その1)に説明がありますので、そちらを参照して下さい。

Notebookの読み込み

　当記事のNotebookは、Watson StudioのJupyter NotebookからURL指定で読み込む場合、次のURLを指定します。
(読み込み手順については上の記事を参照して下さい)
https://github.com/makaishi2/cplex-samples/raw/master/house-building/house-building3.ipynb

また、ブラウザからコードを確認する場合、以下のURLを使って下さい。
https://github.com/makaishi2/cplex-samples/blob/master/house-building/house-building3.ipynb

以下の説明は、下の画面のように、すでにNotebookにコードが読み込み済みであることを前提に行います。

問題設定

コードの解説に入る前に、例によって問題設定を説明します。
今回は、次のような最適化問題を解きます。

• 5件の家を、JoeとJimの2人で建てる
  
• 工程毎に担当者が決まっている
• 家ごとにReleaseDate以降に着手する必要がある
• 家ごとにDueDateが決まっている。遅れてもいいが、Weightで規定されるコストが余分にかかる
• 評価関数は「遅延によるコスト」+「家ごとの建築日数」でこれを最小化する

コード解説

それでは、最適化コードの解説を始めます。解説はセル単位で行います。

条件設定

最初のコードは、変数・配列による条件設定です。

# 工程名
TaskNames = ["masonry", "carpentry", "plumbing", 
             "ceiling", "roofing", "painting", 
             "windows", "facade", "garden", "moving"]

# 工程毎の期間
Duration =  [35, 15, 40, 15, 5, 10, 5, 10, 5, 5]

# 工程間依存関係
Precedences = [("masonry", "carpentry"),("masonry", "plumbing"),
               ("masonry", "ceiling"), ("carpentry", "roofing"),
               ("ceiling", "painting"), ("roofing", "windows"),  
               ("roofing", "facade"), ("plumbing", "facade"),
               ("roofing", "garden"), ("plumbing", "garden"),
               ("windows", "moving"), ("facade", "moving"),  
               ("garden", "moving"), ("painting", "moving")]

# 家の件数
NbHouses = 5

# 作業者
WorkerNames = ["Joe", "Jim"]

# 工程毎の作業者分担
Worker = {"masonry"  : "Joe" , 
          "carpentry": "Joe" , 
          "plumbing" : "Jim" , 
          "ceiling"  : "Jim" , 
          "roofing"  : "Joe" , 
          "painting" : "Jim" , 
          "windows"  : "Jim" , 
          "facade"   : "Joe" , 
          "garden"   : "Joe" , 
          "moving"   : "Jim"}

# 着手日 (必ず守る必要あり)
ReleaseDate = [  0,     0,   151,    59,   243]

# 終了日 (遅れていいがコストが発生)
DueDate     = [120,   212,   304,   181,   425]

# 遅延コスト
Weight      = [100.0, 100.0, 100.0, 200.0, 100.0]

# 家IDのリスト([0, 1, 2, 3, 4])を生成
Houses = range(NbHouses)

工程名TaskNames、工程ごとの作業日数Durationと工程間の依存関係Precedencesは前の例題と同じものです。

これ以降の条件(変数)はすべて今回新しく出てきたものになります。

建てる家の件数 NbHouses
作業者リスト WorkerNames
工程ごとの担当者 Workers
着手日 ReleaseDate
終了日 DueDate
遅延コスト Weight
house_idのリスト Houses

モデルインスタンスの生成

次のセルはモデルインスタンスの生成です。こちらは、前のときと同じなので、説明は省略します。

# ライブラリのインポートとモデルインスタンスの生成

import sys
from docplex.cp.model import *
mdl2 = CpoModel()

家単位のInterval変数

5件の家それぞれに対して工事開始日、工事終了日をInterval変数で表現します。

# housesは家毎のInterval変数の配列
# 開始日はReleaseDateにより個別に指定

houses = [mdl2.interval_var(start=(ReleaseDate[i], INTERVAL_MAX), name="house"+str(i)) for i in Houses]

注意すべき点は、各インスタンスのstartオプション指定の方法です。ReleaseDate以前に着手してはいけないので、ReleaseDate[i]を開始日の先頭として指定しています。

タスク単位のInterval変数

次にタスク単位のInterval変数を定義します。

# <house_id>_<タスク名>でタスクごとのラベルを振る
# TaskNames_ids: タスク毎のラベルからタスクIDを取得する辞書
# itvs: (house_id, task_id)から該当するInterval変数を取得する辞書
# タスク毎の期間(size)はDurationで決められている値を初期設定する

TaskNames_ids = {}
itvs = {}
for h in Houses:
    for i,t in enumerate(TaskNames):
        _name = str(h)+"_"+str(t)
        itvs[(h,t)] = mdl2.interval_var(size=Duration[i], name=_name)
        TaskNames_ids[_name] = i

タスクは、複数の家を建てる今回のケースでは、(対象の家, 対象タスク)の2次元のインデックスを必要とします。
そこで、Interval変数の名称を _の形式で表現することします。
また、Interval変数取得用の辞書itvsを用意して、定義したインタバル変数のインスタンスを順次代入していきます。

次のセルでは、今定義した辞書の中身を確認しています。

# itvsの辞書の内容確認

print(itvs[(3, 'masonry')])

次のような結果になるはずです。

"3_masonry" = intervalVar(size=35)

タスク間の依存関係定義

次にタスク間の依存関係を定義していきます。

# タスク間の依存関係定義

for h in Houses:
    for p in Precedences:
        mdl2.add(mdl2.end_before_start(itvs[(h,p[0])], itvs[(h,p[1])]) )

ここは、「その1」で説明した実装とほぼ同じなので、説明は省略します。

span制約の設定

特定の家の建築期間とは、その家に関するすべてのタスクの開始日の最小値から、すべてのタスクの終了日の最大値までとなります。このような制約はCplexではspan関数で簡潔に表現できます。具体的な実装は下記のとおりです。

# span制約の設定

# 家のInterval変数は、個別タスクの集計で定まる。その関係をspan制約で設定
for h in Houses:
    mdl2.add( mdl2.span(houses[h], [itvs[(h,t)] for t in TaskNames] ) )

workersの定義

ここで、作業者ごとのタスク一覧をworkersという変数で定義します。定義は、Interval変数のリストを引数とするsequence_var関数で行います。結果はsequence変数になります。
sequence変数を次に説明するno_overlap関数にかけると、「重なりのないタスクリスト」という制約を表現できます。同一作業者が、同時に2つのタスクを行うことはできないので、この「sequence変数定義」->「no_overlap関数での制約表現」は、このユースケースの制約実装に便利なことがわかります。
具体的な実装は次のとおりです。

# workersの定義

# 作業者(Joe, Jim)ごとのタスクシーケンス変数として定義
# シーケンス変数は、次のno_overlap関数で条件の定義として用いられる
workers = {w : mdl2.sequence_var([ itvs[(h,t)] for h in Houses for t in TaskNames if Worker[t]==w ], 
        name="workers_"+w) for w in WorkerNames}

今、定義したworkersの中身を確認してみましょう。

# 結果の確認('Jim'に対する結果)

print(workers['Jim'])

次のような結果が返ってくるはずです。

workers_Jim = sequenceVar(["0_plumbing", "0_ceiling", "0_painting", "0_windows", "0_moving", "1_plumbing", "1_ceiling", "1_painting", "1_windows", "1_moving", "2_plumbing", "2_ceiling", "2_painting", "2_windows", "2_moving", "3_plumbing", "3_ceiling", "3_painting", "3_windows", "3_moving", "4_plumbing", "4_ceiling", "4_painting", "4_windows", "4_moving"])

no_overlap関数による制約表現

続いては、no_overlap関数による制約表現です。上のシーケンス変数定義とこの関数呼び出しのセットで、「同時に1つのことしかできない」という条件を表現できます。具体的な実装は次のようになります。

# no_overlap関数による制約表現

# 作業者は同時に一つのタスクしかできないことをno_overlap関数で表現する
# (オリジナルのチュートリアルにあった遷移時間の条件は簡略化のため削除)

for w in WorkerNames:
    mdl2.add( mdl2.no_overlap(workers[w]) )

なお、オリジナルのチュートリアルでは、この制約の定義の中にtransisionという概念を付け加えています。タスク切り替えに際して他の家に取り組む場合は、余分な時間(transision_time)がかかるという概念です。ここでは、説明を簡略化するため、その機能の実装を落としています。こういう高度な使い方もできることだけ、頭に入れておいて下さい。

目的関数の定義

これで、制約条件はすべて表現できました。最後に目的関数の定義を行います。
目的関数は、基本的に「各家の建築日数の最小化」なのですが、冒頭で説明したように、家ごとに、完成目標日があり、この日をすぎると1日ごとに余分なペナルティーがかかることにします。そのように目的関数を定義した後、その最小化を目的とすることを条件とします。
具体的なコード実装は以下のとおりです。

# 目的関数の定義
# 「遅延コスト」+「建築期間」を家ごとに集計した値を最小化する

mdl2.add( 
    mdl2.minimize( 
        mdl2.sum(Weight[h] * mdl2.max([0, mdl2.end_of(houses[h])-DueDate[h]]) +\
                 mdl2.length_of(houses[h]) for h in Houses) 
    ) 
)

モデルを解く

これでモデルに対する定義はすべて終わりました。saigonisolve関数を呼び出してモデルを解いてみます。
実装コードとその結果は、次のようになります。

# モデルを解く

print("\nSolving model....")
msol2 = mdl2.solve(FailLimit=30000)
print("done")

Solving model....
done

結果確認(目的関数値)

無事に問題が解けたので、結果に関して確認していきます。
最初に確認するのは、目的関数値です。get_object_value関数で値の取得ができます。
実装は次の形になります

# 目的関数値の確認

print("Cost will be ", msol2.get_objective_values()[0])

こんな結果になるはずです。

Cost will be  9820

結果確認(タスクごとの詳細スケジュール)

次にタスクごとの詳細スケジュールを確認してみましょう。
作業者ごとのスケジュールを調べる方法を取ることにします。
制約の設定の際に、作業者ごとのタスク一覧をシーケンス変数として定義しました。
まず、その値を取得し、次にその中のインタバル変数をループを回して確認することにします。
具体的な実装は次のとおりです。

# Workerごとの詳細スケジュール確認

for worker in WorkerNames:
    # シーケンス変数取得のためのキー値生成
    w_name = 'workers_' + worker
    # キー値からシーケンス変数を取得
    sequence = msol2.get_var_solution(w_name)
    ＃シーケンス変数から、インタバル変数のリスト取得
    intervals = sequence.get_interval_variables()
    for interval in intervals:
        # 個別のインタバル変数から、名前、開始日、終了日を取得
        name = interval.get_name()
        start = interval.get_start()
        end = interval.get_end()
        # 結果表示
        print(worker, name, start, end)

次のような結果が返ってくるはずです。

Joe 0_masonry 0 35
Joe 0_carpentry 35 50
Joe 0_roofing 50 55
Joe 3_masonry 60 95
Joe 0_facade 95 105
Joe 0_garden 105 110
Joe 3_carpentry 110 125
Joe 3_roofing 125 130
Joe 1_masonry 130 165
Joe 3_facade 165 175
Joe 3_garden 175 180
Joe 1_carpentry 180 195
Joe 1_roofing 195 200
Joe 2_masonry 205 240
Joe 1_garden 240 245
Joe 1_facade 245 255
Joe 2_carpentry 255 270
Joe 2_roofing 270 275
Joe 2_facade 300 310
Joe 2_garden 310 315
Joe 4_masonry 315 350
Joe 4_carpentry 350 365
Joe 4_roofing 365 370
Joe 4_facade 390 400
Joe 4_garden 400 405
Jim 0_plumbing 35 75
Jim 0_windows 75 80
Jim 0_ceiling 80 95
Jim 0_painting 95 105
Jim 3_ceiling 105 120
Jim 0_moving 120 125
Jim 3_plumbing 125 165
Jim 3_painting 165 175
Jim 3_windows 175 180
Jim 3_moving 180 185
Jim 1_ceiling 185 200
Jim 1_plumbing 200 240
Jim 1_painting 240 250
Jim 1_windows 250 255
Jim 1_moving 255 260
Jim 2_plumbing 260 300
Jim 2_ceiling 300 315
Jim 2_painting 315 325
Jim 2_windows 325 330
Jim 2_moving 330 335
Jim 4_plumbing 350 390
Jim 4_ceiling 390 405
Jim 4_painting 405 415
Jim 4_windows 415 420
Jim 4_moving 420 425

Joe、Jimの2人の作業者は、いつからいつまでどの家のどのタスクを作業しているかが確認できました。
最後に、この結果をガントチャートでグラフ化してみることにします。

グラフ表示

事前準備

グラフ表示用のライブラリロードを事前準備として行います。
実装は、次のとおりです。

# グラフ表示

# グラフ表示の準備
import docplex.cp.utils_visu as visu
import matplotlib.pyplot as plt
%matplotlib inline
from pylab import rcParams

家ごとのスケジュール確認

最初に家単位でみたときの建築スケジュールを確認してみます。
制約を定義するときに、house0-house4という変数名で、家の建築期間を示すインタバル変数を定義したので、この値を使えばいいはずです。
具体的な実装は次のようになります。

# 家単位の建築スケジュール

rcParams['figure.figsize'] = 15, 3

for h in Houses:
    hkey = 'house' + str(h)
    wt = msol2.get_var_solution(hkey)
    visu.interval(wt, 'lightblue', hkey)
visu.show()

下のような結果が表示されるはずです。

作業者ごとのスケジュール

次に作業者単位でみた、建築スケジュールを確認してみましょう。
詳細スケジュール確認の時と同様に workers_<Name>という名前で結果のシーケンス変数を取得します。
APIのグラフ表示機能には、シーケンス変数を渡すだけで、丸ごとガントチャートを表示してくれる関数があるので、その関数を利用することにします。
具体的な実装は次のとおりです。

# 作業者単位の建築スケジュール
# ダブルクリックで拡大表示されます

rcParams['figure.figsize'] = 80, 3

for worker in WorkerNames:
    name = 'workers_' + worker
    seq = msol2.get_var_solution(name)
    visu.sequence(name=name, intervals=seq)
visu.show()

結果は、こんな感じになるはずです。

一瞬「小さすぎて字が読めない」と思いますが、心配はいりません。グラフ領域をダブルクリックすると、下のように拡大表示されます。

まとめ

いかがだったでしょうか？今回は前回と違い、より現実に近い問題の定義と、その解決がDecision Optimizationでできたことがわかるかと思います。
また、その際には、Decision Optimization固有の「シーケンス変数」や「no_overlap関数」が、制約定義に重要な役割を果たしたこともわかると思います。
もちろん、これらの機能を使わなくてもこの問題を解くことはできるのですが、この機能を使うことで、制約条件定義の可読性が向上し、また処理速度も速くなることになります。
次回はより複雑な問題にチャレンジしてみます。