层次分析法 - 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)

E= 

    -1.3166

     5.3166

     2.0000

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

[V,D]=eig(A)

V=

    0.8364   -0.6336   -0.2408

   -0.3875   -0.5470   -0.8427

   -0.3875   -0.5470    0.4815

D=

   -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

   0   0   0

   1   0   0

---

3.判断语句

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

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

if a>=85

  disp('成绩优秀')

else if  a>=60

  disp('成绩合格')

else

  disp('成绩挂科')

end

---

层次分析法代码部分

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

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

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

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

---

计算权重：

%%输入判断矩阵

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

0   0   0   0   0

0   0   0   0   0

0   0   0   0   0

r=

    1

c=

    1

%找到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知识点

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

---

代码优化：

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

1.

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

         end

      end

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

     ERROR=2;

  end

end

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

if ERROR==0

    if n>15

       ERROR=3;

    end

end

if ERROR==0

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

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

       ERROR=4;

    end

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

---

感谢观看！