数学建模-评价类模型学习笔记（自存）

评价类问题

各评价类模型的比较： 数学建模中各种评价类模型的优点和缺点总结

一、层次分析法

参考文章：用人话讲明白AHP层次分析法

1、适用场景

• 指标定权 给指标制定权重，打个比方，例如选择旅游地这个决策，可能一般我们由以下5个因素组成，但是每个人（主观）对因素的重视程度不一，ahp可以实现在无需搜集数据的情况下，给这些指标制定权重。

• 量化方案选择 同样是选择旅游地这个决策，可能我们有一些方案，例如苏杭、北戴河、桂林这三个方案，层次分析法可以综合以上5个因素，给这些方案计算得出一个量化得分，例如苏杭0.3分、北戴河0.35分、桂林0.45分，这样根据分值大小，我们就可以选择得到内心或者经验上最心怡的方案了。

2、分析步骤

2.1 构建层次评价模型

确认整个决策事件的目标层、准则层、方案层。

2.2 构造判断矩阵

设现在要比较n个因子$X=x_1,\cdots ,x_n$对某因素$Z$的影响大小,对因子进行两两比较建立成对比较矩阵,每次取两个因子$x_i$和$x_j$，以$a_{ij}$表示 $x_i$和$x_j$对$Z$的影响大小之比，全部比较结果用矩阵$A={(a_{ij})}_{m\times n}$ 表示，称 A 为 $Z-X$之间的成对比较判断矩阵（简称判断矩阵）。判断矩阵中的元素需满足：
$$\begin{gathered} 1)a_{ij}>0 \\ 2)a_{ij}=\frac{1}{a_{ji}} \\ 3)a_{ii}=1 \end{gathered}$$
判断矩阵构造方法：

举个例子：

	A因素	B因素	C因素
A因素	1	1/2
B因素	2	1
C因素			1

若B因素比A因素重要，则以B因素为行首中对应A因素的位置元素为一个大于1的数，相应的对称位置为其倒数

2.3 层次单排序与一致性检验

各个因素之间的重要程度排序，在合理情况下应有如下关系：
$$a_{ij}=\frac{i的重要程度}{j的重要程度}，a_{jk}=\frac{j的重要程度}{k的重要程度}⟹a_{ik}=\frac{i的重要程度}{k的重要程度}=a_{ij}\times a_{jk}$$
若判断矩阵A满足这个关系，则其各行应该成比例，称这种矩阵为一致矩阵，其一个特征值$tr(A)=n$,其余特征值为0。

但由于实际构造判断矩阵元素时含有巨大主观成分，故判断矩阵的一致性（即与一致矩阵之间的误差）需要经过检验，检验步骤如下：

1.计算一致性指标$CI$
$$CI=\frac{{\lambda }_{max}-n}{n-1}$$
其中${\lambda }_{max}$为判断矩阵最大特征值，n为元素的个数。

2.查找对应的随机一致性指标 RI

$n$	1	2	3	4	5	6	7	8	9	10	11
RI	0	0	0.58	0.90	1.12	1.24	1.32	1.41	1.45	1.49	1.51

3.计算一致性比例
$$CR=\frac{CI}{RI}$$
如果$CR<0.1$, 则可认为判断矩阵的一致性可以接受,否则需要对判断矩阵进行修正。

完成一致性检验后，若矩阵一致性满足要求，则可以进行权重计算。

设判断矩阵为$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right]$,三种权重计算方法得到的权重向量如下：

• 算术平均法：${\mathbf{\omega }}_{\mathbf{i}}=\frac{1}{\mathbf{n}}\sum _{j=1}^n\frac{a_{ij}}{\sum _{k=1}^n{a}_{kj}}$

• 几何平均法：${\mathbf{\omega }}_{\mathbf{i}}=\frac{{\left(\prod _{j=1}^n{a}_{ij}\right)}^{\frac{1}{\mathbf{n}}}}{\sum _{k=1}^n{\left(\prod _{j=1}^n{a}_{kj}\right)}^{\frac{1}{n}}}$

• 特征值法：对矩阵$A$的特征向量进行归一化得到权重向量。

2.4 层次总排序

从最高层次到最低层次，依次计算某一层次所有因素对于最高层（总目标）相对重要性的权值。

3、代码

disp('请输入判断矩阵A')
A=input('A=');
[n,n] = size(A);
% % % % % % % % % % % % %方法1： 算术平均法求权重% % % % % % % % % % % % %
Sum_A = sum(A);
SUM_A = repmat(Sum_A,n,1);
Stand_A = A ./ SUM_A;
disp('算术平均法求权重的结果为：');
disp(sum(Stand_A,2)./n)
% % % % % % % % % % % % %方法2： 几何平均法求权重% % % % % % % % % % % % %
Prduct_A = prod(A,2);
Prduct_n_A = Prduct_A .^ (1/n);
disp('几何平均法求权重的结果为：');
disp(Prduct_n_A ./ sum(Prduct_n_A))
% % % % % % % % % % % % %方法3： 特征值法求权重% % % % % % % % % % % % %
[V,D] = eig(A);
Max_eig = max(max(D));
[r,c]=find(D == Max_eig , 1);
disp('特征值法求权重的结果为：');
disp( V(:,c) ./ sum(V(:,c)) )
% % % % % % % % % % % % %下面是计算一致性比例CR的环节% % % % % % % % % % % % %
CI = (Max_eig - n) / (n-1);
RI=[0 0.0001 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];  %注意哦，这里的RI最多支持 n = 15
% 这里n=2时，一定是一致矩阵，所以CI = 0，我们为了避免分母为0，将这里的第二个元素改为了很接近0的正数
CR=CI/RI(n);
disp('一致性指标CI=');disp(CI);
disp('一致性比例CR=');disp(CR);
if CR<0.10
    disp('因为CR<0.10，所以该判断矩阵A的一致性可以接受!');
else
    disp('注意：CR >= 0.10，因此该判断矩阵A需要进行修改!');
end

二、优劣解距离法Topsis

参考文章：TOPSIS(逼近理想解)算法原理详解与代码实现 - 知乎 (zhihu.com)

1、适用场景

TOPSIS法（Technique for Order Preference by Similarity to Ideal Solution）可翻译为逼近理想解排序法，国内常简称为优劣解距离法
TOPSIS 法是一种常用的综合评价方法，其能充分利用原始数据的信息，其结果能精确地反映各评价方案之间的差距。

为了对众多方案给出一个排序，在给出所有方案之后，可以根据这些数据，构造出一个所有方案组成的系统中的理想最优解和最劣解。而TOPSIS的想法就是，通过一定的计算，评估方案系统中任何一个方案距离理想最优解和最劣解的综合距离。如果一个方案距离理想最优解越近，距离最劣解越远，我们就有理由认为这个方案更好。那理想最优解和最劣解又是什么呢？很简单，理想最优解就是该理想最优方案的各指标值都取到系统中评价指标的最优值，最劣解就是该理想最劣方案的各指标值都取到系统中评价指标的最劣值。

理想最优解中的数据都是各方案中的数据，而不要选择方案中没有的数据，理想最劣解同理。

一方面，层次分析法具有一定的局限性，决策层数量即n值不能太大，否则很难满足一致性要求。另一方面，某些分析场景已知决策层中指标的数据，考虑如何利用这些数据来使得评价更加准确，这时候就要用到TOPSIS法。

2、算法实现

2.1 数据预处理

（1）正向化处理

首先进行数据正向化处理。指标类型分为4类，如下图所示：

正向化处理，就是将上述的四种指标数据进行处理，将其全部转化为极大型指标数据，这样我们计算时问题就少一点，码代码时也更加统一。正向化处理得到的矩阵设为X，其中有n个要评价的对象，m个评价指标，则
$$X=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1m}\\ x_{21}& x_{22}& \cdots & x_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ x_{n1}& x_{n2}& \cdots & x_{nm}\end{array}\right]$$

极小型指标：利用 $\widehat{x_i}=max-x_i$ 将其转化为极大型，如果所有元素都为正数，也可以使用 $\widehat{x_i}=\frac{1}{x_i}$

中间型指标：如果其最佳数值为$x_{best}$, 取 M = m a x { ∣ x i − x b e s t ∣ } M=max\{|x_i-x_{best}|\} M=max{∣xi​−xbest​∣},之后按照$\widehat{x_i}=1-\frac{x_i-x_{best}}{M}$ 转化。

区间型指标：如果最佳区间为$[a,b]$,取 M = m a x { a − m i n { x i } , m a x { x i } − b } M=max\{a-min\{x_i\},max\{x_i\}-b \} M=max{a−min{xi​},max{xi​}−b},按照$\widehat{x_i}=\{\begin{array}{ll}1-\frac{a-x_i}{M},& x_i<a\\ 1,& a<x_i<b\\ 1-\frac{x_i-b}{M},& x_i>b\end{array}$ 转化

（2）标准化处理

为了消除不同的数据指标量纲的影响，我们还有必要对已经正向化的矩阵进行标准化。记标准化后的矩阵为Z，其中$z_{ij}=\frac{x_{ij}}{\sqrt{{\sum }_{i=1}^n{x}_{ij}^2}}$，也就是$\frac{每一个元素}{\sqrt{其所在列元素的平方和}}$。对数据进行了相应的处理后，可以用向量$z_i$来表达第$i$个方案。假设有n个待评价的方案，m个指标，此时 $z_i=[z_{i1},z_{i2},...,z_{im}]$ 。由这n个向量构成的矩阵也就是我们的标准化矩阵$Z$了。
$$Z=\left[\begin{array}{cccc}{z}_{11}& z_{12}& \cdots & z_{1m}\\ z_{21}& z_{22}& \cdots & z_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ z_{n1}& z_{n2}& \cdots & z_{mn}\end{array}\right]$$
此外还需判断$Z$矩阵中是否含有负数，若存在负数则需要对X矩阵进行另一种标准化方法，标准化公式为：
z ~ i j = x i j − min ⁡ { x 1 j , x 2 j , ⋯   , x n j } max ⁡ { x 1 j , x 2 j , ⋯   , x n j } − min ⁡ { x 1 j , x 2 j , ⋯   , x n j } \tilde{z}_{i j}=\frac{x_{i j}-\min \left\{x_{1 j}, x_{2 j}, \cdots, x_{n j}\right\}}{\max \left\{x_{1 j}, x_{2 j}, \cdots, x_{n j}\right\}-\min \left\{x_{1 j}, x_{2 j}, \cdots, x_{n j}\right\}} z~ij​=max{x1j​,x2j​,⋯,xnj​}−min{x1j​,x2j​,⋯,xnj​}xij​−min{x1j​,x2j​,⋯,xnj​}​

2.2 最优解与最劣解计算

经过了正向化处理和标准化处理的评分矩阵Z，里面的数据全部是极大型数据。我们就可以从中取出理想最优解和最劣解。因此我们取出每个指标，即每一列中最大的数，构成理想最优解向量，即
z + = [ z 1 + , z 2 + , … , z m + ] = [ max ⁡ { z 11 , z 21 , … , z n 1 } , max ⁡ { z 12 , z 22 , … , z n 2 } , … , max ⁡ { z 1 m , z 2 m , … , z n n } ] z^{+}=\left[z_1^{+}, z_2^{+}, \ldots, z_m^{+}\right]=\left[\max \left\{z_{11}, z_{21}, \ldots, z_{n 1}\right\}, \max \left\{z_{12}, z_{22}, \ldots, z_{n 2}\right\}, \ldots, \max \left\{z_{1 m}, z_{2 m}, \ldots, z_{n n}\right\}\right] z+=[z1+​,z2+​,…,zm+​]=[max{z11​,z21​,…,zn1​},max{z12​,z22​,…,zn2​},…,max{z1m​,z2m​,…,znn​}]
同理，取每一列中最小的数计算理想最劣解向量：
z − = [ z 1 − , z 2 − , … , z m − ] = [ min ⁡ { z 11 , z 21 , … , z n 1 } , min ⁡ { z 12 , z 22 , … , z n 2 } , … , min ⁡ { z 1 m , z 2 m , … , z r m } ] z^{-}=\left[z_1^{-}, z_2^{-}, \ldots, z_m^{-}\right]=\left[\min \left\{z_{11}, z_{21}, \ldots, z_{n 1}\right\}, \min \left\{z_{12}, z_{22}, \ldots, z_{n 2}\right\}, \ldots, \min \left\{z_{1 m}, z_{2 m}, \ldots, z_{r m}\right\}\right] z−=[z1−​,z2−​,…,zm−​]=[min{z11​,z21​,…,zn1​},min{z12​,z22​,…,zn2​},…,min{z1m​,z2m​,…,zrm​}]
$z^{+}\text{就是 }{z}_{\max }，z^{-}\text{就是 }{z}_{\min }\text{ 。 }$

在得到理想最优解和理想最劣解的基础上就能计算每个方案的评分了。根据距离评分公式：$\frac{z_i-z_{\min }}{z_{\max }-z_{\min }}$,对其进行变型，也就是：$\frac{z_i-z_{\min }}{\left(z_{\max }-z_i\right)+\left(z_i-z_{\text{min }}\right)}$ 。变型的目的是为了使用欧几里得距离来衡量两个方案的距离，变形前后分母的计算结果其实是不同的(因为这里$z_i$是一个向量)。这样更能体现出是综合距离。否则所有方案计算得分时分母都是相同的，相当于只衡量了分子，也就是距离最劣解的距离。

2.3 TOPSIS评分计算

对于第i个方案 $z_i$，我们计算它与最优解的距离：$d_i^{+}=\sqrt{\sum _{j=1}^m{\left(z_j^{+}-z_{ij}\right)}^2}$

与最劣解的距离：$d_i^{-}=\sqrt{\sum _{j=1}^m{\left(z_j^{-}-z_{ij}\right)}^2}$

定义第i个方案的评分为$S_i$,根据前面的评分公式有：$S_i=\frac{d_i^{-}}{d_i^{+}+d_i^{-}}$

这个时候我们就有了每个方案的分数了，按分数排排序，就知道哪个方案比较好哪个方案比较差。

3.熵权法的加入

熵权法是一种客观赋权方法，其基本原理在于：指标的变异程度越小，所反映的信息量也越少，其对应的权值也应该越低。

3.1 计算概率

计算第j项指标下第i个样本所占的比重，并将其看作相对熵计算中用到的概率。可以理解为指标归一化，即将标准化矩阵$Z$将每一个数除以这一列数之和。

在标准化处理得到的非负矩阵 $\tilde{Z}=\left[\begin{array}{cccc}{\tilde{z}}_{11}& {\tilde{z}}_{12}& \cdots & {\tilde{z}}_{1m}\\ {\tilde{z}}_{21}& {\tilde{z}}_{22}& \cdots & {\tilde{z}}_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ {\bar{z}}_{n1}& {\tilde{z}}_{n2}& \cdots & {\tilde{z}}_{nm}\end{array}\right]$ 中，计算概率矩阵 $P$ 的公式为：
$$p_{ij}=\frac{{\tilde{z}}_{ij}}{\sum _{i-1}^n{\tilde{z}}_{ij}}$$
容易验证：$\sum _{i=1}^n{p}_{ij}=1$，即保证了每个指标对应的概率和为1。

3.2 计算信息熵

计算每个指标的信息熵，并计算信息效用值，并归一化得到每个指标的熵权。对于第j个指标而言，其信息熵计算公式为：
$$e_j=-\frac{1}{\ln n}\sum _{i=1}^n{p}_{ij}\ln (p_{ij})(j=1,2,\cdots ,m)$$
计算信息效用值的公式为：
$$d_j=1-e_j$$
信息效用值越大，该指标对应的信息越多。将信息效用值归一化，就能得到每个指标的熵权：
$$W_j=\frac{d_j}{\sum _{j=1}^m{d}_j}(j=1,2,\cdots ,m)$$

4.代码

4.1 指标正向化部分

Positivization.m:

% function [输出变量] = 函数名称(输入变量）  
% 函数的中间部分都是函数体
% 函数的最后要用end结尾
% 输出变量和输入变量可以有多个，用逗号隔开
% function [a,b,c]=test(d,e,f)
%     a=d+e;
%     b=e+f;
%     c=f+d;
% end
% 自定义的函数要单独放在一个m文件中，不可以直接放在主函数里面（和其他大多数语言不同）

function [posit_x] = Positivization(x,type,i)
% 输入变量有三个：
% x：需要正向化处理的指标对应的原始列向量
% type： 指标的类型（1：极小型， 2：中间型， 3：区间型）
% i: 正在处理的是原始矩阵中的哪一列
% 输出变量posit_x表示：正向化后的列向量
    if type == 1  %极小型
        disp(['第' num2str(i) '列是极小型，正在正向化'] )
        posit_x = Min2Max(x);  %调用Min2Max函数来正向化
        disp(['第' num2str(i) '列极小型正向化处理完成'] )
        disp('~~~~~~~~~~~~~~~~~~~~分界线~~~~~~~~~~~~~~~~~~~~')
    elseif type == 2  %中间型
        disp(['第' num2str(i) '列是中间型'] )
        best = input('请输入最佳的那一个值： ');
        posit_x = Mid2Max(x,best);
        disp(['第' num2str(i) '列中间型正向化处理完成'] )
        disp('~~~~~~~~~~~~~~~~~~~~分界线~~~~~~~~~~~~~~~~~~~~')
    elseif type == 3  %区间型
        disp(['第' num2str(i) '列是区间型'] )
        a = input('请输入区间的下界： ');
        b = input('请输入区间的上界： '); 
        posit_x = Inter2Max(x,a,b);
        disp(['第' num2str(i) '列区间型正向化处理完成'] )
        disp('~~~~~~~~~~~~~~~~~~~~分界线~~~~~~~~~~~~~~~~~~~~')
    else
        disp('没有这种类型的指标，请检查Type向量中是否有除了1、2、3之外的其他值')
    end
end

Min2Max.m(用于极小型指标正向化)

function [posit_x] = Min2Max(x)
  posit_x = max(x) - x;
   %posit_x = 1 ./ x;   %如果x全部都大于0，也可以这样正向化
end

Mid2Max.m(用于中间型指标正向化)

function [posit_x] = Mid2Max(x,best)
    M = max(abs(x-best));
    posit_x = 1 - abs(x-best) / M;
end

Inter2Max.m(用于区间型指标正向化)

function [posit_x] = Inter2Max(x,a,b)
    r_x = size(x,1);  % row of x 
    M = max([a-min(x),max(x)-b]);
    posit_x = zeros(r_x,1);   %zeros函数用法: zeros(3)  zeros(3,1)  ones(3)
    % 初始化posit_x全为0  初始化的目的是节省处理时间
    for i = 1: r_x
        if x(i) < a
           posit_x(i) = 1-(a-x(i))/M;
        elseif x(i) > b
           posit_x(i) = 1-(x(i)-b)/M;
        else
           posit_x(i) = 1;
        end
    end
end

4.2 熵权法部分

由于MATLAB中，ln（0）的结果返回不是0，所以需要重新定义一个log函数，说明一下这种特殊情况，使得ln（0）的返回结果为0

% 重新定义一个mylog函数，当输入的p中元素为0时，返回0
function [lnp] =  mylog(p)%lnp为输出变量，p为输入变量
n = length(p);   % 向量的长度
lnp = zeros(n,1);   % 初始化最后的结果
    for i = 1:n   % 开始循环
        if p(i) == 0   % 如果第i个元素为0
            lnp(i) = 0;  % 那么返回的第i个结果也为0
        else
            lnp(i) = log(p(i));  
        end
    end
end

熵权法的代码：

function [W] = Entropy_Method(Z)
% 计算有n个样本，m个指标的样本所对应的的熵权
% 输入
% Z ： n*m的矩阵（要经过正向化和标准化处理，且元素中不存在负数）
% 输出
% W：熵权，m*1的行向量

%% 计算熵权
    [n,m] = size(Z);
    D = zeros(1,m);  % 初始化保存信息效用值的行向量
    for i = 1:m
        x = Z(:,i);  % 取出第i列的指标
        p = x / sum(x);
        % 注意，p有可能为0，此时计算ln(p)*p时，Matlab会返回NaN，所以这里我们自己定义一个函数
        e = -sum(p .* mylog(p)) / log(n); % 计算信息熵
        D(i) = 1- e; % 计算信息效用值
    end
    W = D ./ sum(D);  % 将信息效用值归一化，得到权重    
end

4.3 主函数部分

%%  第一步：把数据复制到工作区，并将这个矩阵命名为X
% （1）在工作区右键，点击新建（Ctrl+N)，输入变量名称为X
% （2）在Excel中复制数据，再回到Matlab中右键，点击粘贴Excel数据（Ctrl+Shift+V）
% （3）关掉这个窗口，点击X变量，右键另存为，保存为mat文件（下次就不用复制粘贴了，只需使用load命令即可加载数据）
% （4）注意，代码和数据要放在同一个目录下哦，且Matlab的当前文件夹也要是这个目录。
clear;clc
load X.mat 
%% 注意：如果提示: 错误使用 load，无法读取文件 'data_water_quality.mat'。没有此类文件或目录。
% 那么原因是因为你的Matlab的当前文件夹中不存在这个文件
% 可以使用cd函数修改Matlab的当前文件夹
% 比如说，我的代码和数据放在了: D:第2讲.TOPSIS法（优劣解距离法）\代码和例题数据
% 那么我就可以输入命令：
% cd 'D:第2讲.TOPSIS法（优劣解距离法）\代码和例题数据'


%%  第二步：判断是否需要正向化
[n,m] = size(X);
disp(['共有' num2str(n) '个评价对象, ' num2str(m) '个评价指标']) 
Judge = input(['这' num2str(m) '个指标是否需要经过正向化处理，需要请输入1 ，不需要输入0：  ']);

if Judge == 1
    Position = input('请输入需要正向化处理的指标所在的列，例如第2、3、6三列需要处理，那么你需要输入[2,3,6]： '); %[2,3,4]
    disp('请输入需要处理的这些列的指标类型（1：极小型， 2：中间型， 3：区间型） ')
    Type = input('例如：第2列是极小型，第3列是区间型，第6列是中间型，就输入[1,3,2]：  '); %[2,1,3]
    % 注意，Position和Type是两个同维度的行向量
    for i = 1 : size(Position,2)  %这里需要对这些列分别处理，因此我们需要知道一共要处理的次数，即循环的次数
        X(:,Position(i)) = Positivization(X(:,Position(i)),Type(i),Position(i));
    % Positivization是我们自己定义的函数，其作用是进行正向化，其一共接收三个参数
    % 第一个参数是要正向化处理的那一列向量 X(:,Position(i))   回顾上一讲的知识，X(:,n)表示取第n列的全部元素
    % 第二个参数是对应的这一列的指标类型（1：极小型， 2：中间型， 3：区间型）
    % 第三个参数是告诉函数我们正在处理的是原始矩阵中的哪一列
    % 该函数有一个返回值，它返回正向化之后的指标，我们可以将其直接赋值给我们原始要处理的那一列向量
    end
    disp('正向化后的矩阵 X =  ')
    disp(X)
end

%% 第三步：对正向化后的矩阵进行标准化
Z = X ./ repmat(sum(X.*X) .^ 0.5, n, 1);
disp('标准化矩阵 Z = ')
disp(Z)

%% 第四步：熵权法计算权重
weight = Entropy_Method(Z)

%% 第五步：计算与最大值的距离和最小值的距离，并算出得分
%不使用权重的得分计算
%D_P = sum([(Z - repmat(max(Z),n,1)) .^ 2 ],2) .^ 0.5;   % D+ 与最大值的距离向量
%D_N = sum([(Z - repmat(min(Z),n,1)) .^ 2 ],2) .^ 0.5;   % D- 与最小值的距离向量
%使用权重的得分计算
D_P = sum([(Z - repmat(max(Z),n,1)) .^ 2 ] .* repmat(weight,n,1) ,2) .^ 0.5;   % D+ 与最大值的距离向量
D_N = sum([(Z - repmat(min(Z),n,1)) .^ 2 ] .* repmat(weight,n,1) ,2) .^ 0.5;   % D- 与最小值的距离向量
S = D_N ./ (D_P+D_N);    % 未归一化的得分
disp('最后的得分为：')
stand_S = S / sum(S)
[sorted_S,index] = sort(stand_S ,'descend')

% A = magic(5)  % 幻方矩阵
% M = magic(n)返回由1到n^2的整数构成并且总行数和总列数相等的n×n矩阵。阶次n必须为大于或等于3的标量。
% sort(A)若A是向量不管是列还是行向量，默认都是对A进行升序排列。sort(A)是默认的升序，而sort(A,'descend')是降序排序。
% sort(A)若A是矩阵，默认对A的各列进行升序排列
% sort(A,dim)
% dim=1时等效sort(A)
% dim=2时表示对A中的各行元素升序排列
% A = [2,1,3,8]
% Matlab中给一维向量排序是使用sort函数：sort（A），排序是按升序进行的，其中A为待排序的向量；
% 若欲保留排列前的索引，则可用 [sA,index] = sort(A,'descend') ，排序后，sA是排序好的向量，index是向量sA中对A的索引。
% sA  =  8     3     2     1
% index =  4     3     1     2

三、灰色关联分析法

参考文章：灰色关联分析(GRA)的理论及应用（matlab和python）

​ 【数学建模】灰色关联分析+Matlab代码实现

1.适用场景

灰色关联分析的基本思想是根据序列曲线几何形状的相似程度来判断其联系是否紧密，曲线越接近，相应序列之间的关联度就越大，反之则越小。通常可以运用此方法来分析各个因素对于结果的影响程度，即用于 进行系统分析，也可以运用此方法解决随时间变化的综合评价类问题，其核心是按照一定规则确立随时间变化的母序列，把各个评估对象随时间的变化作为子序列，求各个子序列与母序列的相关程度，依照相关性大小得出结论。

一般的抽象系统，如社会系统、经济系统、农业系统、生态系统、教育系统等都包含有许多种因素,多种因素共同作用的结果决定了该系统的发展态势。人们常常希望知道在众多的因素中，哪些是主要因素，哪些是次要因素；哪些因素对系统发展影响大，哪些因素对系统发展影响小；哪些因素对系统发展起推动作用需强化发展，哪些因素对系统发展起阻碍作用需加以抑制；这些都是系统分析中人们普遍关心的问题。例如，粮食生产系统,人们希望提高粮食总产量，而影响粮食总产量的因素是多方面的，有播种面积以及水利、化肥、土壤、种子、劳力、气候、耕作技术和政策环境等。为了实现少投人多产出，并取得良好的经济效益、社会效益和生态效益，就必须进行系统分析。

2.分析步骤

下面以系统分析为例说明灰色关联分析法的具体步骤。

2.1 确定分析数列

确定反映系统行为特征的参考数列和影响系统行为的比较数列。反映系统行为特征的数据序列，称为参考数列。影响系统行为的因素组成的数据序列，称比较数列。

（1）参考数列（又称母序列）为：$ Y = Y ( k ) ∣ k = 1 , 2,…, n $；

（2）比较数列（又称子序列）为：$ X_i = X_i ( k ) ∣ k = 1 , 2,\cdots,n ;\ i = 1 , 2,\cdots,m$ 。

如上图所示，要分析三种产业对于国内总产值的影响大小，参考序列就是国内总产值，比较序列就是第一产业、第二产业和第三产业。可以理解为参考序列是因变量，比较序列是自变量。

2.2 变量的无量纲化

由于系统中各因素列中的数据可能因量纲不同，不便于比较或在比较时难以得到正确的结论。因此在进行灰色关联度分析时，一般都要进行数据的无量纲化处理。主要有一下两种方法：

（1）初值化处理：$x_i(k)=\frac{x_i(k)}{x_i(1)},k=1,2\dots n;i=0,1,2\dots m$

（2）均值化处理：$x_i(k)=\frac{x_i(k)}{{\bar{x}}_i},k=1,2\dots n;i=0,1,2\dots m$

其中k对应时间段，i对应比较数列中的一行（即一个特征）。

2.3 计算关联系数

$${\xi }_i(k)=\frac{\underset{i}{\min }\underset{k}{\min }|y(k)-x_i(k)|+\mathbf{\rho }\underset{i}{\max }\underset{k}{\max }|y(k)-x_i(k)|}{|y(k)-x_i(k)|+\mathbf{\rho }\underset{i}{\max }\underset{k}{\max }|y(k)-x_i(k)|}$$

设$a$为两级最小差，$b$为两级最大差，即：
$$a=\underset{i}{\min }\underset{k}{\min }|y(k)-x_i(k)|,b=\underset{i}{\max }\underset{k}{\max }|y(k)-x_i(k)|$$
则有：
$${\xi }_i(k)=\frac{a+\mathbf{\rho }b}{|y(k)-x_i(k)|+\mathbf{\rho }b}$$
$\rho \in (0,\infty )$，称为分辨系数。$\rho$越小，分辨力越大，一般$\rho$的取值区间为$(0,1)$,具体取值可视情况而定。当$\rho \le 0.5463$时，分辨力最好，通常取$\mathbf{\rho }=0.5$。

2.4 计算灰色关联度

关联系数是比较数列与参考数列在各个时刻（即曲线中的各点）的关联程度值，所以它的数不止一个，而信息过于分散不便于进行整体性比较。因此有必要将各个时刻（即曲线中的各点）的关联系数集中为一个值，即求其平均值，作为比较数列与参考数列间关联程度的数量表示，灰色关联度 $r_i$ 公式如下：
$$r_i=\frac{1}{n}\sum _{k=1}^n{\xi }_i(k)(k=1,2,\cdots ,n)$$

2.5 关联度排序

将灰色关联度按大小排序，关联度越大，认为对应的比较序列与参考序列的关联度更高。

2.6 综合评价类问题分析步骤

综合评价类问题与系统分析类问题类似，不同之处在于缺少一个明显的参考序列，分析步骤如下：

1. 对指标进行正向化
2. 对正向化矩阵进行预处理
3. 将预处理后的具有真每一行取出最大值构成母序列
4. 计算各个指标与母序列的灰色关联度
5. 计算各个指标的权重
6. 计算每个评价对象的得分

3.代码

clear;clc
load gdp.mat  % 导入数据 以6*4为例
Mean = mean(gdp);  % 求出每一列的均值以供后续的数据预处理
gdp = gdp ./ repmat(Mean,size(gdp,1),1);  %size(gdp,1)=6, repmat(Mean,6,1)可以将矩阵进行复制，复制为和gdp同等大小，
disp('预处理后的矩阵为：'); disp(gdp)
Y = gdp(:,1);  % 母序列
X = gdp(:,2:end); % 子序列
absX0_Xi = abs(X - repmat(Y,1,size(X,2)))  % 计算|Y-Xi|矩阵
a = min(min(absX0_Xi))    % 计算两级最小差a
b = max(max(absX0_Xi))  % 计算两级最大差b
rho = 0.5; % 分辨系数取0.5
gamma = (a+rho*b) ./ (absX0_Xi  + rho*b)  % 计算子序列中各个指标与母序列的关联系数
disp('子序列中各个指标的灰色关联度分别为：')
disp(mean(gamma))

4.补充说明

• 1.选择用标准化回归或者灰色关联分析的具体时机：
  当样本个数 n 较大时，一般使用标准化回归。当样本个数较少时，选择使用灰色关联分析。
• 2.如果母序列中含有多个指标应如何分析？
  例如：Y1 和 Y2 是母序列，X1,X2,…,Xm 是子序列。那么应首先计算 Y1 和 X1,X2,…,Xm 的灰色关联度进行分析，然后再计算 Y2 和 X1,X2,…,Xm 的灰色关联度进行分析。
• 3.灰色关联分析与熵权法都可进行客观的赋权，主要区别在于熵权法是根据数据的变异程度进行赋权，而灰色关联分析是根据序列几何形状的相似程度来进行赋权。这里推荐用熵权法，如果要用灰色关联分析，那么就需要有相应的统计图，这样更具说服力。
• 4.灰色数学的相关概念是国人提出来的，美赛不建议用！

四、模糊综合评价

参考文章：数学建模笔记——评价类模型之模糊综合评价 - 知乎 （原理讲解）

​ 模糊综合评价_面壁十年忘何图的博客-CSDN博客_模糊综合评价 （分析实例）

1.适用场景

模糊性常常是信息浓缩所致，目的是为了提高交换的概率，所以不是毫无用处，而是积极的特性。模糊是指客观事物差异的中间过渡中的“不分明性”或“亦此亦彼性”。模糊数学是用来处理涉及模糊概念的问题，尝试使用某种方法将模糊的概念量化，方便进行处理计算。模糊综合评价是模糊数学在评价类问题的一大应用，也就是处理涉及模糊概念的评价类问题。模糊综合评价法是在模糊环境下,考虑了多因素的影响,为了某种目的对一事物作出综合决策的方法。模糊评价问题是要将论域中的对象对应评语集中一个指定的评语，或将方案作为评语集并选择一个最优的方案。

（1）优点

模糊综合评价法具有结果清晰，系统性强的特点，能较好地解决模糊的、难以量化的问题，适合各种非确定性问题的解决。

（2）缺点

• 计算复杂，对指标权重矢量的确定主观性较强。
• 当指标集U较大，即指标集个数较大时，在权矢量和为1的条件约束下，相对隶属度权系数往往偏小，权矢量与模糊矩阵R不匹配，结果会出现超模糊现象，分辨率很差，无法区分谁的隶属度更高，甚至造成评判失败，此时可用分层模糊评估法加以改进。

2.分析步骤

以一级模糊综合评价为例,引入几个关键概念：

• 因素集（评价指标集）： U = { u 1 , u 2 , u 3 ⋯   , u n } \boldsymbol{U=\{u_1,u_2,u_3\cdots,u_n\}} U={u1​,u2​,u3​⋯,un​}
• 评语集（评价的结果）： V = { v 1 , v 2 , v 3 , ⋯   , v m } \boldsymbol{V=\{v_1,v_2,v_3,\cdots,v_m\}} V={v1​,v2​,v3​,⋯,vm​}
• 权重集（指标的权重）： A = { a 1 , a 2 , a 3 , ⋯   , a n } \boldsymbol{A=\{a_1,a_2,a_3,\cdots,a_n\}} A={a1​,a2​,a3​,⋯,an​}

2.1 建立综合评价的因素集

因素集是以影响评价对象的各种因素为元素所组成的一个普通集合，通常用$U$表示， U = { u 1 , u 2 , u 3 ⋯   , u n } \boldsymbol{U=\{u_1,u_2,u_3\cdots,u_n\}} U={u1​,u2​,u3​⋯,un​}，其中元素$u_i$代表影响评价对象的第 i 个因素。这些因素，通常都具有不同程度的模糊性。

以对员工的表现为例，需要从多个方面进行综合评判，如员工的工作业绩、工作态度、沟通能力、政治表现等。所有这些因素构成了评价指标体系集合，即因素集，记为： U = { 政治表现 u 1 , 工作能力 u 2 , 工作态度 u 3 , 工作成绩 u 4 } U = \{政治表现u_1,工作能力u 2,工作态度u_3,工作成绩u_4 \} U={政治表现u1​,工作能力u2,工作态度u3​,工作成绩u4​}。

2.2 建立综合评价的评价集

评价集是评价者对评价对象可能做出的各种结果所组成的集合，通常用V表示， V = { v 1 , v 2 , v 3 , ⋯   , v m } \boldsymbol{V=\{v_1,v_2,v_3,\cdots,v_m\}} V={v1​,v2​,v3​,⋯,vm​}，代表第 j 种评价结果，可以根据实际情况的需要，用不同的等级、评语或数字来表示。

对企业员工的评价有好、良好、中等、较差、很差等。由各种不同决断构成的集合称为评语集，记为： V = { 优秀 v 1 ，良好 v 2 ，中等 v 3 ，较差 v 4 ，很差 v 5 } V = \{优秀v_1 ，良好v_2 ，中等v_3，较差v_4，很差v_5\} V={优秀v1​，良好v2​，中等v3​，较差v4​，很差v5​}。

2.3 确定各因素的权重

评价工作中各因素的重要程度不同，为此给各因素$u_i$一个权重$a_i$，各因素的权重集合的模糊集用A表示： A = { a 1 , a 2 , a 3 , ⋯   , a n } \boldsymbol{A=\{a_1,a_2,a_3,\cdots,a_n\}} A={a1​,a2​,a3​,⋯,an​}

在没有数据时，我们可以通过层次分析法确定权重；在有数据时，我们可以通过熵权法确定权重。

2.4 进行单因素模糊评价，获得评价矩阵

若因素集U中第 $i$ 个元素对评价集 $V$中第 $k$ 个元素的隶属度为 $r_{ik}$ ，则对第 $i$ 个元素单因素评价的结果用模糊集合表示为： R i = { r i 1 , r i 2 , ⋯   , r i m } R_i=\{r_{i1},r_{i2},\cdots,r_{im}\} Ri​={ri1​,ri2​,⋯,rim​}，以 $m$ 个单因素评价集 $R_1,R_2,\cdots ,R_n$ 为行组成矩阵 $R_{n\times m}$ ，称为模糊综合评价矩阵。
$$R=\left[\begin{array}{cccc}{r}_{11}& r_{12}& \cdots & r_{1m}\\ r_{21}& r_{22}& \cdots & r_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ r_{n1}& r_{n2}& \cdots & r_{nm}\end{array}\right]$$

单因素评价结果需进行归一化，即满足：$\sum _{k=1}^m{r}_{ik}=1(i=1,2,\cdots ,n)$

---

模糊子集：

设U是论域，称映射$u_A=U\to [0,1](x\in U,u_A\in [0,1])$，确定了$U$上的模糊子集$A$。映射 $u_A$ 称为A的隶属函数，$u_A(x)$ 称为$x$对$A$的隶属程度，简称隶属度。模糊子集通常简称模糊集。

隶属函数：

隶属函数是隶属度对各个元素的函数，定义域是我们所研究的元素，函数值就是隶属度。隶属度的范围是$[0，1]$，其值越大，就代表越属于这个模糊集合。

模糊分布：

如果考虑以实数集R作论域，通常把实数集R上模糊集的隶属函数称为模糊分布。当所讨论的客观模糊现象的隶属函数与某种给定的模糊分布相类似时，即可选择这个模糊分布作为所求的隶属函数，然后再通过先验知识或数据实验确定符合实际的参数，从而得到具体的隶属函数。

隶属函数有三种确定方式：

• 模糊统计法 ：找多个人去对同个模糊概念进行描述，用隶属频率去定义隶属度 。

• 借助已有的客观尺度

• 指派法：根据问题的性质直接套⽤某些分布

比赛中通常使用指派法确定隶属函数，如下为一些常见的模糊分布：

---

2.5 建立综合评价模型

确定单因素评判矩阵R和因素权向量A之后，通过模糊变化将U上的模糊向量A变为V上的模糊向量B,即：
$$\mathbf{B}=\mathbf{A}\cdot \mathbf{R}=\left[\begin{array}{ccccc}{a}_1& a_2& a_3& \cdots & a_n\end{array}\right]\cdot \left[\begin{array}{cccc}{r}_{11}& r_{12}& \cdots & r_{1m}\\ r_{21}& r_{22}& \cdots & r_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ r_{n1}& r_{n2}& \cdots & r_{nm}\end{array}\right]=\left[\begin{array}{cccc}{b}_1& b_2& \cdots & b_m\end{array}\right]$$
取数值最大的评判结果作为评语，即若 max ⁡ { b 1 , b 2 , ⋯   , b m } = b k \max\{b_1,b_2,\cdots,b_m\}=b_k max{b1​,b2​,⋯,bm​}=bk​,即第k个评判结果最大，则在评语集 V = { v 1 , v 2 , v 3 , ⋯   , v m } V=\{v_1,v_2,v_3,\cdots,v_m\} V={v1​,v2​,v3​,⋯,vm​}中的元素 $v_k$ 即为评价结果。

2.6 多级模糊综合评价方法

• 1.将因素集 U = { u 1 , u 2 , u 3 ⋯   , u n } U=\{u_1,u_2,u_3\cdots,u_n\} U={u1​,u2​,u3​⋯,un​}划分为若干组，得到 U = { u 1 , u 2 ⋯   , u k } U=\{u_1,u_2\cdots,u_k\} U={u1​,u2​⋯,uk​}，其中$U=\bigcup _k^{i=1}{U}_i,U_i\cap U_j=\varnothing (i\ne j)$。称 U = { u 1 , u 2 ⋯   , u k } U=\{u_1,u_2\cdots,u_k\} U={u1​,u2​⋯,uk​}为第一因素集。

• 2.设评语集 V = { v 1 , v 2 , v 3 , ⋯   , v m } V=\{v_1,v_2,v_3,\cdots,v_m\} V={v1​,v2​,v3​,⋯,vm​}，先对第二级因素集$U_i=\left\{u_1^{(i)},u_2^{(i)},\dots ,u_{n_i}^{(i)}\right\}$的$n_i$个因素进行单因素评价，得到单因素评价矩阵$R_i=\left[\begin{array}{cccc}{r}_{11}^{(i)}& r_{12}^{(i)}& \dots & r_{1m}^{(i)}\\ r_{21}^{(i)}& r_{22}^{(i)}& \dots & r_{2m}^{(i)}\\ \dots & \dots & \dots & \dots \\ r_{n_{i}1}^{(i)}& r_{n_{i}2}^{(i)}& \dots & r_{n_{i}m}^{(i)}\end{array}\right]$，设$U_i=\left\{u_1^{(i)},u_2^{(i)},\dots ,u_{n_i}^{(i)}\right\}$的权重为$A=\left\{a_1^{(i)},a_2^{(i)},\dots ,a_{n_i}^{(i)}\right\}$，求得综合评价为$B_i=A_i\cdot R_i(i=1,2,\cdots ,k)$

• 3.再对第一级因素集 U = { u 1 , u 2 ⋯   , u k } U=\{u_1,u_2\cdots,u_k\} U={u1​,u2​⋯,uk​}进行综合评价，设其权重为 A = { a 1 , a 2 , a 3 , ⋯   , a k } A=\{a_1,a_2,a_3,\cdots,a_k\} A={a1​,a2​,a3​,⋯,ak​}，则总评价矩阵为$R=\left[\begin{array}{c}{B}_1\\ B_2\\ \cdots \\ B_k\end{array}\right]$

从而得到综合评价$\mathbf{B}=\mathbf{A}\cdot \mathbf{R}$，按照最大隶属度原则即得相应评语。