この記事は ウェブクルー Advent Calendar 2017の5日目の記事です。
昨日は@sankaku_qさんの「Emacsで指定範囲の色を変える」でした。
はじめに
業務で最適化を行う案件が出てきたため、PythonのPuLPライブラリを使って実装してみました。
Python自体も学習中ですので、ここが間違っている、こうした方が良いなどございましたらコメントをいただけますと幸いです。
数理最適化とは
ある条件（制約）の中で、目的関数を最大または最小にする解を見つけること。
勤務表の作成や広告配信など色々な分野で使われることがあるようです。
ちなみに、学校で習う2次関数：f(x)が最小となるxを求めよ。も数理最適化。
例題
原料A、B、C、Dからなる製品X、Y、Zがある。
Xを作るのに必要な原料は、Aが3個、Cが8個、Dが5個である。
Yを作るのに必要な原料は、Bが2個、Cが12個である。
Zを作るのに必要な原料は、Aが5個、Bが8個、Dが10個である。
Xを生産したときの利益は70円。
Yを生産したときの利益は100円。
Zを生産したときの利益は20円。
今ある原料はAが74個、Bが62個、Cが99個、Dが80個だとすると、
利益を最大化するにはX、Y、Zを何個ずつ作成すればよいかを求めよ。
問題を定式化すると下記のように表すことができます。
変数 ： $ X, Y, Z >= 0 $（整数）
目的関数 ： $ 70X + 100Y + 20Z $ → 最大化
制約条件 ： $ 3X + 5Z <= 74 $
$ 2Y + 8Z <= 62 $
$ 8X + 12Y <= 99 $
$ 5X + 10Z <= 80 $
PythonのPuLPでモデル化
from pulp import *
m = LpProblem(sense=LpMaximize) # 数理モデル
x = LpVariable('x', lowBound=0, cat=LpInteger) # 変数
y = LpVariable('y', lowBound=0, cat=LpInteger) # 変数
z = LpVariable('z', lowBound=0, cat=LpInteger) # 変数
m += 70 * x + 100 * y + 20 * z # 目的関数
m += 3 * x + 5 * z <= 74 # 材料Aの上限の制約条件
m += 2 * y + 8 * z <= 62 # 材料Bの上限の制約条件
m += 8 * x + 12 * y <= 99 # 材料Cの上限の制約条件
m += 5 * x + 10 * z <= 80 # 材料Dの上限の制約条件
m.solve() # ソルバーの実行
print(value(x), value(y), value(z), value(m.objective)) # 3.0, 6.0, 6.0, 930
Xを3個、Yを6個、Zを6個作ると利益が最大の930円になることがわかりました。
目的関数も制約も式のような形で書けるのでラクです。
以下、順番に説明します。
パッケージのインポート
from pulp import *
数理モデルの作成
m = LpProblem(sense=LpMaximize)
- sense：最小化か最大化を指定する
- LpMinimize：最小化
- LpMaximize：最大化
変数の作成
x = LpVariable(name, lowBound=None, upBound=None, cat='Continuous')
- name：変数名 ※重複NG
- lowBound：下限
- upBound：上限
- cat：変数の種類
- LpContinuous：連続変数
- LpInteger：整数変数
- LpBinary：0-1変数
目的関数
m += 式
制約条件
m += 式 == 式
m += 式 <= 式
m += 式 >= 式
ソルバーの実行
m.solve()
ソルバーとは、複数の変数を含む数式において、目標とする値を得るための最適な値を求めることができる機能のこと。
各値の表示
変数の値：value(変数名)
目的関数の値：value(m.objective)
本題
さて、始めに業務で使いたい案件があると書きましたが、内容は以下のようなものです。
ユーザがX社、Y社、Z社それぞれの想定成約率を持っている。
各社、配信数の下限と上限があり、1人1社しか配信できない。
ユーザをどう配分すれば成約率（成約数）が最大になるかを求めよ。
- 各ユーザに最適な会社を配信したい → 変数
- 成約率（成約数）の最大化 → 目的関数
- 配信は1人1社 → 制約条件①
- 配信の下限はX社/Y社：リストの29%、Z社：なし → 制約条件②
- 配信の上限はX社/Y社：7,000件、Z社：なし → 成約条件③
データは下の表を用意しました。
| ユーザID | X社 | Y社 | Z社 | |
|---|---|---|---|---|
| 0 | 1 | 0.017 | 0.026 | 0.004 |
| 1 | 2 | 0.038 | 0.026 | 0.008 |
| 2 | 3 | 0.009 | 0.026 | 0.011 |
| 3 | 4 | 0.009 | 0.004 | 0.011 |
| : | : | : | : | : |
| 23743 | 23744 | 0.038 | 0.026 | 0.011 |
データを準備
用意したデータは1人1行で各社の成約率が横持ちのものなので、縦持ちに直します。
import pandas as pd
import math
from pulp import *
# データ取得
data = pd.read_csv('user_list.csv')
# 下限算出用に全体のデータ量取得
record_cnt = len(data.index) # 23744
# 各社ごとにデータを取得
data_x = data[['USER_ID', 'X']]
data_x.rename(columns={'X': 'CONTRACT_RATIO'}, inplace=True)
data_x.loc[:, 'COMPANY_NAME'] = 'X'
data_y = data[['USER_ID', 'Y']]
data_y.rename(columns={'Y': 'CONTRACT_RATIO'}, inplace=True)
data_y.loc[:, 'COMPANY_NAME'] = 'Y'
data_z = data[['USER_ID', 'Z']]
data_z.rename(columns={'Z': 'CONTRACT_RATIO'}, inplace=True)
data_z.loc[:, 'COMPANY_NAME'] = 'Z'
# 各社のデータを結合し、indexの振りなおし
data = data_x.append([data_y, data_z])
data = data.reset_index(drop=True)
縦持ちに直した結果のデータがこちら↓
| USER_ID | COMPANY_NAME | CONTRACT_RATIO | |
|---|---|---|---|
| 0 | 1 | X | 0.017 |
| 1 | 2 | X | 0.038 |
| : | : | : | : |
| 23744 | 1 | Y | 0.026 |
| 23745 | 2 | Y | 0.026 |
| : | : | : | : |
| 47488 | 1 | Z | 0.004 |
| 47489 | 2 | Z | 0.008 |
| : | : | : | : |
数理モデル作成
m = LpProblem(sense=LpMaximize) # 数理モデル
data['Var'] = [LpVariable('v%d' % i, cat=LpBinary) for i in data.index] # 変数
m += lpDot(data.CONTRACT_RATIO, data.Var) # 目的関数
for k, v in data.groupby('USER_ID'): # 制約条件
m += lpSum(v.Var) == 1
send_ratio = {"X":0.29, "Y":0.29, "Z":0.0}
for k, v in data.groupby('COMPANY_NAME'): # 制約条件
m += lpSum(v.Var) >= math.floor(record_cnt * send_ratio[k])
max_cnt = {"X":7000, "Y":7000, "Z":999999}
for k, v in data.groupby('COMPANY_NAME'): # 制約条件
m += lpSum(v.Var) <= max_cnt[k]
m.solve() # ソルバーの実行
# 配信するべきとなったもののみ表示
data['Val'] = data.Var.apply(value)
print(data[data.Val == 1])
print(value(m.objective)) # 507.268
最終的に表示されたデータがこちら↓
| USER_ID | COMPANY_NAME | CONTRACT_RATIO | Var | Val | |
|---|---|---|---|---|---|
| 0 | 1 | X | 0.017 | Var0 | 1.0 |
| 1 | 2 | X | 0.038 | Var1 | 1.0 |
| 23746 | 3 | Y | 0.026 | Var23746 | 1.0 |
| 47491 | 4 | Z | 0.011 | Var47491 | 1.0 |
| ： | ： | ： | ： | ： | ： |
全体のユーザ数23,744人に対して、想定成約数が507人なので、想定成約率は2.1%となりました。
現状は1.5%なので、0.6%UPということですね。
ちなみに、データ取得から最適化を完了するまでの所要時間は、大体2分くらいでした。
以下、簡単に説明します。
リストの変数：バイナリ変数(1-0)のリスト
data['Var'] = [LpVariable('v%d' % i, cat=LpBinary) for i in data.index]
目的関数：各社成約率の総和
m += lpDot(data.CONTRACT_RATIO, data.Var)
制約条件①：1人1社
for k, v in data.groupby('USER_ID'):
m += lpSum(v.Var) == 1
制約条件②：配信下限数
send_ratio = {"X":0.29, "Y":0.29, "Z":0.0}
for k, v in data.groupby('COMPANY_NAME'):
m += lpSum(v.Var) >= math.floor(record_cnt * send_ratio[k])
各社の配信率を辞書として用意し、上で求めたユーザ数（record_cnt）と掛け合わせて最小数を算出
制約条件③：配信上限数
max_cnt = {"X":7000, "Y":7000, "Z":999999}
for k, v in data.groupby('COMPANY_NAME'):
m += lpSum(v.Var) <= max_cnt[k]
各社の配信上限数を辞書として用意し、そこから取得
まとめ
pythonで最適化をしてみましたが、数式に近い形で記述できるのでとてもわかりやすいと思いました。
今回は成約しやすい会社への最適化を試しましたが、メルマガ配信や商品の紹介等、色々な場面で活用できそうです。
注意点としては、目的関数と制約を正しく定めなければいけないので、そこはしっかりと考えなければいけませんね。
明日は、@noko_k さんです。よろしくお願いします。