层次分析法 - MATLAB代码详解

1.Matlab中求特征值和特征向量

在Matlab中，计算矩阵A的特征值和特征向量的函数是eig(A),其中最常用的两个用法：

A=[1 2 3 ; 2 2 1 ; 2 0 3]

(1) E=eig(A) :求矩阵A的全部特征值，构成向量E。（几阶方阵就有几个特征值）

E=eig(A)

-1.3166

5.3166

2.0000

(2) [V,D]=eig(A) :求矩阵A的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量。（V的每一列都是D中与之相同列的特征值的特征向量）。

[V,D]=eig(A)

0.8364 -0.6336 -0.2408

-0.3875 -0.5470 -0.8427

-0.3875 -0.5470 0.4815

-1.3166 0 0

0 5.3166 0

0 0 2.0000

2.矩阵与常数的大小判断运算

共有三种运算符：大于 > ; 小于 < ; 等于 == （一个等号表示赋值，两个等号表示判断）

X=[1 -3 0;0 0 8;4 0 6]

1 -3 0

0 0 8

4 0 6

X>0

ans=

3×3 logical 数组（逻辑数组）（大于0，是，1；小于0，否，0）

1 0 0

0 0 1

1 0 1

X==4

ans=

3×3 logical 数组

0 0 0

1 0 0

3.判断语句

Matlab的判断语句，if所在的行不需要冒号，语句的最后一定要以end结尾;中间的语句要注意缩进。

a=input('请输入考试分数:')

if a>=85

disp('成绩优秀')

else if a>=60

disp('成绩合格')

else

disp('成绩挂科')

end

层次分析法代码部分

过程：一致性检验 $\rightarrow$ 判断矩阵计算权重（算术平均值法，几何平均值法，特征值法）

在论文写作中，应该先对判断矩阵进行一致性检验，然后再计算权重，因为只有判断矩阵通过了一致性检验，其权重才是有意义的。

代码顺序：先计算权重 $\rightarrow$ 一致性检验（仅为了顺应计算过程，在逻辑上过不去）

注意：只有非一致矩阵的判断矩阵才需要进行一致性检验，如果你的判断矩阵本身就是一个一致矩阵，则没有必要进行一致性检验。（自己给的时候最好给一个非一致矩阵，再通过一下一致性检验。）

计算权重：

%%输入判断矩阵

clear;clc

disp('请输入判断矩阵A: ')

%A=input('判断矩阵A=')

A=[1 1 4 1/3 3;

1 1 4 1/3 3;

1/4 1/4 1 1/3 1/2;

3 3 3 1 3;

1/3 1/3 2 1/3 1] 只能给1到9和他们的倒数！！

%matlab矩阵有两种写法，可以直接写到一行；

%A=[1 1 4 1/3 3;1 1 4 1/3 3;1/4 1/4 1 1/3 1/2;3 3 3 1 3;1/3 1/3 2 1/3 1]

%也可以写成多行：

A=[1 1 4 1/3 3;

1 1 4 1/3 3;

1/4 1/4 1 1/3 1/2;

3 3 3 1 3;

1/3 1/3 2 1/3 1]

%两行之间以分号结尾（最后一行分号可加可不加），同行元素之间可以用空格（或者逗号）分开

%%方法1：算术平均法求权重

%第一步：将判断矩阵按照列归一化（每一个元素除以其所在列的和）

Sum_A=sum(A)

[n,n]=size(A) %也可以写成n=size(A,1)

%因为我们的判断矩阵A是一个方阵，所以这里的r和c相同，就可以用同一个字母n表示。

SUM_A=repmat(Sum_A,n,1) %repeat martix在（重复）的缩写

%另一种替代的方法如下：

%SUM_A=[];

%for i=1:n %循环，这一行后面不能加冒号（和python不同），这里表示循环n次

%SUM_A=[SUM_A;Sum_A];

%end

SUM_A

Stand_A=A./SUM_A ——归一化完成

%这里直接将两个矩阵对应的元素相除即可

%第二步：将归一化的各列相加（按行求和）

sum(Stand_A,2)

%第三步：将相加后得到的向量中每个元素除以n即可得到权重向量

disp('算术平均法求权重的结果为： ');

disp(sum(sum(Stand_A,2))/n)

%首先对标准化后的矩阵按照行求和，得到一个列向量

%然后再将这个列向量的每个元素同时除以n即可（注意这里也可以用./哦）

%%方法2：几何平均法求权重

%第一步：将A的元素按照行相乘得到一个新的列向量

product_A=prod(A,2)

%prod函数和sum函数类似，一个用于乘，一个用于加dim=2 维度是行（每一行的所有元素相加或相乘，得到的结果是一个列向量）

%第二步：将新的向量的每个分量开n次方

Product_n_A=product_A.^(1/n)

% 这里对每个元素进行乘方操作，因此要加.号。

% ^符号表示乘方，开n次方等价于求1/n次方

%第三步：对该列向量进行归一化即可得到权重向量

%将这个列向量中的每一个元素除以这一个向量的和即可

disp('几何平均法求权重的结果为: ');

disp(Product_n_A./sum(Product_n_A))

%%方法三：特征值法求权重

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

[V,D]=eig(A)

%V是特征向量，D是由特征值构成的对角矩阵（除了对角元素外，其余元素全为0）

Max_eig=max(max(D)) %也可以写成max(D(:)),得到一个行向量,max函数是求每一列的最大值，最后结果是一个行向量，若出现虚数，则比较其模的大小。

%那么怎么找到最大特征值所在的位置呢？

%需要用到find函数，它可以用来返回向量或矩阵中不为0的元素的位置索引。

%那么问题来了，我们要得到最大特征值的位置，就需要将包含所有特征值的这个对角矩阵D中，不等于最大特征值的位置全变为0

%这时候可以用到矩阵与常数的大小判断运算

D==Max_eig

[r,c]=find(D==Max_eig,1)

ans=

5×5 logical 数组

1 0 0 0 0

0 0 0 0 0

%找到D中第一个与最大特征值相等的元素的位置，记录他的行和列。

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

V(:,c) %(提取某一列的所有元素的函数)

disp('特征值法求权重的结果为: ');

disp( V(:,c) ./ sum(V(:,c)) )

%我们先根据上面找到的最大特征值的列数c找到对应的特征向量，然后再进行标准化。

一致性检验

%%计算一致性比例CR

CI=(Max_eig - n) / (n-1); %Max_eig是最大特征值，n是矩阵的大小

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];

%这里的RI最多支持n=15。

CR=CI/RI(n);

disp('一致性指标CI='); disp(CI); %也可用字符串连接函数disp(['一致性指标CI='],num2str(CI)])

disp('一致性比例CR='); disp(CR);

if CR<0.10

disp('因为CR<0.10 ，所以该判断矩阵A的一致性可以接受！');

else

disp('注意：CR>=0.10 ，因此该判断矩阵A需要进行修改！');

end

总结：层次分析法中需要用到的matlab知识点

Matlab基本的小常识：分号的作用，注释的快捷键，clc和clear,disp和input
sum函数
Matlab中如何提取矩阵中指定位置的元素？
size函数
repmat函数
Matlab中矩阵的运算（加点和不加点）
Matlab中求特征值和特征向量
find函数的基本用法
矩阵和常数的大小判断运算
判断和循环语句

代码优化：

请对代码进行优化，例如输入判断矩阵A时，是否能自动检查矩阵A是否为正互反矩阵？
如果我们输入的是一个二阶的判断矩阵，请观察结果有什么问题？怎么改进代码来修正这个问题？

disp('请输入判断矩阵A')

A=input('A=');

%在开始下面正式的步骤之前，我们有必要检验A是否因为粗心而输入有误

ERROR=0; %默认输入是没有错误的

%（1）检查矩阵A的维数是否不大于1或不是方阵

[r,c]=size(A);

if r~=c || r<=1

%~=是不等号； ||表示逻辑运算符 ‘ 或 ’

ERROR=1;

end

%（2）检查矩阵A是否为正互反矩阵 a_ij > 0 且 a_ij*a_ji = 1

if ERROR == 0

[n,n]=size(A)

%判断是否有元素小于0

for i=1:n

for j=1:n

if A(i,j)<=0

ERROR=2;

end

if sum(sum(A<=0))>0

ERROR=2;

end

%检验n是否超过了15，因为RI向量为15维

if ERROR==0

if n>15

ERROR=3;

end

if ERROR==0

%判断a_ij*a_ji = 1是否成立

if sum(sum(A'.*A~= ones(n)))>0

ERROR=4;

end

%A'表示A的转置矩阵，即将a_ij和a_ji互换位置

%ones(n)函数生成一个n*n的全为1的方阵，zero(n)函数生成一个n*n的全为0的方阵

%ones(m,n)函数生成一个m*n的全为1的方阵

if ERROR==0

后面继续计算权重。

...........................................................................................................................................................

else if ERROR==1

disp('请检查矩阵A的维数是否不大于1或不是方阵')

else if ERROR==2

disp('请检查矩阵A中有元素小于等于0')

else if ERROR==3

disp('A的维数n超过了15，请减少准则层的数量')

else if ERROR==4

disp('请检查矩阵A中存在i,j不满足A_ij*A_ji=1')

end

%%代码仅供参考，一定不要用于自己的数模论文中！

2.RI=2，CR>=0.10，结果需要修正。

二阶正互反矩阵一定是一致矩阵，所以不需要进行一致性检验。

代码优化：

CI=(Max_eig - n) / (n-1);

RI=[0 0.001 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);

disp('一致性指标CI='); disp(CI);

disp('一致性比例CR='); disp(CR);

if CR<0.10

disp('因为CR<0.10 ，所以该判断矩阵A的一致性可以接受！');

else

disp('注意：CR>=0.10 ，因此该判断矩阵A需要进行修改！');

end

感谢观看！