数据包络分析DEA

一些好文：https://zhuanlan.zhihu.com/p/390594011

多属性决策

采取加法或者乘法等模式，计算方案之加总价值，假设方案之价值衡量在每个属性上都是越大越好。

好文：https://blog.csdn.net/qq_44431690/article/details/108198631

进行信息集结的方法：加权算术平均（WAA）算法、加权几何平均（WGA）算子、有序加权平均（OWA）算法。

由于不同属性指标的值大小有所不同，为了消除不同物理量纲对决策结果的影响，决策时可对数据进行规范化处理。

多属性决策模型处理步骤：
1.对属性值进行归一化处理

用函数mapminmax
默认的map范围是[-1, 1]，所以如果需要[0, 1]，则按这样的格式提供参数:
MappedData = mapminmax(OriginalData, 0, 1);
只按行归一化，如果是矩阵，则每行各自归一化，如果需要对整个矩阵归一化，用如下方法：
FlattenedData = OriginalData(:)'; % 展开矩阵为一列，然后转置为一行。
MappedFlattened = mapminmax(FlattenedData, 0, 1); % 归一化。
MappedData = reshape(MappedFlattened, size(OriginalData)); % 还原为原始矩阵形式。此处不需转置回去，因为reshape恰好是按列重新排序

2.计算属性权重
2.由属性权重和属性值计算总得分

disp('请输入判断矩阵A(n阶)');
A=input('A=');
[n,n]=size(A);
x=ones(n,100);
y=ones(n,100);
m=zeros(1,100);
m(1)=max(x(:,1));
y(:,1)=x(:,1);
x(:,2)=A*y(:,1);
m(2)=max(x(:,2));
y(:,2)=x(:,2)/m(2);
p=0.0001;i=2;k=abs(m(2)-m(1));
while  k>p
  i=i+1;
  x(:,i)=A*y(:,i-1);
  m(i)=max(x(:,i));
  y(:,i)=x(:,i)/m(i);
  k=abs(m(i)-m(i-1));
end
a=sum(y(:,i));
w=y(:,i)/a;
t=m(i);
disp(w); % w表示了底层对上层的权重
         %以下是一致性检验
CI=(t-n)/(n-1);RI=[0 0 0.52 0.89 1.12 1.26 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59];
CR=CI/RI(n);
if CR<0.10
    disp('此矩阵的一致性可以接受!');
    disp('CI=');disp(CI);
    disp('CR=');disp(CR);
end

数据包络分析

以【效率】为概念作为加总模式，而效率则等于加总产出÷总投入，并以效率最大化为目标。
同时DEA中的权重不是事先预设的，而是在实证资料中推导产生。

好文：https://wenku.baidu.com/link?url=5QTVRpdlhN82xkuKnmBox0csPCggStzVo7SEqakP-ipLwu5757VouBDyQGUl0skclxgGctpirFukXA4pi_bWbwh5KTVmMwFbKw1GSk2SspK

CCR模型

https://zhuanlan.zhihu.com/p/390594011

function [eff,slackX,slackY,Xeff,Yeff ] = dea_ccr( X, Y)
%DEA Data envelopment analysis radial model
%   X表示投入变量矩阵
%   Y表示产出变量矩阵
    
    %1. 提取变量的维度
    [n, m] = size(X);
    s = size(Y,2);
    
    %2. 依据提取的维度信息，创建求解变量的存储矩阵
   
    slackX = nan(n, m);%投入的松弛变量
    slackY = nan(n, s);%产出的松弛变量
    eff = nan(n, 1);%效率值，第一阶段可计算得的
    Eflag = nan(n, 2);%求解的过程信息
    lambda = nan(n, n);% 每个单元的intensity variable
    Xeff = nan(n, m);%投入对应的前沿值
    Yeff = nan(n, s);%产出对应的前沿值
    
    %3. 参数准备及求解线性规划过程
   
        % 对每个 DMU
        for j=1:n
            
            %第一阶段，求解效率θ
                % 目标函数
                f = [zeros(1,n), 1];
                
                % 约束条件
                A = [ X', -X(j,:)';
                     -Y', zeros(s,1)];
                b = [zeros(m,1);
                    -Y(j,:)'];
                Aeq = [];
                beq = [];
                lb = zeros(1, n + 1);
                ub = [];                
                % linprog（）求解
                [z, ~, exitflag] = linprog(f, A, b, Aeq, beq, lb, ub);
                if exitflag ~= 1
                    fprintf(1,'DMU %6d 第一阶段的exitflag: %6d', j, exitflag)
                end
                if isempty(z)
                    fprintf(1,'DMU %6d 第一阶段的解不存在.', j)
                    z = nan(n + 1, 1);
                end
                
                % 提取计算结果
                theta = z(end);
                Eflag(j, 1) = exitflag;
                eff(j) = theta;
    
                
           %第二阶段，计算松弛变量及每个单元各变量对应的前沿值
                
                if ~isnan(theta)
                
                    % 目标函数
                    f = [zeros(1, n), -ones(1, m + s)];

                    % 约束条件
                    Aeq = [ X', eye(m,m)  , zeros(m,s);
                            Y', zeros(s,m), -eye(s,s)];
                    beq = [theta .* Xeval(j,:)';
                            Yeval(j,:)'];
                    lb = zeros(n + s + m, 1);

                    % 求解
                    [z, ~, exitflag] = linprog(f, [], [], Aeq, beq, lb, []);
                    if exitflag ~= 1
                        fprintf(1, 'DMU %6d 第二阶段的exitflag: %6d: %i', j, exitflag)
                    end
                    if isempty(z)
                        fprintf(1,'DMU %6d 第二阶段的解不存在.', j)
                        z = nan(n + m + s, 1);                   
                    end

                    % 提取结果
                    lambda(j,:) = z(1:n);
                    slackX(j,:) = z(n + 1 : n + m);
                    slackY(j,:) = z(n + m + 1 : n + m + s);                
                    Eflag(j, 2) = exitflag;

                    % 计算每个变量的前沿值
                    Xeff(j,:) = repmat(eff(j), 1, m) .* X(j,:) - slackX(j,:);
                    Yeff(j,:) = Y(j,:) + slackY(j,:);

                end
                
        end
end