PuLPを利用した包絡分析モデル

包絡分析法について

包絡分析法について、わかりやすい解説とPythonでのデモが少なかったので、
自分なりの理解をまとめる。
環境としてJupyter Notebookを利用している。
またpulpの利用法は自分の過去のページに記載しているので特に解説などは行わない。
なおPythonで包絡分析をするだけの場合はライブラリが存在する。

1.概要

包絡分析法(DEA:Data Envelopment Analysis)とは、複数の事業の相対的な効率を評価する手法

1.1 例

ある企業のn個の店舗で、それぞれの店舗について支出(人件費、売り場面積、広告代など)を入力、収入(売上、来客数、リピート率など)を出力と捉えれば「出力/入力」が大きいほど事業の効率がいいといえる。

これら各入力と各出力の値に適当な重みを付けることで、それらの和を仮想的な入力と仮想的な出力を考える。このとき各店舗の効率を公平に評価するために、重みを付けて得られる 「仮想的な出力/仮想的な入力」 を比較していく。

1.2 方針

➀ある店舗kに着目して、重さを設定する。
➁他店舗に同じ重さを適用して効率が１を超えないよう➀の重さを調整する。
➂全店舗で➀、➁を繰り返す
以上を行って、「仮想的な出力/仮想的な入力」 を求める。

2.包絡分析(DEA)モデル

2.1 データ

データとして以下を利用する。
2.4 包絡分析法（DEA）モデル
このデータでは入力が店員数と稼働時間、出力が売上とみられる。

import pandas as pd
df = pd.DataFrame()
df["店員数"] = [5,10,20,20,30,50]
df["稼働時間"] = [24,12,12,24,12,12]
df["売上"]= [2,6,10,12,12,20]

店舗	店員数	稼働時間	売上
0	5	24	2
1	10	12	6
2	20	12	10
3	20	24	12
4	30	12	12
5	50	12	20

これをPuLPを利用して線形計画問題として包絡分析していく。

2.2 コード

import pulp

kai_df = pd.DataFrame(columns=["x1","x2","y","目的変数"])
for k in range(df.shape[0]):
    problem  = pulp.LpProblem('DEA', pulp.LpMaximize)

    #変数設定
    x = {}
    y = {}
    for p in range(df.shape[0]):
        #入力の重み
        x[p] = pulp.LpVariable.dicts(f"x_{p}",range(df.shape[1]-1),lowBound=0,cat="Continuous")
        #出力の重み
        y[p] = pulp.LpVariable(f"y_{p}",lowBound=0,cat="Continuous")

    #目的関数
    problem += y[k] * df.T[k][2]

    #制約条件1(他店舗にも重みを適用)
    for i in range(df.shape[0]):
        problem += y[k] * df.T[i][2] <= pulp.lpSum([x[k][j] * df.T[i][j] for j in range(df.shape[1]-1)])
    
    #制約条件2(仮想的な出力を1とする)
    problem += pulp.lpSum([x[k][i] * df.T[k][i] for i in range(df.shape[1]-1)]) == 1
    
    #解
    problem.solve()
    
    #結果を保存
    kai_df.loc[k] = [x[k][0].value(),x[k][1].value(),y[k].value(),y[k].value()*df.T[k][2]]

コードは改善の余地がありそうではある。
包絡分析法は一つの事業に固定(上記コードでいうk)して、他の事業との比較をしつつ最適化していくので店舗数分最適化問題を解く必要がある。

2.3 結果

kai_dfを出力すると下記が得られる。(小数は小数第四位で四捨五入)

店舗	x1	x2	y	目的変数
0	0.200	0	0.333	0.667
1	0.100	0	0.167	1.000
2	0.040	0.017	0.100	1.000
3	0.050	0	0.083	1.000
4	0.025	0.021	0.075	0.900
5	0.000	0.083	0.050	1.000

以上参考ページと同じ値が得られた。
表から、店舗1,2,3,5が他店舗と比較して効率的に売上を出しているのがわかる。
特に店舗0の効率が悪く、他店舗と比較すると店員数を3人にしても同じ売り上げが得られそうという予測が出来る。
ただし、もちろん他の要因も考えるべきではあり、例えば出力として 「従業員満足度」を考慮すると、結果は変わってくるだろう。

3.終わりに

今回はpulpによる包絡分析モデルの作成を行った。
自分の理解のために行ったが、上でも書いた通りPythonで包絡分析法を行う場合はライブラリが存在するのでそちらを利用していくのだろうが、本記事で理論的な理解に繋がる部分があればと思う。

追記

コードの使用が気に入らなかったので改善した。
クラスにしつつ、包絡分析法の定義に基づいた形にして入力と出力の個数を任意に対応できるようにした。

import pandas as pd
import pulp

class DEA:
    def __init__(self, df, m, n):
        self.df = df
        self.m = m
        self.n = n
        
    def solve(self, df):
        list = [f'x{i}' if i <= (self.m-1) else f'y{i-self.m}' for i in range(self.m+self.n)]
        list.append("効率")
        kai_df = pd.DataFrame(columns=list)

        for k in range(df.shape[0]):
            problem  = pulp.LpProblem('DEA', pulp.LpMaximize)

            x = {}
            y = {}
            for p in range(df.shape[0]):
                #入力の重み変数
                x[p] = pulp.LpVariable.dicts(f"x_{p}",range(self.m),lowBound=0,cat="Continuous")
                #出力の重み変数
                y[p] = pulp.LpVariable.dicts(f"y_{p}",range(self.n),lowBound=0,cat="Continuous")


            #目的関数
            problem += pulp.lpSum([y[k][r] * df.T[k][r+self.m] for r in range(self.n)])

            #制約条件
            for i in range(df.shape[0]):
                problem += pulp.lpSum([y[k][r] * df.T[i][r+self.m] for r in range(self.n)]) <= pulp.lpSum([x[k][s] * df.T[i][s] for s in range(self.m)])

            #仮の入力を１とする
            problem += pulp.lpSum([x[k][s] * df.T[k][s] for s in range(self.m)]) == 1

            #解
            problem.solve()

            #結果を保存
            kai_list = [f'{x[k][i].value()}' if i <= (self.m-1) else f'{y[k][i-self.m].value()}' for i in range(self.m+self.n)]
            kai_list.append(sum([y[k][r].value() * df.T[k][r+self.m] for r in range(self.n)]))
            kai_df.loc[k] = kai_list

        return kai_df

参考文献

この記事は以下の情報を参考にして執筆しました。

・「しっかり学ぶ数理最適化　モデルからアルゴリズムまで」梅谷俊治 (https://bookclub.kodansha.co.jp/product?item=0000275620)

・株式会社NTTデータ数理システム　2.4 包絡分析法（DEA）モデル(https://www.msi.co.jp/solution/nuopt/docs/examples/html/02-04-00.html)