DEA模型MATLAB实现（CCR、BCC、超效率）

DEA模型MATLAB实现（CCR、BCC、超效率）

一、模型原理与代码框架

1. CCR模型（规模报酬不变）

MATLAB代码：

function theta = DEA_CCR(X, Y)[m, n] = size(X);s = size(Y,1);theta = zeros(n,1);epsilon = 1e-10; % 非阿基米德无穷小for k = 1:n% 构建参考集（排除当前DMU）X_k = X(:,1:k-1);X_k = [X_k, X(:,k+1:n)];Y_k = Y(:,1:k-1);Y_k = [Y_k, Y(:,k+1:n)];% 目标函数系数f = [zeros(1,m), -epsilon*ones(1,s), 1];% 约束矩阵Aeq = [X_k', Y_k'];beq = Y(:,k);lb = zeros(m+s+1,1);% 求解线性规划options = optimoptions('linprog','Display','off');[sol, ~, exitflag] = linprog(f, [], [], Aeq, beq, lb, [], options);if exitflag == 1theta(k) = 1 / sol(end);elsetheta(k) = NaN; % 无解情况endend
end

2. BCC模型（规模报酬可变）
​​改进点​​：增加凸性约束 $\sum \lambda j=1$
​​MATLAB代码​​：

function theta = DEA_BCC(X, Y)[m, n] = size(X);s = size(Y,1);theta = zeros(n,1);epsilon = 1e-10;for k = 1:nX_k = X(:,1:k-1);X_k = [X_k, X(:,k+1:n)];Y_k = Y(:,1:k-1);Y_k = [Y_k, Y(:,k+1:n)];f = [zeros(1,m), -epsilon*ones(1,s), 1, 0]; % 新增凸性约束变量Aeq = [X_k', Y_k', ones(1,m), -Y(:,k)];beq = zeros(m+1,1);lb = zeros(m+s+2,1);options = optimoptions('linprog','Display','off');[sol, ~, exitflag] = linprog(f, [], [], Aeq, beq, lb, [], options);if exitflag == 1theta(k) = 1 / sol(end-1);elsetheta(k) = NaN;endend
end

3. 超效率DEA模型
​​核心改进​​：排除当前DMU作为参考集，允许效率值>1
​​MATLAB代码​​：

function super_theta = Super_Efficiency_DEA(X, Y)[m, n] = size(X);s = size(Y,1);super_theta = zeros(n,1);for k = 1:n% 构建参考集（排除第k个DMU）X_ref = [X(:,1:k-1), X(:,k+1:n)];Y_ref = [Y(:,1:k-1), Y(:,k+1:n)];% 目标函数（最大化θ）f = [-1, zeros(1,m+s)];% 输入约束：∑λX ≤ θX_kA1 = [X_ref', -X(:,k)];b1 = zeros(m,1);% 输出约束：∑λY ≥ Y_kA2 = [-Y_ref', Y(:,k)];b2 = zeros(s,1);A = [A1; A2];b = [b1; b2];lb = zeros(m+s,1);options = optimoptions('linprog','Display','off');[sol, ~, exitflag] = linprog(f, A, b, [], [], lb, [], options);if exitflag == 1super_theta(k) = 1 / sol(1);elsesuper_theta(k) = NaN;endend
end

二、数据预处理与调用示例

1. 数据准备

% 示例数据（3个DMU，2输入，1输出）
X = [100 80; 80 60; 90 60]; % 输入矩阵
Y = [30; 40; 50];        % 输出矩阵% 数据标准化（可选）
X = zscore(X);
Y = zscore(Y);

2. 模型调用

% CCR模型
CCR_eff = DEA_CCR(X, Y);
disp('CCR效率值:');
disp(CCR_eff);% BCC模型
BCC_eff = DEA_BCC(X, Y);
disp('BCC效率值:');
disp(BCC_eff);% 超效率模型
Super_eff = Super_Efficiency_DEA(X, Y);
disp('超效率值:');
disp(Super_eff);

三、结果解析与可视化

1. 效率值分析

2. 可视化代码

figure;
subplot(1,3,1);
bar(CCR_eff);
title('CCR效率');
xlabel('DMU'); ylabel('效率值');subplot(1,3,2);
bar(BCC_eff);
title('BCC效率');subplot(1,3,3);
bar(Super_eff);
title('超效率');

四、关键改进与优化

1. 非期望产出处理
对于污染排放等非期望指标，采用SBM模型修正：

function [theta, s_minus, s_plus] = DEA_SBM(X, Y)[m, n] = size(X);s = size(Y,1);theta = zeros(n,1);for k = 1:nX_k = X(:,1:k-1);X_k = [X_k, X(:,k+1:n)];Y_k = Y(:,1:k-1);Y_k = [Y_k, Y(:,k+1:n)];f = [zeros(1,m+s), 1];A = [X_k', Y_k', -Y(:,k)];b = zeros(m+1,1);lb = zeros(m+s+1,1);options = optimoptions('linprog','Display','off');[sol, ~, exitflag] = linprog(f, A, b, [], [], lb, [], options);if exitflag == 1theta(k) = sol(end);s_minus(k) = sol(m+1:end-1);s_plus(k) = sol(end);elsetheta(k) = NaN;endend
end

2. 并行计算加速
利用MATLAB并行工具箱加速大规模计算：

parfor k = 1:n% 并行执行各DMU计算X_ref = [X(:,1:k-1), X(:,k+1:n)];Y_ref = [Y(:,1:k-1), Y(:,k+1:n)];% ...（后续计算同上）
end

参考文献
Charnes A, Cooper W W, Rhodes E. Measuring the efficiency of decision making units[J]. European journal of operational research, 1978.

参考代码 求解DEA　CCR BCC 超效率 matlab程序 youwenfan.com/contentcsc/84774.html

超效率DEA模型MATLAB实现代码库. CSDN文库, 202 wenku.csdn.net/answer/7t1s9r0cz7.
范巧. 数据包络分析-Matlab实现. 2018.
精通DEA分析：MATLAB实现CCR、BCC模型及超效率计算. CSDN文库, 2025.