层次分析法（AHP)--代码书写部分

层次分析法（AHP)–代码书写部分

在正常的层次分析法的过程中，如果判断矩阵是一致矩阵，就直接进行权重计算就可以了。但是如果判断矩阵是非一致性矩阵，我们是先进行一致性检验，再进行判断矩阵的权重计算。但是在一致性检验的时候需要计算一致性指标CI，在CI中需要计算最大特征值${\lambda }_{max}$，而${\lambda }_{max}$是在计算判断矩阵权重时计算的。
$$CI=\frac{{\lambda }_{max}-n}{n-1}$$
所以我们在代码中先进行判断矩阵权重的计算，再进行一致性检验。接下来开始我们的代码书写：

输入判断矩阵

%% 输入判断矩阵
clear;clc
disp{'请输入判断矩阵A'} %disp为输出函数，类似c语言中的printf函数 
A=[1 1 4 1/3 3;
   1 1 4 1/3 3;
   3 3 3 1 3;
   1/3 1/3 2 1/3 1]

算数平均法求权重

第一步：因为判断矩阵归一化是判断矩阵的每一个元素除以其列的和，所以我们构造一个和判断矩阵相同大小的各列和的矩阵，再利用“./”来实现归一化。

%% 对判断矩阵进行归一化
Sum_A=sum(A) %对A进行各行相加（按列求和）,得到一个行向量
[n,n]=size(A) %% 因为判断矩阵是一个方阵，所以这里的r（行）和c（列）相同，我们用一个字母n代替
SUM_A=repmat（Sum_A,n,1) % 这里的repmat函数是将Sum_A重复n行一列，得到新的矩阵。
Stand_A=A./SUM_A

第二步：

%% 将归一化的各列相加（按行求和）
sum(Stand_A,2)%%这里的2是各列相加的意思（默认是1，各行相加）

第三步：

disp('算数平均法求权重的结果为：')
disp(sum(Stand_A,2)/n)
%首先对标准化的矩阵按照行求和，得到一个列向量
%然后再将这个列向量的每个元素同时除以n即可

代码总结：

Sum_A=sum(A)
[n,n]=size(A)
SUM_A=repmat（Sum_A,n,1)
Stand_A=A./SUM_A
sum(Stand_A,2)
disp('算数平均法求权重的结果为：')
disp(sum(Stand_A,2)/n)

几何平均法求权重

第一步:将A的元素的各列相乘，得到一个新的列向量

% 将A的元素的各列相乘，得到一个新的列向量
clc;A
Prduct_A=prod(A,2)
% prod函数和sum函数类似，不过是用于乘法，dim=2表示维度，各列相乘。

第二步:将得到的列向量每一个元素开n次方

%将第一步得出的列向量的每个分量开n次方
Prduct_n_A=Prduct_A.^(1/n)
% 这里对每一个元素进行乘方操作，因此需要加.，
% ^号代表乘方，这里是开n次方，所以等价于求1/n次方

第三步:将开过n次方的列向量元素进行归一化

% 将开方后的列向量中的每一个元素除以这个列向量的和
disp('几何平均法求得的权重的结果为：')
disp(Prduct_n_A./sum(prduct_n_A))

代码总结：

clc;A
Prduct_A=prod(A,2)
Prduct_n_A=Prduct_A.^(1/n)
disp('几何平均法求得的权重的结果为：')
disp(Prduct_n_A./sum(prduct_n_A))

特征值法求权重

第一步:求出矩阵A的最大特征值以及其对应的特征向量

clc
[V,D]=eig(A)%V是特征向量，D是由特征值构成的对角矩阵（除了对角线元素，其余位置元素全为0）
Max_dig=max(max(D)) %max函数求每一列的最大值得到一个行向量，再进行一次max操作就可以得到最大特征值了
% 这里也可以写成 max(D(:))
%接下来需要找到最大特征值所对应的特征向量，用到find函数，但是find函数只能找到矩阵中不为0元素的位置
%所以我们先将矩阵中不是最大特征值的位置的元素变为0，再使用find函数
D==Max_eig
[r,c]=find(D==Max_dig,1)
%找到D中与最大特征值相等的位置，记录它的行和列

第二步:对求出的特征向量进行归一化即可得到我们的权重

disp('特征值法求权重的结果为：');
disp(V(:,c)./sum(V(:,c)))
%这里的V(:,c)表示取第c列（最大特征值所在列）的全部元素

代码总结：

clc
[V,D]=eig(A)
Max_dig=max(max(D)) 
D==Max_eig
[r,c]=find(D==Max_dig,1)
disp('特征值法求权重的结果为：');
disp(V(:,c)./sum(V(:,c)))

计算一致性比例CR

计算一致性指标CI和查找平均随机一致性指标RI，最后计算一致性比例CR
$$CI=\frac{{\lambda }_{max}-n}{n-1}$$

n	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
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

clc
CI=(Max_eig-n)/(n-1)；
RI=[0 0 0.52 0.89 0.12 0.26 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59];
%这里RI最多支持n=15
CR=CI/RI(n);

disp(['一致性指标CI=',num2str(CI)]);
disp('一致性比例RI='):disp(CR);
%这里输出使用了两种方法，第一种是以向量的形式输出，第二种是直接使用disp函数输出
if CR<0.01
	disp('因为CR<0.10，所以该判断矩阵A的一致性可以接受！');
else
	disp('注意：CR>=0.10，因此该判断矩阵A需要进行修改！'）

到这里我们层次分析法的代码部分就结束了，下面是总的代码：

Sum_A=sum(A)
[n,n]=size(A)
SUM_A=repmat（Sum_A,n,1)
Stand_A=A./SUM_A
sum(Stand_A,2)
disp('算数平均法求权重的结果为：');
disp(sum(Stand_A,2)/n);

Prduct_A=prod(A,2)
Prduct_n_A=Prduct_A.^(1/n)
disp('几何平均法求得的权重的结果为：');
disp(Prduct_n_A./sum(prduct_n_A));

[V,D]=eig(A)
Max_dig=max(max(D)) 
D==Max_eig
[r,c]=find(D==Max_dig,1)
disp('特征值法求权重的结果为：');
disp(V(:,c)./sum(V(:,c)));

CI=(Max_eig-n)/(n-1)；
RI=[0 0 0.52 0.89 0.12 0.26 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59];
CR=CI/RI(n);
disp(['一致性指标CI=',num2str(CI)]);
disp('一致性比例RI='):disp(CR);
if CR<0.01
	disp('因为CR<0.10，所以该判断矩阵A的一致性可以接受！');
else
	disp('注意：CR>=0.10，因此该判断矩阵A需要进行修改！'）;