こちらはZOZO Advent Calendar 2024 シリーズ8 7日目の記事です。

はじめに

こんにちは、kokoakumaと申します。普段はAndroidアプリを担当していますが、放送大学の科目として数理最適化が提供されていたのを見つけたため、入門してみました。今回はその記録です。

数理最適化とは

数理最適化とは、特定の制約条件の中で、対象となる目的関数の最適解や近似解を求める手法です。
売上の最大化や費用の最小化といった問題が例としてよく出てきます。

使った教材と感想

マンガでわかる数理最適化
放送大学「数理最適化演習」
放送大学「問題解決の数理」

マンガでわかる数理最適化

全体の雰囲気を掴むために使用しました。
問題の定式化や解を見つける過程をイメージしやすくなり、一冊目として良い本でした。ただ、数学的な説明は少ないため、他の本で補う必要があります。

放送大学「数理最適化演習」

放送大学にて提供されている科目です。最小費用問題や最短路問題など、基本的な問題について説明動画と練習問題を通して学習できます。基本的には、問題の性質を把握し、定式化する能力を身につけることに重きが置かれていました。個人的には、後述する「問題解決の数理」だけでは定式化のイメージが掴みづらかったため、取ってよかったと感じる科目でした。また、数式を含むレポート課題の書き方を丁寧に説明してくれている点も、個人的にとても嬉しいポイントでした。

放送大学「問題解決の数理」

こちらも放送大学にて提供されている科目です。上述の演習科目と扱っている問題は共通していますが、こちらは、数式として定式化するまでの過程や解を求める手法の説明にページを割いています。数式や解を求める手法といっても連立一次方程式で表されるものが多いですが、いきなり数式だけ追っていくと、対象の数式が何を表すのかを見失いがちです。そのため、自分は、演習科目やその他の本でイメージを掴んでから、数式を一つずつ追っていきました。正直なところ、すべては理解していませんが、それはおいおい頑張ることとします。

解を求めてみる

講義の感想だけだと少し短いため、以下の経路が存在する場合の最短路を試しに求めてみます。
矢印に振ってある数字は各枝の重みです。

定式化する

まずは問題を定式化します。

目的関数

点の集合Vと枝の集合Eを以下のように定義します。

V = \{1,2,3,4,5\}

E = \{ (1,2), (1,3),(2,3), (2,4),(3,2),(3,4), (3,5), (4,5)\}

また、それぞれの枝を通る場合はXij=1（iは経路の始点、jは経路の終点）、通らない場合はXij=0と表します。

X_{ij} \in \{0,1\}, (i,j) \in E

その場合、目的関数は以下のようになります。

Z = 5X_{12} + X_{13} + 5X_{23} + X_{24} + 2X_{32} + 3X_{34} + 5X_{35} + X_{45}

今回は最短路を求める問題のため、Zの最小化を目指します。

制約条件

　前提として最短路は一つ見つければよいため、最短路は一本であるとします。まず点1について、点1から出る枝の一つは最短路に含まれます。また、点1に入る枝はありません。よって、点1から出る枝の和は1になるため、

X_{12} + X_{13} = 1

となります。

　次に終点である点5の制約を考えます。点5から出る枝はなく、点5に入る枝の一つは最短路に含まれます。よって、入る枝の和は-1になるため、

- (X_{34} + X_{35} ) = -1

となります。
　他の点は、

（最短路に含まれる場合）入る枝と出る枝の内、一つずつ最短路に含まれます
（最短路に含まれない場合）入る枝と出る枝が最短路に含まれません

となるため、出る枝の和から入る枝の和を引くと0になります。
よって、

\displaylines{
X_{23} + X_{24} - (X_{12} + X_{32})= 0 \\
X_{32} + X_{34} + X_{35} - (X_{13}) = 0 \\ 
X_{45} - (X_{24} + X_{34}) = 0 \\
}

と表せます。

ここまでの内容をまとめると、以下の0-1整数最適化問題として定式化できることが分かりました。

\displaylines{
最小化  \\
Z = 5X_{12} + X_{13} + 5X_{23} + X_{24} + 2X_{32} + 3X_{34} + 5X_{35} + X_{45} \\

制約条件 \\ 
V = \{1,2,3,4,5\} \\
E = \{ (1,2), (1,3),(2,3), (2,4),(3,2),(3,4), (3,5), (4,5)\} \\
X_{ij} \in \{0,1\}, (i,j) \in E \\
X_{12} + X_{13} = 1 \\
- (X_{34} + X_{35} ) = -1 \\
X_{23} + X_{24} - (X_{12} + X_{32})= 0 \\
X_{32} + X_{34} + X_{35} - (X_{13}) = 0 \\ 
X_{45} - (X_{24} + X_{34}) = 0 \\
}

ライブラリを使用して解を求める

次に、PythonのライブラリであるPuLPを使用して解を求めてみます。
以下はGoogle Colaboratoryでの使用を想定したコードです。

!pip install pulp
import pulp

vertexList = [i+1 for i in range(5)]
edgeList = [(1,2), (1,3), (2,3), (2,4),(3,2), (3,4), (3,5), (4,5)]
weightList = {(1,2):5, (1,3):1, (2,3):5, (2,4):1,(3,2):2,(3,4):4, (3,5):5, (4,5):1}
start = 1
terminal = 5

x = {edge:pulp.LpVariable(f'x{edge[0]}_{edge[1]}', cat=pulp.LpBinary) for edge in edgeList}

problem = pulp.LpProblem('最短路問題', sense=pulp.LpMinimize)
problem += pulp.lpSum(weightList[edge]*x[edge] for edge in edgeList), '目的関数　経路の長さ'
for vertex in vertexList:
    in_flow = []
    out_flow = []
    for edge in edgeList:
        if edge[1] == vertex:
            in_flow.append(edge)
        if edge[0] == vertex:
            out_flow.append(edge)
    net_outflow = pulp.lpSum(x[edge] for edge in out_flow) - pulp.lpSum(x[edge] for edge in in_flow)
    if vertex == start:
        problem += net_outflow == 1, f'スタート地点{vertex}'
    elif vertex == terminal:
        problem += net_outflow == -1, f'ゴール地点{vertex}'
    else:
        problem += net_outflow == 0, f'点{vertex}'

result = problem.solve()
print(f"経路の長さは{pulp.value(problem.objective)}")
print("通る経路:")
for v in problem.variables():
    if pulp.value(v) > 0:
        print(f'- {v} = {pulp.value(v)}')

出力結果は以下の通りです。

経路の長さは5.0
通る経路:
- x1_3 = 1.0
- x2_4 = 1.0
- x3_2 = 1.0
- x4_5 = 1.0

1→3→2→4→5と進むのが正解のようですね。

アルゴリズムで解を求めてみる

動的計画法の一つであるダイクストラ法を用いて、同じ問題の解を求めてみます。


import math
from enum import Enum

# 到達地点の色分け
class Color(Enum):
    WHITE = 1
    GRAY = 2
    BLACK = 3
    
vertexList = [i+1 for i in range(5)]
edgeList = [(1,2), (1,3), (2,3), (2,4),(3,2), (3,4), (3,5), (4,5)]
weightList = {(1,2):5, (1,3):1, (2,3):5, (2,4):1,(3,2):2,(3,4):4, (3,5):5, (4,5):1}
start = 1
terminal = 5

reachedList = [Color.WHITE for _ in range(terminal)]
routeList = [None for _ in range(terminal)]
routeList[start-1] = -1
distanceList = [math.inf for _ in range(terminal)]
distanceList[start-1] = 0

weightTable = [[None for _ in range(terminal)] for _ in range(terminal)]
for edge in edgeList:
    weightTable[edge[0]-1][edge[1]-1] = weightList[edge]

while True:
    minCost = math.inf
    for i in range(terminal):
        if reachedList[i] != Color.BLACK and distanceList[i] < minCost:
            minCost = distanceList[i]
            u = i

    if minCost == math.inf:
        break
        
    reachedList[u] = Color.BLACK

    for v in range(terminal):
        if reachedList[v] != Color.BLACK and weightTable[u][v] != None:
            if distanceList[u] + weightTable[u][v] < distanceList[v]:
                distanceList[v] = distanceList[u] + weightTable[u][v]
                routeList[v] = u
                reachedList[v] = Color.GRAY
                
for index, route in enumerate(routeList):
    print(f"{route + 1} to {index + 1}")

for index, distance in enumerate(distanceList):
    print(f"{index + 1}'s distance is { distance}")

出力は以下の通りです。

0 to 1
3 to 2
1 to 3
2 to 4
4 to 5
1's distance is 0
2's distance is 3
3's distance is 1
4's distance is 4
5's distance is 5

先ほどと同じ結果になっていそうですね🙌

経路の長さは5
通る経路は1→3→2→4→5

さいごに

今回は数理最適化に入門した話を書きました。普段触っていない分野のため、理解しきれていない箇所も多々ありますが、良い気分転換になりました。
何かしら参考になれば幸いです。