MATLAB实现SBM与超效率SBM模型：DEA效率评价完整指南

跑DEA效率评价这几年我最常被问到的不是“模型怎么选”而是“跑出来的结果为什么没法用”。最典型的一种情况BCC模型跑完三十个决策单元有十八个效率值都是1论文里排名没法写政策含义也没法解释。后来我基本都换成SBM模型再配合超效率SBM做二次排序整个分析逻辑才真正顺起来。这篇文章就是把这条完整的路线拆开讲清楚从原始面板数据怎么整理到SBM模型背后的数学逻辑再到MATLAB用linprog一步步实现最后是我自己在实际跑数据时踩过的坑。适合那些不想被DEAP、MaxDEA界面限制想用MATLAB灵活处理效率评价的科研人员和行业分析师。1. 为什么我弃用传统DEA模型转向SBM非径向框架1.1 径向模型看不见的松弛变量传统DEA里的CCR和BCC模型属于径向模型核心逻辑是“所有投入同时等比例压缩多少倍”或“所有产出同时等比例扩张多少倍”。这个假设在思路上很简洁但放到实际数据里很容易出问题一家工厂的能耗、人力、原料投入怎么可能完全同比例下降现实中往往是某项投入冗余严重、另一项已经压到底线了。径向模型不认这种差异它只会告诉你一个平均的压缩比例。更麻烦的是松弛变量的问题。传统模型求解完的约束条件里会残留一部分投入冗余或产出不足这部分在模型里被称为松弛变量。因为松弛变量不进入目标函数效率值明明是1的决策单元实际可能还藏着大量可改进空间。我记得有一篇早期的数据包络分析文献里就专门讽刺过这种现象“效率值为1只是说你到了前沿面上不代表你没有可压缩的冗余。”这句话我后来写进过好几篇报告里。SBM模型Slack-Based Measure就是冲着这个缺陷去的。它把投入松弛、期望产出松弛、非期望产出松弛全部写进目标函数效率值计算公式变成了“投入平均可缩减比例”除以“产出平均可扩张比例”。这样算出来的效率不仅告诉你谁好谁差还能直接告诉你每个决策单元应该在哪个指标上改进。这恰恰是实证分析最需要的东西。1.2 SBM和超效率SBM的分工先判有效再排座次SBM模型的效率值最大是1达到1说明决策单元已经处于效率前沿上。问题是面板数据经常出现几十个决策单元同时处于前沿的情况这时候就遇到跟传统DEA一样的难题都是1怎么排序超效率SBM模型就是为这个场景设计的。它的核心操作很直白把被评价的决策单元从参考集里剔除让“其他所有单元”构成一个新的生产前沿再看这个被剔除的单元离新前沿有多远。因为无法拿自己当参照物效率值就可以超过1。超效率值越高说明这个单元在效率前沿上的位置越领先。所以在我的分析流程里这两者分工很明确第一轮用常规SBM判断哪些单元有效第二轮用超效率SBM对有效单元做排名。一个解决“是不是”一个解决“排第几”。2. 面板数据的“脏活”从Excel到MATLAB三维数组的规范整理2.1 数据结构怎么摆代码才不拧巴接触过面板数据的都知道Excel里的表格摆放方式五花八门。有人按年拆成多个sheet有人用一列年份字段做长表还有人把年份放在列名里。我踩过几次坑之后给自己定了个规矩优先使用按年拆分的sheet每个sheet内部行是决策单元、列是指标且指标列名在所有sheet里完全一致。这不是洁癖纯粹是为了代码好写。按年拆开后在MATLAB里直接循环sheet名每次取出的都是形状完全相同的矩阵组合成三维数组就是水到渠成的事file panel_data.xlsx; sheets sheetnames(file); for t 1:length(sheets) tbl readtable(file, Sheet, sheets{t}); X(:, :, t) table2array(tbl(:, {labor, capital})); YG(:, :, t) table2array(tbl(:, {gdp})); YB(:, :, t) table2array(tbl(:, {co2})); end这里X、YG、YB分别是投入矩阵、期望产出矩阵和非期望产出矩阵第三维就是年份。这样后面所有模型函数都可以按“某一年份的截面数据”直接传参逻辑非常清晰。如果你的数据是长表格式也不要紧只是多一步拆分循环按年份和DMU编号筛子集后再用reshape或者逐行填充到三维数组里。我个人的经验是长表适合做数据清洗但一旦清洗完尽量转成三维数组再进入模型计算否则每个决策单元都要做逻辑索引效率低不说代码还特别容易出边界错误。2.2 三类指标的边界和方向SBM建模前最重要的一件事是把所有指标归类为投入、期望产出、非期望产出三类而且要反复确认方向。投入指标的逻辑是“越少越好”比如劳动力人数、固定资产投资、能源消耗。期望产出的逻辑是“越多越好”比如GDP、产值、利润。非期望产出则是“越少越好”但同时又是生产过程中必然产生的比如二氧化硫排放、二氧化碳排放、废水排放。这里有一个非常常见的错误有人图省事把非期望产出直接当成投入丢进模型。理论上这确实可行因为两者在SBM里都是“越少越好”但这样做的代价是丧失了经济含义的解释力。投入是可以主动调整的非期望产出往往是伴随期望产出而来的副产物把它跟资本、劳动混在一起效率评价的改进建议就不伦不类。我在这里用的是Tone2003给出的带非期望产出的SBM形式非期望产出单独一组约束松弛独立计算这样结果更贴合碳排放、污染物这类场景。2.3 零值、负值和缺失值处理不好模型直接翻车SBM模型的求解前提是数据严格为正。为什么因为目标函数里要用投入值做分母分母为0数学上没有意义。非期望产出的松弛也要除以对应的产出值同理。所以拿到数据的第一个动作不是建模而是做数据体检assert(all(X(:) 0) all(YG(:) 0) all(YB(:) 0), ... 数据必须全部为正数请先处理零值和负值);我处理这类问题的方式是分情况讨论。零值一般出现在非期望产出里比如某个单元确实没有排放某类污染物。直接把零替换成所有正值最小值的百分之一是最常用的做法但必须清醒地知道替换值会影响效率结果替换得越小这个单元在非期望产出上的表现就被夸得越极端。所以替换后我通常还要做一次敏感性分析看看效率排名是否稳定。负值相对少见一般出现在增长率或利润这类指标里。碰到负值我建议做正平移也就是整列统一加上一个常数使最小值变成某个正数。注意这个平移会改变相对关系解释结论时要额外谨慎。缺失值就更简单了面板数据我习惯于前向插值如果首年缺失就用相邻年份均值做完之后专门写一段文字记录填充规则因为审稿人大概率会问。3. SBM模型线性化与MATLAB求解代码逐段拆解3.1 分式规划怎么变成线性规划带非期望产出的SBM模型基本形式是一个分式规划ρ (1 - (1/m) Σ(s_i⁻ / x_ik)) / (1 (1/(s1s2)) (Σ(s_r^g / y_rk^g) Σ(s_r^b / y_rk^b)))分子衡量投入的平均可缩减程度分母衡量期望产出与非期望产出的平均可扩张或可缩减程度。这个形式本身是非线性的不能直接用linprog求解。解决办法是Charnes-Cooper变换引入辅助变量t把所有变量都乘上t令Λ_j t λ_jS t s然后强行让分母等于1目标函数就只剩分子的线性表达式。我一直觉得这个过程像做“单位换算”原来是个比值不好优化现在把分母固定成1把比值问题变成“分子最小化”问题技术含量不高但很巧。实际实现时变量数会增加约束也更复杂但换来的是全局最优解和极快的求解速度。这也是我坚持用linprog而不是fmincon的原因——非线性求解器既慢又不保证全局最优在效率评价这种要跑几千次循环的场景里完全不可取。3.2 变量排列和约束矩阵在MATLAB里用linprog首先要确定每个决策变量在向量里的位置。我的习惯是固定排列顺序投入松弛、期望产出松弛、非期望产出松弛、所有DMU的Λ、辅助变量t。有了固定顺序构造约束矩阵才不会乱。等式约束来自四条逻辑。第一条每个投入要素的加权参考组合加上投入松弛等于被评价单元的投入乘t。第二条每个期望产出的加权参考组合减去期望产出松弛等于被评价单元的期望产出乘t。第三条每个非期望产出的加权参考组合加上非期望松弛等于被评价单元的非期望产出乘t。第四条是归一化等式也就是前面说的让分母等于1。如果选择VRS可变规模报酬还要额外加一个约束所有Λ之和等于t。这个约束的作用是允许参考前沿包含规模效应很多实证论文都默认用VRS因为CRS假设所有单元都在最优规模下经营这在现实中太苛刻了。3.3 完整代码与调用下面是我实际项目里在用的SBM函数做了一点脱敏和简化处理核心逻辑完整保留function rho sbm_eff(X, YG, YB, k, vrs_flag) % 常规SBM效率值计算含非期望产出 % X : n x m 投入矩阵 % YG : n x s1 期望产出矩阵 % YB : n x s2 非期望产出矩阵 % k : 被评价DMU序号 % vrs_flag : true - VRS, false - CRS [n, m] size(X); [s1, s2] size(YG, 2), size(YB, 2); assert(size(YG, 1) n size(YB, 1) n, 样本数不一致); x0 X(k, :); yg0 YG(k, :); yb0 YB(k, :); nvar m s1 s2 n 1; idx_Si 1 : m; idx_Sg m 1 : m s1; idx_Sb m s1 1 : m s1 s2; idx_lam m s1 s2 1 : m s1 s2 n; idx_t nvar; f zeros(nvar, 1); f(idx_Si) -1 ./ (m .* x0); f(idx_t) 1; Aeq []; beq []; % 1. 投入约束: sum_j lam_j*x_ij S_i t*x_ik for i 1:m row zeros(1, nvar); row(idx_Si(i)) 1; row(idx_lam) X(:, i); row(idx_t) -x0(i); Aeq [Aeq; row]; beq [beq; 0]; end % 2. 期望产出约束: sum_j lam_j*yg_rj - S_g t*yg_rk for r 1:s1 row zeros(1, nvar); row(idx_Sg(r)) -1; row(idx_lam) YG(:, r); row(idx_t) -yg0(r); Aeq [Aeq; row]; beq [beq; 0]; end % 3. 非期望产出约束: sum_j lam_j*yb_rj S_b t*yb_rk for r 1:s2 row zeros(1, nvar); row(idx_Sb(r)) 1; row(idx_lam) YB(:, r); row(idx_t) -yb0(r); Aeq [Aeq; row]; beq [beq; 0]; end % 4. VRS约束: sum lam_j t if vrs_flag row zeros(1, nvar); row(idx_lam) ones(1, n); row(idx_t) -1; Aeq [Aeq; row]; beq [beq; 0]; end % 5. 归一化约束: t (1/(s1s2)) * (sum S_g/y_g sum S_b/y_b) 1 row zeros(1, nvar); row(idx_t) 1; row(idx_Sg) 1 ./ ((s1 s2) .* yg0); row(idx_Sb) 1 ./ ((s1 s2) .* yb0); Aeq [Aeq; row]; beq [beq; 1]; lb zeros(nvar, 1); lb(idx_t) 1e-8; options optimoptions(linprog, Display, off, Algorithm, dual-simplex); [x, fval, exitflag] linprog(f, [], [], Aeq, beq, lb, [], [], options); if exitflag 0 warning(DMU %d SBM求解失败, exitflag%d, k, exitflag); rho NaN; else rho fval; end end调用方式非常直观。比如计算第k个决策单元在VRS假设下的SBM效率rho sbm_eff(X, YG, YB, 5, true);需要注意的是目标函数的最优值fval直接就是SBM效率值不需要再额外换算。我在第一次写的时候也犹豫过后来用一个小数据集和MaxDEA的结果对拍过效率值完全对得上。4. 超效率SBM的MATLAB实现把效率值和排名一起解决4.1 把自己从参考集里剔除之后发生了什么超效率SBM的思想一句话就能说清楚不让被评价单元拿自己当参照物。常规SBM里被评价的单元k出现在参考集里最坏情况下它能选择一个和自己完全相同的组合所以效率值最多到1就封顶了。这是“自我评价”的天然缺陷。超效率模型把k从参考集里拿掉用剩下的n-1个单元构建一个新的生产前沿。如果k是有效的那么这个新前沿一定在k的“下方”或者“内侧”。k需要把自己的投入放大、好产出缩小、坏产出放大到一定程度才会掉进新前沿构成的生产可能集里。这个放大缩小的程度就是超效率值。我常用的类比是跳远比赛运动员跳出的距离必须超过起跳线才算有效成绩。如果你把起跳线往后挪还能跳出同样的距离说明你的成绩含金量更高。超效率本质上就是“把起跳线往后挪”。4.2 超效率SBM的变量与约束超效率SBM的计算里虚拟投入Xbar、虚拟期望产出Ybar_g、虚拟非期望产出Ybar_b都成了决策变量。这些变量描述了一个“如果k变差到还能保持有效最差能差到什么程度”的临界状态。线性化之后的超效率SBM是这样一个线性规划目标函数尽量最小化虚拟投入相对真实投入的比值宇宙靠一堆不等式约束来限定虚拟变量和参考集的关系。我不在这里把全部公式展开成论文格式直接看代码更容易理解。4.3 完整代码function se supersbm_eff(X, YG, YB, k, vrs_flag) % 超效率SBM效率值计算剔除自身参考集 % X : n x m 投入矩阵 % YG : n x s1 期望产出矩阵 % YB : n x s2 非期望产出矩阵 % k : 被评价DMU序号 % vrs_flag : true - VRS, false - CRS [n, m] size(X); [s1, s2] size(YG, 2), size(YB, 2); x0 X(k, :); yg0 YG(k, :); yb0 YB(k, :); J setdiff(1:n, k); nJ length(J); nvar m s1 s2 nJ 1; idx_Xi 1 : m; idx_Yg m 1 : m s1; idx_Yb m s1 1 : m s1 s2; idx_lam m s1 s2 1 : m s1 s2 nJ; idx_t nvar; % 目标: min sum(Xbar_i / x_ik) / m f zeros(nvar, 1); f(idx_Xi) 1 ./ (m .* x0); Aeq []; beq []; % 归一化约束 row zeros(1, nvar); row(idx_Yg) 1 ./ ((s1 s2) .* yg0); row(idx_Yb) 1 ./ ((s1 s2) .* yb0); Aeq [Aeq; row]; beq [beq; 1]; % VRS if vrs_flag row zeros(1, nvar); row(idx_lam) ones(1, nJ); row(idx_t) -1; Aeq [Aeq; row]; beq [beq; 0]; end A []; b []; % 1) Xbar_i sum_j lambda_j * x_ij for i 1:m row zeros(1, nvar); row(idx_Xi(i)) -1; row(idx_lam) X(J, i); A [A; row]; b [b; 0]; end % 2) Ybar_g sum_j lambda_j * yg_rj for r 1:s1 row zeros(1, nvar); row(idx_Yg(r)) 1; row(idx_lam) -YG(J, r); A [A; row]; b [b; 0]; end % 3) Ybar_b sum_j lambda_j * yb_rj for r 1:s2 row zeros(1, nvar); row(idx_Yb(r)) -1; row(idx_lam) YB(J, r); A [A; row]; b [b; 0]; end % 4) Xbar_i t * x_ik for i 1:m row zeros(1, nvar); row(idx_Xi(i)) -1; row(idx_t) x0(i); A [A; row]; b [b; 0]; end % 5) Ybar_g t * yg_rk for r 1:s1 row zeros(1, nvar); row(idx_Yg(r)) 1; row(idx_t) -yg0(r); A [A; row]; b [b; 0]; end % 6) Ybar_b t * yb_rk for r 1:s2 row zeros(1, nvar); row(idx_Yb(r)) -1; row(idx_t) yb0(r); A [A; row]; b [b; 0]; end lb zeros(nvar, 1); lb(idx_t) 1e-8; options optimoptions(linprog, Display, off, Algorithm, dual-simplex); [x, fval, exitflag] linprog(f, A, b, Aeq, beq, lb, [], [], options); if exitflag 0 warning(DMU %d 超效率SBM求解失败, exitflag%d, k, exitflag); se NaN; else se fval; end end说一个容易忽略的细节超效率SBM里参考集大小变成了n-1。如果你的样本量本来就很小比如只有5个单元剔除一个之后参考集只剩4个求出来的前沿会很敏感超效率值可能异常大或者直接无解。我的经验是样本量至少要在指标总数的三倍以上超效率结果才比较稳定。4.4 与常规SBM的衔接实际项目里我从来不对所有单元直接跑超效率而是先做一次常规SBM筛选出效率值在1附近的单元再对这些单元跑超效率。这样有几个好处效率低下的单元在常规SBM阶段就可以给出具体的投入冗余和产出不足分析有效单元进入超效率阶段后输出的是一个可排名的数值报告里能直接写进表格。衔接代码也很短rho sbm_eff(Xt, YGt, YBt, k, true); if ~isnan(rho) abs(rho - 1) 1e-6 eff(k, t) supersbm_eff(Xt, YGt, YBt, k, true); else eff(k, t) rho; end注意判断条件里我用了abs(rho - 1) 1e-6而不是rho 1。浮点数运算里很多效率值其实是0.9999999或者1.0000001直接判等会漏掉一批有效单元。5. 面板数据循环求解、结果导出和踩坑记录5.1 主循环怎么写面板数据的主循环结构我一直严格保证两层外层按年份内层按决策单元。三维数组X、YG、YB按年份切片出来以后传给SBM函数得到的结果存进一个n行T列的矩阵。n size(X, 1); T size(X, 3); eff NaN(n, T); options optimoptions(linprog, Display, off, Algorithm, dual-simplex); for t 1:T Xt X(:, :, t); YGt YG(:, :, t); YBt YB(:, :, t); for k 1:n rho sbm_eff(Xt, YGt, YBt, k, true); if ~isnan(rho) abs(rho - 1) 1e-6 eff(k, t) supersbm_eff(Xt, YGt, YBt, k, true); else eff(k, t) rho; end end end这里有一个性能细节optimoptions不要写在每个DMU循环里。我在早期版本里就犯过这个错误每求解一次就重新建一次options对象导致30个单元跑一年数据都要好几秒。把它提到循环外面效率能提升不少。如果你的数据规模很大内层循环还可以改成parfor但要注意parfor里调用的函数必须都在静态路径上并且每个worker都需要能访问数据矩阵。5.2 结果导出与可视化算出来的eff矩阵行列含义分别是“决策单元”和“年份”直接写回文件就行。out array2table(eff, VariableNames, string(years)); out.Properties.RowNames dmu_names; writetable(out, SBM_results.csv, WriteRowNames, true);我一般还会顺手画一张图把每个决策单元的效率值随着年份变化的趋势画成带点的折线图。要是单元太多就挑前五名和后五名单独画图上一眼就能看出谁在改善、谁在恶化。figure; plot(years, eff(1:5, :), o-, LineWidth, 1.2); xlabel(年份); ylabel(效率值); legend(string(dmu_names(1:5)), Location, best);这张图放到汇报材料里比密密麻麻的表格直观得多。5.3 常见报错和排查链路我在实际项目中遇到过下面几类高频问题排查思路基本固定。第一类是exitflag -2线性规划不可行。出现这种问题的原因通常是数据里有负值或零值导致约束互相矛盾。排查时先跑一遍数据体检脚本看看有没有异常值如果数据全正但依然不可行再检查是不是某个DMU的某一列指标值恒等于0或者恒等于同一个常数这种情况会让约束矩阵退化。第二类是超效率值出现Inf。这个基本上都是因为参考集太小或者某个产出值趋近于0。比如一个单元的某类非期望产出几乎是0超效率公式里分母里出现了极小的数效率值就被拉得很大。处理方式不是直接删掉这个单元而是回到原始数据检查这个指标是否真的有区分度必要时考虑合并指标。第三类是常规SBM效率结果全部接近1。我遇到过一次后来发现是把非期望产出错当成了期望产出放进模型等于模型里鼓励污染增加所有单元自然都显得很高效。修正指标归类后结果马上就正常了。所以每次建模前我都会花十分钟检查指标方向绝不跳过。第四类是循环里出现了NaN但不报错。原因很可能是某个年份的数据缺失切片后矩阵里带空值求解器直接返回了异常结果。建议在主循环开始前加一行检查assert(all(isfinite(X(:))))把问题挡在模型前面。还有一个小技巧如果样本量特别大比如几千个DMUMATLAB里循环求解确实会慢。可以考虑把每个DMU的线性规划独立包装然后用parfor并行。但并行池的启动开销在数据量不够大的时候反而更慢我一般是在DMU数量超过500的时候才启用。这套流程做完剩下的事情就是写报告了。常规SBM给出低效率单元的松弛变量超效率SBM给出有效单元的排序面板循环给出时间趋势三个层次凑齐基本上评审关心的效率评价问题都能覆盖。我自己跑了两年多这套代码最大的感受是它比商业软件灵活太多加约束、改模型、批量处理、和后续计量回归衔接全都自由。如果你也做效率评价相关的工作不妨从SBM开始把MATLAB这套流程搭建起来后面再扩展任何变体都会轻松很多。