GoogleのOR-Toolsで組合せ最適化問題を解く(5)デートコースを決めよう

はじめに

本記事は、数理最適化についての参考テキスト「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」の練習問題を解く第五回目の記事です。

第一回目の記事は以下です。まずはこちらをご覧下さい。

GoogleのOR-Toolsで数理最適化モデルの練習問題を解く(1)一番易しいマス埋め問題
https://qiita.com/ttlabo/private/7e8c93f387814f99931f

デートコースを決めよう

参考テキストに出てくる応用問題です。
以下の問題をさっそく行ってみましょう。

問題

8.4.問題

8つのアトラクションがある遊園地でデートをする(8-4図)。
200分の制限時間の中で総満足度を最大化する。
なお、デートコースは全て回る必要はない。
条件を満たすデートコースを最適化せよ。

考え方

■モデル化の方針
　200分の制限時間を考慮し、全てを回らずとも可能な範囲で満足度を最大化することを考えます。そのために、次のような変数を用意します。

変数定義

\begin{eqnarray*}

S_i\;&\small{:=}&\;\small{i番目のアトラクションを選ぶかどうか}\\
M_{ij}\;&\small{:=}&\;\small{i番目のアトラクションからj番目のアトラクションへ移動するかどうか}\\
satisfy\;&\small{:=}&\;\small{最終的な満足度(の合計)}\\
stayTime\;&\small{:=}&\;\small{最終的な滞在時間(の合計)}\\
moveTime\;&\small{:=}&\;\small{最終的な移動時間(の合計)}\\
totalTime\;&\small{:=}&\;\small{全体時間(最終的な満足度と滞在時間の和)}

\end{eqnarray*}

最適化を解くとこれらの変数の値のうち最適な値が得られ、デートコースが決まります。
また、予めわかっている値や情報を次のように定数として定義とします。

定数定義

\begin{eqnarray*}

満足度定数　satisfyDef\\
滞在時間定数　stayTimeDef\\
移動時間定数　moveTimeDef

\end{eqnarray*}

順序が逆になりますが、まず定数について次のように値を設定します。
①満足度定数について、以下のように値を設定します。

8.4.renshu.py

satisfyDef = [0,50,36,45,79,55,63,71,42]

satisfyDefは配列で、i番目のアトラクションの満足度をsatisfyDef[i]とします。
例えば、2番目のアトラクションの満足度は36です。アトラクションのインデックス番号は0から数えることとします。
②滞在時間定数について、以下のように値を設定します。

8.4.renshu.py

stayTimeDef = [0,20,28,15,35,17,18,14,22]

stayTimeDefは配列で、i番目のアトラクションの滞在時間をstayTimeDef[i]とします。
③移動時間定数について、以下のように値を設定します。

8.4.renshu.py

moveTimeDef = {}
moveTimeDef[0,3] = 7
～

moveTimeDefはDictinary型で、iからjへ移動するときの時間をmoveTimeDef[i,j]とします。
例えば、アトラクション0番目(入口)からアトラクション3番目(コーヒーカップ)へ移動するときの移動時間はmoveTimeDef[0,3]で、7分という意味です。

①満足度定数、②滞在時間定数、③移動時間定数の設定値については、参考テキストで定義されている値です。

次に、変数Si,Mijについて次のように定義します。
①Siについて、以下のように定義します。

8.4.renshu.py

# アトラクションを選ぶかどうかの変数
# 1:選ぶ / 0:選ばない
S = {}
for i in range(9):
    S[i] = solver.BoolVar("S%i" % i)

ここでは、アトラクション0番目から8番目について、「1:選ぶ」か「0:選ばない」のどちらかの変数(BoolVar型)で定義しています。例えば、S[4]=1だった場合、4番目のアトラクションを選ぶ、という意味になります。

②Mijについて、以下のように定義します。

8.4.renshu.py

# 移動に関して、i→jに移動するかどうかの変数
# 1:移動する / 0:移動しない
M = {}
for i in range(9):
    for j in range(9):
        M[i,j] = solver.BoolVar("M%i%i" % (i,j))

ここでは、アトラクション0番目から8番目のそれぞれの組み合わせについて、「1:移動する」か「0:移動しない」のどちらかの変数(BoolVar型)で定義しています。9通り×9通り、全部で81通りの組み合わせを持ちます。例えばM[1,3]=1だった場合、1番目のアトラクションから3番目のアトラクションへ移動する、という意味になります。

以上の変数及び定数を用いて、変数や定数間に様々な条件(制約条件)を設定しながら最適化問題を解きます。

解答

プログラムを検討します。
プログラムの記載内容については、基本的にGoogleのOR-Toolsのガイドに沿っています。(https://developers.google.com/optimization)

プログラムの冒頭におまじないを書きます。

8.4.renshu.py

from __future__ import print_function
from ortools.linear_solver import pywraplp

混合整数計画ソルバーで解きますので、以下に宣言します。

8.4.renshu.py

# Create the mip solver with the CBC backend.
solver = pywraplp.Solver('simple_mip_program',
    pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

次に、定数について定義します。

8.4.renshu.py

# 定数定義 =================================================================

# アトラクション名
attrDef = ['','入口','喫茶','ボート','カップ','レストラン','観覧車',
           'お化け屋敷','コースター','迷路']

# 満足度定数
satisfyDef = [0,50,36,45,79,55,63,71,42]

# 滞在時間定数
stayTimeDef = [0,20,28,15,35,17,18,14,22]

上記、最初の定義attrDefは本プログラムでは使用しません。今回はとりあえず定義しました。
定数定義の説明で書きましたが、satisfayDefとstayTimeDefは配列で、インデックス番号とアトラクション番号をひもづけています。インデックス番号0はアトラクション番号0番目(入口)に対応し、満足度は0、滞在時間も0です。
移動時間については以下のように定義します。

8.4.renshu.py

# 移動時間定数
# moveTimeDef[i,j] ← iからjへ移動するときの時間
moveTimeDef = {}
moveTimeDef[0,0] = 0
moveTimeDef[0,1] = 9
moveTimeDef[0,2] = 0
moveTimeDef[0,3] = 7
〜

ここでは、移動できない経路の移動時間は0にしています。
上記では4経路のみ記載しましたが、後方に掲載するソース全体では全経路分を定義しています。
次に変数を定義します。

8.4.renshu.py

# 変数定義 =================================================================

# アトラクションを選ぶかどうかの変数
# 1:選ぶ / 0:選ばない
S = {}
for i in range(9):
    S[i] = solver.BoolVar("S%i" % i)

# 移動に関して、i→jに移動するかどうかの変数
# 1:移動する / 0:移動しない
M = {}
for i in range(9):
    for j in range(9):
        M[i,j] = solver.BoolVar("M%i%i" % (i,j))

上記、S,Mについては前に説明した通り、0/1を取る変数として宣言します。
どちらの値を取るかはOR-Toolsソルバーが計算して決めます。

8.4.renshu.py

# 満足度を宣言する
satisfy = solver.IntVar(0.0, 1000, 'satisfy')
# 滞在時間を宣言する
stayTime = solver.IntVar(0.0, 1000, 'stayTime')
# 移動時間変数
moveTime = solver.IntVar(0.0, 1000, 'moveTime')
# 全体時間の変数
totalTime = solver.IntVar(0.0, 1000, 'totalTime')

上記、満足度、滞在時間、移動時間、全体時間(滞在時間+移動時間)を、それぞれ0〜1000の整数を取る変数として宣言します。
定数及び変数を用意しましたので、次は条件を設定します(制約式を設定します)。

制約式の設定について、参考テキストに以下の補足がありますので、この条件も含めて組み立てを行います。

8.4.問題

・アトラクションを選んだら、訪れること。
・アトラクションに入ったら出ること。
・入口から接続していること。(途中で終了しないこと)

以下、制約条件のソースです。

8.4.renshu.py

# 制約式 ===================================================================

# 入口は必ず通過
solver.Add(S[0] == 1)

上記、入口に関する制約式を設定します。
入口は必ず選択するので、S[0]の値を1に設定します。

8.4.renshu.py

# 移動について制限を設定
# 移動しないルートについてMを0にする
solver.Add(M[0,0] == 0)
solver.Add(M[0,2] == 0)
solver.Add(M[0,5] == 0)
~

上記、定数定義と同様に、経路について設定を行います。Mは変数なのでソルバーが決定する場合もありますが、ここでは明示的に設定します。
移動しない経路についてM[i,j]を0にします。例えば、アトラクション0番目から0番目には移動できませんし、アトラクション0番目からアトラクション5番目についても同様です。全部で41通りの経路を0に設定します。

8.4.renshu.py

# 重複を避ける設定(アトラクションは1回だけ訪れる。両方訪れないケース可)
solver.Add(M[0,1] + M[1,0] <= 1)
solver.Add(M[0,3] + M[3,0] <= 1)
solver.Add(M[0,4] + M[4,0] <= 1)
~

上記、あるアトラクションiからjへ移動した後、移動経路を逆戻りして再びiへ戻ることは無いための条件を設定します。例えば、アトラクション0番目からアトラクション1番目に移動した場合はM[0,1]=1となるため、その逆経路のM[1,0]は値が0になる(移動しない)ように設定します。逆の移動経路の場合も同様です。
また、M[0,1] + M[1,0] == 1とすると必ず移動することになるので、移動しない場合も含めM[0,1] + M[1,0] <= 1とします。全部で20通りの経路について設定します。

8.4.renshu.py

# Mについての条件
# 一つのアトラクションから複数のアトラクションへ同時に移動しない
for i in range(9):
    solver.Add(solver.Sum([M[i,j] for j in range(9)]) <= 1)

上記、一つのアトラクションから複数のアトラクションへ同時に移動しない条件を設定します。これは、それぞれのアトラクションiにおいて、各jへ移動する経路は一つしか取れないこととして定義します(下図参照)。

次に、経路の連続性条件について次のように定義します。経路の連続性とは、下図のように経路がつながっていないアトラクションを選択しないという意味です。

経路がつながっている(連続している)条件の制約式は以下です。

8.4.renshu.py

# ルートの連続性条件
# 0起点
step = [1,3,4]
solver.Add(solver.Sum([M[0,s] for s in step]) == solver.Sum([M[s,0] for s in step]))
# 1起点
step = [0,2,3,4,5]
solver.Add(solver.Sum([M[1,s] for s in step]) == solver.Sum([M[s,1] for s in step]))
~

上記について、アトラクション1番目を起点に見てみます。
アトラクション1へ入る経路はM[0,1],M[3,1],M[4,1]の3通りあります。
下図(パターン1)の通り、例えばアトラクション1へ何らかの経路で入ったとします。
ここでは、M[0,1],M[3,1],M[4,1]のうち、M[4,1]で入ったと仮定しました(赤矢印)。

アトラクション1へ入る経路は3通りのうち1つしか取れませんので、M[0,1],M[3,1],M[4,1]のうちいずれか一つの値が1で、その他の値は0になります。
制約式として以下を設定します。

　solver.Sum([M[s,1] for s in step]) = 1　…　式①

制約式の設定条件に「アトラクションに入ったら出ること。」とありますので、アトラクション1から別のアトラクションへ移動する経路を考えます。
上図(パターン2)の通り、アトラクション1から移動できるアトラクションは3つありますが、この中から一つを選択することになります。3つの中から一つを選択する制約式は

　solver.Sum([M[1,s] for s in step]) = 1　…　式②

になります。但し、先ほどM[4,1]の経路を取ることを仮定したため、実際に取れる経路はM[1,2]かM[1,5]のいずれかです(上図(パターン3))。どの経路を取るかはソルバーが算出します。

また、アトラクション1へ入らない場合もあり得ますし、この場合は次へ移動しないことになります。
この場合は、入る経路と出る経路のいずれのMの和は0です。
よって、アトラクション1へ出入りするための経路の条件は、式①と式②を用いて以下とすれば良いことになります。

　solver.Sum([M[1,s] for s in step]) == solver.Sum([M[s,1] for s in step])

以上、経路の連続性の制約式をアトラクション1について説明しました。
実際のソース上では、アトラクション0からアトラクション8までそれぞれについて同様の考え方で定義します。残りのアトラクションに関しての制約式は、後方に掲載するソース全体で定義しています。

さらに、訪れるアトラクションと前後の出入りの経路の関係性について考えます。
経路の連続性条件と似ていますが、SとMの関係性を制約式として定義します。

8.4.renshu.py

# SとMの関係性の定義
# 0起点
step = [1,3,4]
solver.Add(S[0] - (solver.Sum([M[0,s] for s in step]) + solver.Sum([M[s,0] for s in step])) / 2 == 0)
# 1起点
step = [0,2,3,4,5]
solver.Add(S[1] - (solver.Sum([M[1,s] for s in step]) + solver.Sum([M[s,1] for s in step])) / 2 == 0)
~

上記について、制約式の内容は「あるアトラクションに訪れる場合は、経路に出入りがある。訪れない場合は、経路も取らない。」という意味です。
これについて、アトラクション1を起点に見てみます。まず以下の経路図を確認します。

アトラクション1を選択する場合は、それに入る経路と出る経路が考えられます。
入る経路についてパターン1を見ると、青経路M[0,1]、緑経路M[3,1]、赤経路M[4,1]のいずれかの経路をたどることになるため、以下が成立します。

　solver.Sum([M[s,1] for s in step]) = 1　…　式①

出る経路についても同様に各色の経路の和の条件は

　solver.Sum([M[1,s] for s in step]) = 1　…　式②

アトラクション1を訪れた場合はS[1]の値は1なので、この値と式①式②は以下のようにつなげることができます。

　S[1]の値　-　(式①　+　式②)　÷　2　=　0　…　式③

また、アトラクション1を訪れない場合についてもこの式③は成立します。S[1]の値は0となり、式①及び式②の値もそれぞれ0となるため、

　0　-　(0　+　0)　÷　2　=　0

よってSとMの関係性は式③を用いて定義することができ、以下とします。

　S[1] - (solver.Sum([M[1,s] for s in step]) + solver.Sum([M[s,1] for s in step])) / 2 == 0

実際のソース上では、アトラクション0からアトラクション8までそれぞれについて同様の考え方で定義します。残りのアトラクションに関しての制約式は、後方に掲載するソース全体で定義しています。

複雑な制約式の説明は以上となります。以下、残りの制約条件を定義します。

8.4.renshu.py

# 満足度の算出
satisfy = solver.Sum([S[i] * satisfyDef[i] for i in range(9)])

# 滞在時間の算出
stayTime = solver.Sum([S[i] * stayTimeDef[i] for i in range(9)])

# 移動時間の算出
moveTime = solver.Sum([M[i,j] * moveTimeDef[i,j] for i in range(9) for j in range(9)])

# 全体時間の算出
totalTime = stayTime + moveTime

上記、満足度及び滞在時間は、選択するアトラクションS[i]値とそれぞれの定数値をかけ、それらの和を取ります。アトラクションを選択するする場合はS[i]値が1、しない場合はS[i]が0なので、和を取ると選択するアトラクションの値だけ算出されます。

移動時間については、各i,jの組み合わせについてM[i,j]値と移動経路の定数値をかけ、それらの和を取ります。
また、全体時間は"滞在時間 + 移動時間"です。

最後に、「200分の制限時間の中で総満足度を最大化する。」を条件にします。

8.4.renshu.py

# 全体時間の制約
solver.Add(totalTime <= 200)

# 満足度を最大化する
solver.Maximize(satisfy)

上記、全体時間が200分以内の条件と、満足度を最大化する条件をそれぞれ分けて定義します。

制約式の設定は以上です。
以下、求解を実行します。

8.4.renshu.py

status = solver.Solve()

以下のソースで、解が見つかった場合、各iにおけるS[i]値と移動する場合のi,jにおけるM[i,j]値を表示するようにします。

8.4.renshu.py

if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Solution:')
    print('ok')
    print('Objective value =', solver.Objective().Value())
    print("culculate Time = ", solver.WallTime(), " milliseconds")

    print('satisfy', satisfy.solution_value())
    print('stayTime', stayTime.solution_value())
    print('moveTime', moveTime.solution_value())
    print('totalTime', totalTime.solution_value())
    for i in range(9):
        print(i,S[i].solution_value())
    for i in range(9):
        for j in range(9):
            if M[i,j].solution_value() >= 0.5:
                print('M',i,j,'=',M[i,j].solution_value())

else:
    print('The problem does not have an optimal solution.')

以上でプログラミングは終了です。
実際の実行結果は以下となりました。

Solution:
ok
Objective value = 405.0
culculate Time = 427 milliseconds
satisfy 405.0
stayTime 141.0
moveTime 59.0
totalTime 200.0
0 1.0
1 1.0
2 0.0
3 1.0
4 1.0
5 1.0
6 1.0
7 1.0
8 1.0
M 0 3 = 1.0
M 1 0 = 1.0
M 3 6 = 1.0
M 4 1 = 1.0
M 5 4 = 1.0
M 6 7 = 1.0
M 7 8 = 1.0
M 8 5 = 1.0

最適化計算の結果、デートコースが決定しました。
ちょうど200分で回る経路で、満足度は405です。

また、滞在時間と移動時間の内訳は、滞在時間が141分、移動時間が51分とわかりました。

選択するアトラクションを見てみると、S[2]以外の値が1なので、
ボート以外のアトラクションを訪れる経路となっています。

また、周回する経路はM[i,j]の値から、
　0 → 3 → 6 → 7 → 8 → 5 → 4 →1→0
という答えが出ました。

以下、プログラム全体を示します。

8.4.renshu.py

from __future__ import print_function
from ortools.linear_solver import pywraplp

# Create the mip solver with the CBC backend.
solver = pywraplp.Solver('simple_mip_program',
    pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)


# 定数定義 =================================================================

# アトラクション名
attrDef = ['','入口','喫茶','ボート','カップ','レストラン','観覧車',
           'お化け屋敷','コースター','迷路']

# 満足度定数
satisfyDef = [0,50,36,45,79,55,63,71,42]

# 滞在時間定数
stayTimeDef = [0,20,28,15,35,17,18,14,22]

# 移動時間定数
# moveTimeDef[i,j] ← iからjへ移動するときの時間
moveTimeDef = {}
moveTimeDef[0,0] = 0
moveTimeDef[0,1] = 9
moveTimeDef[0,2] = 0
moveTimeDef[0,3] = 7
moveTimeDef[0,4] = 12
moveTimeDef[0,5] = 0
moveTimeDef[0,6] = 0
moveTimeDef[0,7] = 0
moveTimeDef[0,8] = 0
moveTimeDef[1,0] = 9
moveTimeDef[1,1] = 0
moveTimeDef[1,2] = 11
moveTimeDef[1,3] = 12
moveTimeDef[1,4] = 7
moveTimeDef[1,5] = 13
moveTimeDef[1,6] = 0
moveTimeDef[1,7] = 0
moveTimeDef[1,8] = 0
moveTimeDef[2,0] = 0
moveTimeDef[2,1] = 11
moveTimeDef[2,2] = 0
moveTimeDef[2,3] = 0
moveTimeDef[2,4] = 14
moveTimeDef[2,5] = 8
moveTimeDef[2,6] = 0
moveTimeDef[2,7] = 0
moveTimeDef[2,8] = 0
moveTimeDef[3,0] = 7
moveTimeDef[3,1] = 12
moveTimeDef[3,2] = 0
moveTimeDef[3,3] = 0
moveTimeDef[3,4] = 11
moveTimeDef[3,5] = 0
moveTimeDef[3,6] = 7
moveTimeDef[3,7] = 12
moveTimeDef[3,8] = 0
moveTimeDef[4,0] = 12
moveTimeDef[4,1] = 7
moveTimeDef[4,2] = 14
moveTimeDef[4,3] = 11
moveTimeDef[4,4] = 0
moveTimeDef[4,5] = 9
moveTimeDef[4,6] = 13
moveTimeDef[4,7] = 9
moveTimeDef[4,8] = 13
moveTimeDef[5,0] = 0
moveTimeDef[5,1] = 13
moveTimeDef[5,2] = 8
moveTimeDef[5,3] = 0
moveTimeDef[5,4] = 9
moveTimeDef[5,5] = 0
moveTimeDef[5,6] = 0
moveTimeDef[5,7] = 13
moveTimeDef[5,8] = 7
moveTimeDef[6,0] = 0
moveTimeDef[6,1] = 0
moveTimeDef[6,2] = 0
moveTimeDef[6,3] = 7
moveTimeDef[6,4] = 13
moveTimeDef[6,5] = 0
moveTimeDef[6,6] = 0
moveTimeDef[6,7] = 7
moveTimeDef[6,8] = 0
moveTimeDef[7,0] = 0
moveTimeDef[7,1] = 0
moveTimeDef[7,2] = 0
moveTimeDef[7,3] = 12
moveTimeDef[7,4] = 9
moveTimeDef[7,5] = 13
moveTimeDef[7,6] = 7
moveTimeDef[7,7] = 0
moveTimeDef[7,8] = 6
moveTimeDef[8,0] = 0
moveTimeDef[8,1] = 0
moveTimeDef[8,2] = 0
moveTimeDef[8,3] = 0
moveTimeDef[8,4] = 13
moveTimeDef[8,5] = 7
moveTimeDef[8,6] = 0
moveTimeDef[8,7] = 6
moveTimeDef[8,8] = 0


# 変数定義 =================================================================

# アトラクションを選ぶかどうかの変数
# 1:選ぶ / 0:選ばない
S = {}
for i in range(9):
    S[i] = solver.BoolVar("S%i" % i)

# 移動に関して、i→jに移動するかどうかの変数
# 1:移動する / 0:移動しない
M = {}
for i in range(9):
    for j in range(9):
        M[i,j] = solver.BoolVar("M%i%i" % (i,j))

# 満足度を宣言する
satisfy = solver.IntVar(0.0, 1000, 'satisfy')
# 滞在時間を宣言する
stayTime = solver.IntVar(0.0, 1000, 'stayTime')
# 移動時間変数
moveTime = solver.IntVar(0.0, 1000, 'moveTime')
# 全体時間の変数
totalTime = solver.IntVar(0.0, 1000, 'totalTime')


# 制約式 ===================================================================

# 入口は必ず通過
solver.Add(S[0] == 1)

# 移動について制限を設定
# 移動しないルートについてMを0にする
solver.Add(M[0,0] == 0)
solver.Add(M[0,2] == 0)
solver.Add(M[0,5] == 0)
solver.Add(M[0,6] == 0)
solver.Add(M[0,7] == 0)
solver.Add(M[0,8] == 0)
solver.Add(M[1,1] == 0)
solver.Add(M[1,6] == 0)
solver.Add(M[1,7] == 0)
solver.Add(M[1,8] == 0)
solver.Add(M[2,0] == 0)
solver.Add(M[2,2] == 0)
solver.Add(M[2,3] == 0)
solver.Add(M[2,6] == 0)
solver.Add(M[2,7] == 0)
solver.Add(M[2,8] == 0)
solver.Add(M[3,2] == 0)
solver.Add(M[3,3] == 0)
solver.Add(M[3,5] == 0)
solver.Add(M[3,8] == 0)
solver.Add(M[4,4] == 0)
solver.Add(M[5,0] == 0)
solver.Add(M[5,3] == 0)
solver.Add(M[5,5] == 0)
solver.Add(M[5,6] == 0)
solver.Add(M[6,0] == 0)
solver.Add(M[6,1] == 0)
solver.Add(M[6,2] == 0)
solver.Add(M[6,5] == 0)
solver.Add(M[6,6] == 0)
solver.Add(M[6,8] == 0)
solver.Add(M[7,0] == 0)
solver.Add(M[7,1] == 0)
solver.Add(M[7,2] == 0)
solver.Add(M[7,7] == 0)
solver.Add(M[8,0] == 0)
solver.Add(M[8,1] == 0)
solver.Add(M[8,2] == 0)
solver.Add(M[8,3] == 0)
solver.Add(M[8,6] == 0)
solver.Add(M[8,8] == 0)

# 重複を避ける設定(アトラクションは1回だけ訪れる。両方訪れないケース可)
solver.Add(M[0,1] + M[1,0] <= 1)
solver.Add(M[0,3] + M[3,0] <= 1)
solver.Add(M[0,4] + M[4,0] <= 1)
solver.Add(M[1,2] + M[2,1] <= 1)
solver.Add(M[1,3] + M[3,1] <= 1)
solver.Add(M[1,4] + M[4,1] <= 1)
solver.Add(M[1,5] + M[5,1] <= 1)
solver.Add(M[2,4] + M[4,2] <= 1)
solver.Add(M[2,5] + M[5,2] <= 1)
solver.Add(M[3,4] + M[4,3] <= 1)
solver.Add(M[3,6] + M[6,3] <= 1)
solver.Add(M[3,7] + M[7,3] <= 1)
solver.Add(M[4,5] + M[5,4] <= 1)
solver.Add(M[4,6] + M[6,4] <= 1)
solver.Add(M[4,7] + M[7,4] <= 1)
solver.Add(M[4,8] + M[8,4] <= 1)
solver.Add(M[5,7] + M[7,5] <= 1)
solver.Add(M[5,8] + M[8,5] <= 1)
solver.Add(M[6,7] + M[7,6] <= 1)
solver.Add(M[7,8] + M[8,7] <= 1)


# Mについての条件
# 一つのアトラクションから複数のアトラクションへ同時に移動しない
for i in range(9):
    solver.Add(solver.Sum([M[i,j] for j in range(9)]) <= 1)

# ルートの連続性条件
# 0起点
step = [1,3,4]
solver.Add(solver.Sum([M[0,s] for s in step]) == solver.Sum([M[s,0] for s in step]))
# 1起点
step = [0,2,3,4,5]
solver.Add(solver.Sum([M[1,s] for s in step]) == solver.Sum([M[s,1] for s in step]))
# 2起点
step = [1,4,5]
solver.Add(solver.Sum([M[2,s] for s in step]) == solver.Sum([M[s,2] for s in step]))
# 3起点
step = [0,1,4,6,7]
solver.Add(solver.Sum([M[3,s] for s in step]) == solver.Sum([M[s,3] for s in step]))
# 4起点
step = [0,1,2,3,5,6,7,8]
solver.Add(solver.Sum([M[4,s] for s in step]) == solver.Sum([M[s,4] for s in step]))
# 5起点
step = [1,2,4,7,8]
solver.Add(solver.Sum([M[5,s] for s in step]) == solver.Sum([M[s,5] for s in step]))
# 6起点
step = [3,4,7]
solver.Add(solver.Sum([M[6,s] for s in step]) == solver.Sum([M[s,6] for s in step]))
# 7起点
step = [3,4,5,6,8]
solver.Add(solver.Sum([M[7,s] for s in step]) == solver.Sum([M[s,7] for s in step]))
# 8起点
step = [4,5,7]
solver.Add(solver.Sum([M[8,s] for s in step]) == solver.Sum([M[s,8] for s in step]))

# SとMの関係性の定義
# 0起点
step = [1,3,4]
solver.Add(S[0] - (solver.Sum([M[0,s] for s in step]) + solver.Sum([M[s,0] for s in step])) / 2 == 0)
# 1起点
step = [0,2,3,4,5]
solver.Add(S[1] - (solver.Sum([M[1,s] for s in step]) + solver.Sum([M[s,1] for s in step])) / 2 == 0)
# 2起点
step = [1,4,5]
solver.Add(S[2] - (solver.Sum([M[2,s] for s in step]) + solver.Sum([M[s,2] for s in step])) / 2 == 0)
# 3起点
step = [0,1,4,6,7]
solver.Add(S[3] - (solver.Sum([M[3,s] for s in step]) + solver.Sum([M[s,3] for s in step])) / 2 == 0)
# 4起点
step = [0,1,2,3,5,6,7,8]
solver.Add(S[4] - (solver.Sum([M[4,s] for s in step]) + solver.Sum([M[s,4] for s in step])) / 2 == 0)
# 5起点
step = [1,2,4,7,8]
solver.Add(S[5] - (solver.Sum([M[5,s] for s in step]) + solver.Sum([M[s,5] for s in step])) / 2 == 0)
# 6起点
step = [3,4,7]
solver.Add(S[6] - (solver.Sum([M[6,s] for s in step]) + solver.Sum([M[s,6] for s in step])) / 2 == 0)
# 7起点
step = [3,4,5,6,8]
solver.Add(S[7] - (solver.Sum([M[7,s] for s in step]) + solver.Sum([M[s,7] for s in step])) / 2 == 0)
# 8起点
step = [4,5,7]
solver.Add(S[8] - (solver.Sum([M[8,s] for s in step]) + solver.Sum([M[s,8] for s in step])) / 2 == 0)

# 満足度の算出
satisfy = solver.Sum([S[i] * satisfyDef[i] for i in range(9)])

# 滞在時間の算出
stayTime = solver.Sum([S[i] * stayTimeDef[i] for i in range(9)])

# 移動時間の算出
moveTime = solver.Sum([M[i,j] * moveTimeDef[i,j] for i in range(9) for j in range(9)])

# 全体時間の算出
totalTime = stayTime + moveTime

# 全体時間の制約
solver.Add(totalTime <= 200)

# 満足度を最大化する
solver.Maximize(satisfy)


# 求解実行
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Solution:')
    print('ok')
    print('Objective value =', solver.Objective().Value())
    print("culculate Time = ", solver.WallTime(), " milliseconds")

    print('satisfy', satisfy.solution_value())
    print('stayTime', stayTime.solution_value())
    print('moveTime', moveTime.solution_value())
    print('totalTime', totalTime.solution_value())
    for i in range(9):
        print(i,S[i].solution_value())
    for i in range(9):
        for j in range(9):
            if M[i,j].solution_value() >= 0.5:
                print('M',i,j,'=',M[i,j].solution_value())

else:
    print('The problem does not have an optimal solution.')

出展

本記事は、数理最適化についての参考テキスト「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」に記載の練習問題を参考にさせて頂いております。

■参考テキスト
「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」
　斉藤努[著]
　近代科学社[出版]

ご意見など

ご意見、間違い訂正などございましたらお寄せ下さい。