用Matlab实现AHP的算法

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

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

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

4、部特征值，构成对角矩阵D，并求A 得特征向量构成V 的列向量。V,D=eig(A,nobablance)：与第 2 种格式中先对 A 作相似变换后求矩阵A 的特征值和特征向量，而格式3 直接求矩阵 A 的特征值和特征向量。例如：A=1,1,;1,1,;,2;V,D=eig(A)V =D =000000求得的 3 个特征值是、和，各特征值对应的特征向量为V 的各列构成的向量。MATLAB的 M 文件用 MATLAB语言编写的程序称为 M 文件。 M 文件是由若干 MATLAB命令组成在一起构成的，它可以完成某些操作，也可以实现某种算法。M 文件可以根据调用方式的不同分为两类：命令文件（Scrip

5、t File）和函数文件（FunctionFile）。它们的扩展名均为 .m。函数文件由 function 语句引导，其基本结构为：function输出形参表 =函数名（输入形参表）注释说明部分函数体语句我们通过举例说明如下：例 2-2 分别建立命令文件和函数文件， 将求矩阵的一致性指标CI：CI=(max-n)/(n-1)程序 1 建立命令文件并以文件名存盘：max=input(please input max:);n=input(please input n:);CI=(max-n)/(n-1)然后在 MATLAB的命令窗口中输入CI 即可。程序 2 建立函数文件。function c=C

6、I(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的计算流程框图根据层次分析法的一般步骤我们得到在图，如图 2 所示 16开 始输入准则层层数输入准则层对目标层aij判断矩阵元素判断矩阵归一化计算R0.1输 入 准 则 层 第 j 个 准 则包 含 的 方 案 数 量 n输 入 方 案 层 第 j 个 准 则层 的 判 断 矩 阵 元 素 bMATLAB 工具上实

7、现的计算程序流程框结 束CR0.1总CR计算Fa权值的计算计算权值矩阵输 入 方 案 层 第 m 个 准则层的相关矩阵元素判断矩阵归一化计算C R 0(2) aij a ji = l(3) aii =l按照事物逻辑要求，该矩阵还应具备一致性，即满足：aij a jk = aik前面已经给出由于客观事物的复杂性与决策者的认识的多样性，实际问题的成对比较矩阵不可能做到严格上的一致性，因而，借助平均随机一致性指标RI 来相对判定其一致性程度。其中表1-4 是已经计算好的 115 阶矩阵的 RI 值表，但未给出其实现过程，且各文献的RI 值表不完全相同。究其原因除没有太大必要介绍外，真正去实现它却有如

8、下三个难度 ：(1)随机两两判断矩阵中的元素要求是19 和它们对应的倒数共17 个整数与小数的均匀分布很难处理。(2)一般高级编程语言实现成对比较矩阵及相关计算，非常复杂，且占用内存巨大，耗时多。随机种子源不能控制。本文使用数学软件包MATLAB对其进行计算。 其设计解决思路为： 先用软件包中随机函数产生数117 的均匀分布的 n 阶矩阵，然后在软件包中采用不同技巧将它转化为成对比较矩阵，最后用循环语句计算出RI 值。结果如下表：表 3-1 计算的 RI 值阶数123456789RI00以 n=4 为例，过程详见如下程序清单， 其中随机成对比较矩阵的实现见相应注释部分。 15MATLAB的程序

9、 M 文件：function rin=4; ri=0; m=100;%计算RI 值的命令文件rand(seed,21)for i=1:m%控制随机发生器a=ceil (17*rand(n);%产生n 阶l17 的随机阵a(find(a=8)=;%消除0 为分母b=1./(a-8);%产生一个辅助阵a(find(a9)=b(find(a9);%借助b，将917 分别转化为1129a(find(a=)=8;e=eye(n);%产生一个4 阶单位阵c=1./a;%将a 中每个元素换成相应倒数c=c;%将c 转置c=tril(c,-1);%抽取c 的下三角(不含主对角线)a=triu(a,1);%抽取

10、a 的上三角 (不含主对角线)a=a+c+e;%实现随机成对比较阵ak=size(a,1);%计算a 的行维数ri=ri+(max(abs(eig(a)-k)/(k-1);%计算100 次RI 值endri/m%计算平均RI 值23 AHP各环节的MATLAB实现1351135以目标矩阵 A= 1 313213,准则层矩阵为 P1=, P2= 13，1 51 311 211 311 512P3=为例，运用 MATLAB进行数据处理。特征向量及其归一化的MATLAB实现MATLAB中求矩阵特征值和特征向量的函数是eig，其调用的格式为 V,D=eig(A),其中， V 为特征向量矩阵， D 为特

11、征值矩阵。层次分析法中需要求得是最大特征值及对应的归一化特征向量，而且考虑到eig 函数在求得的特征值中可能会存在复数。因此，运用直接输入程序代码会产生一定的误差。在此需要对求得的V、D 进行适当选择，定义一个 M-file 来实现。functionmaxeigval,w=maxeigvalvec(A)%求最大特征值及对应的归一化特征向量%A 为判断矩阵eigvec,eigval=eig(A);eigval=diag(eigval);%特征向量eigvalmag=imag(eigval);realind=find(eigvalmag=disp(input( 矩阵没通过一致性检验，请重新调整判断

12、矩阵)elsedisp(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);运行结果如下：矩阵通过一致性检验矩阵通过一致性检验矩阵通过一致性检验矩阵通过一致性检验一致性检验及总排序的MATLAB实现通过层次单排序（权重）计算后，进行层次合成计算，在此本文定义函数计算层次总排序的权重并进

13、行一致性检验。function tw=tolsortvec(utw,dw,CIC,RIC)求层次总排序权重并进行一致性检验utw 为上一层因素的总排序权重行向量dw 为下一层因素相对于上一层各因素的层次单排序权重矩阵CIC为一致性指标列向量RIC为随机一致性指标列向量tw=dw*utwCR=utw*CIC/(utw*RIC);if CR=disp(input( 层次总排序没通过一致性检验，请重新调整判断矩阵);elsedisp(input( 层次总排序通过一致性检验);end在 MATLAB中输入如下指令：dw=zeros(7,3);dw=(1:2,1)=wP1; dw=(3:5,2)=Wp2

14、; dw=(6:7,3)=wP3;CIC=CIP1;CIP2;CIP3;RIC=RIP1;RIP2;RIP3;tw= tolsortvec(wA,dw,CIC,RIC);运行结果如下：tw=层次总排序通过一致性检验其中 tw 是层次总排序结果。因此，根据数据建立如下的层次总排序表。表 3-3层次总排序表（权重）AA层次P的注总排序结果PP10000P0020000P30000注：按概率乘法，P 层次总排序指标的权重值为N P 层次指标的权重与相应上一层次指标A N 层权重的积，且总排序权重值的和为1。选择最优排序计算出层次总排序后， 为了使决策者能迅速得出结果， 本文对层次总排序进行最优排序。

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

16、，并给出最终选择决策。现以 A、B、C 表示选择工作的三个层次，建立如下结构模型：图 3 选择单位层次结构图根据成对比较法，得到相应判断矩阵如下表：表 4-1A-B 判断矩阵AB1B2B3B4B5B6B1111411/2B2112411/2B11/21531/23B41/41/41/511/31/3B5111/3311B6222331表 4-2B1C判断矩阵B1C1C2C3C111/41/2C 2413C 321/31表 4-3B2C判断矩阵B2C1C2C3C111/41/5C2411/2C3521表 4-4B3 判断矩阵CB3C1C2C3C1131/3C21/311/7C3371表 4-5B

17、4C判断矩阵B4C1C2C3B5C1C2C3B6C1C2C3C1C2C311/353171/51/71表 4-6B5判断矩阵CC1C2C31171171/71/71表 4-7B6C判断矩阵C1C2C31791/7111/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 =

18、 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;BS = B1,B2,B3,B4,B5,B6;m = length(B1);n = length(A);%随机一致性指标RIRI = 0,0,;Wa,LA = eig (A)%求A 的特征向量WA 和特征根LAMaxn=input(please input largest eigenvalue:);CIn = (Maxn-

19、n) / (n - 1);%输入最大特征根CRn = CIn / RI(n);%A 的一致性比率CRnWA=Wa(:,1)/sum(Wa(:,1);%特征向量归一化if CRn fprintf(A的 CR %f 通过一致性检验!n,CRn);%控制文本格式elsefprintf(A的 CR %f 未通过一致性检验!n,CRn);endfor k = 1:n%求 B 的特征向量WK 和特征根LKWB,LK = eig( BS(1:3,(k-1)*m+1:(k-1)*m+3) )Max(k)=input(please input largest eigenvalue:);CIm(k) = (Max

20、(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:nifCRm(k)fprintf(B%d 的 CR %f 通过一致性检验!n,k,CRm(1,k);%控制文本格式elsefprintf(B%d 的 CR %f 未通过一致性检验!n,k,CRm(1,k);endenddisp( 准则层对目标层权向量); disp(WA);disp( 方案层对准则层权向量); disp(WK);E=WK*WAdisp( 方案层组合权向量);

21、disp(E);CI = CIm * WA;RI = RIm * WA;CR = CI / RI;%组合一致性比率CRif CR fprintf( 组合一致性比率CR %f 通过一致性检验!n,CRn);elsefprintf( 组合一致性比率CR %f 未通过一致性检验!n,CRn);endMAX,CHOICE = max(E);%最佳选择CHOICEMATLAB运行结果如下：Wa =-+-+-+-+-+-+LA =000000+000000-000000000000+000000-please input largest eigenvalue:A 的 CR通过一致性检验!WB =-+LK

22、=000+000-please input largest eigenvalue:WB =-+LK =000+000-please input largest eigenvalue:WB =-+LK =000+000-please input largest eigenvalue:WB =-+LK =000+000-please input largest eigenvalue:WB =LK =000000please input largest eigenvalue:WB =+-LK =000+000-please input largest eigenvalue:B1 的CR通过一致性检验

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

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

25、CR=disp(input( 层次总排序没通过一致性检验，请重新调整判断矩阵);elsedisp(input( 层次总排序通过一致性检验);end%主程序clear;A = 1,1,1,4,1,1/2;1,1,2,4,1,1/2;1,1/2,1,5,3,12;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,13,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,1/7,1;B5 = 1,

26、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