经济效率计算SBM

资源下载链接
https://download.csdn.net/download/qq_42830971/20586629

可联系邮箱找SBM资源下载：leon_leon@yeah.net

简介

SBM模型最早由Tone K(2001)构想并提出，该模型在目标函数中设定了松弛变量，以此有效弥补因松过变量产生误差的不足，2002年，Tone c基于原模型构建了SuperSBM便型，这一模型在测度不同样本效率值的基础上提供了对比有效样本效率高低的方法，通过允许有效样本的效率值大于1来达到目的。

原理

书上的原理是这样的



这个是不是突然看不懂，有过CCR或者BCC基础的应该能发现，约束由不等式变成了等式，这是因为加入了松弛变量S+/S-,S+/S-仅仅是区分投入松弛与产出松弛的两个符号。听我细细道来。以上是最简单的SBM方法。只需要记住经济效率是在数据中寻找一个标准；如果相同产出，更小的投入即是有效；相同投入，更多产出，也为有效，反之无效。

那就让我们瞧一瞧代码如何写，花了很大的功夫解读了一个代码包。解读完成就可以为所欲为了。

python线性规划基础

考虑如下线性规划问题
$max:Z=2x_1+3x_2-5x_3$
$st:x_1+x_2+x_3=7$
$2x_1-5x_2+x_3>=10$
$x_1+3x_2+x_3<=12$
$x_1,x_2,x_3>=0$

Step1: 导入相关库

import numpy as np 
from scipy import optimize as op 

Step2: 定义决策变量

 #给出变量取值范围
x1=(0,None)  
x2=(0,None)
x3=(0,None)

Step3: 将原问题化为标准形式

$min:-Z=-2x_1+-3x_2+5x_3$
$st:x_1+x_2+x_3=7$
$-2x_1+5x_2-x_3<=-10$
$x_1+3x_2+x_3<=12$
$x_1,x_2,x_3>=0$

注意：编程时默认为最小化目标函数，因此这里改为 ；第二个约束为大于等于约束，这里化为小于等于约束；

Step4: 定义目标函数系数和约束条件系数

c=np.array([-2,-3,5])   # 目标函数系数,3x1列向量
A_ub=np.array([[-2,5,-1],[1,3,1]]) # 不等式约束系数A，2x3维矩阵
B_ub=np.array([-10,12])  # 等式约束系数b, 2x1维列向量
A_eq=np.array([[1,1,1]])  # 等式约束系数Aeq，3x1维列向量
B_eq=np.array([7])   # 等式约束系数beq，1x1数值

Step5: 求解

res=op.linprog(c,A_ub,B_ub,A_eq,B_eq,bounds=(x1,x2,x3)) #调用函数进行求解
res

     con: array([0.])
     fun: -14.571428571428571
 message: 'Optimization terminated successfully.'
     nit: 3
   slack: array([0.        , 3.85714286])
  status: 0
 success: True
       x: array([6.42857143, 0.57142857, 0.        ])

目标函数及约束条件

首先明确变量都有啥

决策单元，即dmu1，dmu2…

t:

投入松弛 S-

期望产出松弛(好的) S+

非期望产出松弛(坏的) Sb

普通SBM模型代码：

def sbmeff2(input_variable, desirable_output, undesirable_output, dmu,data,method = 'revised simplex'):
    """用于求解sbm模型
    
    Parameters:
    -----------
    input_variable:
        投入[v1,v2,v3,...] 
    desirable_output:
        期望产出[v1,v2,v3,...]
    undesirable_output:
        非期望产出[v1,v2,v3,...] 
    dmu:
        决策单元
    data:
        主数据
    method:
        求解方法.默认'revised simplex',可选'interior-point'

	Return:
	------
	res : DataFrame
		结果数据框[dmu	TE	slack...]
    """
 
    res = pd.DataFrame(columns = ['dmu','TE'], index = data.index)
    res['dmu'] = data[dmu]
    ## lambda有dmu个数个，S有变量个数个
    dmu_counts = data.shape[0]
     ## 投入个数
    m = len(input_variable)
    ## 期望产出个数
    s1 = len(desirable_output)
    ## 非期望产出个数
    s2 = len(undesirable_output)
    #以下对应求解后的参数意义
    ## x[:dmu_counts] 为lambda
    ## x[dmu_counts:dmu_counts+1] 为 t
    ## x[dmu_counts+1 :dmu_counts + m + 1] 为投入slack
    ## x[dmu_counts+ 1 + m:dmu_counts + 1 + m + s1] 为期望产出slack
    ## x[dmu_counts + 1 + m + s1 :] 为非期望产出lack
    total = dmu_counts + m + s1 + s2 + 1
    cols = input_variable+desirable_output+ undesirable_output
    newcols = []
    for j in cols:
        newcols.append(j+'_slack')
        res[j+'_slack'] = np.nan
    for i in range(dmu_counts):
        
        ## 优化目标
        c = [0] * dmu_counts + [1] +  list(-1 / (m * data.loc[i, input_variable])) + [0] * (s1 + s2)
        
        ## 约束条件
        A_eq = [[0] * dmu_counts + [1] + [0] * m  + list(1/((s1 + s2) * data.loc[i, desirable_output])) + 
                                                    list(1/((s1 + s2) * data.loc[i, undesirable_output]))]
        
        ## 约束条件（1）
        for j1 in range(m):
            list1 = [0] * m
            list1[j1] = 1
            eq1 = list(data[input_variable[j1]]) + [-data.loc[i ,input_variable[j1]]] + list1 + [0] * (s1 + s2)
            A_eq.append(eq1)
        ## 约束条件（2）
        for j2 in range(s1):
            list2 = [0] * s1
            list2[j2] = -1
            eq2 = list(data[desirable_output[j2]]) + [-data.loc[i, desirable_output[j2]]] + [0] * m + list2 + [0] * s2
            A_eq.append(eq2)
        ## 约束条件（3）
        for j3 in range(s2):
            list3 = [0] * s2
            list3[j3] = 1
            eq3 = list(data[undesirable_output[j3]]) + [-data.loc[i, undesirable_output[j3]]] + [0] * (m + s1) + list3
            A_eq.append(eq3)         
        b_eq = [1] + [0] * (m + s1 + s2)               
        bounds = [(0, None)]*total
        ## 求解
        op1 = op.linprog(c = c,A_eq=A_eq,b_eq=b_eq,bounds=bounds,method = method)
        res.loc[i, 'TE'] = op1.fun
        res.loc[i, newcols] = op1.x[dmu_counts+1 :]
    return res

理解不困难，主要是目标函数与约束函数的理解

目标函数：

约束条件：
1个
$t+1/(s_1+s_2)\sum _{i=1}^{s_1}{s}_i^{-}/x_{ik}+1/(s_1+s_2)\sum _{i=1}^{s_2}{s}_{bi}^{-}/x_{ik}=1$

m个
$a{x}_1+b{x}_2+...+h{x}_h+系数\ast 第j1个投入指标=第i个dmu的第j1个投入指标\ast t$

s1个
$a{x}_1+b{x}_2+...+h{x}_h+系数\ast 第j1个期望产出指标=第i个dmu的第j1个期望产出指标\ast t$

s2个
$a{x}_1+b{x}_2+...+h{x}_h+系数\ast 第j1个非期望产出指标=第i个dmu的第j1个非期望产出指标\ast t$