做效率评价这些年我最大的感受是模型本身并不难真正让人头秃的是“数据从哪来、怎么洗干净、结果怎么让业务方看懂”。SBMSlacks-Based Measure松弛变量测度模型就是这样——原理不复杂但如果你按教材里的公式硬撸代码十有八九会卡在数据预处理和求解器配置上。这篇文章我就把从数据清洗到结果可视化的完整流程拆开讲附带能直接跑的Python代码适合正在做企业效率评价、区域经济绩效分析、供应链绩效评估的读者参考。1. 先搞清楚SBM模型到底在算什么1.1 从DEA家族说起为什么选SBM而非传统CCR/BCC数据包络分析DEA是一大类基于线性规划的效率评价方法最早由Charnes、Cooper和Rhodes在1978年提出CCR模型后来Banker等人扩展出BCC模型。这类方法的核心理念是不需要预设生产函数直接用“投入产出数据”构造生产前沿面——所有决策单元DMU比如企业、医院、城市中在同样投入下产出最高、或者在同样产出下投入最少的那些DMU就构成效率前沿其他DMU与前沿面的距离就是效率损失。传统CCR和BCC模型有一个共同特点它们属于“径向”模型即假设所有投入或产出按同一比例压缩或扩张。但实际生产过程中不可能所有投入都同比例减少——某台设备定了就是定了人员调整有刚性能耗削减空间和原材料削减空间也完全不同。这时候径向模型算出来的效率值就会偏大而且没办法告诉你“具体哪个投入冗余了多少”。SBM模型正是针对这个问题提出来的它把松弛变量直接放进目标函数允许各投入不均等地缩减同时给出每个指标的松弛量信息量比传统DEA模型高出一个量级。1.2 松弛变量的本质SBM到底在优化什么松弛变量这个名字听起来抽象其实理解起来非常简单。假设你有20家制造企业每家投入固定资产、员工人数、能源消耗产出是营业收入和利润。现在把20家企业的数据放到一起找一个“标杆”在同样投入水平下产出最高的那家就是效率前沿。对某个效率不足的企业来说它的投入松弛量就是“相对于标杆多用了多少固定资产、多少人、多少能源”产出松弛量则是“相对于标杆少产出了多少利润”。SBM模型把这些松弛量放进目标函数里求最小化。以投入导向为例目标函数是ρ 1 - (1/m) * Σ(sᵢ⁻ / xᵢ₀)其中m是投入指标个数sᵢ⁻是第i个投入的松弛量xᵢ₀是当前评价DMU的第i个投入值。这个公式的意思是如果这家企业各投入相对标杆都有冗余ρ就会小于1如果没有任何松弛ρ等于1说明它已经在生产前沿面上。由于分母是各投入的实际值松弛量被“相对化”了所以SBM模型是单位不变的——你用“万元”还是“亿元”表示固定资产计算结果完全一样。1.3 非期望产出怎么进模型做效率评价很多场景绕不开一个现实问题非期望产出。比如制造业的CO₂排放、废水排放银行的坏账率医院的再入院率。这些指标“越少越好”但又不能简单地放进“投入”里处理从生产逻辑上讲排放是产出端的伴随物不是先期投入的资源。Tone在2003年提出了处理非期望产出的SBM扩展模型核心思路是把产出再拆成期望产出和非期望产出两类分别计算松弛量分母部分同时考虑“期望产出不足”和“非期望产出过多”两项。不过我得说句实话非期望产出的SBM虽然理论严谨但求解比基础版复杂需要做Charnes-Cooper变换把分式规划转成线性规划。实际项目中如果只是做初步筛选也可以用“把非期望产出当成投入”的近似做法——虽然这个方法在学术上有争议但在业务场景里速度快、结果也基本符合直觉。两种方式我在后文都会给出代码实现。2. 开工前准备环境、依赖与实验数据2.1 Python环境与依赖库清单我做这个流程用的Python 3.10只依赖四个库全是免费开源的库名用途说明numpy数值计算、矩阵操作数据处理的基础pandas数据清洗、表格操作读Excel、去重、缺失值处理都靠它scipy线性规划求解核心的linprog函数matplotlib seaborn可视化画排名图、热力图、散点图安装命令就一条pip install pandas numpy scipy matplotlib seaborn如果你在公司内网可能需要配置国内镜像源比如清华源、阿里源速度会快很多。IDE我建议用VS CodePython插件装好之后配一下解释器路径就能直接用轻量、启动快配置方法网上教程很多这里不展开。2.2 实验数据设计让案例能复现为了把整个流程讲清楚我构造了一份仿真数据。假设要评估20家制造企业的运营效率决策单元DMU编号从D01到D20指标体系如下投入3个固定资产亿元、员工人数百人、能源消耗万吨标准煤期望产出2个营业收入亿元、净利润亿元非期望产出1个CO₂排放量万吨数据用numpy生成设置随机种子保证可复现。我这里有意给部分企业加入了一些“脏数据”——比如D05的固定资产出现负值D12的员工人数是缺失值D17的营业收入异常偏高——这样后面数据清洗环节才有东西可讲。import numpy as np import pandas as pd np.random.seed(42) n_dmu 20 # 生成基础数据 fixed_assets np.random.uniform(5, 50, n_dmu) # 固定资产亿元 employees np.random.uniform(100, 1000, n_dmu) # 员工人数百人 energy np.random.uniform(10, 100, n_dmu) # 能源消耗万吨标准煤 revenue np.random.uniform(20, 200, n_dmu) # 营业收入亿元 profit np.random.uniform(1, 30, n_dmu) # 净利润亿元 co2 np.random.uniform(5, 60, n_dmu) # CO2排放量万吨 df pd.DataFrame({ DMU: [fD{i1:02d} for i in range(n_dmu)], 固定资产: fixed_assets, 员工人数: employees, 能源消耗: energy, 营业收入: revenue, 净利润: profit, CO2排放: co2 }) # 人为制造脏数据负值、缺失值、异常值 df.loc[4, 固定资产] -6.0 # D05 固定资产为负 df.loc[11, 员工人数] np.nan # D12 员工人数缺失 df.loc[16, 营业收入] 350.0 # D17 营业收入异常偏高明显偏离分布关于DMU数量有两条经验法则必须记住一是DMU数量至少要为投入产出指标总数的3倍二是DMU数量不能少于投入数与产出数的乘积。本文3个投入、3个产出含非期望经验法则要求至少18家20家刚好达标。如果DMU数量不够模型很容易出现大量DMU效率为1的情况结果完全没有区分度——这个问题后面会展开讲。3. 数据清洗这一步决定了结果可信度的一半3.1 DEA对数据的硬性要求DEA模型对数据有几个硬性要求任何一个不满足都会导致求解失败或者结果不可信。第一是数据必须为正。投入产出出现0或者负值SBM模型分母就会出问题因为松弛量要用实际投入值做除法。第二是数据不能有缺失。缺失值在线性规划里会直接导致约束条件无法构造要么填充、要么删除。第三是指标之间不能存在完全线性相关。比如你把“总成本”和“固定成本可变成本”同时放进模型这俩一相加就线性相关了求解器要么报错要么给出极其不稳定的解。第四是量级不能差距太离谱。虽然SBM模型理论上是单位不变的但如果你某个指标是0.01级别、另一个是10000级别线性规划求解器在数值层面上容易遇到病态矩阵问题我实测下来还是建议做一次标准化或者缩放。3.2 缺失值和负值怎么处理来看我们这份数据里的人工“脏数据”。D05固定资产变成了-6这是在制造“录入错误”场景D12员工人数是NaN模拟“统计遗漏”D17营业收入是350明显偏离其他19家企业的分布区间属于“异常记录”。负值的处理原则很明确如果投入产出指标为负在DEA框架下没有直接的数学意义首先要判断是数据录入错误还是业务本身真的会出现负值。如果是录入错误直接用合理值替换如果是业务真实值比如某些企业净利润为负那需要重新考虑指标体系设计——比如把净利润替换成“利润总额1”或者改用其他代理变量。在我们这个案例里固定资产为负显然是录入错误我采用该指标其他19家企业的均值来修正。缺失值处理要看缺失比例。缺失比例低于5%可以用均值、中位数填充如果高于20%建议直接删除该DMU因为填充进去的值会严重扭曲效率评价结果。D12只缺失一个值我用中位数填充。这里强调一下填充值会影响最终效率得分所以在报告里一定要注明“哪些DMU的数据经过填充处理”方便评审时追溯。# 1处理负值用同指标非异常均值修正 valid_assets df.loc[(df[固定资产] 0), 固定资产] df.loc[4, 固定资产] valid_assets.mean() # 2处理缺失值用中位数填充 df.loc[11, 员工人数] df[员工人数].median() # 3检测异常值用IQR方法 def detect_outliers(series): q1 series.quantile(0.25) q3 series.quantile(0.75) iqr q3 - q1 lower q1 - 1.5 * iqr upper q3 1.5 * iqr return (series lower) | (series upper) outlier_flags df[[固定资产, 员工人数, 能源消耗, 营业收入, 净利润, CO2排放]].apply(detect_outliers) print(outlier_flags.sum())运行IQR检测会发现D17的营业收入被标为异常值。对于异常值处理逻辑稍微复杂一些如果是单纯的录入错误修正如果是业务真实情况比如某企业当年有特殊一次性收入则保留但要在后续分析中重点关注这个DMU因为它可能是效率前沿的异常点。我这里选择保留D17的数据理由是企业之间本来就存在规模差异营业收入高并不等于说不合理只是需要单独观察它的效率表现是否异常。3.3 指标去冗余相关性检验该看该信谁指标之间的高度相关性是DEA建模中非常隐蔽的坑。很多人把所有能拿到的指标一股脑塞进模型结果求解倒是没问题但效率值变得毫无规律原因就是指标信息重复。判断标准我一般用Spearman相关系数对分布不做正态假设更适合业务数据。如果两个投入指标相关系数超过0.85我只保留其中一个——因为你放两个高度相关的指标进去相当于在无形中给这个维度的投入“加了权重”会扭曲效率结果。反过来投入和产出之间的关系则不怕相关甚至说高度相关才合理能源消耗和CO₂排放相关性强恰恰说明数据质量好。对这份仿真数据我还额外产生了一个“重复指标”把“能源消耗”的1.2倍加上一个小扰动当作“碳排放强度”验证相关性检验能不能发现这种明显的线性关系。实际分析中做完IQR清洗后建议执行这一步# 相关性矩阵Spearman corr df[[固定资产, 员工人数, 能源消耗, 营业收入, 净利润, CO2排放]].corr(methodspearman) print(corr.round(3)) # 找到高相关对 high_corr_pairs [] for i in range(len(corr.columns)): for j in range(i1, len(corr.columns)): if abs(corr.iloc[i, j]) 0.85: high_corr_pairs.append((corr.columns[i], corr.columns[j], round(corr.iloc[i, j], 3))) print(high_corr_pairs)如果发现高度相关的一对指标我的建议是优先删除业务解释力较弱的那一个。比如“员工人数”和“人工成本”高度相关通常保留员工人数因为更直观且口径统一“能源消耗”和“CO₂排放”虽然高度相关但一个是投入一个是非期望产出在SBM模型里分属不同类别即使相关性强也可以同时保留因为它们的经济含义不同。3.4 量纲归一化到底要不要做关于归一化有不同流派。有的文献要求把数据全部换算到[0,1]区间因为在传统DEA模型里量纲会影响径向比例的计算。但SBM模型由于目标函数里的松弛量都除以了对应的投入/产出值理论上具备单位不变性——意思是你用“亿元”还是“万元”表示固定资产算出来的ρ一模一样。不过我在实操中还是建议做一次“缩放到相近量级”的处理不是为了理论正确性而是为了数值稳定性。把固定资产5~50、员工人数100~1000、能源消耗10~100这些指标直接丢进scipy的线性规划矩阵的条件数可能会很大实际运算时容易出现轻微数值误差。稳妥做法是每个指标除以该指标的均值把数据缩放到1附近这样求解器表现更稳定而且结果完全不受影响。from sklearn.preprocessing import StandardScaler # 注意这里用“除以均值”而不是标准化保持DEA的语义 X_raw df[[固定资产, 员工人数, 能源消耗]].values Y_good_raw df[[营业收入, 净利润]].values Y_bad_raw df[[CO2排放]].values x_mean X_raw.mean(axis0) y_good_mean Y_good_raw.mean(axis0) y_bad_mean Y_bad_raw.mean(axis0) X X_raw / x_mean Y_good Y_good_raw / y_good_mean Y_bad Y_bad_raw / y_bad_mean这份“清洗缩放”之后的数据才是真正可以喂给SBM求解器的输入。很多新手直接拿原始数据跑结果效率值对单位特别敏感其实就是数值稳定性问题在作祟。数据清洗做到这一步后面的事情就顺了。4. SBM模型核心代码实现4.1 用scipy.optimize.linprog搭线性规划骨架SBM投入导向模型可以转化为标准线性规划问题来求解。目标函数ρ 1 - (1/m)Σ(sᵢ⁻/xᵢ₀)里xᵢ₀是已知数所以最大化ρ等价于最小化 -Σ(sᵢ⁻/xᵢ₀)。约束条件有投入等式约束和产出不等式约束再根据是否考虑规模报酬决定加不加Σλ1。用scipy.optimize.linprog的求解套路很固定核心是把决策变量组织成λ和松弛变量的大向量然后把所有约束写成矩阵形式。我直接给出完整代码from scipy.optimize import linprog def sbm_io_crs(X, Y, n_dmu, x_index, y_index): 投入导向SBM模型CRS假设 X: 投入矩阵 shape(m, n_dmu) 已归一化 Y: 期望产出矩阵 shape(s, n_dmu) 已归一化 m X.shape[0] s Y.shape[0] n n_dmu # 决策变量顺序: lambda_1..lambda_n, s_minus_1..s_minus_m n_vars n m # 当前被评价DMU的索引 # 这里逐DMU求解返回效率值列表 efficiencies [] for k in range(n): x0 X[:, k] y0 Y[:, k] # 目标系数: 最小化 -sum(s_i^- / x_i0) c np.zeros(n_vars) for i in range(m): c[n i] -1.0 / x0[i] # 等式约束: X lambda s^- x0 A_eq np.hstack([X.T, np.eye(m)]) b_eq x0 # 不等式约束: -Y lambda -y0 A_ub np.hstack([-Y.T, np.zeros((s, m))]) b_ub -y0 # 变量边界: lambda0, s^-0 bounds [(0, None)] * n_vars res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq, boundsbounds, methodhighs) if res.success: rho 1 - np.sum((res.x[n:]) / x0) / m # 重新计算rho efficiencies.append(rho) else: efficiencies.append(np.nan) print(fDMU {k1} 求解失败: {res.message}) return np.array(efficiencies)几个容易踩的细节我先提醒一下。第一个是linprog默认的最小化目标和SBM最大化ρ之间的关系目标函数里只放s⁻的系数λ的系数全设为0这个别搞混。第二个是A_eq矩阵的拼法X.T形状是(n, m)和单位矩阵I_m横向拼接后形状是(n, nm)对应每个DMU的m个约束等式。第三个是产出约束写成不等式-Yλ ≤ -y₀而不是Yλ ≥ y₀因为linprog只接受“小于等于”的约束形式所以要乘负号。这些细节任何一个写错求解结果就会变得莫名其妙。4.2 考虑规模报酬VRS版本只需要改一行上面的实现是CRS规模报酬不变假设。CRS意味着小企业和大企业的效率是按同一标准衡量的——规模扩大一倍产出也应该扩大一倍。但现实中很多行业存在规模效应中小企业可能因为规模不足而效率偏低但这不一定是“经营不善”只是还没到最优规模。这时候就要用VRS规模报酬可变假设。VRS在约束上只多了一个等式Σλ 1。这个约束的直观含义是构造参照生产前沿时只允许用“相似规模”的企业做线性组合不能拿一个大型企业的投入产出比例去评价小企业。加入方式很简单def sbm_io_vrs(X, Y, n_dmu, x_index, y_index): m X.shape[0] s Y.shape[0] n n_dmu n_vars n m efficiencies [] for k in range(n): x0 X[:, k] y0 Y[:, k] c np.zeros(n_vars) for i in range(m): c[n i] -1.0 / x0[i] A_eq np.hstack([X.T, np.eye(m)]) b_eq x0 # 追加 Σλ 1 约束 A_eq np.vstack([A_eq, np.hstack([np.ones(n), np.zeros(m)])]) b_eq np.append(b_eq, 1.0) A_ub np.hstack([-Y.T, np.zeros((s, m))]) b_ub -y0 bounds [(0, None)] * n_vars res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq, boundsbounds, methodhighs) if res.success: rho 1 - np.sum((res.x[n:]) / x0) / m efficiencies.append(rho) else: efficiencies.append(np.nan) print(fDMU {k1} 求解失败: {res.message}) return np.array(efficiencies)CRS和VRS的结果差异往往很大。具体选哪个取决于你的业务判断如果做区域经济评价可以考虑CRS因为区域可以自由竞争和调整规模如果做企业内部车间/部门效率评价VRS更合适因为部门间规模差异大且不能随意改变体量。如果拿不准两个都跑一遍把CRS算出的效率和VRS算出的效率之差解读为“规模效率”也是一个很常用的分析视角——这个差值反映的就是“规模不当带来的效率损失”。4.3 非期望产出的SBM扩展实现现在处理非期望产出。前面提到Tone的SBM-Undesirable模型理论严谨但要处理分式规划实现复杂度较高。我实际项目中更常用一种折中方案把非期望产出变换后纳入模型。具体做法是把CO₂这类非期望产出除以DMU的期望营收得到“单位产出的排放强度”再把这个强度当作投入指标处理。这样做的好处是目前代码框架不用大幅改动直接在投入矩阵里多塞一行即可缺点是理论上不够严谨如果数据存在某企业产出为0强度值会爆炸。如果你需要严格的SBM-Undesirable模型可以考虑用Charnes-Cooper变换做线性化引入标量t令 tλ、ts⁻、ts⁺ 为新的决策变量目标函数变成一个线性表达式约束条件同样变成线性。变换后的代码如下def sbm_undesirable_vrs(X, Yg, Yb, n_dmu): 非期望产出SBM模型VRS假设 X: 投入矩阵 (m, n) Yg: 期望产出矩阵 (s1, n) Yb: 非期望产出矩阵 (s2, n) m X.shape[0] s1 Yg.shape[0] s2 Yb.shape[0] n n_dmu n_vars 1 n m s1 s2 # t, lambda, s-, s, sb efficiencies [] for k in range(n): x0 X[:, k] yg0 Yg[:, k] yb0 Yb[:, k] # 目标: min t - (1/m)*sum(si^-/x_i0) # 这里使用变换后的目标 c np.zeros(n_vars) c[0] 1.0 for i in range(m): c[1 n i] -1.0 / (m * x0[i]) # 约束1: sum(lambda) t A_eq np.hstack([-1.0, np.ones(n), np.zeros(m s1 s2)]) b_eq 0.0 # 约束2: X lambda s^- - t*x0 0 A_eq2 np.hstack([-x0.reshape(-1,1), X.T, np.eye(m), np.zeros((m, s1s2))]) b_eq2 np.zeros(m) # 约束3: Yg lambda - s^ - t*yg0 0 A_eq3 np.hstack([-yg0.reshape(-1,1), Yg.T, np.zeros((s1, m)), -np.eye(s1), np.zeros((s1, s2))]) b_eq3 np.zeros(s1) # 约束4: Yb lambda s^b - t*yb0 0 A_eq4 np.hstack([-yb0.reshape(-1,1), Yb.T, np.zeros((s2, ms1)), np.eye(s2)]) b_eq4 np.zeros(s2) A_eq np.vstack([A_eq, A_eq2, A_eq3, A_eq4]) b_eq np.concatenate([b_eq, b_eq2, b_eq3, b_eq4]) bounds [(0, None)] * n_vars # t 的边界要设为 (0, None) res linprog(c, A_eqA_eq, b_eqb_eq, boundsbounds, methodhighs) if res.success: t_opt res.x[0] s_minus res.x[1n:1nm] efficiency 1 - np.sum(s_minus / (m * x0)) efficiencies.append(efficiency) else: efficiencies.append(np.nan) print(fDMU {k1} 求解失败: {res.message}) return np.array(efficiencies)这个版本代码量略大但每一步都有明确的数学对应关系。如果你只是做个快速分析用前面“把非期望产出强度当投入”的简化版就够用了如果是正式发论文或者给外部评审机构的报告建议用这个严格版本。4.4 效率值出来了怎么解读才对把两份数据跑完你会得到每个DMU的效率值。这里有几个解读时需要特别注意的地方。效率值为1只代表“相对有效”不代表“绝对最优”。DMU效率为1说明它在样本内位于生产前沿面上但如果整个行业的技术水平都很低这个“前沿面”本身并不先进。所以报告里最好叫“相对效率得分”不要写成“绝对效率”。另外投入导向的效率值含义是“在不减少产出的前提下各投入最多能同时缩减到现在的百分之多少”。比如某企业效率0.82意思是产出不变时理论上所有投入可以等比例压缩18%CRS情境下。如果按SBM非径向的结果还需要看具体哪个投入的松弛量最大那个就是效率改进的最优切入点。我习惯把所有DMU的效率得分整理成一张表格按照降序排列同时标注它是否处于前沿面上。这一步也是后面可视化选图表的前提——效率得分的分布形态直接决定你用柱状图还是密度图来展示。5. 结果可视化让业务方一眼看懂效率5.1 效率得分排名条形图先看整体拿到效率值后第一步先画一张排名条形图。横轴是效率得分纵轴是DMU编号按得分从高到低排列。这样一下子就能看到哪些企业领先、哪些企业垫底。用matplotlib实现import matplotlib.pyplot as plt import seaborn as sns plt.rcParams[font.sans-serif] [SimHei] # 支持中文 plt.rcParams[axes.unicode_minus] False eff_df pd.DataFrame({ DMU: df[DMU], 效率值: efficiencies_vrs }).sort_values(效率值, ascendingTrue) fig, ax plt.subplots(figsize(10, 8)) colors [#2e8b57 if v 0.9 else #3a7ca5 for v in eff_df[效率值]] ax.barh(eff_df[DMU], eff_df[效率值], colorcolors) ax.set_xlabel(SBM效率值VRS) ax.set_title(20家制造企业SBM效率得分排名) for idx, v in enumerate(eff_df[效率值]): ax.text(v 0.005, idx, f{v:.3f}, vacenter, fontsize9) plt.tight_layout() plt.show()这里有一个非常实用的小技巧先按效率值升序排序再画水平条形图这样效率最高的DMU会出现在图的最上方符合人的阅读习惯。颜色区分也很重要——效率值≥0.9的用绿色其余用蓝色业务方一眼就能锁定问题企业。5.2 效率分布直方图和密度曲线排名图能看出个体差异但看不出整体分布形态。我通常还会画一个效率得分的直方图叠加核密度估计曲线。如果分布呈现“双峰”——一堆企业挤在0.6附近另一堆企业挤在0.95以上——就说明行业内存在明显的“分化现象”如果接近正态分布说明效率差距是渐进的。这个判断对后续政策建议的措辞影响很大双峰分布适合用“行业分化严重需重点关注低效群体”这类表述正态分布则适合用“整体水平差距有限可在普遍提升方向上发力”。fig, ax plt.subplots(figsize(8, 5)) sns.histplot(efficiencies_vrs, bins15, kdeTrue, color#3a7ca5) ax.axvline(efficiencies_vrs.mean(), colorred, linestyle--, labelf均值 {efficiencies_vrs.mean():.3f}) ax.set_xlabel(SBM效率值) ax.set_title(效率值分布直方图) ax.legend() plt.show()5.3 松弛量热力图找到每个企业的改进方向效率排名只能告诉你“谁差”松弛量热力图才能告诉你“差在哪”。SBM模型的优势之一就是能给出每个投入的具体松弛量我把这些松弛量整理到一张表里用seaborn画热力图。横轴是各投入和产出指标纵轴是DMU颜色深浅表示松弛量标准化后的大小。颜色越深代表这个指标的改进空间越大。# 假设在求解过程中保存了每个DMU的松弛量矩阵 slack_matrix # shape (n_dmu, ms) slack_df pd.DataFrame(slack_matrix, columns[固定资产冗余, 员工冗余, 能耗冗余, 营收不足, 利润不足]) slack_df[DMU] df[DMU] slack_df.set_index(DMU, inplaceTrue) fig, ax plt.subplots(figsize(10, 8)) sns.heatmap(slack_df, cmapOrRd, annotTrue, fmt.3f, linewidths0.5, cbar_kws{label: 松弛量归一化}) plt.title(各DMU投入产出松弛量热力图) plt.tight_layout() plt.show()画完这张图业务方的反馈通常很积极因为从图上可以直接读出“D07员工冗余最严重”“D12能源消耗冗余高”而不是听你空泛地说“效率偏低”。5.4 关联分析图投入松弛与效率的对角线关系一个高价值的辅助分析是把“能源消耗松弛率”和“效率值”画成散点图看是否存在某种模式。举例来说如果能源松弛率高的企业效率值普遍低说明能源管理是当前效率改善的关键抓手如果有些效率值不低的企业也有一定的能耗松弛说明这些企业还“有余力”进一步节能减排。fig, ax plt.subplots(figsize(8, 6)) sc ax.scatter(slack_df[能耗冗余], efficiencies_vrs, cdf[CO2排放], cmapviridis, s80) ax.set_xlabel(能源消耗松弛量归一化) ax.set_ylabel(SBM效率值) ax.set_title(效率值与能耗松弛的关系颜色CO2排放量) plt.colorbar(sc, labelCO2排放) plt.show()颜色映射CO₂排放后这张图还能直观展示“高排放企业是否高效率”——如果高排放企业颜色偏深效率值很高说明环保与效率在当前样本中并不冲突如果高排放企业的效率值反而低那么节能降耗的改进空间就更值得强调。6. 常见问题排查与避坑指南6.1 求解失败或效率全为1先检查这四件事最让人抓狂的情况是运行代码后有的DMU求解失败或者所有DMU效率都等于1。我总结下来80%是以下四个原因。第一数据里有0或负值。SBM目标函数分母出现0直接除零错误即使不报错结果也会异常。检查方法很简单对投入产出矩阵做np.min(X, axis1)任何一个指标的最小值≤0就得回头清洗。第二DMU数量太少。如果你的DMU数量≤投入产出指标数量模型几乎没有区分能力每个DMU都能找到一条“专属”的组合路径使其效率为1。这就是所谓的“维数灾难”。解决办法要么增加采样样本要么减少指标数量。第三约束方向写反了。linprog的A_ub矩阵默认是“≤”约束很多人把Yλ ≥ y₀直接写成Yλ ≤ y₀结果求解器永远返回一个λ0的“退化解”看起来成功了但效率值全为1。调试技巧是随机取一个已知低效的DMU手算验证一下约束是否真的被满足了。第四没有加Σλ1的VRS约束。如果你的模型本意是VRS但忘了加这个约束效率值会比真实值偏低而且规模差异大的样本还会出现规模越大的企业效率越高/越低这类荒谬结论。6.2 松弛量出现“负数”是怎么回事正常SBM模型里松弛变量都是非负的。如果你在结果里看到负的松弛量不用怀疑八成是代码的符号搞错了。比如在A_eq约束里Xλ s⁻ x₀如果写成了Xλ - s⁻ x₀求出来的s⁻就是负值。同样的道理产出端Yλ - s⁺ y₀写错符号s⁺也会变负。建议每次跑完后顺手print几个DMU的松弛量看一眼负值立刻暴露。另外还要注意linprog返回的解可能包含极小的浮点误差比如效率值被算成1.0000000001或者松弛量出现-1e-15这种数值噪音。判断是否有效时要用np.isclose(x, 0, atol1e-8)而不是x 0效率值是否等于1也要用np.isclose(rho, 1.0)来判断。6.3 结果不稳定换求解器也没用如果你把数据换成另一批样本效率排名波动很大这通常不是求解器的问题而是指标体系设计的问题。常见情况是把两个高度相关的指标同时放进了模型比如“营业收入”和“利润”相关性很高导致生产前沿面对指标选择极其敏感。我在3.3节讲的Spearman相关检验做一遍就能筛查掉大部分这类问题。另外一个隐蔽问题是数据质量如果某家企业的某个投入值极端小比如接近0它会在DEA里占据非常特殊的位置直接把自己变成效率前沿点甚至拉偏整个前沿面。实战中遇到这种情况我会把它单独提出来查看投入产出配比是否异常确认不是录入错误后考虑是否要从样本中剔除。6.4 求解速度太慢DMU规模上来之后怎么优化当DMU数量达到上千时逐个DMU循环调linprog会变得很慢。比如1000个DMU、10个指标等于要解1000个线性规划可能要好几分钟。这时候有几个提速手段。首选是用scipy的HiGHS求解器替代默认方法我们代码里已经指定了methodhighs这个求解器比默认的interior-point快不少。其次可以考虑用并行Python的multiprocessing.Pool按4核并行速度几乎线性提升。最后如果DMU数量上万建议换成Gurobi或CPLEX这类商业求解器它们对批量线性规划的支持更好。from multiprocessing import Pool def solve_one(k): # 对第k个DMU求解SBM return k, sbm_single_dmu(X, Y, k) if __name__ __main__: with Pool(processes4) as pool: results pool.map(solve_one, range(n_dmu)) # 汇总效率值不过说实话对绝大多数业务分析场景几十到几百个DMU单线程跑也够了。最重要的还是先把数据清洗和模型定义做对性能优化是后面的事。再分享一个关于指标设计的经验很多人上来就追求指标多把能拿到的数据全塞进模型。但DEA模型是“指标越多效率前沿越宽松效率区分度越差”。我的经验是投入3~5个、产出2~4个含非期望总共不超过8个指标已经能覆盖大多数制造业、服务业和公共部门的评价需求。指标少于这个范围业务方觉得维度不够多于这个范围结果又变成“全员优等生”两头都不讨好。最后聊两句体感。SBM模型用Python实现最大的门槛其实不在模型本身而在于“数据和业务的衔接”。我跑过很多次效率分析算法的坑大多能靠调试解决真正麻烦的是怎么让业务方理解“效率0.82不等于企业很差而是还有18%的投入压缩空间”。这就要靠可视化和报告叙事来解决了——所以千万别把可视化当成最后一步的锦上添花它是整个分析交付的“最后一公里”。数据清洗、模型求解、可视化这三件事每一件都值得花同样的心思。 SBM模型Python实战：从数据清洗到效率可视化完整指南
做效率评价这些年我最大的感受是模型本身并不难真正让人头秃的是“数据从哪来、怎么洗干净、结果怎么让业务方看懂”。SBMSlacks-Based Measure松弛变量测度模型就是这样——原理不复杂但如果你按教材里的公式硬撸代码十有八九会卡在数据预处理和求解器配置上。这篇文章我就把从数据清洗到结果可视化的完整流程拆开讲附带能直接跑的Python代码适合正在做企业效率评价、区域经济绩效分析、供应链绩效评估的读者参考。1. 先搞清楚SBM模型到底在算什么1.1 从DEA家族说起为什么选SBM而非传统CCR/BCC数据包络分析DEA是一大类基于线性规划的效率评价方法最早由Charnes、Cooper和Rhodes在1978年提出CCR模型后来Banker等人扩展出BCC模型。这类方法的核心理念是不需要预设生产函数直接用“投入产出数据”构造生产前沿面——所有决策单元DMU比如企业、医院、城市中在同样投入下产出最高、或者在同样产出下投入最少的那些DMU就构成效率前沿其他DMU与前沿面的距离就是效率损失。传统CCR和BCC模型有一个共同特点它们属于“径向”模型即假设所有投入或产出按同一比例压缩或扩张。但实际生产过程中不可能所有投入都同比例减少——某台设备定了就是定了人员调整有刚性能耗削减空间和原材料削减空间也完全不同。这时候径向模型算出来的效率值就会偏大而且没办法告诉你“具体哪个投入冗余了多少”。SBM模型正是针对这个问题提出来的它把松弛变量直接放进目标函数允许各投入不均等地缩减同时给出每个指标的松弛量信息量比传统DEA模型高出一个量级。1.2 松弛变量的本质SBM到底在优化什么松弛变量这个名字听起来抽象其实理解起来非常简单。假设你有20家制造企业每家投入固定资产、员工人数、能源消耗产出是营业收入和利润。现在把20家企业的数据放到一起找一个“标杆”在同样投入水平下产出最高的那家就是效率前沿。对某个效率不足的企业来说它的投入松弛量就是“相对于标杆多用了多少固定资产、多少人、多少能源”产出松弛量则是“相对于标杆少产出了多少利润”。SBM模型把这些松弛量放进目标函数里求最小化。以投入导向为例目标函数是ρ 1 - (1/m) * Σ(sᵢ⁻ / xᵢ₀)其中m是投入指标个数sᵢ⁻是第i个投入的松弛量xᵢ₀是当前评价DMU的第i个投入值。这个公式的意思是如果这家企业各投入相对标杆都有冗余ρ就会小于1如果没有任何松弛ρ等于1说明它已经在生产前沿面上。由于分母是各投入的实际值松弛量被“相对化”了所以SBM模型是单位不变的——你用“万元”还是“亿元”表示固定资产计算结果完全一样。1.3 非期望产出怎么进模型做效率评价很多场景绕不开一个现实问题非期望产出。比如制造业的CO₂排放、废水排放银行的坏账率医院的再入院率。这些指标“越少越好”但又不能简单地放进“投入”里处理从生产逻辑上讲排放是产出端的伴随物不是先期投入的资源。Tone在2003年提出了处理非期望产出的SBM扩展模型核心思路是把产出再拆成期望产出和非期望产出两类分别计算松弛量分母部分同时考虑“期望产出不足”和“非期望产出过多”两项。不过我得说句实话非期望产出的SBM虽然理论严谨但求解比基础版复杂需要做Charnes-Cooper变换把分式规划转成线性规划。实际项目中如果只是做初步筛选也可以用“把非期望产出当成投入”的近似做法——虽然这个方法在学术上有争议但在业务场景里速度快、结果也基本符合直觉。两种方式我在后文都会给出代码实现。2. 开工前准备环境、依赖与实验数据2.1 Python环境与依赖库清单我做这个流程用的Python 3.10只依赖四个库全是免费开源的库名用途说明numpy数值计算、矩阵操作数据处理的基础pandas数据清洗、表格操作读Excel、去重、缺失值处理都靠它scipy线性规划求解核心的linprog函数matplotlib seaborn可视化画排名图、热力图、散点图安装命令就一条pip install pandas numpy scipy matplotlib seaborn如果你在公司内网可能需要配置国内镜像源比如清华源、阿里源速度会快很多。IDE我建议用VS CodePython插件装好之后配一下解释器路径就能直接用轻量、启动快配置方法网上教程很多这里不展开。2.2 实验数据设计让案例能复现为了把整个流程讲清楚我构造了一份仿真数据。假设要评估20家制造企业的运营效率决策单元DMU编号从D01到D20指标体系如下投入3个固定资产亿元、员工人数百人、能源消耗万吨标准煤期望产出2个营业收入亿元、净利润亿元非期望产出1个CO₂排放量万吨数据用numpy生成设置随机种子保证可复现。我这里有意给部分企业加入了一些“脏数据”——比如D05的固定资产出现负值D12的员工人数是缺失值D17的营业收入异常偏高——这样后面数据清洗环节才有东西可讲。import numpy as np import pandas as pd np.random.seed(42) n_dmu 20 # 生成基础数据 fixed_assets np.random.uniform(5, 50, n_dmu) # 固定资产亿元 employees np.random.uniform(100, 1000, n_dmu) # 员工人数百人 energy np.random.uniform(10, 100, n_dmu) # 能源消耗万吨标准煤 revenue np.random.uniform(20, 200, n_dmu) # 营业收入亿元 profit np.random.uniform(1, 30, n_dmu) # 净利润亿元 co2 np.random.uniform(5, 60, n_dmu) # CO2排放量万吨 df pd.DataFrame({ DMU: [fD{i1:02d} for i in range(n_dmu)], 固定资产: fixed_assets, 员工人数: employees, 能源消耗: energy, 营业收入: revenue, 净利润: profit, CO2排放: co2 }) # 人为制造脏数据负值、缺失值、异常值 df.loc[4, 固定资产] -6.0 # D05 固定资产为负 df.loc[11, 员工人数] np.nan # D12 员工人数缺失 df.loc[16, 营业收入] 350.0 # D17 营业收入异常偏高明显偏离分布关于DMU数量有两条经验法则必须记住一是DMU数量至少要为投入产出指标总数的3倍二是DMU数量不能少于投入数与产出数的乘积。本文3个投入、3个产出含非期望经验法则要求至少18家20家刚好达标。如果DMU数量不够模型很容易出现大量DMU效率为1的情况结果完全没有区分度——这个问题后面会展开讲。3. 数据清洗这一步决定了结果可信度的一半3.1 DEA对数据的硬性要求DEA模型对数据有几个硬性要求任何一个不满足都会导致求解失败或者结果不可信。第一是数据必须为正。投入产出出现0或者负值SBM模型分母就会出问题因为松弛量要用实际投入值做除法。第二是数据不能有缺失。缺失值在线性规划里会直接导致约束条件无法构造要么填充、要么删除。第三是指标之间不能存在完全线性相关。比如你把“总成本”和“固定成本可变成本”同时放进模型这俩一相加就线性相关了求解器要么报错要么给出极其不稳定的解。第四是量级不能差距太离谱。虽然SBM模型理论上是单位不变的但如果你某个指标是0.01级别、另一个是10000级别线性规划求解器在数值层面上容易遇到病态矩阵问题我实测下来还是建议做一次标准化或者缩放。3.2 缺失值和负值怎么处理来看我们这份数据里的人工“脏数据”。D05固定资产变成了-6这是在制造“录入错误”场景D12员工人数是NaN模拟“统计遗漏”D17营业收入是350明显偏离其他19家企业的分布区间属于“异常记录”。负值的处理原则很明确如果投入产出指标为负在DEA框架下没有直接的数学意义首先要判断是数据录入错误还是业务本身真的会出现负值。如果是录入错误直接用合理值替换如果是业务真实值比如某些企业净利润为负那需要重新考虑指标体系设计——比如把净利润替换成“利润总额1”或者改用其他代理变量。在我们这个案例里固定资产为负显然是录入错误我采用该指标其他19家企业的均值来修正。缺失值处理要看缺失比例。缺失比例低于5%可以用均值、中位数填充如果高于20%建议直接删除该DMU因为填充进去的值会严重扭曲效率评价结果。D12只缺失一个值我用中位数填充。这里强调一下填充值会影响最终效率得分所以在报告里一定要注明“哪些DMU的数据经过填充处理”方便评审时追溯。# 1处理负值用同指标非异常均值修正 valid_assets df.loc[(df[固定资产] 0), 固定资产] df.loc[4, 固定资产] valid_assets.mean() # 2处理缺失值用中位数填充 df.loc[11, 员工人数] df[员工人数].median() # 3检测异常值用IQR方法 def detect_outliers(series): q1 series.quantile(0.25) q3 series.quantile(0.75) iqr q3 - q1 lower q1 - 1.5 * iqr upper q3 1.5 * iqr return (series lower) | (series upper) outlier_flags df[[固定资产, 员工人数, 能源消耗, 营业收入, 净利润, CO2排放]].apply(detect_outliers) print(outlier_flags.sum())运行IQR检测会发现D17的营业收入被标为异常值。对于异常值处理逻辑稍微复杂一些如果是单纯的录入错误修正如果是业务真实情况比如某企业当年有特殊一次性收入则保留但要在后续分析中重点关注这个DMU因为它可能是效率前沿的异常点。我这里选择保留D17的数据理由是企业之间本来就存在规模差异营业收入高并不等于说不合理只是需要单独观察它的效率表现是否异常。3.3 指标去冗余相关性检验该看该信谁指标之间的高度相关性是DEA建模中非常隐蔽的坑。很多人把所有能拿到的指标一股脑塞进模型结果求解倒是没问题但效率值变得毫无规律原因就是指标信息重复。判断标准我一般用Spearman相关系数对分布不做正态假设更适合业务数据。如果两个投入指标相关系数超过0.85我只保留其中一个——因为你放两个高度相关的指标进去相当于在无形中给这个维度的投入“加了权重”会扭曲效率结果。反过来投入和产出之间的关系则不怕相关甚至说高度相关才合理能源消耗和CO₂排放相关性强恰恰说明数据质量好。对这份仿真数据我还额外产生了一个“重复指标”把“能源消耗”的1.2倍加上一个小扰动当作“碳排放强度”验证相关性检验能不能发现这种明显的线性关系。实际分析中做完IQR清洗后建议执行这一步# 相关性矩阵Spearman corr df[[固定资产, 员工人数, 能源消耗, 营业收入, 净利润, CO2排放]].corr(methodspearman) print(corr.round(3)) # 找到高相关对 high_corr_pairs [] for i in range(len(corr.columns)): for j in range(i1, len(corr.columns)): if abs(corr.iloc[i, j]) 0.85: high_corr_pairs.append((corr.columns[i], corr.columns[j], round(corr.iloc[i, j], 3))) print(high_corr_pairs)如果发现高度相关的一对指标我的建议是优先删除业务解释力较弱的那一个。比如“员工人数”和“人工成本”高度相关通常保留员工人数因为更直观且口径统一“能源消耗”和“CO₂排放”虽然高度相关但一个是投入一个是非期望产出在SBM模型里分属不同类别即使相关性强也可以同时保留因为它们的经济含义不同。3.4 量纲归一化到底要不要做关于归一化有不同流派。有的文献要求把数据全部换算到[0,1]区间因为在传统DEA模型里量纲会影响径向比例的计算。但SBM模型由于目标函数里的松弛量都除以了对应的投入/产出值理论上具备单位不变性——意思是你用“亿元”还是“万元”表示固定资产算出来的ρ一模一样。不过我在实操中还是建议做一次“缩放到相近量级”的处理不是为了理论正确性而是为了数值稳定性。把固定资产5~50、员工人数100~1000、能源消耗10~100这些指标直接丢进scipy的线性规划矩阵的条件数可能会很大实际运算时容易出现轻微数值误差。稳妥做法是每个指标除以该指标的均值把数据缩放到1附近这样求解器表现更稳定而且结果完全不受影响。from sklearn.preprocessing import StandardScaler # 注意这里用“除以均值”而不是标准化保持DEA的语义 X_raw df[[固定资产, 员工人数, 能源消耗]].values Y_good_raw df[[营业收入, 净利润]].values Y_bad_raw df[[CO2排放]].values x_mean X_raw.mean(axis0) y_good_mean Y_good_raw.mean(axis0) y_bad_mean Y_bad_raw.mean(axis0) X X_raw / x_mean Y_good Y_good_raw / y_good_mean Y_bad Y_bad_raw / y_bad_mean这份“清洗缩放”之后的数据才是真正可以喂给SBM求解器的输入。很多新手直接拿原始数据跑结果效率值对单位特别敏感其实就是数值稳定性问题在作祟。数据清洗做到这一步后面的事情就顺了。4. SBM模型核心代码实现4.1 用scipy.optimize.linprog搭线性规划骨架SBM投入导向模型可以转化为标准线性规划问题来求解。目标函数ρ 1 - (1/m)Σ(sᵢ⁻/xᵢ₀)里xᵢ₀是已知数所以最大化ρ等价于最小化 -Σ(sᵢ⁻/xᵢ₀)。约束条件有投入等式约束和产出不等式约束再根据是否考虑规模报酬决定加不加Σλ1。用scipy.optimize.linprog的求解套路很固定核心是把决策变量组织成λ和松弛变量的大向量然后把所有约束写成矩阵形式。我直接给出完整代码from scipy.optimize import linprog def sbm_io_crs(X, Y, n_dmu, x_index, y_index): 投入导向SBM模型CRS假设 X: 投入矩阵 shape(m, n_dmu) 已归一化 Y: 期望产出矩阵 shape(s, n_dmu) 已归一化 m X.shape[0] s Y.shape[0] n n_dmu # 决策变量顺序: lambda_1..lambda_n, s_minus_1..s_minus_m n_vars n m # 当前被评价DMU的索引 # 这里逐DMU求解返回效率值列表 efficiencies [] for k in range(n): x0 X[:, k] y0 Y[:, k] # 目标系数: 最小化 -sum(s_i^- / x_i0) c np.zeros(n_vars) for i in range(m): c[n i] -1.0 / x0[i] # 等式约束: X lambda s^- x0 A_eq np.hstack([X.T, np.eye(m)]) b_eq x0 # 不等式约束: -Y lambda -y0 A_ub np.hstack([-Y.T, np.zeros((s, m))]) b_ub -y0 # 变量边界: lambda0, s^-0 bounds [(0, None)] * n_vars res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq, boundsbounds, methodhighs) if res.success: rho 1 - np.sum((res.x[n:]) / x0) / m # 重新计算rho efficiencies.append(rho) else: efficiencies.append(np.nan) print(fDMU {k1} 求解失败: {res.message}) return np.array(efficiencies)几个容易踩的细节我先提醒一下。第一个是linprog默认的最小化目标和SBM最大化ρ之间的关系目标函数里只放s⁻的系数λ的系数全设为0这个别搞混。第二个是A_eq矩阵的拼法X.T形状是(n, m)和单位矩阵I_m横向拼接后形状是(n, nm)对应每个DMU的m个约束等式。第三个是产出约束写成不等式-Yλ ≤ -y₀而不是Yλ ≥ y₀因为linprog只接受“小于等于”的约束形式所以要乘负号。这些细节任何一个写错求解结果就会变得莫名其妙。4.2 考虑规模报酬VRS版本只需要改一行上面的实现是CRS规模报酬不变假设。CRS意味着小企业和大企业的效率是按同一标准衡量的——规模扩大一倍产出也应该扩大一倍。但现实中很多行业存在规模效应中小企业可能因为规模不足而效率偏低但这不一定是“经营不善”只是还没到最优规模。这时候就要用VRS规模报酬可变假设。VRS在约束上只多了一个等式Σλ 1。这个约束的直观含义是构造参照生产前沿时只允许用“相似规模”的企业做线性组合不能拿一个大型企业的投入产出比例去评价小企业。加入方式很简单def sbm_io_vrs(X, Y, n_dmu, x_index, y_index): m X.shape[0] s Y.shape[0] n n_dmu n_vars n m efficiencies [] for k in range(n): x0 X[:, k] y0 Y[:, k] c np.zeros(n_vars) for i in range(m): c[n i] -1.0 / x0[i] A_eq np.hstack([X.T, np.eye(m)]) b_eq x0 # 追加 Σλ 1 约束 A_eq np.vstack([A_eq, np.hstack([np.ones(n), np.zeros(m)])]) b_eq np.append(b_eq, 1.0) A_ub np.hstack([-Y.T, np.zeros((s, m))]) b_ub -y0 bounds [(0, None)] * n_vars res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq, boundsbounds, methodhighs) if res.success: rho 1 - np.sum((res.x[n:]) / x0) / m efficiencies.append(rho) else: efficiencies.append(np.nan) print(fDMU {k1} 求解失败: {res.message}) return np.array(efficiencies)CRS和VRS的结果差异往往很大。具体选哪个取决于你的业务判断如果做区域经济评价可以考虑CRS因为区域可以自由竞争和调整规模如果做企业内部车间/部门效率评价VRS更合适因为部门间规模差异大且不能随意改变体量。如果拿不准两个都跑一遍把CRS算出的效率和VRS算出的效率之差解读为“规模效率”也是一个很常用的分析视角——这个差值反映的就是“规模不当带来的效率损失”。4.3 非期望产出的SBM扩展实现现在处理非期望产出。前面提到Tone的SBM-Undesirable模型理论严谨但要处理分式规划实现复杂度较高。我实际项目中更常用一种折中方案把非期望产出变换后纳入模型。具体做法是把CO₂这类非期望产出除以DMU的期望营收得到“单位产出的排放强度”再把这个强度当作投入指标处理。这样做的好处是目前代码框架不用大幅改动直接在投入矩阵里多塞一行即可缺点是理论上不够严谨如果数据存在某企业产出为0强度值会爆炸。如果你需要严格的SBM-Undesirable模型可以考虑用Charnes-Cooper变换做线性化引入标量t令 tλ、ts⁻、ts⁺ 为新的决策变量目标函数变成一个线性表达式约束条件同样变成线性。变换后的代码如下def sbm_undesirable_vrs(X, Yg, Yb, n_dmu): 非期望产出SBM模型VRS假设 X: 投入矩阵 (m, n) Yg: 期望产出矩阵 (s1, n) Yb: 非期望产出矩阵 (s2, n) m X.shape[0] s1 Yg.shape[0] s2 Yb.shape[0] n n_dmu n_vars 1 n m s1 s2 # t, lambda, s-, s, sb efficiencies [] for k in range(n): x0 X[:, k] yg0 Yg[:, k] yb0 Yb[:, k] # 目标: min t - (1/m)*sum(si^-/x_i0) # 这里使用变换后的目标 c np.zeros(n_vars) c[0] 1.0 for i in range(m): c[1 n i] -1.0 / (m * x0[i]) # 约束1: sum(lambda) t A_eq np.hstack([-1.0, np.ones(n), np.zeros(m s1 s2)]) b_eq 0.0 # 约束2: X lambda s^- - t*x0 0 A_eq2 np.hstack([-x0.reshape(-1,1), X.T, np.eye(m), np.zeros((m, s1s2))]) b_eq2 np.zeros(m) # 约束3: Yg lambda - s^ - t*yg0 0 A_eq3 np.hstack([-yg0.reshape(-1,1), Yg.T, np.zeros((s1, m)), -np.eye(s1), np.zeros((s1, s2))]) b_eq3 np.zeros(s1) # 约束4: Yb lambda s^b - t*yb0 0 A_eq4 np.hstack([-yb0.reshape(-1,1), Yb.T, np.zeros((s2, ms1)), np.eye(s2)]) b_eq4 np.zeros(s2) A_eq np.vstack([A_eq, A_eq2, A_eq3, A_eq4]) b_eq np.concatenate([b_eq, b_eq2, b_eq3, b_eq4]) bounds [(0, None)] * n_vars # t 的边界要设为 (0, None) res linprog(c, A_eqA_eq, b_eqb_eq, boundsbounds, methodhighs) if res.success: t_opt res.x[0] s_minus res.x[1n:1nm] efficiency 1 - np.sum(s_minus / (m * x0)) efficiencies.append(efficiency) else: efficiencies.append(np.nan) print(fDMU {k1} 求解失败: {res.message}) return np.array(efficiencies)这个版本代码量略大但每一步都有明确的数学对应关系。如果你只是做个快速分析用前面“把非期望产出强度当投入”的简化版就够用了如果是正式发论文或者给外部评审机构的报告建议用这个严格版本。4.4 效率值出来了怎么解读才对把两份数据跑完你会得到每个DMU的效率值。这里有几个解读时需要特别注意的地方。效率值为1只代表“相对有效”不代表“绝对最优”。DMU效率为1说明它在样本内位于生产前沿面上但如果整个行业的技术水平都很低这个“前沿面”本身并不先进。所以报告里最好叫“相对效率得分”不要写成“绝对效率”。另外投入导向的效率值含义是“在不减少产出的前提下各投入最多能同时缩减到现在的百分之多少”。比如某企业效率0.82意思是产出不变时理论上所有投入可以等比例压缩18%CRS情境下。如果按SBM非径向的结果还需要看具体哪个投入的松弛量最大那个就是效率改进的最优切入点。我习惯把所有DMU的效率得分整理成一张表格按照降序排列同时标注它是否处于前沿面上。这一步也是后面可视化选图表的前提——效率得分的分布形态直接决定你用柱状图还是密度图来展示。5. 结果可视化让业务方一眼看懂效率5.1 效率得分排名条形图先看整体拿到效率值后第一步先画一张排名条形图。横轴是效率得分纵轴是DMU编号按得分从高到低排列。这样一下子就能看到哪些企业领先、哪些企业垫底。用matplotlib实现import matplotlib.pyplot as plt import seaborn as sns plt.rcParams[font.sans-serif] [SimHei] # 支持中文 plt.rcParams[axes.unicode_minus] False eff_df pd.DataFrame({ DMU: df[DMU], 效率值: efficiencies_vrs }).sort_values(效率值, ascendingTrue) fig, ax plt.subplots(figsize(10, 8)) colors [#2e8b57 if v 0.9 else #3a7ca5 for v in eff_df[效率值]] ax.barh(eff_df[DMU], eff_df[效率值], colorcolors) ax.set_xlabel(SBM效率值VRS) ax.set_title(20家制造企业SBM效率得分排名) for idx, v in enumerate(eff_df[效率值]): ax.text(v 0.005, idx, f{v:.3f}, vacenter, fontsize9) plt.tight_layout() plt.show()这里有一个非常实用的小技巧先按效率值升序排序再画水平条形图这样效率最高的DMU会出现在图的最上方符合人的阅读习惯。颜色区分也很重要——效率值≥0.9的用绿色其余用蓝色业务方一眼就能锁定问题企业。5.2 效率分布直方图和密度曲线排名图能看出个体差异但看不出整体分布形态。我通常还会画一个效率得分的直方图叠加核密度估计曲线。如果分布呈现“双峰”——一堆企业挤在0.6附近另一堆企业挤在0.95以上——就说明行业内存在明显的“分化现象”如果接近正态分布说明效率差距是渐进的。这个判断对后续政策建议的措辞影响很大双峰分布适合用“行业分化严重需重点关注低效群体”这类表述正态分布则适合用“整体水平差距有限可在普遍提升方向上发力”。fig, ax plt.subplots(figsize(8, 5)) sns.histplot(efficiencies_vrs, bins15, kdeTrue, color#3a7ca5) ax.axvline(efficiencies_vrs.mean(), colorred, linestyle--, labelf均值 {efficiencies_vrs.mean():.3f}) ax.set_xlabel(SBM效率值) ax.set_title(效率值分布直方图) ax.legend() plt.show()5.3 松弛量热力图找到每个企业的改进方向效率排名只能告诉你“谁差”松弛量热力图才能告诉你“差在哪”。SBM模型的优势之一就是能给出每个投入的具体松弛量我把这些松弛量整理到一张表里用seaborn画热力图。横轴是各投入和产出指标纵轴是DMU颜色深浅表示松弛量标准化后的大小。颜色越深代表这个指标的改进空间越大。# 假设在求解过程中保存了每个DMU的松弛量矩阵 slack_matrix # shape (n_dmu, ms) slack_df pd.DataFrame(slack_matrix, columns[固定资产冗余, 员工冗余, 能耗冗余, 营收不足, 利润不足]) slack_df[DMU] df[DMU] slack_df.set_index(DMU, inplaceTrue) fig, ax plt.subplots(figsize(10, 8)) sns.heatmap(slack_df, cmapOrRd, annotTrue, fmt.3f, linewidths0.5, cbar_kws{label: 松弛量归一化}) plt.title(各DMU投入产出松弛量热力图) plt.tight_layout() plt.show()画完这张图业务方的反馈通常很积极因为从图上可以直接读出“D07员工冗余最严重”“D12能源消耗冗余高”而不是听你空泛地说“效率偏低”。5.4 关联分析图投入松弛与效率的对角线关系一个高价值的辅助分析是把“能源消耗松弛率”和“效率值”画成散点图看是否存在某种模式。举例来说如果能源松弛率高的企业效率值普遍低说明能源管理是当前效率改善的关键抓手如果有些效率值不低的企业也有一定的能耗松弛说明这些企业还“有余力”进一步节能减排。fig, ax plt.subplots(figsize(8, 6)) sc ax.scatter(slack_df[能耗冗余], efficiencies_vrs, cdf[CO2排放], cmapviridis, s80) ax.set_xlabel(能源消耗松弛量归一化) ax.set_ylabel(SBM效率值) ax.set_title(效率值与能耗松弛的关系颜色CO2排放量) plt.colorbar(sc, labelCO2排放) plt.show()颜色映射CO₂排放后这张图还能直观展示“高排放企业是否高效率”——如果高排放企业颜色偏深效率值很高说明环保与效率在当前样本中并不冲突如果高排放企业的效率值反而低那么节能降耗的改进空间就更值得强调。6. 常见问题排查与避坑指南6.1 求解失败或效率全为1先检查这四件事最让人抓狂的情况是运行代码后有的DMU求解失败或者所有DMU效率都等于1。我总结下来80%是以下四个原因。第一数据里有0或负值。SBM目标函数分母出现0直接除零错误即使不报错结果也会异常。检查方法很简单对投入产出矩阵做np.min(X, axis1)任何一个指标的最小值≤0就得回头清洗。第二DMU数量太少。如果你的DMU数量≤投入产出指标数量模型几乎没有区分能力每个DMU都能找到一条“专属”的组合路径使其效率为1。这就是所谓的“维数灾难”。解决办法要么增加采样样本要么减少指标数量。第三约束方向写反了。linprog的A_ub矩阵默认是“≤”约束很多人把Yλ ≥ y₀直接写成Yλ ≤ y₀结果求解器永远返回一个λ0的“退化解”看起来成功了但效率值全为1。调试技巧是随机取一个已知低效的DMU手算验证一下约束是否真的被满足了。第四没有加Σλ1的VRS约束。如果你的模型本意是VRS但忘了加这个约束效率值会比真实值偏低而且规模差异大的样本还会出现规模越大的企业效率越高/越低这类荒谬结论。6.2 松弛量出现“负数”是怎么回事正常SBM模型里松弛变量都是非负的。如果你在结果里看到负的松弛量不用怀疑八成是代码的符号搞错了。比如在A_eq约束里Xλ s⁻ x₀如果写成了Xλ - s⁻ x₀求出来的s⁻就是负值。同样的道理产出端Yλ - s⁺ y₀写错符号s⁺也会变负。建议每次跑完后顺手print几个DMU的松弛量看一眼负值立刻暴露。另外还要注意linprog返回的解可能包含极小的浮点误差比如效率值被算成1.0000000001或者松弛量出现-1e-15这种数值噪音。判断是否有效时要用np.isclose(x, 0, atol1e-8)而不是x 0效率值是否等于1也要用np.isclose(rho, 1.0)来判断。6.3 结果不稳定换求解器也没用如果你把数据换成另一批样本效率排名波动很大这通常不是求解器的问题而是指标体系设计的问题。常见情况是把两个高度相关的指标同时放进了模型比如“营业收入”和“利润”相关性很高导致生产前沿面对指标选择极其敏感。我在3.3节讲的Spearman相关检验做一遍就能筛查掉大部分这类问题。另外一个隐蔽问题是数据质量如果某家企业的某个投入值极端小比如接近0它会在DEA里占据非常特殊的位置直接把自己变成效率前沿点甚至拉偏整个前沿面。实战中遇到这种情况我会把它单独提出来查看投入产出配比是否异常确认不是录入错误后考虑是否要从样本中剔除。6.4 求解速度太慢DMU规模上来之后怎么优化当DMU数量达到上千时逐个DMU循环调linprog会变得很慢。比如1000个DMU、10个指标等于要解1000个线性规划可能要好几分钟。这时候有几个提速手段。首选是用scipy的HiGHS求解器替代默认方法我们代码里已经指定了methodhighs这个求解器比默认的interior-point快不少。其次可以考虑用并行Python的multiprocessing.Pool按4核并行速度几乎线性提升。最后如果DMU数量上万建议换成Gurobi或CPLEX这类商业求解器它们对批量线性规划的支持更好。from multiprocessing import Pool def solve_one(k): # 对第k个DMU求解SBM return k, sbm_single_dmu(X, Y, k) if __name__ __main__: with Pool(processes4) as pool: results pool.map(solve_one, range(n_dmu)) # 汇总效率值不过说实话对绝大多数业务分析场景几十到几百个DMU单线程跑也够了。最重要的还是先把数据清洗和模型定义做对性能优化是后面的事。再分享一个关于指标设计的经验很多人上来就追求指标多把能拿到的数据全塞进模型。但DEA模型是“指标越多效率前沿越宽松效率区分度越差”。我的经验是投入3~5个、产出2~4个含非期望总共不超过8个指标已经能覆盖大多数制造业、服务业和公共部门的评价需求。指标少于这个范围业务方觉得维度不够多于这个范围结果又变成“全员优等生”两头都不讨好。最后聊两句体感。SBM模型用Python实现最大的门槛其实不在模型本身而在于“数据和业务的衔接”。我跑过很多次效率分析算法的坑大多能靠调试解决真正麻烦的是怎么让业务方理解“效率0.82不等于企业很差而是还有18%的投入压缩空间”。这就要靠可视化和报告叙事来解决了——所以千万别把可视化当成最后一步的锦上添花它是整个分析交付的“最后一公里”。数据清洗、模型求解、可视化这三件事每一件都值得花同样的心思。