利用python的pulp库进行CCR、BCC、超效率模型的数学建模

本文介绍了如何使用Python的pulp库建立CCR、BCC和超效率模型的数学模型,提供了数据示例和面向对象的类函数实现。通过实例展示了如何利用pulp库进行线性规划,解决数据包络分析问题,包括最小化目标函数和约束条件的设定。并给出了基于VRS的CCR和超效率模型的计算结果。
Python3.8

Python 是一种高级、解释型、通用的编程语言,以其简洁易读的语法而闻名,适用于广泛的应用,包括Web开发、数据分析、人工智能和自动化脚本


利用python的pulp库进行CCR、BCC、超效率模型的数学建模

本文参考书目为
《数据包络分析方法与MaxDEA软件》


一、CCR、BCC、超效率模型公式?

1. 投入导向的CCR对偶公式

示例:pandas 是基于NumPy 的一种工具,该工具是为了解决数据分析任务而创建的。

2. BCC模型公式

在这里插入图片描述

3. 超效率模型公式

在这里插入图片描述

二、数学模型构建

1.数据示例

import pandas as pd
data = pd.DataFrame([[4,3,1],[7,3,1],[8,1,1],[4,2,1],[2,4,1],[10,1,1],[12,1,1],[10,1.5,1]],columns=['x1','x2','y'],index=['A','B','C','D','E','F','G','H'])
data
>>>
    x1	x2	y
A	4	3.0	1
B	7	3.0	1
C	8	1.0	1
D	4	2.0	1
E	2	4.0	1
F	10	1.0	1
G	12	1.0	1
H	10	1.5	1

1.1 pulp 包使用方法

在这里插入图片描述

d = np.array([2,3,1])
da = np.array([1,2,4])
da1 = np.array([[1,3],[4,2],[2,0]])
r = np.array([8,6])
m = pulp.LpProblem(sense=pulp.LpMinimize)
x1 = pulp.LpVariable.dict("x",range(3),lowBound=0)
x = list(x1.values())
m += pulp.lpSum([d[i0] * x[i0] for i0 in range(3)])
for i in range(2):
    m += pulp.lpSum([x[i0] * da1[i0][i] for i0 in range(3)]) >= r[i]
m += ( pulp.lpDot(da , x) == 101 )
m.solve()
pulp.value(m.objective)
>>>
28.75
>>>
# 以上可能写的太过复杂,可以参考一下简单例子
z = [2,3,1]
a = [[1,4,2],[3,2,0]]
b = [8,6]
c = [1,2,4]
d = 101
m1 = pulp.LpProblem(sense=pulp.LpMinimize)  #定义函数
x = [pulp.LpVariable(f"x_{i}",lowBound=0) for i in range(3)]  #设置变量
m1 += pulp.lpDot(z, x)  # "m1 +=" 是固定符号,后面若无约束,则说明这条语句为目标函数 
m1 += pulp.lpDot(c,x) == d # 约束条件
for i in range(len(a)):    # 约束条件
    m1 += (pulp.lpDot(a[i], x) >= b[i])
m1.solve()
pulp.value(m1.objective)   # 得到的值
>>>
28.75
>>>
m1 #查看定义函数的整个公式情况
>>>
NoName:
MINIMIZE
2*x_0 + 3*x_1 + 1*x_2 + 0
SUBJECT TO
_C1: x_0 + 2 x_1 + 4 x_2 = 101

_C2: x_0 + 4 x_1 + 2 x_2 >= 8

_C3: 3 x_0 + 2 x_1 >= 6

VARIABLES
x_0 Continuous
x_1 Continuous
x_2 Continuous
>>>

1.2 学会基础建模后,套用书上给的CCR公式

# 非超效率CCR模型
import pulp as lp
DMUs = data.index.values
d_var = data.columns.values
# 变量
var = lp.LpVariable.dict("DMU",data.index.values,lowBound=0)
thita = lp.LpVariable('thita')

for DMU in DMUs:
    # 目标函数
    pr = lp.LpProblem('CCR',lp.LpMinimize)
    pr += thita
    # 约束条件
    pr += lp.lpSum(data.to_dict()['x1'][i] * var[i] for i in DMUs) <= thita * data.to_dict()['x1'][DMU]
    pr += lp.lpSum(data.to_dict()['x2'][i] * var[i] for i in DMUs) <= thita * data.to_dict()['x2'][DMU]
    pr += lp.lpSum(data.to_dict()['y'][i] * var[i] for i in DMUs) >=  data.to_dict()['y'][DMU]
    #pr += lp.lpSum(var) == 1
    result = pr.solve()
    print(DMU,":",lp.value(thita))
>>>
A : 0.85714286
B : 0.63157895
C : 1.0
D : 1.0
E : 1.0
F : 1.0
G : 1.0
H : 0.75
#超效率CCR模型
import pulp as lp
DMUs = data.index.values
d_var = data.columns.values
var = lp.LpVariable.dict("DMU",data.index.values,lowBound=0)
thita = lp.LpVariable('thita')
for DMU in DMUs:
    pr = lp.LpProblem('CCR',lp.LpMinimize)
    pr += thita
    pr += lp.lpSum(data.to_dict()['x1'][i] * var[i] for i in DMUs if i != DMU) <= thita * data.to_dict()['x1'][DMU]
    pr += lp.lpSum(data.to_dict()['x2'][i] * var[i] for i in DMUs if i != DMU) <= thita * data.to_dict()['x2'][DMU]
    pr += lp.lpSum(data.to_dict()['y'][i] * var[i] for i in DMUs if i != DMU) >=  data.to_dict()['y'][DMU]
    # pr += lp.lpSum(var) == 1
    result = pr.solve()
    print(DMU,":",lp.value(thita))
>>>
A : 0.85714286
B : 0.63157895
C : 1.1428571
D : 1.25
E : 2.0
F : 1.0
G : 1.0
H : 0.75    

2. 面向对象,写入类函数

代码如下(示例):

#CCR模型
import pandas as pd
import numpy as np
import pulp as lp
class DEA_CCR:
    def __init__(self, data, inputs_numbers, outputs_numbers, returns="CRS"):
        """
        data(数据为pd.DataFrame形式); 
        inputs_numbers(输入投入个数); 
        outputs_numbers(输入产出个数);
        returns (输入 "CRS" or "VRS")

        """
        self.data = data
        self.inputs_numbers = inputs_numbers
        self.outputs_numbers = outputs_numbers
        self.returns = returns
        self.DMUs = data.index.values
        self.d_var = data.columns.values

    def _CCR(self):

        self.var = lp.LpVariable.dict("DMU",self.data.index.values,lowBound=0)
        thita = lp.LpVariable('thita')

        dums_dict = {}
        for DMU in self.DMUs:
            pr = lp.LpProblem('CCR',lp.LpMinimize)
            pr += thita
            if self.returns == "VRS":
                pr += lp.lpSum(self.var) == 1
            for j0 in self.d_var[0:self.inputs_numbers]:
                pr += lp.lpSum(self.data.to_dict()[j0][i] * self.var[i] for i in self.DMUs) <= thita * self.data.to_dict()[j0][DMU]
            for r0 in self.d_var[self.inputs_numbers:self.inputs_numbers+self.outputs_numbers]:   
                pr += lp.lpSum(self.data.to_dict()[r0][i] * self.var[i] for i in self.DMUs) >= self.data.to_dict()[r0][DMU]
            pr.solve()
            dums_dict[DMU] = lp.value(pr.objective)
            dums_dict_data = pd.DataFrame.from_dict(dums_dict, orient='index')   
        return dums_dict_data


    def _CCR_super(self):
      
        self.var = lp.LpVariable.dict("DMU",self.data.index.values,lowBound=0)
        thita = lp.LpVariable('thita')
        
        dums_dict = {}
        for DMU in self.DMUs:
            pr = lp.LpProblem('CCR',lp.LpMinimize)
            pr += thita
            if self.returns == "VRS":
                pr += lp.lpSum([self.var[i] for i in DMUs if i!= DMU]) == 1  #超效率的VRS,和非超效率有些许不同
            for j0 in self.d_var[0:self.inputs_numbers]:
                pr += lp.lpSum([self.data.to_dict()[j0][i] * self.var[i] for i in self.DMUs if i != DMU]) <= thita * self.data.to_dict()[j0][DMU]
            for r0 in self.d_var[self.inputs_numbers:self.inputs_numbers+self.outputs_numbers]:   
                pr += lp.lpSum([self.data.to_dict()[r0][i] * self.var[i] for i in self.DMUs if i != DMU]) >= self.data.to_dict()[r0][DMU]
            result = pr.solve()
            dums_dict[DMU] = lp.value(pr.objective)
            dums_dict_data = pd.DataFrame.from_dict(dums_dict, orient='index')   
        return dums_dict_data  

3. 测试

代码如下(示例):

C_data = DEA_CCR(data,2,2,returns = "VRS")._CCR()
C_data
>>>
A	0.857143
B	0.631579
C	1.000000
D	1.000000
E	1.000000
F	1.000000
G	1.000000
H	0.750000
C_data = DEA_CCR(data,2,2,returns = "VRS")._CCR_super()
C_data
>>>
A	0.857143
B	0.631579
C	1.142857
D	1.250000
E	2.000000
F	1.000000
G	1.000000
H	0.750000

总结

  1. 先对基础公式进行建模测试
  2. 然后放入类打包

您可能感兴趣的与本文相关的镜像

Python3.8

Python3.8

Conda
Python

Python 是一种高级、解释型、通用的编程语言,以其简洁易读的语法而闻名,适用于广泛的应用,包括Web开发、数据分析、人工智能和自动化脚本

对于城市群多年低碳发展效率研究,在选择效率SBM模型的规模报酬假设(BCC vs. CCR)和计算流程时,需结合研究对象的经济学特征和研究目标进行决策。 ### 一、BCCCCR模型的选择分析 选择的核心在于判断城市群内的各城市是否在“最优规模”下运营,以及你是否希望将“规模”因素从效率评价中剥离。下表对比了两种假设在本研究场景下的适用性: | 考量维度 | CCR (CRS) 模型 | BCC (VRS) 模型 | **对本研究的建议** | | :--- | :--- | :--- | :--- | | **规模报酬假设** | 不变规模报酬。假设城市投入增加N倍,低碳产出(如GDP)和非期望产出(如碳排放)同比例增减。 | 可变规模报酬。承认城市规模报酬可能递增、递减或不变,更符合城市发展的现实[ref_2]。 | **更倾向于BCC(VRS)**。城市发展受土地、政策、集聚效应等限制,很难始终处于CRS状态。 | | **效率值内涵** | 测算**综合技术效率**,包含了纯技术效率和规模效率。 | 测算**纯技术效率**,剔除了规模因素的影响,反映城市在既定规模下的管理与技术水平[ref_3]。 | 若关注城市**管理与技术**层面的低碳绩效,选BCC;若关注包含规模因素的**整体**低碳表现,选CCR。 | | **前沿面构建** | 前沿面为锥形,可能导致效率投影点指向一个不现实的“大城市”或“微型城市”[ref_1]。 | 前沿面为实际观测城市的凸包,投影点更贴近现实存在的城市组合。 | BCC的前沿面更贴近城市群内实际存在的规模结构,评估结果更稳健。 | | **结果解读** | 效率值低可能源于管理不善或规模不当,难以区分。 | 可结合CCR结果进行分解:`规模效率 = CCR效率 / BCC效率`,从而识别效率损失的根源[ref_4]。 | **强烈建议同时计算CCRBCC效率**。通过分解,可以清晰判断一个城市效率低下,究竟是因为“管理水平差”(纯技术效率低)还是“规模不匹配”(规模效率低)。 | **结论与建议**: 1. **首选方案:同时计算CCRBCC效率值**。这是最严谨的做法。你可以先使用BCC模型计算各城市的**纯技术效率**,评估其低碳管理与技术水平;再用CCR模型计算**综合技术效率**,并通过公式 `规模效率 = CCR效率 / BCC效率` 计算出**规模效率**。这样,你能得到一个三维度的效率分析:综合表现如何?技术管理水平如何?规模是否最优? 2. **若仅选其一,推荐BCC(VRS)模型**。因为城市发展普遍存在规模报酬可变的情况(如初期集聚效应带来报酬递增,后期可能出现拥堵等报酬递减)。BCC模型能更公平地比较不同规模城市在低碳技术和管理上的优劣,避免因城市规模固有差异导致评价偏误[ref_3]。 3. **研究问题导向**:如果你的核心问题是“哪个城市的低碳**发展模式**更优?”,应侧重BCC结果。如果是“哪个城市的低碳**综合绩效**更优?”,则应侧重CCR结果。 ### 二、逐年计算与综合分析流程 “先逐年计算,再综合分析”是处理面板数据的标准且推荐的方法。具体流程如下: **步骤一:逐年截面效率计算** 将每一年的数据视为一个独立的截面数据集,分别应用效率SBM模型(基于你选择的BCCCCR假设)进行计算。这能得到每个城市在各年份的效率值序列($\theta_{it}$,其中 $i$ 代表城市,$t$ 代表年份)。 **代码示例(Python思路 - 基于`pulp`或`cvxopt`)**: ```python import pandas as pd import numpy as np # 假设使用 pulp 求解线性规划 import pulp def calculate_super_sbm_vrs(inputs, good_outputs, bad_outputs, year_data): """ 计算某一年所有DMU的效率SBM-VRS效率值。 参数 year_data: DataFrame, 列包括['city', 'input1', ..., 'good_output1', ..., 'bad_output1', ...] """ num_dmu = len(year_data) efficiency_scores = np.zeros(num_dmu) for k in range(num_dmu): # 对每个被评价城市k # 1. 定义决策变量 # t, s_i_minus, s_r_plus, s_l_bad, Lambda_j (j != k) prob = pulp.LpProblem(f'SuperSBM_VRS_{k}', pulp.LpMinimize) t = pulp.LpVariable('t', lowBound=0) s_minus = [pulp.LpVariable(f's_minus_{i}', lowBound=0) for i in range(m)] s_plus = [pulp.LpVariable(f's_plus_{r}', lowBound=0) for r in range(s1)] s_bad = [pulp.LpVariable(f's_bad_{l}', lowBound=0) for l in range(s2)] lambdas = [pulp.LpVariable(f'lambda_{j}', lowBound=0) for j in range(num_dmu) if j != k] # 2. 构建目标函数 (效率SBM目标函数,此处为示例简化形式) # 目标函数通常为 Min (1/m) * sum(s_i_minus / x_ik) / [1 + (1/(s1+s2)) * (sum(s_r_plus/y_rk) + sum(s_l_bad/b_lk))] # 需线性化处理,此处省略具体转换步骤 # ... # 3. 添加约束条件 # a) 投入约束: sum(Lambda_j * x_ij) + s_i_minus = t * x_ik for i in range(m): prob += pulp.lpSum([lambdas[j_idx] * year_data.iloc[other_idx, input_cols[i]] for j_idx, other_idx in enumerate([idx for idx in range(num_dmu) if idx != k])]) + s_minus[i] == t * year_data.iloc[k, input_cols[i]] # b) 期望产出约束: sum(Lambda_j * y_rj) - s_r_plus = t * y_rk for r in range(s1): prob += pulp.lpSum([lambdas[j_idx] * year_data.iloc[other_idx, good_output_cols[r]] for j_idx, other_idx in enumerate([idx for idx in range(num_dmu) if idx != k])]) - s_plus[r] == t * year_data.iloc[k, good_output_cols[r]] # c) 非期望产出约束: sum(Lambda_j * b_lj) + s_l_bad = t * b_lk for l in range(s2): prob += pulp.lpSum([lambdas[j_idx] * year_data.iloc[other_idx, bad_output_cols[l]] for j_idx, other_idx in enumerate([idx for idx in range(num_dmu) if idx != k])]) + s_bad[l] == t * year_data.iloc[k, bad_output_cols[l]] # d) VRS凸性约束: sum(Lambda_j) = t [这是BCCCCR的关键区别] prob += pulp.lpSum(lambdas) == t # 4. 求解并计算效率值 prob.solve(pulp.PULP_CBC_CMD(msg=False)) # 从求解结果中计算最终的rho (效率值) # efficiency_scores[k] = 1 / (t.varValue) * ... (根据目标函数线性化后的形式计算) # ... return efficiency_scores # 假设df是包含多年份的面板数据DataFrame efficiency_panel = {} for year in df['year'].unique(): year_data = df[df['year'] == year].reset_index(drop=True) efficiency_panel[year] = calculate_super_sbm_vrs(inputs=['资本', '劳动力', '能源'], good_outputs=['GDP'], bad_outputs=['CO2排放'], year_data=year_data) # 将结果整合为一个面板DataFrame result_df = pd.DataFrame(efficiency_panel).T # 行为年份,列为城市 ``` *注意:以上代码为目标函数和约束的框架性展示,效率SBM目标函数的线性化过程较为复杂,需参考完整模型公式实现[ref_6]。* **步骤二:跨期综合分析** 在获得面板效率值后,可采用以下方法进行深入分析: 1. **效率趋势分析**:绘制每个城市效率值随时间变化的折线图,直观观察其低碳效率的演进趋势(上升、下降或波动)。 2. **效率水平比较**:计算整个研究期内各城市效率的均值、标准差,进行排序和分类(如高效、中效、低效城市群)。 3. **收敛性分析**:使用σ收敛或β收敛模型,检验城市群内各城市的低碳效率差异是否随时间推移而缩小(即是否存在“追赶效应”)[ref_1]。 * **σ收敛**:计算各年份所有城市效率值的标准差或变异系数,观察其随时间的变化趋势。 * **β收敛**:建立模型:`ln(θ_i,t / θ_i,t-1) = α + β ln(θ_i,t-1) + ε_i,t`。若β显著为负,则存在绝对β收敛,即初始效率低的城市有更快的增长率。 4. **影响因素分析(第二阶段回归)**:将计算得到的效率值作为被解释变量,建立如Tobit回归模型,探究其与经济发展水平、产业结构、能源结构、环境规制强度等外部因素的关系[ref_5]。 ```stata * 假设 efficiency.dta 中包含变量:city, year, eff_score (效率值), gdp_pc, ind_struc, energy_struc, reg_strength use "efficiency.dta", clear * 使用Tobit模型进行回归,因变量eff_score理论上限可大于1(效率),下限为0 tobit eff_score gdp_pc ind_struc energy_struc reg_strength i.year i.city, ll(0) margins, dydx(*) ``` **为什么不直接使用包含时间维度的全局前沿面(如将多年数据混合)?** * **优点**:将所有年份数据置于同一个前沿面下计算,效率值完全可比,便于直接进行跨期比较。 * **缺点**:忽略了技术进步。假设2010年的最佳实践技术到2020年仍然是最优的,这不符合现实,会高估早期年份的效率,低估后期年份的效率。 * **折中方法**:可考虑使用**窗口DEA**或**序列DEA**,但复杂度更高。对于多数研究,“逐年计算”已能很好地反映效率动态,且更易于理解和实施。 **最终建议操作流程**: 1. **数据准备**:整理成“城市-年份”面板数据格式。 2. **模型选择与计算**: * **方案A(推荐)**:分别运行**效率SBM-BCC**和**效率SBM-CCR**模型,逐年计算,得到两套效率值。 * **方案B**:若只选一个,运行**效率SBM-BCC**模型,逐年计算。 3. **效率分解(如果采用方案A)**:利用两套结果计算`规模效率 = CCR效率 / BCC效率`。 4. **描述性统计与可视化**:分析效率的时序变化和截面差异。 5. **深入分析**:根据研究问题,选择进行收敛性检验或影响因素回归分析。 通过“逐年计算”获得动态效率序列,再通过“综合分析”揭示时空演变规律和驱动因素,这一流程能系统、深入地回答关于城市群低碳发展效率的各类研究问题[ref_4][ref_5]。
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值