PuLP による数理最適化超入門

●数理最適化とは？

「数理最適化 (Mathematical Optimization)」あるいは「数理計画法 (Mathematical Programing)」は、変数に関する不等式や等式で表される制約の条件下で、目的の関数を最小 (あるいは最大) にする変数の値を求める問題です。

特に、制約条件と目的関数が線形方程式で表される問題を「線形計画法 (Liner Programming, LP)」とか「線形最適化」といいます。この分野には高速な解法アルゴリズム (単体法 : simplex method や内点法 : internal point method など) があって、非常に広範囲な分野で使用されています。

また、変数の値が連続的な実数ではなく、整数値しかとらない場合を「整数計画法 (integer programing)」といいます。特に、変数の値が 0 と 1 しかとらない場合を「0-1 整数計画法 (0-1 integer programming)」といいます。ただし、整数計画法はいわゆる NP 問題になるので、変数の数や変数の領域 (domain) が大きくなると、現実的な時間で厳密解を求めるのは困難になる場合があります。

数理最適化は、問題を解くプログラム (最適化ソルバー) が商用非商用とわず開発されていて、それを使って解くのが一般的です。もちろん商用のソルバーのほうが高性能でしょうが、近年ソルバーの性能は著しく向上していて、非商用 (無料) のソルバーでもかなりの問題を解くことができるようです。無料のソルバーであれば、私達でも気軽に試してみることができますね。

ソルバーは数理モデル (変数、制約条件、目的関数など) を入力とし、それを解いた結果 (最適値) を出力します。数理モデルの記述には標準的なフォーマット (LP ファイルや MPS ファイルなど) があるようですが、それらを自分で書くのは大変なので、モデリング言語を使うのが一般的です。

たとえば、Python のモジュール PuLP を使うと、Python で変数や数式を記述することができます。PuLP は COIN-OR プロジェクトで開発されたモジュールです。同じく COIN-OR プロジェクトで開発された CBC (COIN-OR Branch and Cut) というオープンソースのソルバーが同梱されているので、PuLP をインストールするだけで数理最適化の問題を解くことができます。

PuLP は次のコマンドでインストールすることができます。

python3 -m pip install pulp

$ python3 -m pip freeze | grep PuLP
PuLP==2.7.0

M.Hiroi がインストールしたのは ver 2.7.0 (2023 年 5 月時点) です。PuLP はインストールされているがバージョンが古い場合、次のコマンドでパッケージをアップグレードすることができます。

python3 -m pip install --upgrade pulp

PuLP が扱うことができるのは「混合整数計画法」といって、実数となる変数と整数となる変数が混在している問題です。もちろん、線形計画法や整数計画法 (0-1 整数計画法) でも解くことができます。

本ページでは PuLP (CBC) を使って数理最適化の問題にチャレンジしてみようと思います。なお、M.Hiroi は最適化ソルバーを使うのは初めてなので、何か勘違いや間違いがあると思います。お気づきの点がありましたら、メールでご指摘いただけると助かります。とりあえず、拙作のページ制約論理プログラミング超入門で取り上げた問題を PuLP で解くことができればいいな、と思っています。たいしたことはできませんが、よろしければお付き合いくださいませ。

●参考 URL

GitHub - coin-or/pulp, (本家)
Optimization with PuLP - PuLP 2.7.0 Documentation, (本家, 英語)
PuLPによるモデル作成方法, (PyQ ドキュメント)
組合せ最適化入門：線形計画から整数計画まで - SlideShare, (梅谷俊治さん)
Python を用いた最適化ソルバー Gurobi 入門 - SlideShare, (久保幹雄さん)
整数計画法による定式化入門 (PDF), (藤江哲也さん)
巡回セールスマン問題から始まる数理最適化 - Qiita, (@panchovie さん)

●参考文献

奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991

●PuLP の基礎知識

PuLP の基本的な使い方は簡単です。手順は次のようになります。

prob = LpProblem(name='NoName', sense=LpMinimize) で数理モデルを格納するオブジェクトを生成する

name : 名前
sense : LpMinimize (最小値を求める), LpMaximize (最大値を求める)

var = LpVariable(name, lowBound=None, upBound=None, cat='Continuous') で変数を定義する

name: 変数名

同じ名前の変数があるとソルバー側でエラーになる, PuLP 側では NO CHECK のようだ

lowBound: 変数の下限値, None で -∞
upBound: 変数の上限値, None で ∞
cat: 変数の種類

Continuous, 実数
Integer, 整数
Binary, 0 または 1

prob += 式で目的関数と制約条件を指定する

式には LpVariable() で生成した変数を使用する
制約条件は左辺式演算子右辺式の形式で指定する
使用できる演算子は <=. ==, >= の 3 つだけ (それ以外の演算子は NG)
制約条件はいくつでも指定することができる
演算子を含まない式が目的関数となる
目的関数を複数指定することはできない (省略することは可能)
式には以下の関数が使用できる

lpSum(list), list の要素の合計値を表す
lpDot(list1, list2), list1 と list2 の内積を表す

print(prob) で設定した数理モデルを表示できる

status = prob.solve() で問題を解く (ソルバーを起動して解を求める)

引数には起動するソルバーを渡す
無引数または PULP_CBC_CMD() の返り値を渡すとパッケージに同梱されている CBC が起動される
最適値が得られたならば 1 を返す
文字列で表示したい場合は LpStatus[status] とする
1 は Optimal と表示される
変数 var の値は var.value() で取得する
目的関数の値は prob.objective.value() で取得する

それでは実際に PuLP を使って線形計画法を解いてみましょう。問題は参考文献の線形計画法の例題をお借りしました。

目的関数:
    z = x + y + 1
制約条件:
    3 * x + 5 * y <= 15
    2 * x + y >= 4
    x - y == 1
    x >= 0
    y >= 0

上記の条件で z の最小値を求めます。プログラムは次のようになります。

リスト : 線形計画法 (test1.py)

import pulp

prob = pulp.LpProblem('sample')

# 変数の定義
x = pulp.LpVariable('x', lowBound = 0)
y = pulp.LpVariable('y', lowBound = 0)

# 目的関数
prob += x + y + 1  

# 制約条件
prob += 3 * x + 5 * y <= 15
prob += 2 * x + y >= 4
prob += x - y == 1

print(prob)

# 実行
status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])

# 結果の表示
print("Result")
print("x", x.value())
print("y", y.value())
print("z", prob.objective.value())

LpVariable() で生成した変数を Python の変数 x, y にセットします。あとは x と y を使って prob に数式を追加していくだけです。数式の書き方は Python と同じなので難しいところはないでしょう。

制約条件 x > = 0, y > = 0 は変数の値が非負であることを表しています。今回は LpVariable() の lowBound に 0 にすることで、変数の範囲を [0, ∞] に設定しています。もちろん、prob += x >= 0 のように数式を prob に追加してもかまいません。

あとは solve() でソルバーを実行します。引数には呼び出すソルバーを指定します。無指定または PULP_CBC_CMD() の返り値を渡すと、パッケージに同梱されているソルバー (CBC) を呼び出します。CBC の動作 (オプション指定) は PULP_CBC_CMD() の引数で指定することができます。msg=0 を指定すると CBC が出力するログを表示しません。

実行結果は次のようになります。

$ python test1.py
sample:
MINIMIZE
1*x + 1*y + 1
SUBJECT TO
_C1: 3 x + 5 y <= 15

_C2: 2 x + y >= 4

_C3: x - y = 1

VARIABLES
x Continuous
y Continuous

Status Optimal
Result
x 1.6666667
y 0.66666667
z 3.33333337

実行環境 : Python 3.10.6, PuLP 2.7.0, Ubunts 22.04 LTS (WSL2, Windows 10), Intel Core i5-6200U 2.30GHz

数式を入力するだけで、簡単に解を求めることができました。

●線形計画法の典型的な例題

次は、線形計画法の典型的な例題を 3 つ取り上げます。

[栄養問題]

下表に示す 3 種類の食品 (a, b, c) を使って、2 種類の栄養素 (x, y) の摂取量を満たす一番安い組み合わせを求めてください。

表 : 食品の栄養素
	食品a	食品b	食品c	摂取量
栄養素x	3	1	2	15
栄養素y	1	2	4	10
単価	4	2	5

食品の購入量を変数 a, b, c で表すことにすると、目的関数は次式のようになります。

目的関数:
    z = 4 * a + 2 * b + 5 * c

この式が最小となる各変数の値を求めます。一日に摂取すべき栄養素には下限値があるので、これが制約条件となります。

制約条件:
    3 * a + 1 * b + 2 * c >= 15,   栄養素 x の制約条件
    1 * a + 2 * b + 4 * c >= 10,   栄養素 y の制約条件
    a, b, c >= 0,                  非負条件

あとはこの式を PuLP でプログラムするだけです。

リスト : 栄養問題

import pulp

prob = pulp.LpProblem('sample1')

# 変数の定義
a = pulp.LpVariable('a', lowBound = 0)
b = pulp.LpVariable('b', lowBound = 0)
c = pulp.LpVariable('c', lowBound = 0)

# 目的関数
prob += 4 * a + 2 * b + 5 * c

# 制約条件
prob += 3 * a + 1 * b + 2 * c >= 15
prob += 1 * a + 2 * b + 4 * c >= 10

print(prob)

# 実行
status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])

# 結果の表示
print("Result")
print("a", a.value())
print("b", b.value())
print("c", c.value())
print("z", prob.objective.value())

プログラムは簡単なので説明は割愛させていただきます。実行結果を示します。

sample1:
MINIMIZE
4*a + 2*b + 5*c + 0
SUBJECT TO
_C1: 3 a + b + 2 c >= 15

_C2: a + 2 b + 4 c >= 10

VARIABLES
a Continuous
b Continuous
c Continuous

Status Optimal
Result
a 4.0
b 3.0
c 0.0
z 22.0

食品 a と b だけを購入したほうが安くなりました。もし、食品 b の単価が 2 から 3 に値上がりしたとすると、結果は次のようになります。

sample1:
MINIMIZE
4*a + 3*b + 5*c + 0
SUBJECT TO
_C1: 3 a + b + 2 c >= 15

_C2: a + 2 b + 4 c >= 10

VARIABLES
a Continuous
b Continuous
c Continuous

Status Optimal
Result
a 4.0
b 0.0
c 1.5
z 23.5

この場合は食品 a と c だけを購入したほうが安くなります。ですが、経費は 22.0 から 23.5 に増えてしまいました。

[生産計画問題]

4 つの原料 (w, x, y, z) を使って 4 つの製品 (a, b, c, d) を生産します。製品 1 kg を作るのに必要な原材料、一日に使用できる原料の量 (kg)、製品の利益 (万円) が下表のように定義されています。一日の利益が最大となるように、各製品の生産量を決めてください。

表 : 原料の使用料と利益
	原料w	原料x	原料y	原料z	利益
製品a	2	1	0	0	5
製品b	0	2	1	0	3
製品c	0	0	1	2	2
製品d	1	0	0	2	4
使用量	4	8	6	10

各製品の生産量を変数 a, b, c, d で表すことにすると、目的関数は次式のようになります。

目的関数:
    g = 5 * a + 3 * b + 2 * c + 4 * d

この式が最大となる各変数の値を求めます。一日に使用できる原料は上限値があるので、これが制約条件となります。

制約条件:
    2 * a + 1 * d <= 4,   原料 w の制約
    1 * a + 2 * b <= 8,   原料 x の制約
    1 * b + 1 * c <= 6,   原料 y の制約
    2 * c + 2 * d <= 10,  原料 z の制約
    a, b, c, d >= 0,      非負条件

あとはこの式を PuLP でプログラムするだけです。

リスト : 生産計画問題

import pulp

prob = pulp.LpProblem('sample2', sense = pulp.LpMaximize)

# 変数の定義
a = pulp.LpVariable('a', lowBound = 0)
b = pulp.LpVariable('b', lowBound = 0)
c = pulp.LpVariable('c', lowBound = 0)
d = pulp.LpVariable('d', lowBound = 0)

# 目的関数
prob += 5 * a + 3 * b + 2 * c + 4 * d

# 制約条件
prob += 2 * a + d <= 4
prob += a + 2 * b <= 8
prob += b + c <= 6
prob += 2 * c + 2 * d <= 10

print(prob)

# 実行
status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])

# 結果の表示
print("Result")
print("a", a.value())
print("b", b.value())
print("c", c.value())
print("d", d.value())
print("g", prob.objective.value())

実行結果を示します。

sample2:
MAXIMIZE
5*a + 3*b + 2*c + 4*d + 0
SUBJECT TO
_C1: 2 a + d <= 4

_C2: a + 2 b <= 8

_C3: b + c <= 6

_C4: 2 c + 2 d <= 10

VARIABLES
a Continuous
b Continuous
c Continuous
d Continuous

Status Optimal
Result
a 0.0
b 4.0
c 1.0
d 4.0
g 30.0

製品 a の生産量が 0 になってしまいました。そこで、変数の制約を a, b, c, d >= 1 に変更して実行してみましょう。

sample2:
MAXIMIZE
5*a + 3*b + 2*c + 4*d + 0
SUBJECT TO
_C1: 2 a + d <= 4

_C2: a + 2 b <= 8

_C3: b + c <= 6

_C4: 2 c + 2 d <= 10

VARIABLES
1 <= a Continuous
1 <= b Continuous
1 <= c Continuous
1 <= d Continuous

Status Optimal
Result
a 1.0
b 3.5
c 2.5
d 2.0
g 28.5

すべての製品を生産しますが、利益 g は 30.0 から 28.5 に低下しました。製品 a の利益が高いと違った結果になるでしょう。興味のある方は試してみてください。

[輸送問題]

工場 (x, y) から商品を店 (a, b, c) に配送します。供給量、需要量、輸送コストが下表で与えられているとき、総輸送コストが最小となる配送の仕方を求めてください。

表 : 輸送コスト
	店 a	店 b	店 c	供給量
工場 x	10	6	16	8
工場 y	8	8	4	16
需要量	12	4	8

工場 x から店に配送する商品の量を変数 xa, xb, xc とし、工場 y から店に配送する商品の量を変数 ya, yb, yc とします。目的関数は次式のようになります。

目的関数:
    z = 10 * xa + 6 * xb + 16 * xc + 8 * ya + 8 * yb + 4 * yc

この式が最小となる各変数の値を求めます。工場の供給量と店の需要には制限があるので、これが制約条件となります。

制約条件:
    xa + xb + xc == 8,  工場 x の供給量
    ya + yb + yc == 16, 工場 y の供給量
    xa + ya == 12,      店 a の需要
    xb + yb == 4,       店 b の需要
    xc + yc == 8,       店 c の需要
    xa, xb, xc >= 0,    非負条件
    ya, yb, yc >= 0

あとはこの式を PuLP でプログラムするだけです。

リスト : 輸送問題

import pulp

prob = pulp.LpProblem('sample3')

# 変数の定義
xa = pulp.LpVariable('xa', lowBound = 0)
xb = pulp.LpVariable('xb', lowBound = 0)
xc = pulp.LpVariable('xc', lowBound = 0)
ya = pulp.LpVariable('ya', lowBound = 0)
yb = pulp.LpVariable('yb', lowBound = 0)
yc = pulp.LpVariable('yc', lowBound = 0)

# 目的関数
prob += 10 * xa + 6 * xb + 16 * xc + 8 * ya + 8 * yb + 4 * yc 

# 制約条件
prob += xa + xb + xc == 8
prob += ya + yb + yc == 16
prob += xa + ya == 12
prob += xb + yb == 4
prob += xc + yc == 8

print(prob)

# 実行
status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])

# 結果の表示
print("Result")
print("xa", xa.value())
print("xb", xb.value())
print("xc", xc.value())
print("ya", ya.value())
print("yb", yb.value())
print("yc", yc.value())

print("z", prob.objective.value())

実行結果を示します。

sample3:
MINIMIZE
10*xa + 6*xb + 16*xc + 8*ya + 8*yb + 4*yc + 0
SUBJECT TO
_C1: xa + xb + xc = 8

_C2: ya + yb + yc = 16

_C3: xa + ya = 12

_C4: xb + yb = 4

_C5: xc + yc = 8

VARIABLES
xa Continuous
xb Continuous
xc Continuous
ya Continuous
yb Continuous
yc Continuous

Status Optimal
Result
xa 4.0
xb 4.0
xc 0.0
ya 8.0
yb 0.0
yc 8.0
z 160.0

この例では、工場 x, y ともに半分ずつの量を店に配送するのが一番安くなりました。輸送コストが変化すると、当然ですが違った結果になります。興味のある方はいろいろ試してみてください。

●鶴亀算

PuLP を使って連立方程式を解くこともできます。簡単な例として「鶴亀算」を解いてみましょう。

鶴と亀、合わせて 100 匹いる。足の合計が 272 本のとき、鶴と亀はそれぞれ何匹ずついるか。
鶴と亀とトンボが合わせて 10 匹いる。足の合計が 38 本で羽の合計が 14 枚であるとき、鶴と亀とトンボはそれぞれ何匹ずついるか。(トンボの足は 6 本で羽は 4 枚)
鶏と犬とタコ、合わせて 24 匹が台所にいる。足の合計が 102 本のとき、鶏、犬、タコはそれぞれ何匹ずついるか。

問題 1 は次の連立方程式を解けば求めることができます。

x + y = 100
2x + 4y = 272

問題 2 は次の連立方程式を解けば求めることができます。

 x +  y +  z = 10
2x + 4y + 6z = 38
2x +      4z = 14

問題 3 は答えを 1 つに決めることはできませんが、適切な目的関数と制約条件を設定すれば、それを満たす解を求めることができます。

プログラムは次のようになります。

リスト : 鶴亀算

import pulp

# 変数の定義
x = pulp.LpVariable('x', cat='Integer')
y = pulp.LpVariable('y', cat='Integer')
z = pulp.LpVariable('z', cat='Integer')

# 問題１
p1 = pulp.LpProblem('turukame1')
p1 += x + y == 100
p1 += 2 * x + 4 * y == 272
print(p1)
st1 = p1.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[st1])
print("x", x.value())
print("y", y.value())

# 問題２
p2 = pulp.LpProblem('turukame2')
p2 += x + y + z == 10
p2 += 2 * x + 4 * y + 6 * z == 38
p2 += 2 * x + 4 * z == 14
print(p2)
st2 = p2.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[st2])
print("x", x.value())
print("y", y.value())
print("z", z.value())

# 問題３
p3 = pulp.LpProblem('turukame3', sense=pulp.LpMaximize)
p3 += z
p3 += x + y + z == 24
p3 += 2 * x + 4 * y + 8 * z == 102
p3 += x >= 1
p3 += y >= 1
p3 += z >= 1
print(p3)
st3 = p3.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[st3])
print("x", x.value())
print("y", y.value())
print("z", z.value())

変数は整数 (cat='Integer') に設定します。問題１と２の場合、目的関数の設定は不要です。問題３の場合、変数の値を正、目的関数を z とし、z が最大となる解を求めます。結果は次のようになりました。

turukame1:
MINIMIZE
None
SUBJECT TO
_C1: x + y = 100

_C2: 2 x + 4 y = 272

VARIABLES
x free Integer
y free Integer

Status Optimal
x 64.0
y 36.0
turukame2:
MINIMIZE
None
SUBJECT TO
_C1: x + y + z = 10

_C2: 2 x + 4 y + 6 z = 38

_C3: 2 x + 4 z = 14

VARIABLES
x free Integer
y free Integer
z free Integer

Status Optimal
x 3.0
y 5.0
z 2.0
turukame3:
MAXIMIZE
1*z + 0
SUBJECT TO
_C1: x + y + z = 24

_C2: 2 x + 4 y + 8 z = 102

_C3: x >= 1

_C4: y >= 1

_C5: z >= 1

VARIABLES
x free Integer
y free Integer
z free Integer

Status Optimal
x 13.0
y 3.0
z 8.0

一瞬で答えを求めることができます。

●部分和問題

部分和問題 (subset-sum problem) は、要素が数値の集合 S において、要素の総和が M となる部分集合があるか判定する問題です。たとえば、集合 {2, 3, 5, 8} の場合、総和が 10 となる部分集合は {2, 3, 5} と {2, 8} がありますが、14 となる部分集合はありません。部分集合の総数は、要素数を n とすると 2ⁿ 個になるので、n が大きくなるとナイーブな方法では時間がかかってしまいます。

実際には、動的計画法を使って現実的な時間で部分和問題を解くことができるといわれています。最適化ソルバーでも高速に解くことができるか試してみましょう。

部分和問題は集合 X の要素を x_i とし、その係数を a_iすると、次の等式を満たすか判定する問題になります。

\( \displaystyle \sum_{i=1}^{n} a_i x_i = M, \quad (a_i = 0 \ or \ 1) \)

係数 a_i の値が 0 か 1 かを決める問題になります。このような問題を「0-1 整数計画問題」といいます。PuLP を使うと、とても簡単にプログラムすることができます。次のリストを見てください。

リスト : 部分和問題 (subsetsum.py)

import pulp

nums = [  1,   2,   3,   5,   8,   13,   21,   34,   55,    89,
        144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946]

def solver(xs, n):
    prob = pulp.LpProblem('subset-sum')

    # 変数
    vs = [pulp.LpVariable('x{}'.format(x), cat='Binary') for x in range(len(xs))]
    # 制約条件
    prob += pulp.lpDot(xs, vs) == n

    print(prob)
    status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
    print("Status", pulp.LpStatus[status])
    if status == 1:
        print([xs[i] for i in range(len(xs)) if vs[i].value() == 1])

関数 solver の引数 xs が数の集合で、引数 n が求める部分集合の総和を表します。PuLP の変数は xs の要素と同じ数だけ生成して、リスト vs に格納します。つまり、xs[i] に対応する変数が vs[i] になります。これは Python の内包表記を使えば簡単ですね。このとき、cat に 'Binary' を指定して変数の値を 0 or 1 にします。

制約条件の指定も簡単です。xs と vs の内積を計算し、その値が n になればいいので、prob += pulp.lpDot(xs, vs) == n とするだけです。総和が n となる部分集合があるか判定するだけでよければ、目的関数を設定する必要はありません。解があれば status は Optimal (1) になり、解がない場合は Infeasible (-1) になります。

リスト nums はフィボナッチ数列になっています。要素の総和を M とすると、1 から M までの整数は、要素を組み合わせて必ず作ることができます。これはフィボナッチ数列の面白い特徴です。これをテストに使います。それでは実際に試してみましょう。

>>> from subsetsum import *
>>> solver(nums, sum(nums) // 2)
subset-sum:
MINIMIZE
None
SUBJECT TO
_C1: x0 + 2 x1 + 144 x10 + 233 x11 + 377 x12 + 610 x13 + 987 x14 + 1597 x15
 + 2584 x16 + 4181 x17 + 6765 x18 + 10946 x19 + 3 x2 + 5 x3 + 8 x4 + 13 x5
 + 21 x6 + 34 x7 + 55 x8 + 89 x9 = 14327

VARIABLES
0 <= x0 <= 1 Integer
0 <= x1 <= 1 Integer
0 <= x10 <= 1 Integer
0 <= x11 <= 1 Integer
0 <= x12 <= 1 Integer
0 <= x13 <= 1 Integer
0 <= x14 <= 1 Integer
0 <= x15 <= 1 Integer
0 <= x16 <= 1 Integer
0 <= x17 <= 1 Integer
0 <= x18 <= 1 Integer
0 <= x19 <= 1 Integer
0 <= x2 <= 1 Integer
0 <= x3 <= 1 Integer
0 <= x4 <= 1 Integer
0 <= x5 <= 1 Integer
0 <= x6 <= 1 Integer
0 <= x7 <= 1 Integer
0 <= x8 <= 1 Integer
0 <= x9 <= 1 Integer

Status Optimal
[1, 8, 34, 144, 610, 2584, 10946]
>>> solver(nums, sum(nums) + 1)
subset-sum:
MINIMIZE
None
SUBJECT TO
_C1: x0 + 2 x1 + 144 x10 + 233 x11 + 377 x12 + 610 x13 + 987 x14 + 1597 x15
 + 2584 x16 + 4181 x17 + 6765 x18 + 10946 x19 + 3 x2 + 5 x3 + 8 x4 + 13 x5
 + 21 x6 + 34 x7 + 55 x8 + 89 x9 = 28656

VARIABLES
・・・省略・・・

Status Infeasible

変数の個数が少ないので一瞬で答えを求めることができます。ところで、解が無くてもそれに近い値がほしい場合があるでしょう。この場合、目的関数を lpDot(xs, vs)、制約条件を lpDot(xs, vs) <= n として最大値を探索すれば、n 以下でもっとも近い値を求めることができます。プログラムは次のようになります。

リスト : 部分和問題 (2)

def solver1(xs, n):
    prob = pulp.LpProblem('subset-sum', sense = pulp.LpMaximize)

    # 変数
    vs = [pulp.LpVariable('x{}'.format(x), cat='Binary') for x in range(len(xs))]
    # 目的関数
    prob += pulp.lpDot(xs, vs)
    # 制約条件
    prob += pulp.lpDot(xs, vs) <= n

    print(prob)

    status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
    print("Status", pulp.LpStatus[status])
    if status == 1:
        print([xs[i] for i in range(len(xs)) if vs[i].value() == 1])
        print(prob.objective.value())

それでは実際に試してみましょう。

>>> solver1(nums, 10945)
subset-sum:
MAXIMIZE
1*x0 + 2*x1 + 144*x10 + 233*x11 + 377*x12 + 610*x13 + 987*x14 + 1597*x15 + 2584*x16 + 4181*x17 + 6765*x18 + 
10946*x19 + 3*x2 + 5*x3 + 8*x4 + 13*x5 + 21*x6 + 34*x7 + 55*x8 + 89*x9 + 0
SUBJECT TO
_C1: x0 + 2 x1 + 144 x10 + 233 x11 + 377 x12 + 610 x13 + 987 x14 + 1597 x15
 + 2584 x16 + 4181 x17 + 6765 x18 + 10946 x19 + 3 x2 + 5 x3 + 8 x4 + 13 x5
 + 21 x6 + 34 x7 + 55 x8 + 89 x9 <= 10945

VARIABLES
・・・省略・・・

Status Optimal
[1, 3, 8, 21, 55, 144, 377, 987, 2584, 6765]
10945.0
>>> solver1(nums, 30000)
subset-sum:
MAXIMIZE
1*x0 + 2*x1 + 144*x10 + 233*x11 + 377*x12 + 610*x13 + 987*x14 + 1597*x15 + 2584*x16 + 4181*x17 + 6765*x18 + 
10946*x19 + 3*x2 + 5*x3 + 8*x4 + 13*x5 + 21*x6 + 34*x7 + 55*x8 + 89*x9 + 0
SUBJECT TO
_C1: x0 + 2 x1 + 144 x10 + 233 x11 + 377 x12 + 610 x13 + 987 x14 + 1597 x15
 + 2584 x16 + 4181 x17 + 6765 x18 + 10946 x19 + 3 x2 + 5 x3 + 8 x4 + 13 x5
 + 21 x6 + 34 x7 + 55 x8 + 89 x9 <= 30000

VARIABLES
・・・省略・・・

Status Optimal
[1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946]
28655.0

引数 n と求めた解が等しければ、n と等しくなる部分集合を求めることができました。そうでなければ、n と等しい部分集合はありませんが、それに近い値と部分集合を得ることができます。

●ナップザック問題

次は「ナップザック問題」を取り上げます。ナップザック (knapsack) とは辞書を引いてみると、ランドセルのような背中にせおう四角形の袋や箱のことを意味します。ここでは物を入れる袋と簡単に考えてください。

ここで、ナップザックの中に品物を詰め込むことを考えてみます。一つのナップザックと複数の品物が与えられたとき、袋に詰めた品物の合計金額が最大になるような選び方を求めることが「ナップザック問題」です。ナップザック問題にはバリエーションがあって、同じ品物をいくつも入れて良い場合と、一つしか入れてはいけない場合があります。後者の場合を「0-1 ナップザック問題」といいます。

ナップザック問題は、部分和問題と同様に NP 問題になります。これは厳密に解を求めようとすると、全ての場合について総当たりで調べるしか方法がなく、データ数が多くなると時間がかかるため、現実的な時間では解答を出すことができないというものです。品物の詰め方が難問の一つ、といわれてもピンとこないと思いますが、ナップザック問題は品物の種類が増えるにしたがって、その組み合わせ方が爆発的に増えるのです。

ところが、幸いなことにナップザック問題は実用的には解決済みの問題と考えられています。とくに有名なのが「動的計画法」を用いた解法です。ナップザックと品物の大きさを整数値に限定すれば、動的計画法を用いることで厳密解を求めることが可能です。興味のある方は拙作のページ Algorithms with Python 動的計画法をお読みください。

まず最初に、0-1 ナップザック問題を解いてみましょう。0-1 ナップザック問題を数式で表すと次のようになります。

W はナップザックの大きさで、V の最大値を求める

\(\begin{array}{l} \displaystyle \sum_{i=1}^{n} a_i w_i \leq W \\ \displaystyle \sum_{i=1}^{n} a_i p_i = V \end{array}\)

\(w_i\) : 重さ
\(p_i\) : 金額
\(a_i\) : 0 or 1
\(n\) : 品物の個数

基本的には、部分和問題と同じようなプログラムになりますが、ナップザックに入る範囲で、最大の金額になるような入れ方を求める必要があります。それでは実際に簡単な問題を解いてみましょう。

[問題１]

下表に示す品物をサイズ 15 のナップザックに入れるとき、金額が最大となる入れ方を求めてください。

品物	金額	サイズ
A	4	3
B	5	4
C	6	5
D	8	7
E	10	9

出典 : Coprisによる制約プログラミング入門, (田村直之さん)

プログラムは次のようになります。

リスト : 問題１の解法 (knapsack1.py)

import pulp

prob = pulp.LpProblem('knapsack1', sense = pulp.LpMaximize)

# 品物
goods = ['a', 'b', 'c', 'd', 'e']
# 金額
price = [4, 5, 6, 8, 10]
# 重さ
weight = [3, 4, 5, 7, 9]

# 変数の定義
xs = [pulp.LpVariable('{}'.format(x), cat='Binary') for x in goods]
# 目的関数
prob += pulp.lpDot(price, xs)
# 制約条件
prob += pulp.lpDot(weight, xs) <= 15

print(prob)

status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
print([x.value() for x in xs])
print(prob.objective.value())

knapsack1:
MAXIMIZE
4*a + 5*b + 6*c + 8*d + 10*e + 0
SUBJECT TO
_C1: 3 a + 4 b + 5 c + 7 d + 9 e <= 15

VARIABLES
0 <= a <= 1 Integer
0 <= b <= 1 Integer
0 <= c <= 1 Integer
0 <= d <= 1 Integer
0 <= e <= 1 Integer

Status Optimal
[1.0, 0.0, 1.0, 1.0, 0.0]
18.0

金額の最大値は 18 で、選択した品物は A, C, D の 3 つです。

次は同じ品物をいくつ選んでもよい問題を解いてみましょう。この場合、変数は Binary ではなく Integer になります。

[問題２]

下表に示す品物をサイズ 10 のナップザックに入れるとき、金額が最大となる入れ方を求めてください。なお、同じ品物を何個入れてもかまいません。

品物	金額	サイズ
A	6	4
B	4	3
C	1	1

リスト : 問題２の解法 (knapsack2.py)

import pulp

prob = pulp.LpProblem('knapsack2', sense = pulp.LpMaximize)

# 品物
goods = ['a', 'b', 'c']
# 金額
price = [6, 4, 1]
# 重さ
weight = [4, 3, 1]

# 変数の定義
xs = [pulp.LpVariable('{}'.format(x), cat='Integer', lowBound=0) for x in goods]
# 目的関数
prob += pulp.lpDot(price, xs)
# 制約条件
prob += pulp.lpDot(weight, xs) <= 10

print(prob)

status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
print([x.value() for x in xs])
print(prob.objective.value())

knapsack2:
MAXIMIZE
6*a + 4*b + 1*c + 0
SUBJECT TO
_C1: 4 a + 3 b + c <= 10

VARIABLES
0 <= a Integer
0 <= b Integer
0 <= c Integer

Status Optimal
[2.0, 0.0, 2.0]
14.0

金額の最大値は 14 で、A を 2 個、C を 2 個選びます。なお、A を 1 個、B を 2 個選んでも金額は 14 になります。

もう一つ、簡単な例を示しましょう。

[問題３]

下表に示す品物をサイズ w のナップザックに入れるとき、金額が最大となる入れ方を求めてください。なお、同じ品物を何個入れてもかまいません。

品物	金額	サイズ
A	91	3
B	120	4
C	610	20
D	930	30

リスト : 問題３の解法 (knapsack3.py)

import pulp

def solver(w):
    # 品物
    goods = ['a', 'b', 'c', 'd']
    # 金額
    price = [91, 120, 610, 930]
    # 重さ
    weight = [3, 4, 20, 30]

    prob = pulp.LpProblem('knapsack3', sense = pulp.LpMaximize)
    # 変数の定義
    xs = [pulp.LpVariable('{}'.format(x), cat='Integer', lowBound=0) for x in goods]
    # 目的関数
    prob += pulp.lpDot(price, xs)
    # 制約条件
    prob += pulp.lpDot(weight, xs) <= w

    print(prob)

    status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
    print("Status", pulp.LpStatus[status])
    print([x.value() for x in xs])
    print(prob.objective.value())

>>> from knapsack3 import *
>>> solver(500)
knapsack3:
MAXIMIZE
91*a + 120*b + 610*c + 930*d + 0
SUBJECT TO
_C1: 3 a + 4 b + 20 c + 30 d <= 500

VARIABLES
0 <= a Integer
0 <= b Integer
0 <= c Integer
0 <= d Integer

Status Optimal
[0.0, 0.0, 1.0, 16.0]
15490.0
>>> solver(1000)
knapsack3:
MAXIMIZE
91*a + 120*b + 610*c + 930*d + 0
SUBJECT TO
_C1: 3 a + 4 b + 20 c + 30 d <= 1000

VARIABLES
0 <= a Integer
0 <= b Integer
0 <= c Integer
0 <= d Integer

Status Optimal
[2.0, 1.0, 0.0, 33.0]
30992.0

サイズ w が大きな値でも一瞬で解を求めることができました。SWI-Prolog のライブラリ clpfd を使って問題３を解く場合、サイズ w を増やすと少し時間がかかるようになりますが、PuLP ではあまり変わらないようです。興味のある方はいろいろ試してみてください。

初版 2019 年 6 月
改訂 2023 年 5 月 27 日

Python3 Programming