简介这份MATLAB资源围绕数据包络分析DEA设计内置DEA-Solver工具箱面向经管、运筹及系统工程领域的研究者、高校师生和效率评价实践者。压缩包共21个文件涵盖.m脚本、.mat与.xls数据文件、PDF文献和说明文档总大小2.54MB结构紧凑但覆盖了模型定义、数据处理、求解和结果输出等完整环节。已有443人学习下载。借助DEA-Solver使用者能够快速构建CCR、BCC等主流DEA模型通过示例数据与参考论文对照学习效率评估流程还可依据实际需求修改输入输出指标进行定制化分析。资源内含聚类DEA方法、区域创新政策排序等PDF文档及配套Excel数据可辅助深入理解前沿应用场景教育工作者也能直接用于课堂教学演示。需注意运行前请确认MATLAB优化工具箱已正确安装以避免大型数据集求解时出现计算问题。1. 内含DEA Solver的MATLAB程序先看建模前提再谈运行拿到一套“内含DEA Solver”的MATLAB程序你第一时间要做的不是点运行而是先想清楚你手里的数据适不适合跑DEA。数据包络分析DEA解决的是多投入、多产出情境下的效率测算问题它不需要预先设定权重也不需要假设生产函数直接用线性规划构造出一条效率边界。这个特性让它在银行网点绩效、医院运行效率、高校院系评估等场景里被广泛使用。但反直觉的结论是越不需要参数设定的方法对数据质量越敏感。指标选错、量纲失控、样本量不足都会让Solver输出一份漂亮的效率表却完全不可信。这篇文章会把DEA Solver的模型选型、最小MATLAB实现、结果解读和常见坑一次讲透。2. DEA建模三选规模报酬、优化导向与投入产出指标的取舍DEA Solver本身只是一套线性规划求解流程真正决定结果价值的是建模前的三个选择用CCR还是BCC、选投入导向还是产出导向、投入产出指标怎么定。这三个问题没有想清楚后面无论代码多漂亮结论都会失真。2.1 CCR还是BCC规模报酬假设直接决定效率口径CCR模型假设规模报酬不变意思是投入翻倍产出也应该翻倍被评价的所有决策单元被放在同一个规模水平上比较。BCC模型则允许规模报酬可变每个单元可以和与自己规模相近的单元比较。两者的核心区别在于CCR算出的效率包含了规模因素的影响BCC算出的效率剔除了规模影响。我一般会建议先跑BCC再拿CCR结果除以BCC结果得到规模效率。这样一次分析能同时看到三件事综合技术效率、纯技术效率、规模效率。如果只报CCR一个规模偏小但运营方式很好的网点效率值会被拉低看起来像管理问题实际是规模问题。某高校课题组用CCR评了十几个院系效率普遍只有0.6左右换成BCC后大半院系纯技术效率在0.9以上结论完全反转。这就是选错模型的代价。具体选择上没有绝对标准实践里的经验是如果评价对象是同类生产线、同类设备这类规模差异小的场景CCR够用如果是分行、网点、医院、院系这类规模差异大的场景以BCC为主并补充规模效率报告。最稳妥的做法是用代码同时计算两种模型输出效率分解表让读者自己看差异在哪里。2.2 投入导向还是产出导向谁可控谁就是优化方向投入导向回答的问题是保持产出不变投入最多能压缩多少。产出导向回答的问题是保持投入不变产出还能扩张多少。对于银行网点柜员人数、场地面积容易调整优先考虑投入导向对于科研机构经费投入相对固定课题产出增长空间大产出导向更贴合实际。从求解结果看效率值为1的有效单元在两种导向下结果一致无效单元则不同。容易踩坑的是把两种导向的结果混在一起比较。一个在投入导向下效率0.85的经营单元在产出导向下可能是0.91这不是出错了只是优化视角不同。我习惯在结果表的备注列里写明模型和导向避免自己和后来接手数据的人产生误解。如果发现某个DMU在两种导向下效率差异非常大说明这个单元在数据里比较“极端”要回去检查原始指标。2.3 投入产出指标设计的三个硬性要求指标设计是DEA里最接近“玄学”的部分。Solver来者不拒你给什么指标它就算什么结果所以指标质量必须自己把关。三条硬性要求一是同向性。投入增加时产出不能下降。可以在建模前用corr函数做一次相关性检查如果出现投入与产出显著负相关要么指标定义反了要么这组数据本身不适合DEA。二是样本量要够。经验法则是决策单元数量不少于投入产出指标总数之和的2倍3倍更稳。4个指标配8个DMU几乎所有单元都会变成有效解释力接近于零。三是指标不要冗余。两个产出指标相关系数超过0.9比如“总产值”和“销售收入”本质上在测同一件事。把它们同时放进去相当于变相给这个产出维度更高的权重会虚高效率值。冗余指标宁可删掉一个也不要保留。零值和负值也需要处理。标准DEA模型要求投入产出为正数。某个产出指标为0时常见做法是用一个很小的正数如1e-6替换或者重新设计指标口径比如把“利润”换成“利润率”加常数平移。平移会改变效率值所以要在结果报告里注明数据的预处理方式否则别人无法复现你的结果。3. 用linprog实现最小DEA Solver从Excel数据到效率值的完整路径这一章进入代码。核心思路并不复杂DEA对每个决策单元都要单独求解一个线性规划MATLAB自带的linprog函数可以完成这一步。你把投入产出矩阵准备好对每个DMU调用一次linprog记录目标函数值就是效率值。3.1 数据准备把Excel表格整理成X矩阵和Y矩阵我通常用Excel管理原始数据。格式按行放DMU第一列是DMU名称中间列依次放投入和产出指标表头用英文不包含空格。读取时用readtable配合ReadRowNames选项把第一列识别成行名。% 读取数据假设Excel中第一列为DMU名称 T readtable(dea_data.xlsx, ReadRowNames, true); % 按列名取数避免按位置取数时列顺序错位 % X矩阵m行 n列每一列是一个DMUY矩阵同理 X table2array(T(:, {labor, capital})); % 两项投入转置为 2×n Y table2array(T(:, {output, quality})); % 两项产出转置为 2×n % 快速检查维度和缺失值 disp(size(X)); disp(size(Y)); assert(~any(isnan(X(:))) ~any(isnan(Y(:))), 存在NaN请先清洗数据);这里做了两次准备动作用列名而不是列位置取数据防止Excel列顺序变动导致投入产出错位取完数据后立即转置因为Solver接口设计成“每一列是一个DMU”而Excel习惯“每一行是一个DMU”。检查维度和NaN本质上是给数据质量上保险DEA对异常值敏感缺失值会让linprog直接退出。之前某开发者把Excel里一个空单元格读成NaN求解器返回全NaN排查了半天才发现是数据问题。3.2 CCR投入导向的核心求解代码θ、λ与约束矩阵的构造这里给出一个最小可用的DEA Solver核心函数。变量向量z的第一个分量是效率值θ后面n个分量是λ表示每个参考DMU被赋予的权重。约束条件分两组上半部分是投入约束含义是所有被参考DMU的投入组合不能超过θ倍当前DMU的投入下半部分是产出约束含义是参考DMU的产出组合不能低于当前DMU的产出。function res dea_solver(X, Y, model) % 最小DEA Solver投入导向 % X: 投入矩阵m×n每列一个DMU % Y: 产出矩阵s×n每列一个DMU % model: CCR 或 BCC [m, n] size(X); [s, ~] size(Y); res struct(theta, {}, lambda, {}, exitflag, {}); for k 1:n x0 X(:, k); % 当前DMU的投入 y0 Y(:, k); % 当前DMU的产出 % 变量 z [theta; lambda_1; ...; lambda_n] f [1; zeros(n, 1)]; % 目标函数最小化 theta % 投入约束 -x0 * theta X * lambda 0 A1 [-x0, X]; % m 行每行对应一个投入指标 b1 zeros(m, 1); % 产出约束 -Y * lambda -y0 A2 [zeros(s, 1), -Y]; % s 行每行对应一个产出指标 b2 -y0; A [A1; A2]; b [b1; b2]; % 变量下界theta 和 lambda 均非负 lb zeros(n 1, 1); ub []; % BCC模型增加凸性约束sum(lambda) 1 if strcmpi(model, BCC) Aeq [0, ones(1, n)]; beq 1; else Aeq []; beq []; end options optimoptions(linprog, Display, off); [z, ~, exitflag] linprog(f, A, b, Aeq, beq, lb, ub, options); res(k).theta z(1); res(k).lambda z(2:end); res(k).exitflag exitflag; end end代码里的几个参数需要解释清楚。f向量首元素是1因为目标函数只是最小化θλ在目标里没有系数A矩阵的行数是m加s正好对应每个投入和产出拆出来的行lb全零是因为θ和λ在经济学含义上都应为非负θ为0意味着理论上可以完全不消耗投入这在现实中几乎不会出现但对求解器来说是一个安全下界。linprog默认算法在较新版本的MATLAB中输入参数是optimoptions如果你的版本不识别把它换成optimset(Display,off)即可其他部分不用动。跑通CCR后你可以这样调用并查看效率向量res dea_solver(X, Y, CCR); theta [res.theta]; % 1×n 的效率向量 disp(theta);每个DMU得到一组λ权重权重和不一定等于1这正是CCR的特点。如果你想同时支持产出导向需要把投入产出约束做对偶变换但新手阶段建议先把投入导向跑熟产出导向的差异留到结果解读时再研究。3.3 从CCR扩展到BCC加一个约束就够了BCC和CCR的代码差异小到只有两行Aeq [0, ones(1, n)]beq 1。这个约束要求所有λ之和等于1表示参考组合的规模与当前DMU的规模保持可比。加了它效率值就不再包含规模因素。调用方式变成res_bcc dea_solver(X, Y, BCC); theta_bcc [res_bcc.theta];从CCR和BCC的效率值还可以推导规模效率theta_ccr [dea_solver(X, Y, CCR).theta]; theta_bcc [dea_solver(X, Y, BCC).theta]; se theta_ccr ./ theta_bcc; % 规模效率数值越接近1说明规模越合适规模效率小于1且BCC效率高说明这个DMU的问题是规模不当而非运营低效。很多报告把CCR和BCC同时列出再附一个se列读者一眼就能看出效率损失来自规模还是来自管理。4. 读懂DEA Solver的输出有效判定、参考集与改进方向的计算求解完成后你手里有三个关键输出效率值θ、权重向量λ、求解状态exitflag。这三样东西组合起来能回答一个完整的问题谁有效、谁无效、无效的单元向谁学习、要改进多少。只报一个θ等于浪费了大半个模型。4.1 从效率和松弛判断DEA有效性判断一个DMU是否有效不能只看θ是否为1。严格的有效要满足两个条件θ等于1且所有松弛变量为0。θ等于1但存在松弛称为弱有效意思是虽然已经在前沿面上但某个指标还有压缩空间。用程序判断时需要根据λ计算投入和产出的投影值再拿投影值与原始值对比% 对第k个DMU计算投入冗余和产出不足 projX X * res_bcc(k).lambda; % 投入方向的参考组合 projY Y * res_bcc(k).lambda; % 产出方向的参考组合 slackX X(:, k) - projX; % 大于0表示投入冗余 slackY projY - Y(:, k); % 大于0表示产出不足这里没有做第二阶段线性规划属于简化版松弛分析但在Solver原型阶段足够用。真正精细的分析需要在第一阶段结果基础上再解一次线性规划来单独最大化松弛变量大部分商业软件的“松弛改进量”就是那个口径。4.2 用lambda找参考集无效单元的追赶标杆λ向量里非零分量对应的DMU就是当前DMU的参考集。比如某分行B的效率是0.84λ在A分行和C分行上分别是0.3和0.7说明B的参考组合是A和C的加权混合体。参考集的实践价值在于它能告诉被评价者“你应该向谁看齐”而不是笼统地说“你效率低”。下面是一份虚构的结果示例。H2的效率0.742参考集是H1和H4说明H2的前沿位置由这两个分行决定它改进的方向就是向H1和H4的组合看齐。DMUBCC效率Σλ规模报酬参考集人工费冗余H11.0001.000不变H10H20.7420.618递增H1、H431.5万H31.0001.000不变H30H40.9031.112递减H1、H312.2万规模报酬一列在BCC投入导向下这样判断Σλ略小于1是规模报酬递增说明扩大规模有利Σλ大于1是规模报酬递减说明规模偏大。这也是为什么每次运行后都要把λ求和的直接原因。4.3 把效率结果批量导出成报告散点式的结果不利于汇报我习惯把结果汇总成一个表格直接写回ExceloutTable table(T.Properties.RowNames, theta_bcc, se, VariableNames, ... {DMU, TE_BCC, SE}); writetable(outTable, dea_result.xlsx);写回文件之前建议把模型设定也写进表里或附页里。行业里复现DEA结果失败十有八九是因为对方报告里没写清楚用的是CCR还是BCC、投入导向还是产出导向。这不是Solver的锅是记录习惯的锅。记住一个格式模型名、导向、指标数、样本量每一项都不能省。5. 跑DEA的常见问题排查全为1、NaN结果和与文献对不上DEA求解器一旦跑出异常结果绝大多数问题不是出在线性规划算法而是出在数据准备或者维度配置上。这里把最常撞上的五个问题按现象、原因、解决三个步骤写清楚按顺序排查通常几分钟就能定位。5.1 效率结果全为1不是模型强大而是数据构造失误现象跑完CCR所有DMU的效率全是1.0000。看起来像大团圆但做效率分析的人都明白这种情况下结果没有区分度等于白做。原因三个高频来源。一是X矩阵里混入了DMU编号列readtable读进来的第一列如果你没有用ReadRowNames序号会被当成一个投入指标二是DMU数量太少比如8个DMU却放了5个指标线性规划的自由度太大每个DMU都能找到一套只对自己有利的权重三是数据里有大量重复行相当于把同一个DMU复制多份。解决先打印X和Y的size确认维度是否和设计一致。再用数组可视化快速看一眼数据分布比如imagesc(X)如果出现整行整列的规律性图案就要回看数据表。最后执行指标数量控制手册DMU数量应至少是指标总数的2倍。某开发者遇到过12个DMU、6个指标结果9个有效缩到4个指标后正常了。5.2 linprog返回不可行或NaN先查数据再查矩阵维度现象exitflag小于0theta是NaN控制台报错“Linprog stopped because no feasible point was found”。原因DEA模型理论上一定有可行解因为取θ等于1、当前DMU自己权重为1、其他权重为0就满足约束。所以一旦报不可行基本可以断定是数据里有NaN或Inf或者是负的产出值让约束方向失去了可行域。解决在调用linprog之前加一段强力检查assert(~any(isnan(X(:))) ~any(isinf(X(:))), X中包含NaN或Inf); assert(~any(isnan(Y(:))) ~any(isinf(Y(:))), Y中包含NaN或Inf); assert(min(Y(:)) 0, 产出指标必须为正数需要预处理);把零值产出用1e-6这种小正数替代或者重构指标。别用NaN填补DEA没有缺失值插补的概念。5.3 投入产出颠倒造成结果错乱列名比列位置可靠现象算出的效率值分布完全不符合直觉比如公认表现好的DMU效率反而低。把X、Y的尺寸打出来看数值范围也不对劲。原因table2array按列位置取值如果Excel的列顺序调整过比如新加了一列你原来写死的第2到第4列就全错了投入产出交换位置后DEA的约束方向失效结果自然离谱。解决按列名取列放弃列位置。代码切到3.1节那种T(:, {labor, capital})写法。列名保持英文无空格避免readtable在中文列名上做变量名转换时引入麻烦。这个习惯能省下大量复查时间。5.4 结果和文献或商业软件对不上先对齐设置再比数值现象同一份数据你算出来的CCR效率是0.8511文献里写0.8372。数值差了一点但排序基本一致。原因三种常见差异。第一是模型导向不一致对方用产出导向你用了投入导向无效单元的效率值本来就不同第二是对方做了两阶段松弛分析第一阶段目标函数一样但第二阶段重新算过松弛效率值虽然没有变投影量和冗余量会不同第三是求解器精度差异linprog的收敛容差和部分商业求解器不同会让结果在小数点后两位附近浮动。解决先核对模型名、导向、数据预处理方式再比较排序一致性。效率值的零点零几差异通常不改变管理结论不用过度纠结。在报告里写明自己的求解器设置别人才能复现。5.5 指标总数一多有效单元比例飙升识别能力被稀释现象指标从3个加到6个有效DMU数量从3个变成8个。原因每个指标都在给DMU提供新的“维度自由”DMU可以在某个指标上突出从而获得更高的效率值。指标数接近DMU数时DEA的判别力消失结果趋于全有效。这是维度灾难的典型表现不是模型bug。解决严格执行指标数量约束n至少大于2倍的m加s。如果业务上确实需要保留很多指标先做相关性分析和主成分压缩把强相关的指标合并再进入DEA。不要为了“全面”牺牲区分度。6. 进阶扩展超效率DEA的MATLAB实现与有效单元排序标准DEA有个让很多人头疼的局限效率值为1的DMU无法继续比较等于一群满分学生并列第一。超效率模型解决这个问题思路是把被评价的DMU从参考集里剔除让它不能再“自己和自己比”然后用剩下的DMU组成前沿面。此时有效单元的θ会大于1数值越大说明它在保持有效的条件下越有冗余空间排序自然出来了。在dea_solver上做扩展核心改动是构造参考矩阵时去掉当前列。下面是关键修改逻辑function theta_super dea_super(X, Y, k) % 超效率DEA只评价第k个有效DMU [m, n] size(X); [s, ~] size(Y); refIdx true(1, n); refIdx(k) false; % 关键一步把自己从参考集剔除 Xr X(:, refIdx); Yr Y(:, refIdx); nr sum(refIdx); x0 X(:, k); y0 Y(:, k); f [1; zeros(nr, 1)]; A1 [-x0, Xr]; % 注意列数变成1nr A2 [zeros(s, 1), -Yr]; b [zeros(m, 1); -y0]; lb zeros(nr 1, 1); options optimoptions(linprog, Display, off); [z, ~, ~] linprog(f, [A1; A2], b, [], [], lb, [], options); theta_super z(1); end这段代码的前提是被评价的DMU在标准DEA下已经有效。如果无效DMU也放进超效率模型它的效率值不会大于1没有额外信息。所以正确流程是先跑一遍标准DEA筛出θ等于1的单元再逐个调用dea_super。得到超效率值之后配合原有的BCC效率就能给所有DMU一个完整的排名。我在实际项目里习惯把有效单元的超效率值单独列一列并在旁边注明“大于1表示在保持DEA有效前提下的投入容忍度”。数值大的单元不只是排名靠前它还说明这个单元的抗干扰能力强即使投入小幅增加也未必跌出前沿面。使用超效率模型的一个重要提醒是当DMU数量很少时剔除自己后的参考集可能退化导致超效率值异常大。我的习惯是先看参考集大小至少保留3个DMU再做超效率否则结果不稳定。这个扩展虽然只有十几行代码但让DEA分析从“给结论”变成了“给排名”是目前我在效率评价项目里最常用的升级动作。无论你是跑通标准DEA还是准备做超效率都建议先对数据做同向性检查保存好原始指标列表和预处理记录。这些工作看起来琐碎却是整个分析里最不玄学、也最救命的习惯。希望帮到你。本文还有配套的精品资源点击获取