最適化モデリング入門＠新人研修

はじめに

2023年10月にシニアリサーチャーとして入社したアドバンスドテクノロジーラボ（ATL）の梅谷俊治です。 2023年9月まで、大阪大学大学院情報科学研究科にて数理最適化寄附講座教授を務めていました。本記事では、データ推進室に配属される新人を対象とした「最適化モデリング入門」をご紹介します。

データ推進室の新人は、部署に配属される前に「BootCamp」と呼ばれる研修を受けます。この研修は約3ヶ月にわたり、エンジニアリング技術の向上だけでなく、ヒューマンスキルの育成や事業理解を深めるための幅広いプログラムが用意されています。具体的には、データサイエンスの基礎、プログラミング、データベース管理の分野で徹底的なトレーニングが行われます。例えば、PythonやRを使用したデータ解析の手法を学び、SQLを用いたデータベースの操作方法を習得します。また、実際のデータ分析の問題を解決するプロジェクトにチームで取り組む機会もあり、理論だけではなく実践的なスキルを磨くことができます。研修の一部をブログにて紹介していますのでこちらもご覧ください。

本記事は、データ推進室 Advent Calendar 2024 17日目の記事です

講義の概要

私の担当は4月10日と11日の2日間でした。前職の大学では数理最適化の授業を担当し、線形計画問題、非線形計画問題、整数計画問題など代表的な最適化問題に対する基本的なアルゴリズムの考え方を半年間かけて教えていました。しかし、企業の研修は時間が限られますし、受講者がすぐに実務に活用できる知識を教授することが望ましいです。数理最適化を実務に活用する際には、(1)最適化問題の定式化、(2)アルゴリズムによる求解、(3)計算結果の分析・検証、(4)最適化問題の修正という一連の手続きを踏みます。

現在では、多くの数理最適化ソルバー(最適化問題を解くソフトウェア)が公開されており、自身でアルゴリズムを設計・実装しなくても、数理最適化を実務に活用できるようになりました。そこで、本講義では、データサイエンティストが最速で数理最適化を実務に活用することを目指して、最適化問題の定式化と数理最適化ソルバーの利用法に焦点を当てた以下のカリキュラムを準備しました。

数理最適化の概要
組合せ最適化問題とその応用
現実問題への数理最適化の適用
線形計画問題の定式化
整数計画問題の定式化
整数計画ソルバーの利用法

1日目に上記の内容を講義した後に、2日目には、受講者の皆さんに、Pythonライブラリ Python-MIP を用いて、組合せ最適化問題を解くという演習に取り組んでもらいました。はじめに、ナップサック問題、ビンパッキング問題、機械スケジューリング問題を例題に、サンプルコードを見せながらPython-MIPの使い方を説明しました。その後に、受講者の皆さんにナーススケジューリング問題に取り組んでもらいました。最適化問題の定式化からプログラムの作成と、半日で取り組むにはタフな課題だと思いましたが、新人研修を企画している人事の皆さんから遠慮なくやっちゃって下さいとのお言葉をいただいたので、難易度は下げずに、つまずくようなことがあれば、一人で解決しようとせずに受講者の皆さんで一緒に相談しながら進めて下さいと声をかけるのみに留めました。一方で、最適化問題に限らず、プログラミング演習では、受講者が入力データの読込みや出力データの表示の実装に手間取ることが少なくありません。特に、ナーススケジューリング問題で計算結果をシフト表に出力するのは簡単ではありません。そこで、こちらで入出力部分のサンプルコードを用意して、受講者が整数計画問題の定式化に集中できるように努めました。

当日の様子ですが、受講者には、大学で数理最適化を研究していた人、競技プログラミングやデータ分析コンペティションに取り組んでいた人が何人かいたようで、1日目の講義では専門家でなければ思い付かないような鋭い質問がポンポンと飛び出し、2日目の演習では30分足らずで課題を済ませ、私が焦ってしまう場面もありました。

以降では、1日目後半の「線形計画問題の定式化」「整数計画問題の定式化」と2日目の演習内容の一部をかいつまんで紹介します。

線形計画問題の定式化

現実問題を最適化問題に定式化する際には、売上げを最大化したい、費用を最小化したいなど全体のKPIの改善が常に求められると思われがちです。しかし、実際には、予算の配分額や従業員の勤務時間などのバラつきを小さくしたい、平準化したいという事例が少なくありません。例えば、あれも良くしたい、これも良くしたいと複数の目的関数を同時に最小化する多目的最適化問題を考えたくなることがよくあります。このような事例では、各目的関数に重み係数 $w_k$ を掛けた重み付き和を最小化する定式化がよく使われます。

$$ \begin{array}{ll} \textnormal{最小化} & \displaystyle w_1 \sum_{j=1}^n c_{1j} x_j + w_2 \sum_{j=1}^n c_{2j} x_j + \cdots + w_k \sum_{j=1}^n c_{kj} x_j \cr \textnormal{条件} & \displaystyle\sum_{j=1}^n a_{ij} x_j \ge b_i, \quad i=1,\dots,m, \cr & x_j \ge 0, \quad j=1,\dots,n. \end{array} $$

しかし、実際には、こちらを立てればあちらが立たないということが方々で生じて、全ての目的関数をバランス良く最小化するにはほど遠く頭を抱えることが少なくありません。頑張って重み係数 $w_k$ の値を色々と変えて試行錯誤してもなかなか良い感じの結果は得られません。バラつきを小さくするなら分散いわゆる目的関数の2乗和を最小化すれば良い気もしますが、目的関数が2次関数の最適化問題は求解が難しくなる上に、実際には思ったほどにバラつきは小さくなりません。

そこで、方針を切り替えて、以下の通り、最も結果が悪いすなわち値が最大の目的関数を最小化することを考えます。

$$ \begin{array}{ll} \textnormal{最小化} & \displaystyle \max \{ \sum_{j=1}^n c_{1j} x_j, \sum_{j=1}^n c_{2j} x_j, \dots, \sum_{j=1}^n c_{kj} x_j \} \end{array} $$

ただし、このままでは線形関数にならないので、トリックを導入して線形計画問題に書き換えます。具体的には、全ての目的関数の最大値を表す人工変数 $z$ を導入して、以下の通りに書き換えます。要するに、全ての目的関数の値が $z$ 以下となる制約条件を追加した上で、 $z$ の値を最小化する定式化に書き換えたわけです。

$$ \begin{array}{ll} \textnormal{最小化} & z \cr \textnormal{条件} & \displaystyle\sum_{j=1}^n c_{1j} x_j \le z,\cr & \displaystyle\sum_{j=1}^n c_{2j} x_j \le z, \cr & \cdots \cr & \displaystyle\sum_{j=1}^n c_{kj} x_j \le z, \cr \end{array} $$

図書館における雑誌の購読計画にこの定式化を適用した結果を紹介します。図書館では、幅広い分野の雑誌の需要がある一方で購読に使える予算は限られています。下図は雑誌の需要に対してどれだけ購読できているかの充足率をグラフで表しています。全体の充足率を最大化してしまうと左図のように分野ごとに大きなバラつきが生じます。そこで、最小の充足率を最大化するという定式化を導入すると、右図のように分野ごとのバラつきを見事に抑えることができます。

このように発想の柔軟さが求められることが、最適化問題の定式化の難しくかつ面白い部分です。

ビンパッキング問題

2日目の演習で例題に取りあげたビンパッキング問題の定式化とPython-MIPによる実装を紹介します。ビンパッキング問題は、 $n$ 個の荷物を詰め込むのに必要な箱の数を最小化する問題です。ビン(bin)は英語で保管用の箱を指す言葉で瓶ではないことに注意して下さい。箱はいずれも同じ大きさで、詰め込める荷物の重さ合計を $C$ とします。また、荷物 $j$ の重さを $w_j (< C)$ とします。

変数 $x_{ij}$ と $y_i$ を用意して、荷物 $j$ が箱 $i$ に入っていれば $x_{ij} = 1$ 、そうでなければ $x_{ij}=0$ 、箱 $i$ に1つでも荷物が入っていれば $y_i=1$ 、そうでなければ $y_i=0$ とします。このとき、ビンパッキング問題は以下の整数計画問題に定式化できます。

$$ \begin{array}{ll} \textnormal{最小化} & \displaystyle\sum_{i=1}^n y_i \cr \textnormal{条件} & \displaystyle\sum_{i=1}^n x_{ij} = 1, \quad j = 1,\dots,n, \cr & \displaystyle\sum_{j=1}^n w_j x_{ij} \le C y_i, \quad i=1,\dots,n, \cr & x_{ij} \in \{ 0,1 \}, \quad i=1,\dots,n, ; j=1,\dots,n, \cr & y_i \in \{ 0,1 \}, \quad i=1,\dots,n. \end{array} $$

1番目の制約条件は、各荷物 $j$ をちょうど1つの箱に割り当てることを表します。 2番目の制約条件は、箱 $i$ に詰め込んだ荷物の重さ合計が上限 $C$ を超えないことを表します。この制約条件の右辺が $C y_i$ となっているところがビンパッキング問題の特徴です。箱 $i$ に1つでも荷物を詰め込むと変数 $y_i=1$ とする必要が生じます。この制約条件によって箱の利用の有無も併せてチェックしています。

現実問題では、コストが必ずしも取り扱い量に比例するわけではなく、ほんのわずかでも取り扱い量があれば段取り替えや設備費などの固定費用が発生することが少なくありません。ビンパッキング問題は、製造や物流など幅広い分野に現れる固定費用を考慮した最適化問題のひな形としてとらえることができます。

では、Python-MIPを使ってビンパッキング問題を解くプログラムを実装しましょう。 Python-MIPでは、(1)入力データの読込み、(2)変数の定義、(3)制約条件・目的関数の定義、(4)数理最適化ソルバーの実行、(5)計算結果の取得と出力という手順で処理します。これは、Python-MIP独自ではなく、他の数理最適化ソルバーのPython APIでも同じ手順を踏みます。

まず、以下の通りcsv形式の入力データを用意しました。

num,cap
20,100
id,wt
A,36.1
B,25.7
...
S,46
T,26.7

入力データの読込みは以下の通りです。ここで、numは荷物の数 $n$ 、capは箱に詰め込める荷物の重さ合計の上限 $C$ 、idは荷物の識別名(上記の入力データではAからTの文字)、itemsは荷物の集合、wtは各荷物の重さを表します。実務では、荷物に $1,\dots,n$ と都合の良い番号が割り振られている訳ではないので、識別名をそのまま扱えるようにしました。

input = open(sys.argv[1], 'r', encoding='uft-8-sig')
data = [_ for _ in csv.reader(input)]
input.close()
num, cap = int(data[1][0]), float(data[1][1])
items = set()
wt = {}
for line in data[3:]:
    items.add(line[0])
    wt[line[0]] = float(line[1])

変数・制約条件・目的関数の定義は以下の通りです。個人的には、Pythonでは多次元配列を使わず、添字の組をtupleでまとめたものをkeyとして辞書で扱う方が楽だと思っています。数理最適化ソルバーは、同じ定式化であっても変数が現れる順番が変わると異なる結果が出力されることが少なくないので、繰り返し処理では常にsortedで整列して走査順を固定しています。シフトスケジューリング問題ではしばしば複数の最適解が生じるため、実行のたびに異なる結果が出力されて頭を抱えたこともありました。 xsumは $\sum$ を表すPython-MIP独自の関数ですが、他のPython APIでも同じ機能を持つ関数が用意されています。

# initialize model
model = mip.Model(solver_name=CBC)
# variables
x = {}
loop = ((i,id) for i in range(num) for id in items)
for i,j in sorted(loop):
    x[(i,j)]= model.add_var(var_type='B', name=f'x({i},{j})')
y = {}
for i in range(num):
    y[i] = model.add_var(var_type='B', name=f'y({i})')
# constraints
for id in sorted(items):
    model.add_constr(xsum(x[(i,id)] for i in range(num)) == 1, name=f'ITEM({id})')
for i in range(num):
    model.add_constr(xsum(wt[id] * x[(i,id)] for id in sorted(items)) - cap * y[i] <= 0, name=f'BIN({i})')
# objective function
model.objective = minimize(xsum(y[i] for i in range(num)))

ソルバーの実行と計算結果の取得は以下の通りです。ここでは、ソルバーを実行する前に最適化問題をLP形式と呼ばれるフォーマットで書かれたテキストファイルに書き出しています。このファイルは最適化問題が思った通りに定式化できているかどうか確認するのに便利で、よくお世話になっています。ソルバーの実行はoptimizeです。整数計画問題はしばしば最適解を求めるのに非常に長い時間かかるので、max_sedondsで実行時間の上限を設定しています。最適解が求まる前に実行が終了した場合には、その時点での暫定解が得られます。各変数 $x_{ij}$ の値はx[(i,id)].xに格納されます。ここで注意ですが、変数の定義でvar_type='B'と{0,1}の2値であると宣言しているにも関わらず、数理最適化ソルバーが「0.99999…」という値を出力することがあります。このとき、変数値が1.0と等しいかどうかで判定してしまうと、どの箱にも荷物が1つも詰め込まれていないという誤った結果が出力されることがあります。そこで、下記のコードでは値が0.5よりも大きければ1と判断して、各荷物が詰め込まれた箱の番号をsolに格納しています。

# LP format file
model.write('bpp.lp')
# run solver
model.max_seconds = 60
model.optimize()
# get solution
sol = {}
if model.num_solutions:
    obj = model.objective_value
    for i,id in x:
        if x[(i,id)].x > 0.5: sol[id] = i

ナーススケジューリング問題

このようにPython-MIPの使い方を一通り紹介した後で、受講者の皆さんにナーススケジューリング問題の定式化と求解に取り組んでもらいました。ナーススケジューリング問題は多様なバリエーションがある上に現場ごとに特有の追加要件が生じることが多いです。演習では、各シフトの必要人数および各スタッフの勤務日数に関する制約条件のみを課した最も簡単なナーススケジューリング問題を扱うことにしました。対象とするナーススケジューリング問題のシフト表を下図に示します。昼勤(TD)と夜勤(TN)の2種類のシフトがあり、各勤務日に対して必要人数が設定されています。また、各スタッフは計画期間内の各シフトの勤務日数の上下限が設定されています。

詳細の説明は省きますが、整数計画問題に定式化すると以下の通りになります。ナーススケジューリング問題は厳密には最適化問題ではないので、各スタッフの各シフトの勤務日数の上下限に対する違反量の合計を最小化する定式化にしています。演習では、受講者の皆さんに整数計画問題の定式化から取り組んでもらいました。

$$ \begin{array}{ll} \textnormal{最小化} & \displaystyle \sum_{i \in I} \sum_{j \in J} (z_{ij}^- + z_{ij}^+) \cr \textnormal{条件} & \displaystyle \sum_{j \in J_i} x_{ijt} \le 1, \quad i \in I, ; t \in T, \cr & \displaystyle \sum_{i \in I_j} x_{ijt} \ge d_{jt}, \quad j \in J, ; t \in T, \cr & \displaystyle \sum_{t \in T} x_{ijt} + z_{ij}^- \ge L_{ij}, \quad i \in I, ; j \in J, \cr & \displaystyle \sum_{t \in T} x_{ijt} - z_{ij}^- \le U_{ij}, \quad i \in I, ; j \in J, \cr & x_{ijt} \in \{ 0,1 \}, \quad i \in I, ; j \in J, ; t \in T, \cr & z_{ij}^-, z_{ij}^+ \ge 0, \quad i \in I, ; j \in J. \end{array} $$

Pythonは手軽にデータを可視化できるライブラリが充実していますので、今回はstreamlitを用いたシフト表を出力するプログラムを用意しました。計算結果を表示するだけの簡単なプログラムですが数行の実装でデータを可視化できるのはありがたいです。

研修を終えて

最後に、受講者の皆さんからいただいた感想をご紹介します。

一見複雑で解き方の見当もつかない問題をどうスマートに定式化するか考えるのは、パズルを解くような面白さがありました。制約を守りながら何らかの最大化を目指す場面は実務でしばしば登場し、配属されてからより一層、本講義の内容の大切さを実感しています。(K.B.)
大学の講義ではシンプレックス法や内点法を使って手計算で解を求めていたのですが、本講義でライブラリを使った演習ができたことで、業務レベルのより複雑な問題に対応するイメージが湧きました。(高橋駿一)
この最適化モデリング入門は1日目で「最適化の理論が実務にどのように使えるのか」「最適化理論によってどれだけのインパクトをおよぼせるのか」を学び、2日目で実装と演習まで行う非常に贅沢な研修です。最適化理論についてこれほど深く学べる場は他になく、この講義のおかげで、実務で課題を解く手法の一つとして最適化を挙げることができています。また、見返すたびに新たな発見があり、普段の勉強にも活きています。(大塚空来)
一般的に最適化を勉強するツールはたくさんありますが、「最適化を勉強する」ことと「最適化を実装する」ことが別になっており、実装はできても中身は分からない、中身は分かっても実装ができないことがありました。しかし、本講義ではその両方を同時に学ぶことができ、中身を理解するだけではなく自分の手を動かすことができたため、業務で利用するイメージが湧きました。(東将己)

おわりに

本記事では、データ推進室の新人研修BootCampにおいて、私が担当した「最適化モデリング入門」の講義の一部をご紹介しました。数理最適化はデータ分析の後に続く問題解決の出口を担う技術と位置付けられます。データサイエンティストが、個人としても組織としても差別化を図る上で重要な技術であり、受講者の皆さんがそれぞれの配属先でサービスやプロダクトの開発に広く数理最適化を活用してくれることを期待しています。数理最適化を実務に活用するためのノウハウは実務につなげる数理最適化で紹介しましたので、ぜひ、そちらの記事もご参照下さい。

リクルートでは、人生の節目となる「ライフイベント」領域、日常消費の「ライフスタイル」領域に加えて、近年では、業務・経営支援（クラウドを活用したSaaSソリューション）にも力を入れており、多岐に渡る業界で様々なサービスを提供しています。これらのサービスにおける問題解決に関わる中で、数理最適化の新たな活用事例を開拓する多くの機会に恵まれたことを嬉しく思っています。数理最適化を軸とする新しいサービスを立ち上げて世の中を変える、そういう仕事ができるよう今後も努めたいと思います。