用Matlab实现AHP的算法

1、1MATLAB的基本内容MATLAB（MATrix LABoratory，矩阵实验室的缩写）是一种特殊用途的计算机程序优化执行工程和科学计算。它开始为旨在执行矩阵数学程式的生活，但多年来它已发展成为一个灵活的计算系统基本上能够解决任何技术问题。MATLAB具有编程语言的基本特征，使用MATLAB也可以使用像BASIC、FORTRAN、C等传统编程语言一样，进行程序设计，而且简单易学、编程效率高。正因为MATLAB的强大的功能，使得它在许多领域得到广泛应用。在科研与工程应用领域，MATLAB已被广泛地用于科学研究和解决各种具体的实际问题。许多科技工作者选用MATLAB做为计算工具，避免了繁琐的底

2、层编程，从而可以把主要精力和时间花在科学研究和解决实际问题是上，提高了工作效率。1.1 MATLAB矩阵矩阵是MATLAB的基本处理对象，因此根据本文所需，简单介绍所涉及MATLAB矩阵内容。1.1.1 MATLAB矩阵的建立1、直接输入法最简单的建立矩阵的方法是从键盘直接输入矩阵的元素。例如：A=1 2 3;4 5 6;7 8 9A = 1 2 3 4 5 6 7 8 9也可以用回车键代替分号，按下列方式输入：A=1 2 3 4 5 6 7 8 9 2、利用M文件建立矩阵比较大且复杂的矩阵，可以为它专门建立一个M文件，如同下例。利用M文件建立矩阵。启动有关编辑程序或MATLAB文本编辑器，并

3、输入待建矩阵：MYMAT= 111 ,112,113,114,115,116,117,118,119; 211,212,213,214,215,216,217,218,219;把输入的内容以纯文本方式存盘（设文件名为mymatrix.m）。在MATLAB命令窗口中输入mymatrix，即运行该M文件，就会自动建立一个名为MYMAT的矩阵，可供以后使用。1.1.2 矩阵的特征值与特征向量特征值和特征向量在科学研究和工程计算中都有非常广泛地应用。在MATLAB中，计算矩阵A的特征值和特征向量的函数是eig(A),常用的调用格式有3种14：E = eig( A ) ：求矩阵A的全部特征值，构成向量E

4、。V,D=eig(A)：求矩阵A的全部特征值，构成对角矩阵D，并求A得特征向量构成V的列向量。V,D=eig(A,nobablance)：与第2种格式中先对A作相似变换后求矩阵A的特征值和特征向量，而格式3直接求矩阵A的特征值和特征向量。例如：A=1,1,0.5;1,1,0.25;0.5,0.25,2; V,D=eig(A)V = 0.7212 0.4443 0.5315 -0.6863 0.5621 0.4615 -0.0937 -0.6976 0.7103D = -0.0166 0 0 0 1.4801 0 0 0 2.5365求得的3个特征值是-0.0166、1.4801和2.5365，

5、各特征值对应的特征向量为V的各列构成的向量。1.2 MATLAB的M文件用MATLAB语言编写的程序称为M文件。M文件是由若干MATLAB命令组成在一起构成的，它可以完成某些操作，也可以实现某种算法。M文件可以根据调用方式的不同分为两类：命令文件（Script File）和函数文件（Function File）。它们的扩展名均为.m。 函数文件由function语句引导，其基本结构为：function 输出形参表=函数名（输入形参表）注释说明部分函数体语句我们通过举例说明如下：例2-2 分别建立命令文件和函数文件，将求矩阵的一致性指标CI：CI=(max-n)/(n-1)程序1 建立命令文件并

6、以文件名CI.m存盘：max=input(please input max:);n=input(please input n:);CI=(max-n)/(n-1)然后在MATLAB的命令窗口中输入CI即可。程序 2 建立函数文件CI.m。function c=CI(max,n)c=(max-n)/(n-1)然后在MATLAB的命令窗口调用该函数文件。max=input(please input max:);n=input(please input n:);c=CI(max,n)2基于MATLAB的AHP实现21 AHP的MATLAB的计算流程框图根据层次分析法的一般步骤我们得到在MATLAB工

7、具上实现的计算程序流程框图，如图2所示16输入准则层层数输入准则层对目标层判断矩阵元素判断矩阵归一化计算CR0.1输入准则层第j个准则包含的方案数量n判断矩阵归一化计算输入方案层第j个准则层的判断矩阵元素bijCR0.1另m=1;n=1输入方案层第m个准则层的相关矩阵元素计算权值矩阵权值的计算总CR计算开始结束CR0 (2) = l (3) =l按照事物逻辑要求，该矩阵还应具备一致性，即满足：=前面已经给出由于客观事物的复杂性与决策者的认识的多样性，实际问题的成对比较矩阵不可能做到严格上的一致性，因而，借助平均随机一致性指标RI来相对判定其一致性程度。其中表1-4是已经计算好的115阶矩阵的R

8、I值表，但未给出其实现过程，且各文献的RI值表不完全相同。究其原因除没有太大必要介绍外，真正去实现它却有如下三个难度：(1)随机两两判断矩阵中的元素要求是19和它们对应的倒数共17个整数与小数的均匀分布很难处理。(2)一般高级编程语言实现成对比较矩阵及相关计算，非常复杂，且占用内存巨大，耗时多。(3) 随机种子源不能控制。本文使用数学软件包MATLAB对其进行计算。其设计解决思路为：先用软件包中随机函数产生数117的均匀分布的n阶矩阵，然后在软件包中采用不同技巧将它转化为成对比较矩阵，最后用循环语句计算出RI值。结果如下表：表3-1 计算的RI值阶数123456789RI000.54400.8

9、9801.13131.25151.34951.41901.4542以n=4为例，过程详见如下程序清单，其中随机成对比较矩阵的实现见相应注释部分。15MATLAB的程序M文件：function ri %计算RI值的命令文件n=4; ri=0; m=100;rand(seed,21) %控制随机发生器for i=1:m a=ceil (17*rand(n); %产生n阶l17的随机阵 a(find(a=8)=8.1; %消除0为分母 b=1./(a-8); %产生一个辅助阵 a(find(a9)=b(find(a9); %借助b，将917分别转化为 a(find(a=8.1)=8; e=eye(n

10、); %产生一个4阶单位阵 c=1./a; %将a中每个元素换成相应倒数 c=c; %将c转置 c=tril(c,-1); %抽取c的下三角(不含主对角线) a=triu(a,1); %抽取a的上三角(不含主对角线) a=a+c+e; %实现随机成对比较阵a k=size(a,1); %计算a的行维数 ri=ri+(max(abs(eig(a)-k)/(k-1); %计算100次RI值endri/m %计算平均RI值23 AHP各环节的 MATLAB实现以目标矩阵A=,准则层矩阵为P1=, P2=，P3=为例，运用MATLAB进行数据处理。2.3.1 特征向量及其归一化的MATLAB实现MAT

11、LAB中求矩阵特征值和特征向量的函数是eig，其调用的格式为V,D=eig(A),其中，V为特征向量矩阵，D为特征值矩阵。层次分析法中需要求得是最大特征值及对应的归一化特征向量，而且考虑到eig函数在求得的特征值中可能会存在复数。因此，运用直接输入程序代码会产生一定的误差。在此需要对求得的V、D进行适当选择，定义一个M-file maxeigvalvec.m来实现。functionmaxeigval,w=maxeigvalvec(A) %求最大特征值及对应的归一化特征向量 %A为判断矩阵eigvec,eigval=eig(A);eigval=diag(eigval); %特征向量eigvalm

12、ag=imag(eigval);realind=find(eigvalmag=0.10 disp(input(矩阵没通过一致性检验，请重新调整判断矩阵)else disp(input(矩阵通过一致性检验);end在MATLAB中键入如下指令：RIA,CIA= sglsortexamine(max(1),A);RIP1,CIP1= sglsortexamine(max(2),P1);RIP2,CIP2= sglsortexamine(max(3),P2);RIP3,CIP3= sglsortexamine(max(4),P3);运行结果如下：矩阵通过一致性检验矩阵通过一致性检验矩阵通过一致性检验

13、矩阵通过一致性检验2.3.3 一致性检验及总排序的MATLAB实现通过层次单排序（权重）计算后，进行层次合成计算，在此本文定义tolsortvec.m函数计算层次总排序的权重并进行一致性检验。function tw=tolsortvec(utw,dw,CIC,RIC)% 求层次总排序权重并进行一致性检验% utw为上一层因素的总排序权重行向量% dw为下一层因素相对于上一层各因素的层次单排序权重矩阵% CIC为一致性指标列向量% RIC为随机一致性指标列向量tw=dw*utwCR=utw*CIC/(utw*RIC);if CR=0.10 disp(input(层次总排序没通过一致性检验，请重新

14、调整判断矩阵);else disp(input(层次总排序通过一致性检验);end在MATLAB中输入如下指令：dw=zeros(7,3);dw=(1:2,1)=wP1; dw=(3:5,2)=Wp2; dw=(6:7,3)=wP3; CIC=CIP1;CIP2;CIP3;RIC=RIP1;RIP2;RIP3;tw= tolsortvec(wA,dw,CIC,RIC);运行结果如下：tw=0.42470.21230.06670.16460.02700.06980.0349层次总排序通过一致性检验其中tw 是层次总排序结果。因此，根据数据建立如下的层次总排序表。表3-3 层次总排序表（权重）AP

15、A层次P的总排序结果注0.6370 0.25830.1047P10.6667000.42470.3333000.2123P200.258300.066700.637000.164600.104700.0270P3000.66670.0698000.33330.0349注：按概率乘法，P层次总排序指标的权重值为NP层次指标的权重与相应上一层次指标AN层权重的积，且总排序权重值的和为1。2.3.4 选择最优排序计算出层次总排序后，为了使决策者能迅速得出结果，本文对层次总排序进行最优排序。运用MATLAB键入如下指令：n=length(tw);for i=1:nt=max(tw);b(i)=t;m

16、n=find(a=t);tw(n)=;endb运行结果如下： b=0.42470.21230.16460.06980.06670.03490.0276利用MATLAB大大缩短了计算复杂矩阵的时间，为决策者节省了宝贵的时间，从而有更多的精力投入其他事务。3基于MATLAB的AHP应用31 挑选合适工作问题某毕业生选择工作，经双方恳谈，假设已有三个单位C1,C2,C3表示愿意录用他。该生对三个单位进行了解后，选取了一些中间指标进行考察，例如单位的研究课题，发展前途，待遇，同事情况，地理位置，单位名气等。根据层次分析法，试求该生工作优先排序（给出权值、计算程序），并给出最终选择决策。现以A、B、C表

17、示选择工作的三个层次，建立如下结构模型：图3 选择单位层次结构图根据成对比较法，得到相应判断矩阵如下表：表4-1 A-B判断矩阵B1B2B3B4B5B6B1111411/2B2112411/2B311/21531/2B41/41/41/511/31/3B5111/3311B6222331表4-2 B1C判断矩阵11/41/241321/31表4-3 B2C判断矩阵C1C2C3C111/41/5C2411/2C3521表4-4 B3C判断矩阵C1C2C3C1131/3C21/311/7C3371表4-5 B4C判断矩阵C1C2C3C111/35C2317C31/51/71表4-6 B5C判断矩阵

18、C1C2C3C1117C2117C31/71/71表4-7 B6C判断矩阵C1C2C3C1179C21/711C31/911现在在MATLAB中分别用直接输入程序法和M文件方法求解。1）、直接输入代码法：在MATLAB中输入如下程序：A = 1,1,1,4,1,1/2; 1,1,2,4,1,1/2; 1,1/2,1,5,3,1/2; 1/4,1/4,1/5,1,1/3,1/3; 1,1,1/3,3,1,1; 2,2,2,3,3,1; B1 = 1,1/4,1/2;4,1,3;2,1/3,1; B2 = 1,1/4,1/5;4,1,1/2;5,2,1; B3 = 1,3,1/3;1/3,1,1/

19、7;3,7,1; B4 = 1,1/3,5;3,1,7;1/5,1/7,1; B5 = 1,1,7;1,1,7;1/7,1/7,1; B6 = 1,7,9;1/7,1,1;1/9,1,1; BS = B1,B2,B3,B4,B5,B6; m = length(B1);n = length(A); %随机一致性指标RI RI = 0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51; Wa,LA = eig (A) %求A的特征向量WA和特征根LAMaxn=input(please input largest eigenvalue:); %输入最大特征

20、根 CIn = (Maxn- n) / (n - 1); CRn = CIn / RI(n); %A的一致性比率CRnWA=Wa(:,1)/sum(Wa(:,1); %特征向量归一化 if CRn 0.10 fprintf(A 的CR %f 通过一致性检验!n,CRn); %控制文本格式 else fprintf(A 的CR %f 未通过一致性检验!n,CRn); end for k = 1:n %求B的特征向量WK和特征根LK WB,LK = eig( BS(1:3,(k-1)*m+1:(k-1)*m+3) )Max(k)=input(please input largest eigenva

21、lue:); CIm(k) = (Max(k)- m) / (m - 1);RIm(k) = RI(m); CRm(k) = CIm(k) / RIm(k); %B的一致性比率CRm WK(:,k)= WB(:,1)/sum(WB(:,1); end for k = 1:n if CRm(k) 0.10 fprintf(B%d的CR %f 通过一致性检验!n,k,CRm(1,k); %控制文本格式 else fprintf(B%d的CR %f 未通过一致性检验!n,k,CRm(1,k); end enddisp(准则层对目标层权向量); disp(WA);disp(方案层对准则层权向量); d

22、isp(WK);E = WK * WAdisp(方案层组合权向量);disp(E);CI = CIm * WA;RI = RIm * WA;CR = CI / RI; %组合一致性比率CR if CR 0.10 fprintf(组合一致性比率CR %f 通过一致性检验!n,CRn); else fprintf(组合一致性比率CR %f 未通过一致性检验!n,CRn); endMAX,CHOICE = max(E); %最佳选择CHOICE MATLAB运行结果如下：Wa = -0.3396 -0.1255 - 0.0563i -0.1255 + 0.0563i 0.7354 -0.1896 +

23、 0.3838i -0.1896 - 0.3838i -0.4038 -0.1884 - 0.5736i -0.1884 + 0.5736i -0.6464 -0.4476 - 0.2693i -0.4476 + 0.2693i -0.4249 0.6724 0.6724 0.0834 0.3884 - 0.0605i 0.3884 + 0.0605i -0.1063 -0.0138 + 0.0429i -0.0138 - 0.0429i -0.0405 -0.0592 - 0.0922i -0.0592 + 0.0922i -0.3298 -0.1384 + 0.3417i -0.1384

24、- 0.3417i -0.1217 0.0035 + 0.1640i 0.0035 - 0.1640i -0.6488 -0.1467 - 0.0710i -0.1467 + 0.0710i 0.1337 0.5920 0.5920 LA = 6.6178 0 0 0 0 0 0 -0.1557 + 1.2808i 0 0 0 0 0 0 -0.1557 - 1.2808i 0 0 0 0 0 0 -0.0603 0 0 0 0 0 0 -0.1230 + 0.5461i 0 0 0 0 0 0 -0.1230 - 0.5461i please input largest eigenvalue

25、:6.6178A 的CR 0. 通过一致性检验!WB = 0.1999 0.1000 + 0.1731i 0.1000 - 0.1731i 0.9154 -0.9154 -0.9154 0.3493 0.1747 - 0.3025i 0.1747 + 0.3025iLK = 3.0183 0 0 0 -0.0091 + 0.2348i 0 0 0 -0.0091 - 0.2348iplease input largest eigenvalue: 0.1999 WB = 0.1460 0.0730 + 0.1265i 0.0730 - 0.1265i 0.4994 0.2497 - 0.4325

26、i 0.2497 + 0.4325i 0.8540 -0.8540 -0.8540 LK = 3.0246 0 0 0 -0.0123 + 0.2725i 0 0 0 -0.0123 - 0.2725iplease input largest eigenvalue:3.0246WB = 0.3382 -0.1691 + 0.2929i -0.1691 - 0.2929i 0.1226 -0.0613 - 0.1062i -0.0613 + 0.1062i 0.9331 0.9331 0.9331 LK = 3.0070 0 0 0 -0.0035 + 0.1453i 0 0 0 -0.0035

27、 - 0.1453iplease input largest eigenvalue: 3.0070 WB = 0.3928 -0.1964 + 0.3402i -0.1964 - 0.3402i 0.9140 0.9140 0.9140 0.1013 -0.0506 - 0.0877i -0.0506 + 0.0877iLK = 3.0649 0 0 0 -0.0324 + 0.4448i 0 0 0 -0.0324 - 0.4448iplease input largest eigenvalue:3.0649 WB = -0.7035 -0.3392 0.6619 -0.7035 -0.92

28、33 -0.7495 -0.1005 0.1804 0.0125LK = 3.0000 0 0 0 0.0000 0 0 0 -0.0000please input largest eigenvalue: 3.0000 WB = -0.9844 0.9844 0.9844 -0.1293 -0.0647 - 0.1120i -0.0647 + 0.1120i -0.1189 -0.0595 + 0.1030i -0.0595 - 0.1030iLK = 3.0070 0 0 0 -0.0035 + 0.1453i 0 0 0 -0.0035 - 0.1453iplease input larg

29、est eigenvalue: 3.0070 B1的CR -2. 通过一致性检验!B2的CR 0. 通过一致性检验!B3的CR 0. 通过一致性检验!B4的CR 0. 通过一致性检验!B5的CR 0. 通过一致性检验!B6的CR 0. 通过一致性检验!准则层对目标层权向量 0.1507 0.1792 0.1886 0.0472 0.1464 0.2879方案层对准则层权向量 0.1365 0.0974 0.2426 0.2790 0.4667 0.7986 0.6250 0.3331 0.0879 0.6491 0.4667 0.1049 0.2385 0.5695 0.6694 0.0719

30、 0.0667 0.0965E = 0.3952 0.2996 0.3052方案层组合权向量 0.3952 0.2996 0.3052组合一致性比率CR 0. 通过一致性检验!CHOICE = 12）、M文件法：根据第三章节的层次分析法各环节的MATLAB实现定义如下三个函数：（1） 定义maxeigvalvec.m来实现最大特征根及对应的归一化特征向量。functionmaxeigval,w=maxeigvalvec(A) %求最大特征值及对应的归一化特征向量 %A为判断矩阵eigvec,eigval=eig(A);eigval=diag(eigval); %特征向量eigvalmag=im

31、ag(eigval);realind=find(eigvalmag=0.10 disp(input(矩阵没通过一致性检验，请重新调整判断矩阵)else disp(input(矩阵通过一致性检验);end（3）定义tolsortvec.m函数计算层次总排序的权重并进行一致性检验。function tw=tolsortvec(utw,dw,CIC,RIC)% 求层次总排序权重并进行一致性检验% utw为上一层因素的总排序权重行向量% dw为下一层因素相对于上一层各因素的层次单排序权重矩阵% CIC为一致性指标列向量% RIC为随机一致性指标列向量tw=dw*utw;CR=utw*CIC/(utw*

32、RIC);if CR=0.10 disp(input(层次总排序没通过一致性检验，请重新调整判断矩阵);else disp(input(层次总排序通过一致性检验);end%main.m主程序clear;A = 1,1,1,4,1,1/2;1,1,2,4,1,1/2;1,1/2,1,5,3,1/2;1/4,1/4,1/5,1,1/3,1/3;1,1,1/3,3,1,1;2,2,2,3,3,1;B1 = 1,1/4,1/2;4,1,3;2,1/3,1;B2 = 1,1/4,1/5;4,1,1/2;5,2,1;B3 = 1,3,1/3;1/3,1,1/7;3,7,1;B4 = 1,1/3,5;3,1

33、,7;1/5,1/7,1;B5 = 1,1,7;1,1,7;1/7,1/7,1;B6 = 1,7,9;1/7,1,1;1/9,1,1;max(1),wA=maxeigvalvec(A);max(2),wB1=maxeigvalvec(B1);max(3),wB2=maxeigvalvec(B2);max(4),wB4=maxeigvalvec(B4);max(5),wB3=maxeigvalvec(B3);max(6),wB5=maxeigvalvec(B5);max(7),wB6=maxeigvalvec(B6);RIA,CIA=sglsortexamine(max(1),A);RIB1,CIB1=sglsortexamine(max(2),B1);RIB2,CIB2=sglsortexamine(max(3),B2);RIB3,CIB3=sglsortexamine(max(4),B3);RIB4,CIB4=sglsortexamine(max(5),B4);RIB5,CIB5=sglsortexamine(max(6),B5);RIB6,CIB6=sglsortexamine(max(7),B6);dw=zeros(18,6);dw(1:3,1)=wB1;dw(4:6,2)=wB2;dw(7:9,3)=wB3;dw(10:12,4)=wB4;dw(1