DEA基础模型python代码

最新推荐文章于 2024-08-04 11:40:08 发布

小p127

最新推荐文章于 2024-08-04 11:40:08 发布

阅读量2.2k

点赞数 3

文章标签： python 开发语言

本文链接：https://blog.csdn.net/weixin_63461935/article/details/131501335

版权

DEA简介

DEA是数据包络分析（ Data envelopment analysis）的简称，是通过数学规划模型得出效率的非参数方法，不需要对生产函数的形式进行假定。以上文字是截取自吴杰老师的《行业分析视角下中国区域环境效率研究——基于数据包络分析（DEA）方法》中的描述。（其实我的理解就是初中就学过的不等式组，只不过没有具体数字了，求最佳值也跟以前在可行域内求最大最小值的思想一样。）吴杰老师这本书对于刚接触DEA方法的人来说真的很友好，各种基础模型之间的联系和区别都写的很清楚，不过要是想了解最新的模型，还是建议去看最新的英文文献哦。（本来是想详细写简介的，可是觉得字太多会不想看）

CCR与BCC模型

这两个模型算是DEA模型中基础中的基础了，从包络型来看区别主要就是增加了一个 $\sum \lambda = 1$ ，乘子型二者也有一定区别。

CCR乘子型(已归一化)

$\begin{aligned} max \quad &\sum_{r=1}^s u_{r}{y_{rd}} \\ s.t. \quad &\sum_{i=1}^m {v_{i}x_{id}} = 1;\\ &\sum_{r=1}^s{u_{r}y_{rj}} - \sum_{i=1}^m{v_{i}x_{ij}} \leq 0,\quad j=1,...,n;\\ &\forall v_{i},u_{r} \ge 0. \end{aligned}$

用的gurobi线性规划求解器哦，按道理来说，用linprog也可以，matlab里面就用的它，不过我因为师姐用的这个，我们就传承下来了。后面就把def ccr一整块换到，然后引用的dea.ccr()换成定义的名称就好了，比如dea.bcc().

from gurobipy import *
import pandas as pd
import numpy as np


class DEA(object):
    def __init__(self, DMU_name, X, Y):
        self.m1 = X.shape[1]
        self.m2 = Y.shape[1]
        self.DMUs, self.X, self.Y = multidict({DMU: [X.loc[DMU].tolist(),
                                                     Y.loc[DMU].tolist()] for DMU in DMU_name})

    def ccr(self):
        """计算theta dd 的值"""
        lst = []
        a = []
        for d in self.DMUs:
            m = Model()
            v = m.addVars(self.m1)
            u = m.addVars(self.m2)
            m.update()

            m.setObjective(quicksum(u[i] * self.Y[d][i] for i in range(self.m2)), sense=GRB.MAXIMIZE)
            m.addConstr(quicksum(v[i] * self.X[d][i] for i in range(self.m1)) == 1)
            m.addConstrs(quicksum(u[i] * self.Y[j][i] for i in range(self.m2))
                         - quicksum(v[i] * self.X[j][i] for i in range(self.m1)) <= 0 for j in self.DMUs)
            m.setParam('OutputFlag', 0)
            m.setParam('NonConvex', 2)
            m.optimize()
            lst.append(m.objVal)
        return lst
if __name__ == '__main__':
    data = pd.read_excel('C:/Users/wzu/Desktop/例1.xlsx', index_col=0, header=0)
    i1 = 2  # 表示表格中X有几列
    i2 = 1
    X = data[data.columns[:i1]]     # 原始数据
    Y = data[data.columns[i1:i1 + i2]]
    dea = DEA(DMU_name=data.index, X=X, Y=Y)
    print(dea.ccr()[0])

CCR包络型

$\begin{aligned} min\quad\ &\theta\\ s.t.\quad-&\sum_{j=1}^n \lambda_{j}x_{ij} + \theta x_{id} \ge 0 \quad i=1,...,m;\\ &\sum_{j=1}^n \lambda_{j}y_{rj} \ge y_{rd} \quad r=1,...,s; \\ &\forall \lambda_{j} \ge 0,\quad j=1,...,n. \end{aligned}$

def ccr(self):
        for k in self.DMUs:
            m = Model()
            theta = m.addVar()  #单个变量不加s
            lambdas = m.addVars(self.DMUs)
            m.update()  #以上都是变量
            m.setObjective(theta, sense=GRB.MAXIMIZE)
            #添加约束
            m.addConstrs(self.X[k][j] >= quicksum(lambdas[i] * self.X[i][j] for i in self.DMUs) for j in range(self.m1))
            m.addConstrs(theta*self.Y[k][j] <= quicksum(lambdas[i] * self.Y[i][j] for i in self.DMUs) for j in range(self.m2))

            #setParam(paramname,newvalue)是修改Gurobi参数，使运行结果更快，参考第三章P75
            m.setParam('OutputFlag', 0) #去掉计算过程，只取最后结果;
            #m.setParam('NonConvex', 2) #凸与非凸再查查
            m.optimize()
            print(m.objVal)

还可以在模型中引入松弛变量（其实起不了太大作用），模型为

$\begin{aligned} min\quad\ &\theta\\ s.t.\quad&\sum_{j=1}^n \lambda_{j}x_{ij} + s^+ =\theta x_{id} \quad i=1,...,m;\\ &\sum_{j=1}^n \lambda_{j}y_{rj}-s^- = y_{rd} \quad r=1,...,s; \\ &\forall s^+,s^-,\forall \lambda_{j} \ge 0,\quad j=1,...,n. \end{aligned}$

    def ccr2(self):
        for k in self.DMUs:
            lst = []
            m = Model()
            theta = m.addVar()
            lambdas = m.addVars(self.DMUs)
            sx = m.addVars(self.m1)
            sy = m.addVars(self.m2)
            m.update()
            m.setObjective(theta, sense=GRB.MINIMIZE)
            # 添加约束
            m.addConstrs(theta*self.X[k][j] - sx[j] == quicksum(lambdas[i] * self.X[i][j] for i in self.DMUs) for j in range(self.m1))
            m.addConstrs(self.Y[k][j] + sy[j] == quicksum(lambdas[i] * self.Y[i][j] for i in self.DMUs) for j in range(self.m2))
            m.setParam('OutputFlag', 0)
            m.setParam('NonConvex', 2)
            m.optimize()
            for j in range(self.m1):
                lst.append(sx[j].x)
            for j in range(self.m2):
                lst.append(sy[j].x)
            print(m.objVal)
            print(lst) #导出的为乘数

上述径向处理 $\theta$ 是在投入中，但也可以将径向处理放在产出中，并引入松弛变量模型为：

$\begin{aligned} min\quad\ &\theta\\ s.t.\quad &\sum_{j=1}^n \lambda_{j}y_{rj} -s^- - \theta y_{rd} = 0 \quad r=1,...,s;\\ &\sum_{j=1}^n \lambda_{j}x_{ij}+s^+ = x_{id} \quad i=1,...,m; \\ &\forall s^-,s^+,\lambda_{j} \ge 0,\quad j=1,...,n. \end{aligned}$

    def ccr3(self):
        for k in self.DMUs:
            m = Model()
            theta = m.addVar()
            lambdas = m.addVars(self.DMUs)
            sx = m.addVars(self.m1)
            sy = m.addVars(self.m2)
            m.update()
            m.setObjective(theta, sense=GRB.MAXIMIZE)
            #添加约束
            m.addConstrs(self.X[k][j] - sx[j] == quicksum(lambdas[i] * self.X[i][j] for i in self.DMUs) for j in range(self.m1))
            m.addConstrs(theta*self.Y[k][j] + sy[j] == quicksum(lambdas[i]* self.Y[i][j] for i in self.DMUs) for j in range(self.m2))
            m.setParam('OutputFlag', 0)
            #m.setParam('NonConvex', 2)
            m.optimize()
            print(m.objVal)

BCC乘子型

$\begin{aligned} max \quad &\sum_{r=1}^s u_{r}{y_{rd}}-u_{d} \\ s.t. \quad &\sum_{i=1}^m {v_{i}x_{id}} = 1;\\ &\sum_{r=1}^s{u_{r}y_{rj}} - \sum_{i=1}^m{v_{i}x_{ij}}-u_{d} \leq 0,\quad j=1,...,n;\\ &\forall v_{i},u_{r} \ge 0. \end{aligned}$

BCC包络型

$\begin{aligned} min\quad\ &\theta\\ s.t.\quad-&\sum_{j=1}^n \lambda_{j}x_{ij} + \theta x_{id} \ge 0 \quad i=1,...,m;\\ &\sum_{j=1}^n \lambda_{j}y_{rj} \ge y_{rd} \quad r=1,...,s; \\ &\sum_{j=1}^n \lambda_{j}=1;\\ &\forall \lambda_{j} \ge 0,\quad j=1,...,n. \end{aligned}$

    def bcc(self):
        for k in self.DMUs:
            m = Model()
            theta = m.addVar()  #单个变量不加s
            lambdas = m.addVars(self.DMUs)
            m.update()  #以上都是变量
            m.setObjective(theta, sense=GRB.MAXIMIZE)
            #添加约束
            m.addConstrs(self.X[k][j] >= quicksum(lambdas[i] * self.X[i][j] for i in self.DMUs) for j in range(self.m1))
            m.addConstrs(theta*self.Y[k][j] <= quicksum(lambdas[i] * self.Y[i][j] for i in self.DMUs) for j in range(self.m2))
            m.addConstr(quicksum(lambdas[i] for i in self.DMUs) == 1)

            #setParam(paramname,newvalue)是修改Gurobi参数，使运行结果更快，参考第三章P75
            m.setParam('OutputFlag', 0) #去掉计算过程，只取最后结果;
            #m.setParam('NonConvex', 2) #凸与非凸再查查
            m.optimize()
            print(m.objVal)

SBM模型

模型和内容好难打啊，我下次补上，这次直接放代码吧！

    def SBM(self):
        for k in self.DMUs:
            lst = []
            m = Model()
            t = m.addVar()
            lambdas = m.addVars(self.DMUs)
            Sx = m.addVars(self.m1)
            Sy = m.addVars(self.m2)
            m.update()
            m.setObjective(t - 1/self.m1 * quicksum(Sx[j] / self.X[k][j] for j in range(self.m1)), sense=GRB.MINIMIZE)
            # 添加约束
            m.addConstrs(
                t * self.X[k][j] - Sx[j] == quicksum(lambdas[i] * self.X[i][j] for i in self.DMUs) for j in range(self.m1))
            m.addConstrs(
                t * self.Y[k][j] + Sy[j] == quicksum(lambdas[i] * self.Y[i][j] for i in self.DMUs) for j in range(self.m2))
            m.addConstr(t + 1/self.m2*quicksum(Sy[j] / self.Y[k][j] for j in range(self.m2)) == 1)
            m.setParam('OutputFlag', 0)
            m.setParam('NonConvex', 2)
            m.optimize()
            print(m.objVal)
            for j in range(self.m1):
                lst.append(Sx[j].x)
            for j in range(self.m2):
                lst.append(Sy[j].x)
            print(lst)

哦，对，关于数据，我随便举得数据，形式大概是这样，