python解決DEA模型–CCR

最近本人潛心學習python來解決DEA模型，DEA模型說難其實并不難，就是各種線性規劃，但是它作為一個方向必然有其難點，我認為就是模型太多了，不斷被創新，不斷有新模型，

近來，了解到有gurobi這個庫的存在，對于各種各樣的模型都方便處理，太好用了！

gurobi的下載請自行去官網下載https://www.gurobi.com/，

那么本篇文章就直接附上代碼，

import gurobipy
import pandas as pd

# 分頁顯示資料, 設定為 False 不允許分頁
pd.set_option('display.expand_frame_repr', False)

# 最多顯示的列數, 設定為 None 顯示全部列
pd.set_option('display.max_columns', None)

# 最多顯示的行數, 設定為 None 顯示全部行
pd.set_option('display.max_rows', None)

class DEA(object):
    def __init__(self, DMUs_Name, X, Y, AP=False):
        self.m1, self.m1_name, self.m2, self.m2_name, self.AP = X.shape[1], X.columns.tolist(), Y.shape[1], Y.columns.tolist(), AP
        self.DMUs, self.X, self.Y = gurobipy.multidict({DMU: [X.loc[DMU].tolist(), Y.loc[DMU].tolist()] for DMU in DMUs_Name})
        print(f'DEA(AP={AP}) MODEL RUNING...')
        #multidict擴展字典，便于處理同一個物件的不同屬性約束

    def __CCR(self):
        for k in self.DMUs:
            MODEL = gurobipy.Model()
            OE, lambdas, s_negitive, s_positive = MODEL.addVar(), MODEL.addVars(self.DMUs),  MODEL.addVars(self.m1),\
                                                  MODEL.addVars(self.m2)
            #決策變數 Model.addVar() 和 Model.addVars()，比如說是x = MODEL.addVar(lb=0.0, ub=gurobipy.GRB.INFINITY,
            # vtype=gurobipy.GRB.CONTINUOUS, name="")
            MODEL.update()
            ## 更新變數環境
            MODEL.setObjectiveN(OE, index=0, priority=1)
            #添加目標函式 Model.setObjective() 和 Model.setObjectiveN()，并且index: 目標函式對應的序號 (默認 0，1，2，…),
            # 以 index=0 作為目標函式的值, 其余值需要另外設定引數
            MODEL.setObjectiveN(-(sum(s_negitive) + sum(s_positive)), index=1, priority=0)
            #priority大就先算
            MODEL.addConstrs(gurobipy.quicksum(lambdas[i] * self.X[i][j] for i in self.DMUs if i != k or not self.AP)
                             + s_negitive[j] == OE * self.X[k][j] for j in range(self.m1))
            MODEL.addConstrs(gurobipy.quicksum(lambdas[i] * self.Y[i][j] for i in self.DMUs if i != k or not self.AP)
                             - s_positive[j] == self.Y[k][j] for j in range(self.m2))
            MODEL.setParam('OutputFlag', 0)
            MODEL.optimize()
            self.Result.at[k, ('效益分析', '綜合技術效益(CCR)')] = MODEL.objVal
            self.Result.at[k, ('規模報酬分析', '有效性')] = '非 DEA 有效' if MODEL.objVal < 1 else 'DEA 弱有效' \
                if s_negitive.sum().getValue() + s_positive.sum().getValue() else 'DEA 強有效'
            self.Result.at[k, ('規模報酬分析', '型別')] = '規模報酬固定' if lambdas.sum().getValue() == 1 \
                else '規模報酬遞增' if lambdas.sum().getValue() < 1 else '規模報酬遞減'
            for m in range(self.m1):
                self.Result.at[k, ('差額變數分析', f'{self.m1_name[m]}')] = s_negitive[m].X
                self.Result.at[k, ('投入冗余率',  f'{self.m1_name[m]}')] = 'N/A' if self.X[k][m] == 0 \
                    else s_negitive[m].X / self.X[k][m]
            for m in range(self.m2):
                self.Result.at[k, ('差額變數分析', f'{self.m2_name[m]}')] = s_positive[m].X
                self.Result.at[k, ('產出不足率', f'{self.m2_name[m]}')] = 'N/A' if self.Y[k][m] == 0 \
                    else s_positive[m].X / self.Y[k][m]
        return self.Result
    def dea(self):
        columns_Page = ['效益分析'] + ['規模報酬分析'] *2 + ['差額變數分析'] * (self.m1 + self.m2) + \
                       ['投入冗余率'] * self.m1 + ['產出不足率'] * self.m2
        columns_Group = ['綜合技術效益(CCR)','有效性', '型別'] + (self.m1_name + self.m2_name) * 2
        self.Result = pd.DataFrame(index=self.DMUs, columns=[columns_Page, columns_Group])
        self.__CCR()
        return self.Result

    def analysis(self, file_name=None):
        Result = self.dea()
        file_name = 'DEA 資料包絡分析報告.xlsx' if file_name is None else f'\\{file_name}.xlsx'
        Result.to_excel(file_name, 'DEA 資料包絡分析報告')
if __name__=='__main__':
    innum, outnum = 2, 1
    file = "C://Users//admin//Desktop//資料.xlsx"
    data = pd.read_excel(file, header=0, index_col=0)
    print(data)
    X = data[data.columns[:innum]]
    Y = data[data.columns[innum:innum+outnum]]
    print(X)
    print(Y)
    dea = DEA(DMUs_Name=data.index, X=X, Y=Y)
    #dea.analysis()  # dea 分析并輸出表格
    print(dea.dea())  # dea 分析，不輸出結果