Ryotaro Sato

Data Scientist

2024/07/09

線形計画問題・混合整数計画問題をソルバーで解く

自己紹介

データサイエンティストの佐藤です。エムシーデジタルでは、最適化アルゴリズムの設計・実装や、機械学習モデルの構築、バックエンド開発などの業務に主に携わってきました。最適化をはじめとするアルゴリズムの実装・適用に興味があり、プログラミングコンテストサイト AtCoder では Algorithm 部門と Heuristic 部門の両方でレーティング 3000 に到達しています。

概要

本稿は、最適化問題の一種である線形計画問題や混合整数計画問題を紹介し、一般に入手可能なソルバーを利用してこれらの問題を解く Python プログラムの実装・実行方法を説明します。

導入

現実世界における課題は、しばしば「最適化問題」として定式化されます。例えば、物流のルート決定、製造計画の策定、資源の配分など、さまざまな分野で最適化問題が登場します。これらの問題は、変数に関する制約条件と目的関数を持ちます。制約条件を満たしつつ目的関数の値を最小化または最大化するような変数の値の組合せ（解）を求めることが我々の目標です。

こうした最適化問題は、線形計画問題や混合整数計画問題としてモデル化できることがあります。線形計画問題は、目的関数と制約条件がすべて変数に関して線形な式で表せる問題です。一方、混合整数計画問題は、一部の変数に関して「その値が整数でなければならない」という制約条件が追加で課される問題で、より複雑な状況をモデル化できるため、多くの実務的な問題を表現できます。

幸いなことに、線形計画問題や混合整数計画問題を解くための汎用的なソルバーが多数開発されています。高度にチューニングされたソルバーを適切に活用することで、数百や数千個の変数・制約条件からなる問題の厳密解が 1 秒以内に得られてしまうことも少なくありません。本稿では、線形計画問題や混合整数計画問題の基本的な概念から、具体的なソルバーの使用方法・注意点までを解説します。これにより、適切なソルバーを用いて最適化問題を解くプログラムを作れるようになることを目指します。

本稿は以下の内容を扱います。

線形計画問題と混合整数計画問題の概要。
数理的に定式化された課題を線形計画問題や混合整数計画問題として表す方法の例。
線形計画問題や混合整数計画問題をソルバーを利用して解くプログラムの実装。

一方で、例えば以下の内容は扱いません。

Python 言語に関する基本的な解説。
線形計画問題や混合整数計画問題を解くアルゴリズムの説明 ¹²。
「線形計画問題や混合整数計画問題の主問題を単にソルバーに解かせること」の理解に必要な範囲を超える数理的な話題。
実運用に適した最適化プログラムのコード設計。
現実世界における課題を最適化問題に落とし込む方法 ³。

想定される読者は、方程式や不等式などの数式に抵抗がなく、また Python 言語のプログラムの実装にある程度慣れている方です。

なお、紹介しているソフトウェアのライセンスなどの情報は全て執筆時点のものです。今後変化する可能性があるので、実際に試される際は最新の状況を確認するようにしてください。

線形計画問題・混合整数計画問題とは

本稿で扱う最適化問題である線形計画問題や混合整数計画問題について順に簡単に説明します。

線形計画問題

本稿では、 $n$ 次元の連続変数 $x = (x_1, \dots, x_n)$ に関する以下のような最適化問題を 線形計画問題 (linear programming problem, LP) と呼びます：

$$ \begin{align} \text{minimize} \quad & \sum_{i=1}^n c_i x_i \notag \
\text{subject to} \quad & \sum_{j=1}^n A_{ij} x_j \le b_i \quad (i = 1, \dots, m). \notag \end{align} $$

ここで、 $c$ は $n$ 次元実ベクトル、 $A$ は $m \times n$ 実行列、 $b$ は $m$ 次元実ベクトルです。“$\text{subject to}$” 以下の制約条件は、各 $i = 1, \dots, m$ について「$A$ の $i$ 行目による $x$ の重み付き和が $b_i$ 以下である」という形になっています。このような 線形不等式制約 がいくつか（この式では $m$ 個）課された状況で、 $c = (c_1, \dots, c_n)^\top$ による $x$ の重み付き和（目的関数）を最小化せよ、という問題が LP です。

LP において、全ての制約条件を満たすような $x$ の値を 実行可能解 と呼びます。特に目的関数の最小値が存在する場合、それを実現する実行可能解を 最適解 と呼びます。実行可能解をうまく選ぶことで目的関数をいくらでも小さくできる場合は、実行可能解は存在しても最適解が存在しません。逆に制約条件が厳しすぎると、そもそも実行可能解すら存在しないこともあります。

▼Tips
目的関数を最小化する問題のほか、最大化する問題も LP に含まれます。目的関数 $\sum_i c_i x_i$ に関する 最大化 問題を解くのは、目的関数 $\sum_i -c_i x_i$ に関する 最小化 問題を解くのと同じことだからです。そこで、本稿では以後最小化問題と最大化問題を特に区別せずカジュアルに扱うことにします。

▼Tips
制約条件として線形不等式のほか、不等号が等号になったもの（等式制約）を考えることも可能です。例えば、 $\sum_i a_i x_i = b$ という制約を課したい場合、 $\sum_i a_i x_i \le b$ と $\sum_i -a_i x_i \le -b$ の両方を課していると解釈すると納得できるかと思います。そこで、本稿では以後等式制約もカジュアルに扱うことにします。

線形計画問題の例

ここまでで抽象的な LP の定義を示しましたが、ここでは様々な具体的な問題が LP として表現できることを見ていきます。

例（最大流問題）

最大流問題として知られる以下の最適化問題は LP で表現できます。

$n$ 頂点 $m$ 辺の重み付き無向グラフ $G = (V, E)$, $V = \{1, \dots, n \}$, $E = \{1, \dots, m \}$ が与えられる。辺 $i$ ($1 \le i \le m$) は頂点 $u_i$ と $v_i$ を結び、容量 $c_i$ である。2 頂点 $s, t$ $(s \neq t)$ 間の最大フローを求めよ。

各辺 $i$ の $u_i$ から $v_i$ への符号つき流量を変数 $x_i$ とおきます $(i = 1, \dots, m)$。この問題の制約条件は以下のように列挙できます。

各辺の流量は辺の容量を超えてはいけません。よって、各辺 $i$ について、制約 $-c_i \le x_i \le c_i$ が課されます。
$s$ と $t$ を除く各頂点について、その頂点への流入量とその頂点からの流出量は等しくならなければなりません。よって、 $s$, $t$ 以外の各頂点 $j$ について、制約 $\sum_{i : u_i = j} x_i - \sum_{i : v_i = j} x_i = 0$ が課されます。

また、目的関数は頂点 $s$ からの流出量なので、 $\sum_{i : u_i = s} x_i - \sum_{i : v_i = s} x_i$ と表現できます。目的関数が $x_i$ たちの重み付き和で、制約条件は全て線形不等式や線形方程式として表せるので、これは LP の一種であることが確かめられました。

例（天秤とおもり）

以下のような問題も LP を通して解くことができます。これはプログラミングコンテスト AtCoder Heuristic Contest 025 で出題された問題に着想を得たものです。

手元に $n$ 個のおもりがある。それぞれのおもりの重量は $1$ グラム以上であることは分かっているが、詳しい重量は不明である。これらのおもりと一つの天秤を使って、重量の比較を $m$ 回行ったという。それぞれの比較では、いくつかのおもりを天秤の左右の皿に載せ、 左の皿に載せたおもりの重量の和は右の皿に載せたおもりの重量の和より $1$ グラム以上大きかった という報告を得た。これらの報告に矛盾はないだろうか？

それぞれのおもりのグラム単位の重量を変数 $x_i$ ($i = 1, \dots, n$) とおきます。制約条件は以下のように書けます。

それぞれのおもりの重量は $1$ グラム以上なので、制約 $x_i \ge 1$ が課されます。
$j$ 回目の比較 ($j = 1, \dots, m$) で左の皿に載せたおもりの集合を $L_j$, 右の皿に載せたおもりの集合を $R_j$ と表すと、制約 $\sum_{i \in L_j} x_i - \sum_{i \in R_j} x_i \ge 1$ が課されます。

以上の条件のもと、目的関数を任意に設定した LP を解き、実行可能解が存在しなければ報告のいずれかに矛盾があったことになります。逆に実行可能解が存在する場合、少なくとも各おもりの重量がその解に一致していた場合は報告内容どおりの結果になるはずなので、矛盾はありません。

混合整数計画問題

混合整数計画問題 ・混合整数線形計画問題 (mixed integer (linear) programming problem, MIP(MILP)) とは、 LP の一部の変数に対して その値が整数でなければならない という制約条件を加えたものです。

▼Tips
全ての変数が整数制約を課される場合、特に整数計画問題（integer programming problem, IP）と呼ばれることもあります。

整数制約が使えるようになったおかげで、ナップサック問題や巡回セールスマン問題などを含む様々な最適化問題が MIP として表現可能です。以下では MIP の例を見ていきます。

例（ナップサック問題）

以下のような最適化問題（ナップサック問題）は MIP として表現できます。

店で $n$ 種類の商品を扱っている。商品 $i$ ($i = 1, \dots, n$) の重量は $w_i$, 価値は $v_i$ で、各商品の在庫は十分あるものとする。カバンに合計重量 $W$ 以下の範囲で商品を詰め込み、カバン内の商品の価値の合計を最大化せよ。

カバンに詰め込む各商品 $i$ の個数を $x_i$ $(i = 1, \dots, n)$ とおきます。この問題の制約条件を列挙すると以下のようになります。

各 $i$ について $x_i$ は個数を表すので、$x_i \ge 0$ や「 $x_i$ は整数」という制約が課されます。
合計重量の制約として、 $\sum_{i=1}^n w_i x_i \le W$ が課されます。

また、目的関数は $\sum_{i=1}^n v_i x_i$ と表現でき、これを最大化したいです。目的関数が $x_i$ の重み付き和で、制約条件は全て線形不等式と整数制約なので、これは MIP の一種です。

例（グラフの最大独立集合問題）

以下の最適化問題（最大独立集合問題）も MIP として表現できます。

$n$ 頂点 $m$ 辺の無向グラフ $G = (V, E)$, $V = \{1, \dots, n \}$, $E = \{1, \dots, m \}$ が与えられる。辺 $i$ ($1 \le i \le m$) は頂点 $u_i$ と $v_i$ を結ぶ。 $V$ の部分集合 $V'$ であって、 $V'$ に含まれる頂点間に辺が存在しないもので大きさが最大のものを一つ求めよ。

各頂点 $i$ について、変数 $x_i$ を「$V'$ に $i$ が含まれるとき $1$、含まれないとき $0$ をとる」ものとして定めます。この問題の制約条件は以下のように列挙されます。

各 $x_i$ は　$0 \le x_i \le 1$ という線形不等式を満たします。また $x_i$ は整数です。
各辺 $i$ について、 $x_{u_i} + x_{v_i} \le 1$ という制約が課されます。

あとは目的関数を $\sum_i x_i$ としてこれを最大化すればよく、これも MIP （特に IP）の一種です。

ある選択肢を選ぶか選ばないか（この例だと各頂点を採用するかどうか）を $0$ 以上 $1$ 以下の整数変数で表現する手法は、最適化問題の MIP への帰着で多用される印象です。

実装してみよう

この章では LP や MIP として記述される具体的な最適化問題を解く Python プログラムの作成方法を簡単に説明します。

ソルバー・モデラーの活用

LP や MIP は汎用的な枠組みで、具体的な問題を数値的に解くソルバーが数多く開発されており、我々もこれらを利用するのが望ましいです。一般の LP や MIP を解く基礎的なアルゴリズム・数学的な原理自体は広く知られているものの、現実の問題を計算機上で安定して高速に解くために多くの工夫や調整が施された既存のソルバーを性能面で上回るものを自力で作り上げるのは困難です。特に理由がなければ、既存の高性能なソルバーを活用すべきです。ユーザーが問題の定式化さえ行えばそれを解く部分をソルバー任せにできる点は、課題を素早く解決する上で大きなメリットになるでしょう（独自の効率的なアルゴリズムを自分で考えるのが楽しい方にはちょっと物足りないかもしれませんが……）。

また、自作のプログラムとソルバー間のやり取りを自前で実装するのは手間がかかります。この負担を軽減し、より直感的な実装を可能とするライブラリ（本稿ではモデラーと呼ぶことにします）として様々なものが公開・提供されているので、本稿ではこれも利用することにします。

必要なライブラリの導入

本稿では以下 Python を使い、モデラーとして歴史が長く日本語の資料も多数存在する PuLP を採用します。Python のパッケージ管理に pip を利用している場合、以下のコマンドで PuLP が導入できます。

$ pip install PuLP

PuLP は様々なソルバーに対応しており、設定を切り替えるだけで様々なソルバーに問題を解かせることが可能です。ソルバーは原則別途インストールする必要がありますが、PuLP には CBC というソルバーが付属し、特別な設定をしなければこの CBC が自動的に使用されます。本稿では以下 CBC を使うため、追加のソルバーのインストール作業は不要です。

最適化問題を解く Python プログラムの実装例

「混合整数計画問題の例」で例示したナップサック問題を解くためのプログラムは、 PuLP を利用すると以下のように実装できます。本稿で紹介した他の LP・MIP に帰着可能な問題例も、下記のコードを応用することで実装できるはずです。

from pulp import (
    LpBinary,
    LpContinuous,
    LpInteger,
    LpMaximize,
    LpMinimize,
    LpProblem,
    LpStatusOptimal,
    LpVariable,
    PULP_CBC_CMD,
    lpSum,
)


# ナップサック問題を解く
# 入力：
# - weights: weights[i] = (商品 i の重量)
# - values: values[i] = (商品 i の価値)
# - W: (重量上限)
# 出力：
# - solution: solution[i] = (商品 i の個数)
def solve(weights: list[int], values: list[int], W: int):

    n = len(weights)

    problem = LpProblem(
        name="knapsack",  # 問題の名前。任意に定める
        sense=LpMaximize,  # 最大化問題。最小化問題の場合は LpMinimize
    )

    # 変数の定義
    # NOTE: LpVariable.dicts() を使うと変数を辞書で定義できる
    x = [
        LpVariable(
            name=f"x{i}",
            lowBound=0,  # 変数の下限制約。 None で下限なし
            upBound=None,  # 変数の上限制約。 None で上限なし
            cat=LpInteger,  # 変数の種類。 LpContinuous で連続変数。 LpInteger で整数変数。 LpBinary で 0/1 変数
        )
        for i in range(n)
    ]

    # 制約条件の追加
    problem.addConstraint(
        # NOTE: Python 組み込みの `sum()` よりも `pulp.lpSum()` の方が高速
        lpSum(x[i] * weights[i] for i in range(n)) <= W,
    )

    # 目的関数の定義
    problem.setObjective(
        lpSum(x[i] * values[i] for i in range(n)),
    )

    # 最適化
    problem.solve(
        # ソルバーの指定。 `PULP_CBC_CMD` を与えるとソルバーとして CBC を使う
        PULP_CBC_CMD(msg=False),  # msg=False でログを非表示にできる
    )

    # 結果の出力
    if problem.status == LpStatusOptimal:
        # 解は小数として得られるため整数に丸めて利用する
        solution = [int(round(x[i].value())) for i in range(n)]
        return solution
    else:
        return None


print(solve(weights=[2, 3, 5, 7], values=[2, 4, 7, 11], W=10))  # => [0, 1, 0, 1]
print(solve(weights=[2, 3, 5, 7], values=[2, 4, 7, 11], W=100))  # => [1, 0, 0, 14]
print(solve(weights=[2, 3, 5, 7], values=[2, 4, 7, 11], W=1))  # => [0, 0, 0, 0]

その他特記事項

ソルバーの使用にあたって注意すべき事項や細かいトピックについて列挙します。ソルバーの実用に関してはインターネット上に既に様々な資料が存在するので ⁴ 、それらも併せてご参照ください。

数値誤差の考慮

多くの LP・MIP ソルバーでは内部的には各種数値を浮動小数点数として扱っています。この場合、最終的にソルバーが出力した「最適解」に数値誤差が含まれていたり、整数制約を課されていたはずの変数の値がわずかに整数からずれていたり、制約条件がシビアな場合は本来実行可能解が存在するのに解なしという結果になることがあり得ます。誤差が生じることを想定した上でソルバーを利用し、必要に応じて後段の処理等で吸収するのが重要です。なお、ソルバーの設定パラメータを調整することでこの問題が軽減できることもあります。

ソルバーの挙動の再現性

LP・MIP ソルバーが実行される際、数学的に等価な最適化問題が与えられたとしても、変数を宣言する順序や制約を与える順序に依存して挙動が変わり得ます。特に最適解を与えるような変数 $x$ が複数存在する場合、最適化問題の与え方によって得られる解が変わることはごく普通です。

ソルバーの選定

LP や MIP を解くソルバーは多数公開されていますが、同一の最適化問題を解くのに要する時間（やそもそも解けるかどうか）はソルバーによって大きく異なることがあります。ソルバーを適切に選択することで、より高速に問題が解けるようになるかもしれません。

特定の問題群に対する一部のソルバーのベンチマーク結果は ⁵ などのサイトで公開されています。ただしソルバーの優劣は常に明確に判断できるものではなく、特に解かせる問題の性質や用途によっても大きく変わりうる点に注意が必要です。

サイト ⁵ の “LPopt Benchmark (find optimal basic solution)” の検証内容では、本稿執筆時点で無償で商用利用可能な LP ソルバーのなかでは HiGHS が比較的良好な結果を示しているようです（筆者の感覚とも一致します）。例えば SciPy は scipy.optimize.linprog() や scipy.optimize.milp() 関数を通して LP や MIP を解く機能を提供していますが、最近のバージョンの SciPy ではこれらの関数もデフォルトで HiGHS を使用しています。

おわりに

本稿では LP や MIP の定式化と実例を示したのち、具体的な問題をソルバーを利用して解く方法を簡単に説明しました。我々が解きたい問題が適切に LP や MIP に帰着できるならば、ソルバーは効率的に課題を解決する上で強力な武器になるでしょう。

しかしながら、我々が解きたい全ての問題が汎用的なソルバーと「相性が良い」わけではありません。問題に応じて適切な解法を選ぶ能力はソルバーを使いこなす能力以上に重要でしょう（加えて述べると、実業務ではそもそもの課題設定が妥当かを問える能力がしばしば更に重要だと考えています）。新たな問題に相対した際の解法の選択肢にソルバーの活用という手段を加えるための資料として本稿を活用いただけると幸いです。

エムシーデジタルでは、適切な手段を選んで現実世界の課題を解決していくことに興味のある方を募集しています。弊社に興味を持ってくださった方は、ぜひ採用ページもチェックしてみてください。
▶エムシーデジタルの採用ページはこちら