1. 项目概述为什么数据包络分析是数学建模的“效率裁判”如果你正在备战数学建模竞赛尤其是涉及到评价类、效率分析类的题目比如评价多个银行的经营效率、比较不同医院的资源利用水平或者分析各省份的能源消耗绩效那么“数据包络分析”这个工具你绝对不能错过。它不像回归分析那样需要一个预设的函数形式也不像主成分分析那样需要人为赋权DEA的魅力在于它能让数据自己“说话”客观地找出谁是那个“优等生”以及其他人差在哪里。简单来说DEA就是一个“效率裁判”。它通过构建一个“生产前沿面”把所有被评价的决策单元比如那些银行、医院、省份都放在同一个坐标系里。那些落在前沿面上的单元就是相对有效的效率值为1而那些落在前沿面内部的单元就是相对无效的效率值介于0和1之间。更重要的是DEA不仅能给你一个效率分数还能告诉你无效的单元具体差在哪儿是投入太多了还是产出太少了应该向哪个有效的“榜样”学习这些信息对于提出改进建议至关重要非常契合数学建模论文需要“问题分析-模型构建-结果解读-政策建议”的逻辑链条。在实战中无论是国赛、美赛还是其他建模赛事DEA都是解决多输入多输出效率评价问题的利器。而MATLAB凭借其强大的矩阵运算能力和丰富的优化工具箱是实现DEA模型编程求解的不二之选。接下来我就结合自己多次带队和评审的经验拆解一下如何用MATLAB从零开始实现DEA并分享一些让论文脱颖而出的关键技巧。2. DEA核心模型选择与数学原理拆解DEA模型家族庞大但最常用、最适合入门数学建模的主要是两个CCR模型和BCC模型。选对模型是正确分析的第一步。2.1 CCR模型规模收益不变的“全能基准”CCR模型由Charnes, Cooper和Rhodes提出它假设生产过程是规模收益不变的。这意味着如果你把投入翻倍产出也预期能翻倍。这个假设比较强适用于那些规模大小不太影响效率的场景或者在你不太确定规模影响时作为一个基础分析。它的核心是为一个特定的决策单元假设是第k个求解一个线性规划问题。我们假设有n个决策单元每个单元有m种投入和s种产出。投入记为 $X_{ij}$ (i1,...,m; j1,...,n) 比如员工数、资金、能耗。产出记为 $Y_{rj}$ (r1,...,s; j1,...,n) 比如利润、服务人数、论文数量。对于第k个单元CCR模型投入导向的线性规划形式如下目标函数最大化第k个单元的效率值 $\theta_k$。约束条件产出约束所有单元的线性组合权重为 $\lambda_j$的产出必须不少于第k个单元的产出。这保证了比较的基准至少不比被评价单元差。投入约束所有单元的线性组合的投入必须不大于第k个单元投入的 $\theta_k$ 倍。这里的 $\theta_k$ 就是我们要算的效率值它表示在保持产出不减少的情况下投入最多能按比例缩减多少。权重非负$\lambda_j \geq 0$。用公式表达更清晰 $$ \begin{align*} \text{Max } \theta_k \ \text{s.t. } \sum_{j1}^{n} \lambda_j Y_{rj} \geq Y_{rk}, \quad r1,...,s \ \sum_{j1}^{n} \lambda_j X_{ij} \leq \theta_k X_{ik}, \quad i1,...,m \ \lambda_j \geq 0, \quad j1,...,n \ \theta_k \text{ free} \end{align*} $$求解这个线性规划得到的 $\theta_k^$ 就是第k个单元的相对效率。如果 $\theta_k^ 1$且所有松弛变量实际约束与边界之间的差值为0则该单元是DEA有效的如果 $\theta_k^* 1$ 但存在松弛变量不为0则是弱DEA有效如果 $\theta_k^* 1$则是DEA无效。注意CCR模型求出的效率是“技术效率”包含了规模效率的影响。也就是说一个单元效率低可能是因为技术不行也可能是因为规模不合适。2.2 BCC模型规模收益可变的“现实修正”BCC模型由Banker, Charnes和Cooper提出它放松了CCR的假设允许规模收益可变递增、不变或递减。这更符合现实因为很多机构在规模扩大时效率并不会同比例变化。BCC模型只是在CCR模型的约束条件中增加了一个凸性约束$\sum_{j1}^{n} \lambda_j 1$。这个等式约束保证了前沿面是凸的从而可以将技术效率进一步分解为“纯技术效率”和“规模效率”。纯技术效率反映的是在给定规模下管理和技术水平的优劣。由BCC模型求得。规模效率反映的是实际规模与最优生产规模之间的差距。计算公式为规模效率 CCR效率 / BCC效率。模型选择心法优先使用BCC模型在大多数数学建模场景下由于决策单元规模差异往往很大比如从小县城医院到省级三甲医院BCC模型更贴合实际分析结果也更细腻。CCR模型作为补充可以用CCR模型的结果与BCC结果对比计算规模效率分析无效是源于管理不善还是规模不当。看题目暗示如果题目描述中明显暗示“不考虑规模影响”或强调“同比例缩放”可选用CCR如果提到“规模各异”、“挖掘管理潜力”则BCC更佳。3. MATLAB编程实现全流程与核心代码解析理论懂了关键在实现。下面我将以BCC模型投入导向为例展示完整的MATLAB编程步骤。我们会自己编写核心求解函数而不是依赖黑箱工具箱这样更能体现建模能力也方便自定义和调试。3.1 数据准备与标准化处理数据是模型的基石。首先我们需要将数据整理成两个矩阵投入矩阵X(m×n) 和产出矩阵Y(s×n)其中每一列代表一个决策单元。% 假设我们有20个决策单元(DMU)3种投入2种产出 n 20; % DMU数量 m 3; % 投入指标数 s 2; % 产出指标数 % 示例随机数据实际中应从Excel或CSV读取 % X: 投入矩阵 (m x n) 每一列是一个DMU的投入数据 X [randn(m, n)*10 50; randn(m, n)*5 20; randn(m, n)*2 10]; % Y: 产出矩阵 (s x n) 每一列是一个DMU的产出数据 Y [randn(s, n)*15 100; randn(s, n)*8 60]; % 数据标准化非常重要 % DEA对数据的量纲敏感如果投入产出单位差异大如资金“亿元”和员工“人”必须标准化。 % 常用方法均值标准化或极差标准化 for i 1:m X(i, :) X(i, :) / mean(X(i, :)); % 投入除以均值 end for r 1:s Y(r, :) Y(r, :) / mean(Y(r, :)); % 产出除以均值 end实操心得数据标准化是DEA分析前必须做的步骤否则量纲大的指标会主导模型导致结果失真。除了均值标准化也可以使用“除以最大值”或“极差标准化”方法。在论文中一定要写明你采用了哪种标准化方法及原因。3.2 BCC模型线性规划求解函数编写我们将为每一个DMU循环求解一个线性规划问题。MATLAB自带的linprog函数可以很好地完成这个任务。function [theta, lambda, slack_x, slack_y] dea_bcc_input(X, Y, k) % DEA_BCC_INPUT 计算第k个DMU在BCC模型投入导向下的效率 % 输入: % X: 投入矩阵 (m x n) % Y: 产出矩阵 (s x n) % k: 待评价的DMU索引 % 输出: % theta: 效率值 (0 theta 1) % lambda: 权重向量 (n x 1) % slack_x: 投入松弛量 (m x 1) % slack_y: 产出松弛量 (s x 1) [m, n] size(X); [s, ~] size(Y); % 1. 构建线性规划的目标函数系数 f % 目标最大化 theta 在linprog中默认是最小化所以 f [-1; zeros(n,1); zeros(ms,1)] % 变量顺序: [theta; lambda_1,...,lambda_n; s_x^-; s_y^] f [-1; zeros(n, 1); zeros(ms, 1)]; % 2. 构建不等式约束矩阵 A 和向量 b % 约束1: sum(lambda_j * X_ij) s_x^- theta * X_ik (投入约束) % 移项: -theta * X_ik sum(lambda_j * X_ij) s_x^- 0 写成 A*x b 形式 % 即: [-X(:,k), X, eye(m), zeros(m,s)] * [theta; lambda; s_x^-; s_y^] 0 A1 [-X(:,k), X, eye(m), zeros(m, s)]; b1 zeros(m, 1); % 约束2: sum(lambda_j * Y_rj) - s_y^ Y_rk (产出约束) % 移项: [zeros(s,1), Y, zeros(s,m), -eye(s)] * [theta; lambda; s_x^-; s_y^] Y(:,k) % 为了统一成 A*x b 两边乘以-1: -[zeros(s,1), Y, zeros(s,m), -eye(s)] * x -Y(:,k) A2 [zeros(s,1), -Y, zeros(s,m), eye(s)]; b2 -Y(:,k); % 合并不等式约束 A [A1; A2]; b [b1; b2]; % 3. 构建等式约束BCC特有的凸性约束 % sum(lambda_j) 1 Aeq [0, ones(1, n), zeros(1, ms)]; beq 1; % 4. 变量边界约束 % theta 无界 lambda 0, 松弛变量 0 lb [-inf; zeros(nms, 1)]; % theta下界为负无穷其他变量下界为0 ub []; % 无上界 % 5. 调用linprog求解线性规划 options optimoptions(linprog, Display, off); % 关闭求解过程显示 [x, fval, exitflag] linprog(f, A, b, Aeq, beq, lb, ub, options); if exitflag ~ 1 warning(DMU %d 的线性规划求解可能未收敛或失败。, k); theta NaN; lambda NaN(n,1); slack_x NaN(m,1); slack_y NaN(s,1); return; end % 6. 解析结果 theta x(1); lambda x(2:n1); slack_x x(n2:nm1); slack_y x(nm2:end); end3.3 批量计算与结果整合有了单个DMU的计算函数我们就可以循环计算所有单元的效率。% 初始化结果存储 num_dmu n; theta_scores zeros(num_dmu, 1); % 效率值 lambda_weights zeros(num_dmu, num_dmu); % 权重矩阵第j列是DMU j的lambda slack_x_all zeros(m, num_dmu); % 投入松弛 slack_y_all zeros(s, num_dmu); % 产出松弛 reference_set cell(num_dmu, 1); % 参考集存储每个无效DMU的“榜样” for k 1:num_dmu [theta, lambda, slack_x, slack_y] dea_bcc_input(X, Y, k); theta_scores(k) theta; lambda_weights(:, k) lambda; slack_x_all(:, k) slack_x; slack_y_all(:, k) slack_y; % 找出参考集lambda权重显著大于0的DMU ref_idx find(lambda 1e-5); % 设置一个小的阈值避免数值误差 reference_set{k} ref_idx; end % 计算规模效率 (Scale Efficiency) CCR效率 / BCC效率 % 需要先计算CCR效率这里假设我们也有一个dea_ccr_input函数 % theta_ccr ... (调用CCR模型计算) % scale_eff theta_ccr ./ theta_scores; % 将结果整理成表格便于查看和分析 DMU_ID (1:num_dmu); Efficiency_BCC theta_scores; % Efficiency_CCR theta_ccr; % Scale_Efficiency scale_eff; result_table table(DMU_ID, Efficiency_BCC); %, Efficiency_CCR, Scale_Efficiency); disp(BCC模型效率评价结果); disp(result_table); % 找出有效和无效的DMU effective_dmu find(abs(Efficiency_BCC - 1) 1e-6); % 效率值接近1的视为有效 ineffective_dmu setdiff(1:num_dmu, effective_dmu); fprintf(\nDEA有效的DMU有%s\n, num2str(effective_dmu)); fprintf(DEA无效的DMU有%s\n, num2str(ineffective_dmu));3.4 结果可视化让数据说话图表是论文的加分项。我们可以用几种方式可视化DEA结果。% 1. 效率值排序条形图 figure(Position, [100, 100, 800, 400]) subplot(1,2,1) [~, sorted_idx] sort(Efficiency_BCC, descend); barh(1:num_dmu, Efficiency_BCC(sorted_idx), FaceColor, [0.2 0.6 0.8]); set(gca, YTick, 1:num_dmu, YTickLabel, arrayfun((x) sprintf(DMU%d, x), sorted_idx, UniformOutput, false)); xlabel(效率值 (θ)); ylabel(决策单元 (排序后)); title(BCC模型效率值排序); grid on; axis tight; line([1 1], ylim, Color, r, LineStyle, --, LineWidth, 1.5); % 标出效率前沿线 text(1.02, 1, 效率前沿 (θ1), Color, r); % 2. 投入/产出松弛分析热图针对无效DMU subplot(1,2,2) if ~isempty(ineffective_dmu) % 选取前几个无效DMU的投入松弛进行分析 num_to_show min(5, length(ineffective_dmu)); sample_ineff ineffective_dmu(1:num_to_show); slack_data slack_x_all(:, sample_ineff); % 转置以便热图显示 imagesc(slack_data); colorbar; set(gca, XTick, 1:m, XTickLabel, arrayfun((x) sprintf(投入%d, x), 1:m, UniformOutput, false)); set(gca, YTick, 1:num_to_show, YTickLabel, arrayfun((x) sprintf(DMU%d, x), sample_ineff, UniformOutput, false)); xlabel(投入指标); ylabel(无效决策单元); title(无效DMU的投入松弛量热图); % 为每个格子添加数值 [text_x, text_y] meshgrid(1:m, 1:num_to_show); text_labels arrayfun((v) sprintf(%.3f, v), slack_data, UniformOutput, false); text(text_x(:), text_y(:), text_labels(:), HorizontalAlignment, center, Color, w, FontWeight, bold); else text(0.5, 0.5, 所有DMU均有效无松弛量, HorizontalAlignment, center, FontSize, 12); axis off; end sgtitle(数据包络分析(DEA)结果可视化);4. 数学建模实战技巧与论文写作要点在竞赛中仅仅跑出结果是不够的如何将DEA分析过程清晰地呈现在论文里并挖掘出深度才是取胜关键。4.1 指标体系的科学构建这是DEA应用中最容易失分也最能体现水平的地方。投入和产出指标的选择直接决定了评价结论的合理性。构建原则相关性指标必须与评价目标紧密相关。评价医院效率投入可以是“医生数”、“床位数”、“医疗设备价值”产出可以是“年门诊量”、“年手术量”、“患者满意度得分”。同向性投入指标应越小越好或至少不应越大越好产出指标应越大越好。如果某个指标方向不对如“医疗纠纷数”作为产出需要进行正向化处理例如取其倒数或负数。避免内部强相关性投入指标之间、产出指标之间不应有高度的线性相关性如“员工总数”和“技术人员数”可能高度相关否则会导致模型自由度下降评价结果不稳定。可以用相关系数矩阵检查并考虑删除或合并高度相关的指标。数据可获得性理想指标必须有可靠的数据来源支撑。避坑指南切勿堆砌指标。我曾见过有论文为评价区域创新效率列出了十几项投入产出。这会导致大量DMU的效率值都为1即“有效”单元过多失去区分度这种现象称为“维度诅咒”。一般经验是DMU数量至少应为投入与产出指标数量之和的2到3倍。4.2 模型结果的深度解读与政策建议算出效率值只是开始解读才是灵魂。对于有效单元θ1在论文中可以称其为“标杆”或“前沿单元”。分析其共性特征它们是在哪些投入产出配置上做到了最优这能为其他单元提供发展方向。对于无效单元θ1效率值θ表示该单元有多“无效”。例如θ0.8意味着在产出不变的情况下所有投入可以同时减少20%。松弛变量Slacks这是提出具体改进建议的核心依据。投入松弛表示即使按比例缩减后仍有多余的、可进一步减少的投入。例如对某个医院床位数的松弛量为10意味着在现有产出水平下即使按θ比例缩减了总投入仍可再减少10张床位。产出松弛表示在投入不变的情况下可以额外增加的产出。例如患者满意度得分的产出松弛为5分意味着该医院的管理或服务还有提升空间能使满意度增加5分。参考集λ权重0的单元指明了学习对象。例如DMU5的参考集是{DMU2 DMU8}那么就可以建议DMU5研究DMU2和DMU8在人员配置、设备利用等方面的具体做法。在论文中的呈现方式制作综合结果表除了效率值一定要包含主要松弛变量和参考集。DMUBCC效率规模收益投入1松弛投入2松弛产出1松弛参考标杆A1.000不变000自身B0.856递增15.203.5A, CC1.000不变000自身结合松弛量提出定量建议“针对效率较低的B单元建议其一参照A单元的做法将‘投入1’如行政人员减少15个单位其二借鉴C单元的经验将‘产出1’如门诊量提升3.5个单位。同时其规模收益处于递增阶段说明扩大规模有助于提升效率可考虑适度增加资源投入。”4.3 模型的稳健性检验与扩展分析高级的论文不会只用一个模型就下结论需要进行稳健性检验。超效率DEA模型当有效单元过多时可以用超效率模型对有效单元进行进一步排序。其原理是在评价某个有效单元时将其自身从参考集中排除这样其效率值就可能大于1从而区分出“更有效”的单元。窗口DEA或Malmquist指数如果你的数据是面板数据多年份、多截面可以用这些方法分析效率的动态变化。Malmquist指数能分解出效率变化是源于技术进步还是技术效率改善非常适合用于评价政策实施效果或发展趋势。敏感性分析通过增加或删除一个指标观察效率排名是否发生剧烈变化。如果变化很大说明指标体系不够稳健结论需要谨慎对待。% 一个简单的敏感性分析示例移除一个投入指标看效率值变化 original_eff Efficiency_BCC; sensitivity_result zeros(num_dmu, m); % 存储每次移除一个投入后的平均效率变化 for remove_idx 1:m X_reduced X; X_reduced(remove_idx, :) []; % 移除第remove_idx个投入 m_reduced m - 1; eff_reduced zeros(num_dmu, 1); % ... 使用X_reduced和Y重新计算BCC效率 (需要重写函数或调整此处简化) % 假设我们调用了一个支持可变输入维度的dea函数 % eff_reduced my_DEA_function(X_reduced, Y, bcc); % sensitivity_result(:, remove_idx) abs(eff_reduced - original_eff); end % 分析sensitivity_result如果某列对应某个被移除的指标的变化均值很大说明该指标影响显著。5. 常见问题排查与MATLAB编程避坑指南在实际编程和调试过程中你肯定会遇到各种问题。这里我总结几个最常见的“坑”及其解决方法。5.1 问题linprog求解失败或无解可能原因与解决方案数据未标准化量纲差异导致数值问题使优化问题病态。务必先标准化数据。数据存在零或负值DEA要求投入产出数据一般为正。如果有零可以加一个极小的正数如1e-10避免除零错误如果有负值如利润亏损需要先进行数据平移等正向化处理。线性规划问题无可行解检查约束条件是否矛盾。在BCC模型中凸性约束sum(lambda)1与权重非负结合保证了始终有可行解。如果仍报错检查输入矩阵X和Y的维度是否正确。linprog选项设置尝试调整求解器选项提高求解精度。options optimoptions(linprog, Display, iter, ConstraintTolerance, 1e-9, OptimalityTolerance, 1e-9);5.2 问题所有DMU的效率值都是1可能原因指标过多DMU过少违反了“DMU数量 ≥ 2*(投入数产出数)”的经验法则。解决方案是使用指标筛选方法如相关系数法、主成分分析法降维或尝试收集更多DMU数据。指标间存在线性关系某个投入指标可能是其他几个投入指标的线性组合。检查投入/产出矩阵的秩是否亏缺。数据变异太小所有DMU的投入产出比例高度相似。检查数据这可能是数据本身的特点也可能是数据处理错误。5.3 问题松弛变量数值极小但不为零如何判断有效性这是数值计算中的常见问题。由于浮点数精度理论上应为0的松弛可能得到如1e-12这样极小的值。解决方案设定一个合理的容差阈值。tol 1e-6; % 根据数据规模设定通常1e-6到1e-8 is_effective (abs(theta - 1) tol) all(abs(slack_x) tol) all(abs(slack_y) tol);在判断有效性和筛选参考集lambda 1e-5时都要使用这样的阈值。5.4 问题如何提高大量DMU的计算速度循环调用linprog求解n个线性规划当n很大时如上千个可能会较慢。优化技巧向量化/并行计算如果MATLAB版本支持可以使用parfor循环替代for循环利用多核并行计算。if isempty(gcp(nocreate)) parpool; % 开启并行池 end parfor k 1:num_dmu % 调用dea计算函数 end预分配内存如我们之前所做在循环前用zeros预分配所有结果矩阵避免在循环中动态增长数组这会极大提升速度。使用更专业的优化求解器对于超大规模问题可以考虑CVX、YALMIP等建模工具搭配Gurobi、MOSEK等商业求解器但竞赛中MATLAB自带linprog通常足够。5.5 问题结果与文献或软件如DEAP不一致排查步骤检查模型导向确认一致。是投入导向还是产出导向同样的数据两种导向结果不同。检查数据标准化方法确保完全一致。别人用均值标准化你用极差标准化结果必然不同。检查松弛变量的处理有些软件或文献在计算效率时使用两阶段法第一阶段先求最优θ第二阶段在θ固定的情况下求最大松弛。而我们编写的模型是同时求解的。虽然理论上最终前沿一致但路径可能不同对于弱有效的单元判断可能有细微差别。验证核心算法用一个非常小的、手工可验证的示例数据如3个DMU2个投入1个产出测试你的代码确保逻辑正确。最后在论文中一定要写明你使用的软件是MATLAB并简要说明实现原理如“基于线性规划方法”将核心代码以整洁的格式放入附录。清晰的逻辑、完整的分析、可靠的实现和深度的解读才是让DEA模型在数学建模竞赛中为你赢得高分的核心。