用Matlab实现AHP的算法

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

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

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

4、阵D,并求A得特征向量构成V的列向量。V,D=eig(A, ' nobablanee '):与第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.72120.44430.5315-0.68630.56210.4615-0.0937-0.69760.7103D =-0.01660001.4801000 2.5365求得的3个特征值是-0.0166、1.4801和2.5365，各特征值对应的特征向量为 V 的各列构成的向量。1.2 M

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

6、please in put max:');n=in put('please in put n:');Cl=(max-n)/(n-1)然后在MATLA的命令窗口中输入Cl即可。程序2建立函数文件Cl.m。fun etio n e=CI(max ,n)e=(max- n)/(n-1)然后在MATLA的命令窗口调用该函数文件。 max=in put('please in put max:');n=in put('please in put n:');c=CI(max ,n)2.基于MATLAB勺AHP实现2.1 AHP的MATLAB的计算流程框

7、图根据层次分析法的一般步骤我们得到在 MATLA工具上实现的计算程序流程框图,16如图2所示输 入 准则层层数“判准矩阵元素标层判断矩阵归化计算输入准则层第j个准n则 包含的方案数量层入勺判断矩阵j个素则判断矩阵归一化计算图2 以MATLA实现的层次分析法的计算流程框图通过流程框图，层次分析的基本步骤如下：第一步：准则层对目标层的判断矩阵归一化且判断是否满足一致性；第二步：第一步满足时，将方案层对准则层的判断矩阵归一化并判断其一致性； 第三步：当第一、二步满足时，求方案层的总排序权值与总CR并判断一致性。2. 2平均随机一致性指标的 MATLAB实现运用层次分析法决策者需要通过反复地解决决策问

8、题，将同一层次的各元素与上一层次中某一准则的重要性进行比较，从而构造出两两判断比较矩阵 A=(aij ) nn(称为成对比较矩阵)。前面已经描述了九级标度法，此处运用其他描述，则这些成对比较矩阵 应满足如下条件：(I) aij >0(2) aij aji = l (3)=1按照事物逻辑要求，该矩阵还应具备一致性，即满足：aij a jk =aik前面已经给出由于客观事物的复杂性与决策者的认识的多样性，实际问题的成对比较 矩阵不可能做到严格上的一致性，因而，借助平均随机一致性指标 RI来相对判定其一 致性程度。其中表1-4是已经计算好的115阶矩阵的RI值表，但未给出其实现过程， 且各文献

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

10、.3491.4191.4540035502以n=4为例，过程详见如下程序清单，其中随机成对比较矩阵的实现见相应注释 部分。15MATLAB勺程序M文件：fun ctio n ri%计算RI值的命令文件n=4; ri=0; m=100;ran d('seed',21)%控制随机发生器for i=1:ma=ceil (17*ra nd(n);%产生n阶117的随机阵a(fi nd(a=8)=8.1;%消除0为分母b=1./(a-8);%产生一个辅助阵a(fi nd(a>9)=b(fi nd(a>9);%昔助b,将917分别转化为丄2a(fi nd(a=8.1)=8;e=

11、eye( n);%产生一个4阶单位阵c=1./a;%将a中每个元素换成相应倒数c=c'%将c转置c=tril(c,-1);%抽取c的下三角(不含主对角线a=triu(a,1);%抽取a的上三角(不含主对角线a=a+c+e;%实现随机成对比较阵a)19k=size(a,1);%计算a的行维数ri=ri+(max(abs(eig(a)-k)/(k-1);%计算100次RI值endri/m%计算平均RI值2. 3 AHP各环节的 MATLAB实现135P1=12135以目标矩阵 A= 1313 ,准则层矩阵为,P 2= 1313，1513112 1151311 2P3= 1/2 1为例，运用

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

13、gval);reali nd=fi nd(eigvalmagveps); realeigval=eigval(reali nd);%maxeigval=max(realeigval)%in dex=fi nd(eigval=maxeigval);veci nit=eigvec(:,i ndex);%w=veci nit./sum(veci nit)%在MATLA中键入如下指令：在此需要对求得的V、D进行适当选择，定求最大特征值及对应的归一化特征向为判断矩阵特征向量实特征根最大特征值最大特征值对应的特征向量特征向量归一化A= 1,3,5;1/3,1,3;1/5,1/3,1;P1=1,2;1/2,

14、1;P2= 1,3,5;1/3,1,3;1/5,1/3,1;P3=1,2;1/2,1;max(1),wA= maxeigvalvec(A); max(2)wP1= maxeigvalvec(P1);max (3) ,wP2= maxeigvalvec(P2); max ,wP3= maxeigvalvec(P3);MATLA运行结果如下：maxeigval=3.0385w =0.63700.25830.1047maxeigval=2w=0.66670.3333maxeigval=3.0385w =0.63700.25830.1047maxeigval=2w=0.66670.3333一致性检验及

15、单排序的MATLA实现由AHP勺MATLA的计算流程图知，必须对各层次间的判断矩阵进行层次单排序和 一致性检验。因此，定义sglsortexami ne.m函数来实现层次单排序的一致性检验。function RI,CI=sglsortexamine(maxeigval,A)%层次分析法单排序一致性检验%maxeigval为最大特征值，A为判断矩阵n=size(A,1);RIT= 0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51;RI=RIT( n);CI=(maxeigval-n)/(n-1);CR=CI/RI;if CR>=0.10dis

16、p(i nput('矩阵没通过一致性检验，请重新调整判断矩阵')elsedisp(input('矩阵通过一致性检验');end在MATLA中键入如下指令：RIA,CIA= sglsortexami ne(max(1),A);RIP1,CIP1= sglsortexami ne(max(2),P1);RIP2,CIP2= sglsortexami ne(max (3),P2);RIP3,CIP3= sglsortexami ne(max ,P3);运行结果如下：矩阵通过一致性检验 矩阵通过一致性检验 矩阵通过一致性检验 矩阵通过一致性检验一致性检验及总排序的MAT

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

18、调整判断矩阵elsedisp(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 层次总

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

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

21、、B、C表示选择工作的三个层次，建立如下结构模型：图3选择单位层次结构图根据成对比较法，得到相应判断矩阵如下表：表4-1 A-B 判断矩阵BiB2B3B5B6B111411/2112411/2R11/21531/21/41/41/511/31/3111/3311B6222331表4-2 B iC判断矩阵11/41/241321/31表4-3 B2C判断矩阵C1C2C3C111/41/5C2411/2C3521表4-4 B 3C判断矩阵C1C2C3G131/3C21/311/7C3371表4-5 B 4C判断矩阵GC2G11/35C2317Q1/51/71表4-6 B5C判断矩阵C1aG117C

22、2117Q1/71/71表4-7B 6C判断矩阵C1aQG179Q1/711Q1/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/331,1;222,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,7;1/5,

23、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 = len gth(B1);n = len gth(A);%随机一致性指标RIRI = 0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51;Wa,LA = eig (A)%征根LAMaxn二in put('please in put largest eige nvalue:');%CIn = (Max n- n) / (n - 1);CRn 二 CIn / RI(

24、n);%A求A的特征向量WA和特输入最大特征根的一致性比率CRnWA=Wa(:,1)/sum(Wa(:,1);%特征向量归一化if CRn < 0.10fprintf('A的CR %f通过一致性检验！n',CRn); % 控制文本格式elsefprintf('A end for k = 1:n LK的CR %f未通过一致性检验！n',CRn);%求B的特征向量WK和特征根WB,LK = eig( BS(1:3,(k-1)*m+1:(k-1)*m+3)Max(k)=i nput('please in put largest eige nvalue:&

25、#39;);CIm(k) = (Max(k)- m) / (m - 1);RIm(k) = RI(m);CRm(k) = CIm(k) / RIm(k);%B的一致性比率CRmWK(:,k)= WB(:,1)/sum(WB(:,1);endfor k = 1:nif CRm(k)< 0.10fprin tf('B%dCR %f 通验!n',k,CRm(1,k);控制文本格式elsefprin tf('B%d的CR %f未通过一致性检验!n',k,CRm(1,k);end enddisp('准则层对目标层权向量');disp(WA);disp

26、('方案层对准则层权向量');disp(WK);E = WK * WAdisp（'方案层组合权向量'）;disp（E）;CI = CIm * WA;RI = Rim * WA;CR = CI / RI;%组合一致性比率CRif CR < 0.10fprintf（'组合一致性比率CR %f通过一致性检验！n',CRn）;elsefprintf（'组合一致性比率CR %f未通过一致性检验!n',CRn）;endMAX,CHOICE = max（E）; %最佳选择CHOICEMATLAB运行结果如下：Wa =-0.3396 -0.

27、1255 - 0.0563i -0.1255 + 0.0563i 0.7354 -0.1896 +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.42490.67240.67240.08340.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.09

28、22i-0.3298 -0.1384 + 0.3417i -0.1384 - 0.3417i -0.12170.1640i 0.0035 - 0.1640i-0.6488 -0.1467 - 0.0710i -0.1467 + 0.0710i 0.13370.5920LA =0.0035 +0.59206.617800000-0.1557 + 1.2808i0 0000 -0.1557 -1.2808i000000-0.060300000 -0.1230+ 0.5461i0000-0.1230-0.5461iplease in put largest eige nvalue:6.6178A的C

29、R 0.099645通过一致性检验！WB =0.19990.1000 + 0.1731i 0.1000 - 0.1731i0.9154-0.9154-0.91540.34930.1747 - 0.3025i 0.1747 + 0.3025iLK3.01830 00-0.0091 + 0.2348i000-0.0091 - 0.2348iplease in put largest eige nvalue: 0.1999WB =0.14600.0730 + 0.1265i 0.0730 - 0.1265i0.49940.2497 - 0.4325i 0.2497 + 0.4325i0.8540-0

30、.8540-0.8540LK3.02460 00-0.0123 + 0.2725i000-0.0123 - 0.2725iplease in put largest eige nvalue:3.0246WB =0.3382-0.1691 + 0.2929i -0.1691 - 0.2929i0.1226-0.0613 - 0.1062i-0.0613 + 0.1062i0.93310.93310.9331LK3.0070-0.0035 + 0.1453i00-0.0035 - 0.1453iplease in put largest eige nvalue: 3.0070WB =0.3928-

31、0.1964 + 0.3402i-0.1964 - 0.3402i0.91400.91400.91400.1013-0.0506 - 0.0877i-0.0506 + 0.0877iLK =3.0649000-0.0324 + 0.4448i000-0.0324 - 0.4448iplease in put largest eige nvalue:3.0649WB =-0.7035-0.33920.6619-0.7035-0.9233-0.7495-0.10050.18040.0125LK =3.0000000 0.0000 00 0 -0.0000 please in put largest

32、 eige nvalue: 3.0000WB =-0.98440.98440.9844-0.1293-0.0647 - 0.1120i -0.0647 + 0.1120i-0.1189-0.0595 + 0.1030i -0.0595 - 0.1030iLK =3.0070000-0.0035 + 0.1453i000-0.0035 - 0.1453i致性检验 通过一致性检验 通过一致性检验 通过一致性检验 通过一致性检验 通过一致性检验please in put largest eige nvalue: 3.0070B1 的 CR -2.413879 通过B2 的 CR 0.021207B3

33、 的 CR 0.006034B4 的 CR 0.055948B5 的 CR 0.000000B6 的 CR 0.006034准则层对目标层权向量0.15070.17920.18860.04720.14640.2879方案层对准则层权向量0.13650.09740.24260.27900.46670.79860.62500.33310.08790.64910.46670.10490.23850.56950.66940.07190.06670.09650.39520.29960.3052方案层组合权向量0.39520.29960.3052组合一致性比率CR 0.099645通过一致性检验CHOIC

34、E =12）、M文件法：根据第三章节的层次分析法各环节的MATLABS现定义如下三个函数：（1）定义maxeigvalvec.m来实现最大特征根及对应的归一化特征向量。fun ctio nmaxeigval,w=maxeigvalvec（A） %求最大特征值及对应的归一化特征向%A为判断矩阵eigvec,eigval=eig(A);eigval=diag(eigval);*%特征向量eigvalmag=imag(eigval);reali nd=fi nd(eigvalmagveps);realeigval=eigval(reali nd);%实特征根maxeigval=max(realeig

35、val)%最大特征值in dex=fi nd(eigval=maxeigval);veci nit=eigvec(:,i ndex);%最大特征值对应的特征向量w=veci nit./sum(veci nit)%特征向量归一化(2、定义 sglsortexamine.m函数来实现层次单排序的一致性检验。function RI,CI=sglsortexam in e(maxeigval,A)%层次分析法单排序一致性检验%maxeigval为最大特征值，A为判断矩阵n=size(A,1);RIT= 0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51;R

36、I=RIT( n);CI=(maxeigval-n)/(n-1);CR=CI/RI;if CR>=0.10disp(i nput('矩阵没通过一致性检验，请重新调整判断矩阵')elsedisp(input('矩阵通过一致性检验');end(3)定义tolsortvec.m函数计算层次总排序的权重并进行一致性检验。fun ctio n tw=tolsortvec(utw,dw,CIC,RIC)% 求层次总排序权重并进行一致性检验% utw为上一层因素的总排序权重行向量% dw为下一层因素相对于上一层各因素的层次单排序权重矩阵% CIC为一致性指标列向量% R

37、IC为随机一致性指标列向量tw=dw*utw;CR=utw'*CIC/(utw'*RIC);if CR>=0.10');disp(i nput('层次总排序没通过一致性检验，请重新调整判断矩阵elsedisp(input('层次总排序通过一致性检验');end%mai n.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;B

38、2 = 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,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 ,wB6=maxeigvalvec(B6);RIA,CIA=sglsortexami ne(max(1),A);RIB1,CIB1=sglsortexami ne(max(2),B1);RIB2,CIB2=sglsortexami ne(max (3),B2);RIB3,CIB3=sglsortexami ne(max ,B3);RIB4,CIB4=s