遅ればせながら PuLP の使い方について勉強しました。
案件で使えるレベルとは程遠いですが、まぁ経験が大事ということで。
基本 (Basis)
最初に PuLP の使い方を 勉強します。
大まかにいうと モデルを作って解決するだけです。 モデルには目的関数や制約条件関数を細かく設定できます。
- LpVariable
- 一意な変数名 と 定義域、変数の型 を指定します
- LpProblem
-
解決すべき問題の種類を指定することでモデルを定義します。
モデルの制約条件がある場合、条件式を加算代入で注ぎ足していきます。 条件式は前述した LpVariable で定義した変数を使って作成します。
問題の種類は LpMaximize (最大), LpMinimize (最小) のいずれかです。
最初は 適当な線形の目的関数に 適当な制約条件を追加し、 その前後で最大値、最小値 がどのように変化するか見ていきましょう。
%autosave 0
%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pulp
# 変数の定義には 「名前(一意)」 「定義域(最小,最大)」 「型」 を指定する
# 連続した値(浮動小数点数) は 'Continuous' (default)
a = pulp.LpVariable('a', -3, 3, cat='Continuous')
# 整数は 'Integer'
b = pulp.LpVariable('b', -5, 5, cat='Integer')
# 真偽値 は 'Binary'
c = pulp.LpVariable('c', cat='Binary')
# モデルの初期化時に 「名前」 と 「最大か最小か(default:最小)」 を指定する
# model_max で最大値, model_min で最小値を 平行して求めていく
model_max = pulp.LpProblem('Maximizing', pulp.LpMaximize)
model_min = pulp.LpProblem('Maximizing', pulp.LpMinimize)
def objective_function(a, b):
return -2 * a + 3 * b + 1
# モデルに目的関数と制約条件を追加していく
# 適当な目的関数
model_max += objective_function(a, b)
model_min += objective_function(a, b)
# model.solve() で 求解 (Optimal がでれば OK)
assert pulp.LpStatus[model_max.solve()] == 'Optimal'
solved_max1 = {'x': a.value(), 'y': b.value(), 'z': objective_function(a.value(), b.value())}
assert pulp.LpStatus[model_min.solve()] == 'Optimal'
solved_min1 = {'x': a.value(), 'y': b.value(), 'z': objective_function(a.value(), b.value())}
solved_max1, solved_min1
# 適当な制約条件を追加
model_max += a >= 1
model_min += a >= 1
assert pulp.LpStatus[model_max.solve()] == 'Optimal'
solved_max2 = {'x': a.value(), 'y': b.value(), 'z': -2 * a.value() + 3 * b.value() + 1}
# 制約条件を加えても変わらないケースもある
assert pulp.LpStatus[model_min.solve()] == 'Optimal'
solved_min2 = {'x': a.value(), 'y': b.value(), 'z': -2 * a.value() + 3 * b.value() + 1}
solved_max2, solved_min2
# 適当な制約条件を追加(2)
model_max += (a >= b) == c
model_min += (a >= b) == c
assert pulp.LpStatus[model_max.solve()] == 'Optimal'
solved_max3 = {'x': a.value(), 'y': b.value(), 'z': objective_function(a.value(), b.value())}
assert pulp.LpStatus[model_min.solve()] == 'Optimal'
solved_min3 = {'x': a.value(), 'y': b.value(), 'z': objective_function(a.value(), b.value())}
solved_max3, solved_min3
# 座標をプロットしてみる
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
x = np.arange(a.lowBound, a.upBound + 1)
y = np.arange(b.lowBound, b.upBound + 1)
xs, ys = np.meshgrid(x, y)
zs = -2 * xs + 3 * ys + 1
ax.plot_wireframe(xs, ys, zs)
# 最大値(赤)
ax.plot(
[solved_max1['x'], solved_max2['x'], solved_max3['x']],
[solved_max1['y'], solved_max2['y'], solved_max3['y']],
[solved_max1['z'], solved_max2['z'], solved_max3['z']]
, 'o', color='red')
# 最小値(緑)
ax.plot(
[solved_min1['x'], solved_min2['x'], solved_min3['x']],
[solved_min1['y'], solved_min2['y'], solved_min3['y']],
[solved_min1['z'], solved_min2['z'], solved_min3['z']]
, 'o', color='green')
座標プロットのところはちょっとわかりにくかったですかね?
- objective_function は x, y 座標を受け取って係数をかけて返却する単純な 目的関数です
- 格子状の平面が 目的関数です
- この格子は私が自力で座標を指定して描画してるもので勝手に描画されたりはしません
- 格子は matplotlib の plot_wireframe という関数で描画できます。
- 最初は プロットが 変数自体の定義域の制限により格子上下端に配置されていて、
制約条件を 与えると動ける範囲で 最大、あるいは最小の点に移動します。
- 意味のない制約条件を与えると、座標が変わらないこともあります。
- 最小値はの点は1つカブらせたので2つしかないです
- 意味のない制約条件を与えると、座標が変わらないこともあります。
線形計画法 (Linear programming)
続いて、応用としてつとむさんが作成した以下の問題を解いてみます。斎藤さんだぞ。 (最適化におけるPython)
材料AとBから合成できる化学製品XとYをたくさん作成したい。
Xを1kg作るのに、Aが1kg、Bが3kg必要である。
Yを1kg作るのに、Aが2kg、Bが1kg必要である。
また、XもYも1kg当りの価格は100円である。
材料Aは16kg、Bは18kgしかないときに、XとYの価格の合計が最大になるようにするには、 XとYをどれだけ作成すればよいか求めよ。
PuLP を使って解いてみます。
%autosave 0
%matplotlib inline
import pulp
# 金額が最大になるようにしたいので LpMaximize
model = pulp.LpProblem('Maximizing', pulp.LpMaximize)
x = pulp.LpVariable('x', lowBound=0, cat='Integer')
y = pulp.LpVariable('y', lowBound=0, cat='Integer')
# 目的関数
model += 100 * x + 100 * y
# Aの制約条件: Xなら1kg, Yなら2kg を必要とし 16kg まで利用可能
model += x + 2 * y <= 16
# Bの制約条件: Xなら3kg, Yなら1kg を必要とし 18kg まで利用可能
model += 3 * x + y <= 18
assert pulp.LpStatus[model.solve()] == 'Optimal'
x.value(), y.value()
とても簡単ですね。
PuLP を使わずに解くとこんな感じ
%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from sympy import symbols, Eq, Le, init_printing, solve
init_printing()
a, b, x, y, k = symbols('a b x y k')
変数の定義
- a: A の使用量合計
- b: B の使用量合計
- x: X の作成量
- y: Y の作成量
- k: 価格合計
# x, y を使って 目的関数 k を定義
Eq(k, 100 * x + 100 * y)
# x, y を使って a を定義
Eq(a, x + 2 * y)
# x, y を使って b を定義
Eq(b, 3 * x + y)
# a の制約条件(a <= 16) に代入して
x + 2 * y <= 16
# b の制約条件(b <= 18) に代入して
3 * x + y <= 18
# 交点の座標
solve([x + 2 * y - 16, 3 * x + y - 18])
だいぶ情報は集まってきましたが、この時点では効率の良い(材料をちょうど使い切る)組み合わせを求めただけで金額が最大になるとは言えませんね。
ここで 目的関数 k を使います。しかし k は 2変数関数なので 普通に考えると 3次元のグラフになってしまいそうです。
そこで k を定数として y の式に変形することで y = -x + k/100
のように 傾き:-1
, 切片:k/100
の 1次関数とみなすことができます。
制約条件のグラフを使って 逆像法的に 以下の条件を満たす x, y を探してみます。
- k/100 (切片) が最大になる
- 制約条件と 共有点 を持つ
# 制約条件を図示してみる (塗りつぶしが取りうる範囲)
XLIM = YLIM = 10
fig = plt.figure()
ax = fig.add_subplot(111)
ax.grid(color='gray', linestyle='-')
ax.set_xlabel('X')
ax.set_ylabel('Y')
xs = np.arange(0, XLIM + 1, 1)
ys1 = xs * -1/2 + 8
ys2 = xs * -3 + 18
p1 = ax.plot(xs, ys1, color='red')
p2 = ax.plot(xs, ys2, color='blue')
ax.fill_between(xs, ys1, where=xs <= 4, color="orange")
ax.fill_between(xs, ys2, where= xs >= 4, color="orange")
ax.set_xlim([0, XLIM])
ax.set_ylim([0, YLIM])
ims = [
ax.plot(xs, -1 * xs + seg, 'g')
for seg in range(11)
]
ax.legend(['[a] constraint cond boundary', '[b] constraint cond boundary', 'objective function'])
ani = animation.ArtistAnimation(fig, ims)
plt.show()
k/100
が 制約条件を満たす範囲(塗りつぶし領域と共有点を持ち)、切片が最大となる点は グラフを見て分かる通り 10
で、その共有点は、やはり (4, 6)
でした。
今回は最大になる組み合わせを聞かれているだけなので重要なのは価格の比率だけで、実際にいくらになるかを求める必要はありませんが、 k=10 を 100倍して 1000円
が最大価格と言えそうです。
目的関数に代入してもわかります。
本当はアニメーションで 目的関数が上下移動するようにしたんですが、 静止画でしか保存できないようなので 制約条件下で切片が最大になったときの目的関数を表示しています。
魔方陣 (Magic square)
最後に魔方陣を解きます。 魔方陣のルールは Wikipedia を参照。
魔方陣(まほうじん、英:Magic square)とは、n×n 個の正方形の方陣に数字を配置し、縦・横・対角線のいずれの列についても、その列の数字の合計が同じになるもののことである。 特に1から方陣のマスの総数 n2 までの数字を1つずつ過不足なく使ったものを言う。
魔方陣を通して組合せ最適化を学ぶ を参考にしました。 斉藤さんだぞ。(毎回言うやつ)
重要なのは 必要な条件を漏れなく制約条件として追加してあげることです。 一つでも抜けると正しい解は得られません。
備考
今回はセルの抽出を用途として pandas も使っています。
%autosave 0
%matplotlib inline
import pandas as pd
import pulp
df = pd.DataFrame([
(i, j , k , pulp.LpVariable(f'Var_{i}_{j}_{k}', cat=pulp.LpBinary))
for i in range(3)
for j in range(3)
for k in range(1, 10)
], columns=['row', 'col', 'num', 'Var'])
# 各セル(座標)に 1 ~ 9 までの数字を割り当て、適用可否を示すフラグを変数として持つ
df
# 3(行) * 3(列) * 9(1~9) = 81
len(df)
model = pulp.LpProblem()
# 各セルの値は必ず1回だけ使う
for key, group in df.groupby(['row', 'col']):
model += pulp.lpSum(group.Var) == 1
# 数字は必ず一回だけ使う
for key, group in df.groupby('num'):
model += pulp.lpSum(group.Var) == 1
# 行方向の合計が 15 になるような組み合わせ
for key, group in df.groupby('row'):
model += pulp.lpDot(group.num, group.Var) == 15
# 列方向の合計が 15 になるような組み合わせ
for key, group in df.groupby('col'):
model += pulp.lpDot(group.num, group.Var) == 15
# 対角線上の合計値 が 15 になるような組み合わせ
model += pulp.lpDot(*df.query('row==col')[['num','Var']].T.values) == 15
# 対角線(逆向き)上の合計値 が 15 になるような組み合わせ
model += pulp.lpDot(*df.query('row + col == 2')[['num','Var']].T.values) == 15
assert model.solve() == 1
上記は Var
が 0
or 1
の 2値をとるという性質を活用して制約条件を組んでいる
- 一回だけ使う というのは 合計値が 1
- 組み合わせの抽出は
使う
(1),使わない
(0) との内積
df['Val'] = df.Var.apply(pulp.value)
pd.DataFrame(df[df.Val>0.5].num.values.reshape(3, 3))
# 関数化してみる
def solve_square(num):
df = pd.DataFrame([
(i, j , k , pulp.LpVariable(f'Var_{i}_{j}_{k}', cat=pulp.LpBinary))
for i in range(num)
for j in range(num)
for k in range(1, num ** 2 + 1)
], columns=['row', 'col', 'num', 'Var'])
model = pulp.LpProblem()
expected_sum = sum(range(1, num ** 2 + 1)) / num
# 行, 列 につき一つの値しかとらない
for key, group in df.groupby(['row', 'col']):
model += pulp.lpSum(group.Var) == 1
# 値は 1 回ずつしか使わない
for key, group in df.groupby('num'):
model += pulp.lpSum(group.Var) == 1
# 行方向の内積合計値
for key, group in df.groupby('row'):
model += pulp.lpDot(group.num, group.Var) == expected_sum
# 列方向の内積合計値
for key, group in df.groupby('col'):
model += pulp.lpDot(group.num, group.Var) == expected_sum
# 対角線上の内積合計値
model += pulp.lpDot(*df.query('row == col')[['num','Var']].T.values) == expected_sum
# 対角線(逆向き)上の内積合計値
model += pulp.lpDot(*df.query('row + col == {}'.format(num-1))[['num','Var']].T.values) == expected_sum
assert pulp.LpStatus[model.solve()] == 'Optimal'
df['Val'] = df.Var.apply(pulp.value)
return pd.DataFrame(df[df.Val > 0.5].num.values.reshape(num, num))
solve_square(3)
solve_square(4)
solve_square(5)
solve_square(6)
ただの写経では面白くないので、応用として 違う大きさの 魔方陣も作っています。
6x6 の魔方陣は 多分処理に数時間かかってます。放置してたので正確な時間はわかりませんが。
やっぱり魔法ではないんですね〜
DataFrameの値がどのように抽出されるか見たい方はこちら
import pandas as pd
import pulp
df = pd.DataFrame([
(i, j , k , pulp.LpVariable(f'Var_{i}_{j}_{k}', cat=pulp.LpBinary))
for i in range(3)
for j in range(3)
for k in range(1, 10)
], columns=['row', 'col', 'num', 'Var'])
for key, group in df.groupby(['row', 'col']):
print(f'---------- {key} ----------')
print(group)
for key, group in df.groupby('num'):
print(f'---------- {key} ----------')
print(group)
for key, group in df.groupby('row'):
print(f'---------- {key} ----------')
print(group)
for key, group in df.groupby('col'):
print(f'---------- {key} ----------')
print(group)
# 対角セル (横長だと見にくいので転地した)
pd.DataFrame(df.query('row == col')[['num','Var']].T.values.T)
# 反対角セル (横長だと見にくいので転地した)
pd.DataFrame(df.query('row + col == 2')[['num','Var']].T.values.T)