はじめに

数理最適化問題を定式化する際に、例えば$x_{ij}$といった$i$とか$j$とかの添え字のついた変数や定数を定義して、$\sum$で和をとったりします。ここでは、これらに対応するGurobi OptimizerのPythonインターフェースを使ったコードを紹介します。

復習記事

想定環境

Windows 10 64bit
Gurobi Optimizer 8.0.1
Python 3.6

サンプル2

解きたい最適化問題

社員を仕事に割り当てるのだが、その際にかかる費用を最小化したい。

記法

記号	意味	このサンプルでの値
$I$	社員の集合	$I := \{Aさん, Bさん, Cさん\}$
$J$	仕事の集合	$J := \{東京で内勤, 大阪で出張\}$
$c_{ij}$	社員$i \in I$を仕事$j \in J$に割り当てた場合の費用を表す定数	$c_{Aさん, 東京で内勤} = 7$、など
$x_{ij}$	社員$i \in I$を仕事$j \in J$に割り当てれば$1$、そうでなければ$0$をとる変数

図

最適化問題

\begin{alignat}{3}
 &\mbox{最小化} 
 &\qquad \sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij} & 
 &\qquad & \\
 &\mbox{制約} 
 & \sum_{j \in J} x_{ij} &\leq 1 
 & &\forall i \in I \\
 &
 &\sum_{i \in I} x_{ij} &\geq 1
 & &\forall j \in J \\
 &
 & x_{ij} &\in \{0,1\}
 & &\forall i \in I, \quad \forall j \in J.
\end{alignat}

式の意味

1つめ | 人を仕事に割り当てたときの費用の合計
2つめ | 各人について、割り当てられる仕事は最大で1つまで
3つめ | 各仕事について、割り当てる人は最低で1人
4つめ | 各変数がとり得る値は0か1のどちらか

ソース一式

データをExcelファイルで定義し、Pythonでそれを読み込んで最適化の計算をするコードになっています。一式を↓に置いたので、ダウンロードしてください。
https://github.com/keisukesato-ac/gurobi-sample-2

Pythonコードと実行結果

gurobi_sample_2.py

# coding: utf-8

# Dummy Japanese character: あ（意味はないが確実に日本語を含むファイルにする）



# Gurobi Optimizerのパッケージをインポートし、
# "gp"という短縮名を設定（違う短縮名にしてもよい）
import gurobipy as gp

# Excelファイルを読み込むためにpandasをインポート
import pandas as pd



# クラスを定義
# 社員クラス
class Shain():
    def __init__(self, row):
        self.code = row[0]
        self.name = row[1]

    def __repr__(self):
        return f"{self.code}_{self.name}"

# 仕事クラス
class Shigoto():
    def __init__(self, row):
        self.code = row[0]
        self.name = row[1]

    def __repr__(self):
        return f"{self.code}_{self.name}"



# データを読み込み
excel_file = pd.ExcelFile("./Data.xlsx")

# 「社員」シートをDataFrameとして読み込み
df_shain = excel_file.parse("社員")
# 「社員コード」「社員名」の列を抽出
df_shain_part = df_shain[["社員コード", "社員名"]]
# listに変換
list_shain = df_shain_part.values.tolist()
# 社員オブジェクト集合を定義（※便宜上、setの代わりにlistを使用）
I = list()
# 「社員コード : 社員オブジェクト」の辞書を定義
dict_I = dict()
# 社員オブジェクト集合にShainクラスのオブジェクトを追加し、
# 「社員コード : 社員オブジェクト」の辞書にも追加
for row in list_shain:
    shain = Shain(row)
    I.append(shain)
    dict_I[row[0]] = shain
print(f"[集合・定数]")
print(f"    社員集合: {I}")

# 「仕事」シートについて同様の処理
J = list(
    [Shigoto(row) 
     for row 
     in excel_file.parse("仕事")[["仕事コード", "仕事名"]].values.tolist()])
dict_J = {shigoto.code : shigoto for shigoto in J}
print(f"    仕事集合: {J}")


# 「費用」シートをDataFrameとして読み込み、
# 「社員コード」「仕事コード」「費用」の列を抽出し、listに変換
list_hiyo = excel_file.parse("費用")[
    ["社員コード", "仕事コード", "費用"]].values.tolist()
# c[社員オブジェクト, 仕事オブジェクト] を 費用 を表す辞書として定義
c = {}
# c[社員オブジェクト, 仕事オブジェクト] に 費用 を設定
for (shain_code, shigoto_code, hiyo) in list_hiyo:
    shain = dict_I[shain_code]
    shigoto = dict_J[shigoto_code]
    c[shain, shigoto] = hiyo
print(f"    費用定数: {c}")
print()



# 問題を設定
model_2 = gp.Model(name = "GurobiSample2")


# 変数を設定（変数単体にかかる制約を含む）
# x[社員オブジェクト, 仕事オブジェクト] を gp.Var クラスのオブジェクト を
# 表す辞書として定義
x = {}
for i in I:
    for j in J:
        x[i, j] = model_2.addVar(vtype = gp.GRB.BINARY, name = f"x({i},{j})")


# 目的関数を設定
model_2.setObjective(gp.quicksum(c[i, j] * x[i, j] for i in I for j in J), 
                     sense = gp.GRB.MINIMIZE)


# 制約を設定
con_1 = {}
for i in I:
    con_1[i] = model_2.addConstr(gp.quicksum(x[i, j] for j in J) <= 1, 
                                 name = f"con_1({i})")
con_2 = {}
for j in J:
    con_2[j] = model_2.addConstr(gp.quicksum(x[i, j] for i in I) >= 1, 
                                 name = f"con_2({j})")


# 解を求める計算
print("[Gurobi Optimizerログ]")

model_2.optimize()

print()


# 最適解が得られた場合、結果を出力
print("[解]")
if model_2.Status == gp.GRB.OPTIMAL:
    print("    最適解: ")
    for i in I:
        for j in J:
            if x[i, j].X > 0.98:
                print(f"        {i.name}に{j.name}を割り当て")
    val_opt = model_2.ObjVal
    print(f"    最適値: {val_opt}")

実行結果

[集合・定数]
    社員集合: [11_Aさん, 12_Bさん, 13_Cさん]
    仕事集合: [201_東京で内勤, 202_大阪に出張]
    費用定数: {(11_Aさん, 201_東京で内勤): 7, (11_Aさん, 202_大阪に出張): 2, (12_Bさん, 201_東京で内勤): 10, (12_Bさん, 202_大阪に出張): 13, (13_Cさん, 201_東京で内勤): 3, (13_Cさん, 202_大阪に出張): 6}

[Gurobi Optimizerログ]
Optimize a model with 5 rows, 6 columns and 12 nonzeros
Variable types: 0 continuous, 6 integer (6 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [2e+00, 1e+01]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Presolve time: 0.00s
Presolved: 5 rows, 6 columns, 12 nonzeros
Variable types: 0 continuous, 6 integer (6 binary)
Found heuristic solution: objective 5.0000000

Root relaxation: cutoff, 0 iterations, 0.00 seconds

Explored 0 nodes (0 simplex iterations) in 0.00 seconds
Thread count was 4 (of 4 available processors)

Solution count 1: 5 

Optimal solution found (tolerance 1.00e-04)
Best objective 5.000000000000e+00, best bound 5.000000000000e+00, gap 0.0000%

[解]
    最適解: 
        Aさんに大阪に出張を割り当て
        Cさんに東京で内勤を割り当て
    最適値: 5.0

解説

Excelファイルを読む部分は、回りくどいコードになっており、すみません。意図としては、集合があって、その各要素が同じようなタイプのデータであるので、データをオブジェクトとして扱ったほうがよいかなあというものです。

model_2.setObjective(gp.quicksum(c[i, j] * x[i, j] for i in I for j in J), 
                     sense = gp.GRB.MINIMIZE)

数式の$\sum$を実現するのに、gp.quicksum()と書きます。()の中はPythonの内包表記っぽい書き方になります。
応用として、例えばgp.quicksum(c[i, j] * x[i, j] for i in I for j in J if i.name != "Bさん")と書くと、Bさんが$\sum$の範囲から外れます。数式で書くと$\sum_{i \in I \setminus \{ Bさん \}}$となるということです。

            if x[i, j].X > 0.98:
                print(f"        {i.name}に{j.name}を割り当て")

$x_{ij}$を確かに0-1変数として定義しましたが、計算結果を取得すると、もろもろの関係で、0.00001とかいう値や0.99998とかいう値が返ってくることがあります。整数値のはずの計算結果は少々の誤差を含んだ小数として返ってくると思って対処するのがよいでしょう。
この仕様は、私が以前に使用していた別の数理最適化ソルバーでも同じでした。

おわりに

読んでいただいた方がご自分の数理最適化問題を実装して解が求められるようになれば幸いです。なにかあればコメントください。

Gurobi Optimizer + Pythonのサンプルコード（その2）