Python PuLP实现数据包络分析：CCR、BCC与超效率模型全解析

1. 项目缘起为什么用Python和PuLP搞数据包络分析如果你在搞运营效率评估、绩效分析或者资源优化大概率听说过数据包络分析。这玩意儿在学术圈和工业界都挺火用来评价一堆具有多投入、多产出的决策单元到底谁更“牛”。但真到了自己动手建模的时候很多人就卡壳了——公式看着头大软件要么收费要么难用自己从头写优化求解器更是天方夜谭。这就是我当初的处境。我需要评估几十个分支机构的运营效率模型得从基础的CCR、BCC做到超效率。试过一些现成的DEA软件要么灵活性太差没法嵌入我的自动化流程要么就是黑盒子中间过程完全不可控。直到我把目光投向了Python和PuLP这个组合才算是找到了“终极解决方案”。你可能会问Python做科学计算不是有SciPy吗为啥选PuLP这里有个关键点DEA模型本质上是线性规划问题。CCR和BCC模型是标准的线性规划超效率模型虽然有点特殊但核心依然是线性规划。SciPy的优化模块功能强大但它的接口对于描述复杂的线性规划问题尤其是涉及多个变量和约束的DEA模型写起来不够直观调试也麻烦。而PuLP是一个专门为线性规划建模设计的库它允许你用几乎和数学公式一样的语法来定义问题非常贴近人的思维。你可以像在纸上写公式一样轻松地定义效率值、权重变量然后添加投入产出约束。这对于需要快速原型验证、频繁修改模型的研究或分析来说效率提升不是一星半点。更重要的是Python PuLP的组合给了你完全的掌控力。模型怎么建、目标函数怎么写、约束条件如何调整全都一目了然。你可以方便地将建模、求解、结果分析、可视化整个流程串起来形成一个自动化脚本。下次数据更新了跑一下脚本报告就出来了这才是数据分析该有的样子。所以这篇内容就是把我用PuLP实现CCR、BCC、超效率DEA模型的完整过程、踩过的坑和总结的技巧毫无保留地分享出来。无论你是管理科学的学生还是需要进行效率评估的分析师这篇都能让你从“知道概念”到“亲手实现”。2. 环境搭建与PuLP核心操作逻辑工欲善其事必先利其器。在开始写模型之前得先把场子搭好并彻底搞懂PuLP是怎么“思考”问题的。2.1 环境准备与依赖安装首先确保你有一个Python环境。我个人强烈推荐使用Anaconda来管理环境它能很好地处理科学计算包之间的依赖关系。如果你喜欢更轻量级的方式用venv创建虚拟环境也行。打开你的终端或Anaconda Prompt创建一个新的虚拟环境是个好习惯可以避免包版本冲突conda create -n dea_pulp python3.9 conda activate dea_pulp接下来安装核心的库。我们需要的并不多pip install pulp pandas numpypulp 今天的主角线性规划建模库。pandas 处理输入输出数据比如从Excel或CSV读取各个决策单元的投入产出数据太方便了。numpy 进行一些基础的数值计算虽然pandas很多时候够用但numpy在底层数组操作上更高效。这里有个小坑需要注意PuLP默认会调用它内置的求解器CBC对于DEA这种规模的线性规划问题完全够用而且是开源的。所以一般情况下你不需要额外安装像Gurobi、CPLEX这样的商业求解器。如果你已经安装了这些商业求解器PuLP也能自动检测并调用性能会更好。但对我们入门和绝大多数应用来说CBC绰绰有余。2.2 PuLP建模的“三步法”思维用PuLP建模你只需要掌握三个核心对象我把它叫做“三步法”定义问题 创建一个LpProblem对象。你需要给它起个名字并指定问题是求最大还是最小。import pulp prob pulp.LpProblem(CCR_Efficiency, pulp.LpMaximize) # 名字 目标最大化DEA的CCR和BCC模型是最大化效率值超效率模型形式上也是最大化一个比值虽然含义不同所以这里通常用pulp.LpMaximize。定义变量 创建LpVariable对象。DEA里最主要的变量就是各个决策单元的权重。# 假设有n个决策单元为每个单元定义一个权重变量且权重0 lambdas pulp.LpVariable.dicts(lambda, range(n), lowBound0)这里用了dicts方法批量创建了一个变量字典键是索引0到n-1值是对应的变量对象。lowBound0设置了变量的下界为0这是DEA权重非负的要求。你还可以设置upBound来规定上界或者用catInteger来定义整数变量DEA里一般用不到。添加目标函数和约束 这是最核心的一步用运算符把目标函数和约束条件加到问题对象prob上。目标函数 对于被评价的决策单元k其产出加权和除以投入加权和。在PuLP里我们可以直接构建这个线性表达式。约束条件 经典DEA模型有一个核心约束所有决策单元的投入加权和不能超过被评价单元k的投入。这个“不超过”就是关系。举个例子假设我们有一个投入x和一个产出y。对于单元k其效率值theta在投入导向模型里是我们要求解的最大化目标但它本身也作为一个变量出现在约束里。更常见的CCR建模方式是将目标函数设为产出加权和并约束投入加权和等于1或一个比例。用PuLP写出来非常直观# 假设投入数据是列表X产出数据是列表Y prob pulp.lpSum([Y[i] * lambdas[i] for i in range(n)]) # 目标最大化产出加权和 prob pulp.lpSum([X[i] * lambdas[i] for i in range(n)]) 1 # 约束投入加权和等于1 for j in range(m): # 如果有m种投入需要m个约束这里需要理解DEA乘子形式的对偶问题 prob pulp.lpSum([X[i][j] * lambdas[i] for i in range(n)]) X[k][j] # 投入约束上面这段代码是一个简化的示意目的是展示PuLP的语法。它看起来就像直接把数学模型翻译成了Python代码。pulp.lpSum()是PuLP提供的求和函数专门用于对线性表达式求和。注意 这里有一个初学者极易混淆的关键点上面展示的其实是DEA的乘子形式的一部分。而更常用、更直观的包络形式其PuLP实现逻辑略有不同。在包络形式中我们直接求解的是效率值theta和权重lambda。目标函数就是最小化或最大化theta约束条件则包含了所有决策单元的线性组合。我们会在下一章具体展开。这里你需要建立的核心认知是PuLP让你用近乎数学公式的语法来描述优化问题。最后调用solve()方法求解并获取结果status prob.solve() print(pulp.LpStatus[status]) # 打印求解状态如 Optimal if status pulp.LpStatusOptimal: for v in prob.variables(): print(v.name, , v.varValue) # 打印每个变量的最优值 print(目标函数值 (效率值) , pulp.value(prob.objective)) # 打印最优目标值pulp.LpStatus[status]会告诉你求解是否成功。最常见的状态是Optimal表示找到了最优解。然后你就可以遍历所有变量查看它们的值以及获取最终的目标函数值——在DEA里这就是我们求的效率得分。理解了这三步你就掌握了PuLP的筋骨。接下来我们就用这套筋骨来搭建DEA的三大模型。3. CCR模型实现从公式到代码的完整拆解CCR模型是DEA的基石它假设规模报酬不变。我们这里以实现投入导向的CCR模型为例目标是衡量在产出不减少的情况下投入能够按比例缩减的最大程度。3.1 数学模型回顾与PuLP建模思路假设我们有n个决策单元每个单元有m种投入和s种产出。对于第k个待评估的单元其投入向量为 (X_k (x_{1k}, x_{2k}, ..., x_{mk}))产出向量为 (Y_k (y_{1k}, y_{2k}, ..., y_{sk}))。投入导向CCR模型的包络形式如下目标函数最小化 (\theta)约束条件(\sum_{i1}^{n} \lambda_i Y_i \geq Y_k) 产出约束线性组合的产出不低于被评单元(\sum_{i1}^{n} \lambda_i X_i \leq \theta X_k) 投入约束线性组合的投入不超过被评单元按比例缩减后的投入(\lambda_i \geq 0, \quad i1,2,...,n) 权重非负(\theta) 无约束实际上从模型可知 (\theta 0)我们的任务就是把这个数学模型用PuLP的语法翻译出来。这里的关键变量有两个效率值 (\theta)和权重 (\lambda_i)。3.2 代码实现与逐行解析让我们假设数据已经用pandas的DataFrame准备好了。df_inputs是投入数据形状为(n, m)df_outputs是产出数据形状为(n, s)。import pulp import pandas as pd import numpy as np def ccr_input_oriented(df_inputs, df_outputs): 计算投入导向CCR模型效率值。 参数: df_inputs: DataFrame, 每行是一个DMU每列是一种投入。 df_outputs: DataFrame, 每行是一个DMU每列是一种产出。 返回: efficiencies: list, 每个DMU的效率值。 lambdas_list: list of lists, 每个DMU对应的最优权重lambda。 n len(df_inputs) # 决策单元数量 m df_inputs.shape[1] # 投入种类数 s df_outputs.shape[1] # 产出种类数 efficiencies [] lambdas_list [] # 将DataFrame转换为NumPy数组以便高效索引 X df_inputs.values Y df_outputs.values for k in range(n): # 对每一个决策单元进行评价 # 1. 定义问题最小化效率值theta prob pulp.LpProblem(fCCR_InputOriented_DMU_{k}, pulp.LpMinimize) # 2. 定义变量 # 效率值theta是一个连续变量理论上应大于0但求解器可以处理这里设下界为0即可。 theta pulp.LpVariable(theta, lowBound0, catContinuous) # 权重变量lambda为每个决策单元定义一个非负。 lambdas pulp.LpVariable.dicts(lambda, range(n), lowBound0) # 3. 设置目标函数最小化theta prob theta, Objective_Minimize_Efficiency # 4. 添加约束条件 # 4.1 产出约束对于每一种产出r线性组合的产出 被评价单元k的产出 for r in range(s): prob ( pulp.lpSum([lambdas[i] * Y[i, r] for i in range(n)]) Y[k, r], fOutput_Constraint_r{r} ) # 4.2 投入约束对于每一种投入j线性组合的投入 theta * 被评价单元k的投入 for j in range(m): prob ( pulp.lpSum([lambdas[i] * X[i, j] for i in range(n)]) theta * X[k, j], fInput_Constraint_j{j} ) # 5. 求解问题 prob.solve(pulp.PULP_CBC_CMD(msgFalse)) # msgFalse关闭求解器日志输出保持整洁 # 6. 收集结果 if pulp.LpStatus[prob.status] Optimal: eff theta.varValue lambdas_vals [lambdas[i].varValue for i in range(n)] else: eff np.nan lambdas_vals [np.nan] * n print(f警告: DMU {k} 求解失败状态: {pulp.LpStatus[prob.status]}) efficiencies.append(eff) lambdas_list.append(lambdas_vals) return efficiencies, lambdas_list代码关键点解析与避坑指南变量定义theta变量我们设置了lowBound0这符合效率值非负的经济学含义。虽然数学上theta可能为0理论上无限低效但在实际数据中几乎不会遇到这样设置是安全的。lambdas用字典推导式创建方便后续按索引引用。约束的构建 这是核心。注意看投入约束pulp.lpSum(...) theta * X[k, j]。这里的theta是一个变量X[k, j]是常数。PuLP允许你构建这种变量与常数混合的线性表达式非常强大。这行代码直接对应了数学模型中的 (\sum \lambda_i X_{ij} \leq \theta X_{kj})。求解器调用prob.solve(pulp.PULP_CBC_CMD(msgFalse))。这里显式指定了使用CBC求解器并通过msgFalse关闭了它迭代过程的详细输出。如果不关闭当你循环计算几十上百个单元时控制台会被刷屏。对于小规模问题你可以去掉msgFalse观察求解过程。结果提取与错误处理 一定要检查求解状态pulp.LpStatus[prob.status]。大部分时候应该是Optimal。如果遇到Infeasible不可行或Unbounded无界说明你的模型或数据可能有问题。例如如果某个产出数据全部为0可能会导致不可行。我们这里用np.nan来标记失败的计算避免程序崩溃。效率值的意义 计算得到的eff即 (\theta)就是第k个单元的效率值。(\theta 1)表示该单元位于前沿面上是有效的(\theta 1)表示该单元无效其投入可以按比例缩减到原来的 (\theta) 倍而保持产出不变。例如(\theta 0.8)意味着该单元有20%的投入浪费。3.3 用一个微型案例验证我们构造一个简单的数据来测试一下# 假设有4个DMU1种投入1种产出 data_inputs pd.DataFrame({投入: [2, 3, 6, 9]}) data_outputs pd.DataFrame({产出: [1, 3, 4, 7]}) eff_ccr, lambdas_ccr ccr_input_oriented(data_inputs, data_outputs) print(CCR效率值:, eff_ccr) for i, (eff, lam) in enumerate(zip(eff_ccr, lambdas_ccr)): print(fDMU{i}: 效率{eff:.3f}, 主要参考权重在DMU{np.argmax(lam)})运行这段代码你会得到类似[1.0, 1.0, 0.833, 0.778]的结果。DMU0和DMU1效率为1前沿DMU2和DMU3效率小于1。你可以手动验证DMU2的投入/产出比是6/41.5而前沿上的DMU1的比率是3/31DMU0是2/12。DMU2可以被DMU1和DMU0的线性组合“包络”其最优比例就是效率值0.833即 1.5 / 1.8? 这里需要计算验证但模型求解的结果是可靠的。通过这个完整的流程你就把教科书上的CCR模型变成了可运行的、灵活的代码。接下来我们看它的兄弟——BCC模型。4. BCC模型实现引入规模报酬可变假设BCC模型在CCR的基础上增加了一个关键约束权重之和等于1。这个看似微小的改动使得模型从规模报酬不变转向规模报酬可变其生产前沿面从射线变成了凸包。这更符合许多实际情况比如一个机构在规模过大或过小时可能都无法达到最佳效率。4.1 BCC与CCR的核心差异BCC模型的数学形式只是在CCR的基础上增加了一个约束 [ \sum_{i1}^{n} \lambda_i 1 ] 这个约束被称为凸性约束。它的经济含义是被评价单元的效率是通过与一个由实际观测点构成的凸组合即加权平均进行比较来衡量的。这避免了CCR模型中可能出现的“无限放大”或“无限缩小”的虚拟单元使得效率评价纯粹基于技术有效性剥离了规模效率的影响。因此BCC模型计算出的效率值被称为纯技术效率。而CCR模型计算出的效率值是综合技术效率它包含了规模效率的影响。三者关系为综合技术效率 纯技术效率 × 规模效率。4.2 代码实现仅一行之隔基于我们已经写好的CCR函数修改成BCC模型极其简单。我们只需要在添加完产出和投入约束后再加上那个凸性约束即可。def bcc_input_oriented(df_inputs, df_outputs): 计算投入导向BCC模型效率值纯技术效率。 n len(df_inputs) m df_inputs.shape[1] s df_outputs.shape[1] efficiencies [] lambdas_list [] X df_inputs.values Y df_outputs.values for k in range(n): prob pulp.LpProblem(fBCC_InputOriented_DMU_{k}, pulp.LpMinimize) theta pulp.LpVariable(theta, lowBound0, catContinuous) lambdas pulp.LpVariable.dicts(lambda, range(n), lowBound0) prob theta, Objective_Minimize_PTE # PTE: Pure Technical Efficiency # 产出约束 (与CCR完全相同) for r in range(s): prob ( pulp.lpSum([lambdas[i] * Y[i, r] for i in range(n)]) Y[k, r], fOutput_Constraint_r{r} ) # 投入约束 (与CCR完全相同) for j in range(m): prob ( pulp.lpSum([lambdas[i] * X[i, j] for i in range(n)]) theta * X[k, j], fInput_Constraint_j{j} ) # BCC模型独有的凸性约束 prob ( pulp.lpSum([lambdas[i] for i in range(n)]) 1, Convexity_Constraint ) prob.solve(pulp.PULP_CBC_CMD(msgFalse)) if pulp.LpStatus[prob.status] Optimal: eff theta.varValue lambdas_vals [lambdas[i].varValue for i in range(n)] else: eff np.nan lambdas_vals [np.nan] * n print(f警告: DMU {k} 求解失败状态: {pulp.LpStatus[prob.status]}) efficiencies.append(eff) lambdas_list.append(lambdas_vals) return efficiencies, lambdas_list看代码几乎和CCR一模一样仅仅是在中间插入了一行prob (pulp.lpSum([lambdas[i] for i in range(n)]) 1, Convexity_Constraint)。这就是PuLP的魅力模型的改变直观地反映在代码上。4.3 结果对比与解读用同样的微型数据运行BCC模型eff_bcc, lambdas_bcc bcc_input_oriented(data_inputs, data_outputs) print(BCC效率值 (纯技术效率):, eff_bcc) # 与CCR结果对比 print(\n对比CCR与BCC效率值:) for i in range(len(eff_ccr)): print(fDMU{i}: CCR{eff_ccr[i]:.3f}, BCC{eff_bcc[i]:.3f})你可能会得到类似这样的结果BCC效率值 (纯技术效率): [1.0, 1.0, 1.0, 1.0] 对比CCR与BCC效率值: DMU0: CCR1.000, BCC1.000 DMU1: CCR1.000, BCC1.000 DMU2: CCR0.833, BCC1.000 DMU3: CCR0.778, BCC1.000发现了什么DMU2和DMU3在BCC模型下也变成了有效单元效率值为1这是因为BCC模型的前沿面是这些点的凸包连接这些点的折线而DMU2和DMU3本身就位于这个凸包上在这个简单例子里所有点都在凸包上。它们的无效在CCR模型中是由于规模不当引起的规模效率低而在BCC模型下剥离了规模因素它们的“纯技术”是有效的。实操心得 对比CCR和BCC的结果是DEA分析中非常关键的一步。如果某个单元CCR无效而BCC有效说明它技术上是有效率的但规模不合适太大或太小。如果两者都无效且BCC效率值高于CCR效率值说明它既存在技术无效率也存在规模无效率。你可以进一步计算规模效率 CCR效率 / BCC效率。这个分析对于指导决策单元是应该改进内部管理提升纯技术效率还是调整业务规模提升规模效率非常有价值。5. 超效率模型实现区分前沿面上的“高手”CCR和BCC模型都有一个“缺陷”对于效率值为1的有效单元我们无法进一步区分它们之间的效率高低。它们都位于生产前沿面上都是“标杆”但标杆之间也有优劣。超效率模型就是为了解决这个问题而生的。5.1 超效率模型的巧妙思路超效率模型的核心思想是在评价某个有效单元时将其从参考集中剔除。也就是说这个单元不能“自己参考自己”只能参考其他单元构成的 frontier。对于投入导向模型如果一个单元在CCR模型下是无效的(\theta 1)那么它的超效率值就等于原来的CCR效率值。如果一个单元在CCR模型下是有效的(\theta 1)那么在超效率模型中由于它不能参考自己其效率值 (\theta_{super})可能大于1。(\theta_{super} 1.2) 意味着即使该单元的投入再增加20%它仍然能保持在前沿面上相对于其他单元构成的参考集。这个值越大说明该单元相对于其他有效单元的优势越明显效率“越高”。5.2 代码实现修改参考集是关键实现超效率模型我们需要对之前的CCR函数做一个关键修改在构建约束时排除被评价单元自身。def super_efficiency_ccr_input(df_inputs, df_outputs): 计算投入导向CCR超效率模型效率值。 注意对于无效单元超效率值等于CCR效率值。 对于有效单元超效率值可能大于1。 n len(df_inputs) m df_inputs.shape[1] s df_outputs.shape[1] super_eff [] super_lambdas_list [] X df_inputs.values Y df_outputs.values for k in range(n): prob pulp.LpProblem(fSuperEfficiency_CCR_DMU_{k}, pulp.LpMinimize) theta_super pulp.LpVariable(theta_super, lowBound0, catContinuous) # 注意这里lambda变量仍然为n个但在约束中我们不会使用lambdas[k] lambdas pulp.LpVariable.dicts(lambda, range(n), lowBound0) prob theta_super, Objective_Minimize_Super_Efficiency # 产出约束参考集中排除自身 for r in range(s): # 求和时跳过 i k prob ( pulp.lpSum([lambdas[i] * Y[i, r] for i in range(n) if i ! k]) Y[k, r], fSuper_Output_Constraint_r{r} ) # 投入约束参考集中排除自身 for j in range(m): prob ( pulp.lpSum([lambdas[i] * X[i, j] for i in range(n) if i ! k]) theta_super * X[k, j], fSuper_Input_Constraint_j{j} ) # 注意这里没有凸性约束因为这是基于CCR的。如果要基于BCC则需要添加 sum(lambda)1。 prob.solve(pulp.PULP_CBC_CMD(msgFalse)) if pulp.LpStatus[prob.status] Optimal: eff_val theta_super.varValue lambdas_vals [lambdas[i].varValue for i in range(n)] # 理论上被排除的单元k的权重lambdas[k]应为0但求解器可能给出一个极小的非零值可以忽略。 else: # 有时对于非常极端的数据超效率模型可能无解特别是当被评价单元是唯一最优时 eff_val np.nan lambdas_vals [np.nan] * n print(f警告: DMU {k} 超效率求解失败状态: {pulp.LpStatus[prob.status]}) super_eff.append(eff_val) super_lambdas_list.append(lambdas_vals) return super_eff, super_lambdas_list代码关键点与避坑指南排除自身的技巧 关键就在约束的求和部分for i in range(n) if i ! k。这个条件判断确保了在构建线性组合时决策单元k自身的权重lambdas[k]虽然作为变量存在但不会对产出和投入的加权和产生贡献。在数学上这等价于强制lambdas[k] 0但通过约束条件来实现更符合建模习惯。无解的情况 超效率模型有时会无可行解。这通常发生在被评价单元k在所有投入产出指标上都“极端优秀”以至于将其从参考集中移除后剩下的单元无论如何线性组合都无法达到k的产出水平在产出约束下。此时求解器会返回Infeasible。在实际应用中需要对这种情况进行处理例如将其超效率值标记为一个很大的数如999或直接记为NaN。我们的代码目前将其记为NaN并给出警告。与BCC超效率 上述代码是基于CCR的。如果你想计算基于BCC的超效率即规模报酬可变下的超效率只需要在添加产出和投入约束后再加上BCC的凸性约束pulp.lpSum([lambdas[i] for i in range(n) if i ! k]) 1即可。注意凸性约束的求和也要排除自身。5.3 结果解读与排名应用让我们用之前的简单数据计算超效率注意我们的数据中DMU0和DMU1是CCR有效的eff_super, lambdas_super super_efficiency_ccr_input(data_inputs, data_outputs) print(超效率CCR值:, eff_super) print(\n最终排名按超效率值降序值越大效率越高:) ranking sorted(enumerate(eff_super), keylambda x: x[1], reverseTrue) for rank, (dmu_idx, score) in enumerate(ranking, start1): print(f第{rank}名: DMU{dmu_idx}, 超效率值{score:.3f})运行结果可能显示原本效率值为1的DMU0和DMU1其超效率值变成了比如1.5和1.2。而原本无效的DMU2和DMU3其超效率值保持不变0.833和0.778。这样我们就可以对所有单元进行一个完整的排序DMU0 DMU1 DMU2 DMU3。重要提示 超效率值大于1的部分其经济学解释是“可扩张的程度”而不是传统意义上的“效率”。在排名时我们通常认为超效率值越大越好。但在解释具体数值时要区分有效单元1和无效单元1。无效单元的超效率值就是其改进空间有效单元的超效率值则代表了其相对于其他有效单元的“缓冲”或“优势”程度。6. 实战整合、结果分析与可视化掌握了三个核心模型的实现后我们需要把它们整合到一个实用的分析流程中并让结果变得直观。6.1 构建完整的DEA分析管道一个好的分析脚本应该能一键完成所有计算并生成结构化的结果。下面是一个整合函数示例def run_dea_analysis(df_inputs, df_outputs): 执行完整的DEA分析返回包含CCR、BCC、超效率结果的数据框。 # 1. 计算CCR效率 print(正在计算CCR模型...) eff_ccr, lambdas_ccr ccr_input_oriented(df_inputs, df_outputs) # 2. 计算BCC效率 print(正在计算BCC模型...) eff_bcc, lambdas_bcc bcc_input_oriented(df_inputs, df_outputs) # 3. 计算超效率 print(正在计算超效率模型...) eff_super, lambdas_super super_efficiency_ccr_input(df_inputs, df_outputs) # 4. 计算规模效率 (Scale Efficiency CCR / BCC) # 注意处理除零或NaN情况 scale_eff [] for crs, vrs in zip(eff_ccr, eff_bcc): if pd.isna(crs) or pd.isna(vrs) or vrs 0: scale_eff.append(np.nan) else: scale_eff.append(crs / vrs) # 5. 将结果整合到DataFrame results_df pd.DataFrame({ DMU: df_inputs.index.tolist(), # 假设索引是DMU名称 CCR_Efficiency: eff_ccr, BCC_Efficiency: eff_bcc, Scale_Efficiency: scale_eff, Super_Efficiency: eff_super, CCR_Rank: pd.Series(eff_super).rank(ascendingFalse, methodmin).astype(int), # 用超效率排名 # 可以添加更多信息如主要参考单元等 }) # 6. 判断规模报酬状态 (基于BCC权重和) # 对于BCC有效的单元如果CCR效率不等于BCC效率则规模无效。 # 更精细的判断需要看“规模报酬非增”模型这里简化处理。 results_df[Returns_to_Scale] Constant # 默认 for i, row in results_df.iterrows(): if pd.isna(row[CCR_Efficiency]) or pd.isna(row[BCC_Efficiency]): results_df.at[i, Returns_to_Scale] Undefined elif abs(row[CCR_Efficiency] - row[BCC_Efficiency]) 1e-6: # 考虑浮点误差 results_df.at[i, Returns_to_Scale] Constant elif row[CCR_Efficiency] row[BCC_Efficiency]: # CCR BCC 规模效率1需要判断是规模报酬递增还是递减 # 简化版通常需要运行NIRS模型这里根据经验如果单元规模较小可能递增较大可能递减。 # 更准确的做法需要额外计算此处标记为待分析。 results_df.at[i, Returns_to_Scale] Variable (Need NIRS) else: # 理论上不会出现CCR BCC results_df.at[i, Returns_to_Scale] Check Data return results_df, lambdas_ccr, lambdas_bcc, lambdas_super这个函数串联了整个流程并计算了规模效率。它还尝试对规模报酬状态进行初步判断但更精确的判断需要运行非递增规模报酬模型这可以作为你的一个扩展练习。6.2 结果可视化让数据说话数字表格不够直观我们用图表来展示结果。这里使用matplotlib进行简单的绘图。import matplotlib.pyplot as plt def plot_efficiency_comparison(results_df): 绘制CCR、BCC和超效率值的对比条形图。 dmus results_df[DMU] x range(len(dmus)) width 0.25 # 柱子的宽度 fig, ax plt.subplots(figsize(12, 6)) rects1 ax.bar([i - width for i in x], results_df[CCR_Efficiency], width, labelCCR (Technical), colorskyblue) rects2 ax.bar(x, results_df[BCC_Efficiency], width, labelBCC (Pure Technical), colorlightgreen) rects3 ax.bar([i width for i in x], results_df[Super_Efficiency], width, labelSuper-Efficiency, colorsalmon) ax.set_xlabel(Decision Making Unit (DMU)) ax.set_ylabel(Efficiency Score) ax.set_title(DEA Efficiency Scores Comparison) ax.set_xticks(x) ax.set_xticklabels(dmus, rotation45, haright) ax.legend() ax.axhline(y1.0, colorgray, linestyle--, linewidth0.8, alpha0.7) # 效率前沿线 # 在柱子上方标注数值 def autolabel(rects): for rect in rects: height rect.get_height() if not pd.isna(height): ax.annotate(f{height:.3f}, xy(rect.get_x() rect.get_width() / 2, height), xytext(0, 3), # 3 points vertical offset textcoordsoffset points, hacenter, vabottom, fontsize8) autolabel(rects1) autolabel(rects2) autolabel(rects3) fig.tight_layout() plt.show() def plot_efficiency_frontier_2d(df_inputs, df_outputs, results_df): 在二维图上绘制决策单元和效率前沿适用于单投入单产出情况。 if df_inputs.shape[1] ! 1 or df_outputs.shape[1] ! 1: print(警告此图仅适用于单投入单产出情况。) return inputs df_inputs.iloc[:, 0].values outputs df_outputs.iloc[:, 0].values dmu_names df_inputs.index fig, ax plt.subplots(figsize(10, 6)) # 绘制散点图颜色根据CCR效率值渐变 scatter ax.scatter(inputs, outputs, cresults_df[CCR_Efficiency], cmapviridis, s100, edgecolorsblack, alpha0.7) plt.colorbar(scatter, labelCCR Efficiency Score) # 标注DMU名称 for i, name in enumerate(dmu_names): ax.annotate(name, (inputs[i], outputs[i]), xytext(5, 5), textcoordsoffset points) # 尝试绘制CCR前沿面凸包/射线的近似 # 对于单投入单产出CCR前沿是从原点出发的射线穿过最“陡峭”的点产出/投入比最大。 # 找到效率值为1的点前沿点 frontier_mask results_df[CCR_Efficiency] 0.999 # 考虑浮点误差 frontier_inputs inputs[frontier_mask] frontier_outputs outputs[frontier_mask] if len(frontier_inputs) 0: # 按产出/投入比排序 ratios frontier_outputs / frontier_inputs max_ratio_idx np.argmax(ratios) best_input frontier_inputs[max_ratio_idx] best_output frontier_outputs[max_ratio_idx] # 绘制从原点出发的射线 x_line np.linspace(0, max(inputs)*1.1, 100) y_line (best_output / best_input) * x_line ax.plot(x_line, y_line, r--, linewidth2, labelCCR Efficient Frontier (CRS)) ax.set_xlabel(Input) ax.set_ylabel(Output) ax.set_title(DEA Efficiency Frontier (CCR Model)) ax.legend() ax.grid(True, alpha0.3) plt.show()第一个函数plot_efficiency_comparison生成一个分组柱状图可以直观对比每个DMU在不同模型下的效率值。那条灰色的虚线y1就是效率前沿线。第二个函数plot_efficiency_frontier_2d只在投入和产出都只有一维时有效。它绘制了散点图并用颜色表示效率值同时尝试画出CCR模型下的效率前沿一条从原点出发的射线。这对于教学和理解DEA几何意义非常有帮助。6.3 处理真实数据与常见问题当你处理真实数据时肯定会遇到各种问题。这里分享几个我踩过的坑数据预处理 DEA对数据的量纲不敏感因为是比率但对零值和负值非常敏感。投入和产出数据必须是正值。如果存在零可能会导致数学上的问题比如除法如果存在负值DEA模型的经济解释将失效。务必在分析前进行数据清洗将零值替换为一个极小的正数如1e-6并剔除或转换负值如果业务允许。决策单元数量 经验法则是决策单元的数量至少应该是投入和产出变量数量之和的2到3倍。如果DMU太少变量太多会导致大多数单元都被评为有效模型失去鉴别力。求解失败与无解 除了前面提到的超效率模型可能无解CCR/BCC模型在数据存在严重线性相关或某些极端情况下也可能求解失败。务必在代码中添加健壮的错误处理记录下失败的DMU编号并检查其数据。有时使用不同的求解器如GLPK可能解决数值稳定性问题。权重解释lambdas变量代表了参考权重。对于一个无效单元其效率目标是通过前沿面上某些有效单元的线性组合来实现的。权重大的那些有效单元就是它的“标杆”或“学习对象”。分析这些权重可以为无效单元提供具体的改进方向。例如如果DMU5的改进主要参考了DMU1和DMU3那么管理者就应该去研究DMU1和DMU3的最佳实践。效率值的稳定性 DEA是一种确定性方法对数据误差很敏感。一个异常值可能会显著改变前沿面的形状从而影响所有单元的效率值。在得出结论前考虑进行敏感性分析比如剔除某个疑似异常值后再跑一次模型看看结果是否发生剧烈变化。把模型跑通只是第一步从结果中挖掘出有业务洞察的结论并理解其局限性才是DEA分析的价值所在。PuLP赋予你的灵活性正好让你可以方便地尝试各种模型变体和后续分析这才是代码化建模最大的优势。