Advertisement

DEA基础模型python代码

阅读量：

DEA简介

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

CCR与BCC模型

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

CCR乘子型(已归一化)
egin{aligned} max uad &um_{r=1}^s u_{r}{y_{rd}} s.t. uad &um_{i=1}^m {v_{i}x_{id}} = 1; &um_{r=1}^s{u_{r}y_{rj}} - um_{i=1}^m{v_{i}x_{ij}} eq 0,uad j=1,...,n; &orall v_{i},u_{r} e 0. nd{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])
    
    
    
    
    
![](https://ad.itadn.com/c/weblog/blog-img/images/2025-08-18/ZnMy0BQkcjlPfU8G6Jdtz53AwFea.png)

CCR包络型
egin{aligned} minuad&heta s.t.uad-&um_{j=1}^n ambda_{j}x_{ij} + heta x_{id} e 0 uad i=1,...,m; &um_{j=1}^n ambda_{j}y_{rj} e y_{rd} uad r=1,...,s; &orall ambda_{j} e 0,uad j=1,...,n. nd{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)
    
    
    
    
    
![](https://ad.itadn.com/c/weblog/blog-img/images/2025-08-18/c4ESdhW5FMeuUplkCgAztronsPfG.png)

还可以在模型中引入松弛变量（其实起不了太大作用），模型为
egin{aligned} minuad&heta s.t.uad&um_{j=1}^n ambda_{j}x_{ij} + s^+ =heta x_{id} uad i=1,...,m; &um_{j=1}^n ambda_{j}y_{rj}-s^- = y_{rd} uad r=1,...,s; &orall s+,s-,orall ambda_{j} e 0,uad j=1,...,n. nd{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) #导出的为乘数
    
    
    
    
    
![](https://ad.itadn.com/c/weblog/blog-img/images/2025-08-18/ZFrwckpHxLMzNEhB3VJyqKjTugft.png)

上述径向处理heta是在投入中，但也可以将径向处理放在产出中，并引入松弛变量模型为：
egin{aligned} minuad&heta s.t.uad &um_{j=1}^n ambda_{j}y_{rj} -s^- - heta y_{rd} = 0 uad r=1,...,s; &um_{j=1}^n ambda_{j}x_{ij}+s^+ = x_{id} uad i=1,...,m; &orall s-,s+,ambda_{j} e 0,uad j=1,...,n. nd{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)
    
    
    
    
    
![](https://ad.itadn.com/c/weblog/blog-img/images/2025-08-18/0ORvwed5S7oaZUY92cEh3mJfBlky.png)

BCC乘子型
egin{aligned} max uad &um_{r=1}^s u_{r}{y_{rd}}-u_{d} s.t. uad &um_{i=1}^m {v_{i}x_{id}} = 1; &um_{r=1}^s{u_{r}y_{rj}} - um_{i=1}^m{v_{i}x_{ij}}-u_{d} eq 0,uad j=1,...,n; &orall v_{i},u_{r} e 0. nd{aligned}

BCC包络型
egin{aligned} minuad&heta s.t.uad-&um_{j=1}^n ambda_{j}x_{ij} + heta x_{id} e 0 uad i=1,...,m; &um_{j=1}^n ambda_{j}y_{rj} e y_{rd} uad r=1,...,s; &um_{j=1}^n ambda_{j}=1; &orall ambda_{j} e 0,uad j=1,...,n. nd{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)
    
    
    
    
    
![](https://ad.itadn.com/c/weblog/blog-img/images/2025-08-18/asmloy2I7XUbP4Y1e3r6ASBt9pzx.png)

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)
    
    
    
    
    
![](https://ad.itadn.com/c/weblog/blog-img/images/2025-08-18/MQoa1scwj6zLNWdFTfVE9Bn3rKHu.png)

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

当一个人开始学习时会发现除了学习，其他的所有事物都好趣，所以有了这几篇博客的诞生。

全部评论 (0)

还没有任何评论哟~

相关文章推荐

DEA基础模型python代码

DEA简介 DEA是数据包络分析（Dataenvelopmentanalysis）的简称，是通过数学规划模型得出效率的非参数方法，不需要对生产函数的形式进行假定。以上文字是截取自吴杰老师的《行业分析视...

用python处理DEA模型--CCR

python解决DEA模型–CCR 本人所用软件为gurobi，在python中运算。 gurobi的下载请自行前往官网：https://www.gurobi.com/。 本篇文章针对基础CCR模型，...

DEA各种模型原理及stata代码实现

DEA各种模型原理及stata代码实现 DEA各种模型原理及stata代码实现 一、CCR和BCC 1.原理 2.效率测算stata代码 3.Malmquist指数 3.1M指数 3.2GlobalM...

DEA模型大全（DEAp2.1、DEAslover、MAXDEA、三阶段DEA模型））

ta作DEA的do文件及案例演示数据 5.三阶段DEA模型案例案例演示 、do文件、数据 6.DEA相关经典参考文献 以上方法我 们都搜集了的都比较详细,每份都配备了相关案例,大家都可以通俗易懂的学会...

零和收益DEA模型（ZSG-DEA）的求解

今天要向大家推介的是DEA中的ZSG模型。 Lins和Gomes等人提出了零和收益DEA模型（ZSGDEA），该模型可以根据决策单元的DEA效率值对非期望产出的分配方案进行调整，给出改进的分配方案，该...

基于零和收益的DEA模型研究python实现

零和收益DEA模型（ZSG-DEA）

DEA模型及MATLAB应用2：BCC模型

DEA模型及MATLAB应用2：BCC模型 在上一篇博客中，我们介绍了CCR模型作为DEA模型的一种最早形式。CCR模型假设所有的DMUs都具有相同的规模效率。

python基础代码大全-python基础代码大全

【实例简介】 python代码大全，适用于基础python学习者，里面的代码基本上是基础学习者必经过程。 Python学习入门很快，但学习之路任重道远 【实例截图】 【核心代码】 python代码大全...

matlab 三阶段dea模型,三阶段DEA模型SFA二阶段剔除过程

楼主你好，跑出来的ce值都小于1，啥子情况？对各个值也不清楚如何分析？头疼死了，可以推荐相关资料么 thefinalmleestimatesare: coefficientstandarderrort...

python基础代码

一晃两个月了，好长时间不写了。有太多的借口去推脱，但静下来却也发觉，每个推脱都显的苍白。 今天写python，为什么选择py。因为python属于编程语言，而且也是开源的，在人工智能以后的发展上，目前...