用Matlab实现AHP的算法

第1页1．MATLAB的基本内容MATLAB（MATrixLABoratory，矩阵实验室的缩写）是一种特殊用途的计算机程序优化执行工程和科学计算。它开始为旨在执行矩阵数学程式的生活，但多年来它已发展成为一个灵活的计算系统基本上能够解决任何技术问题。MATLAB具有编程语言的基本特征，使用MATLAB也可以使用像BASIC、FORTRAN、C等传统编程语言一样，进行程序设计，而且简单易学、编程效率高。正因为MATLAB的强大的功能，使得它在许多领域得到广泛应用。在科研与工程应用领域，MATLAB已被广泛地用于科学研究和解决各种具体的实际问题。许多科技工作者选用MATLAB做为计算工具，避免了繁琐的底层编程，从而可以把主要精力和时间花在科学研究和解决实际问题是上，提高了工作效率。1.1MATLAB矩阵矩阵是MATLAB的基本处理对象，因此根据本文所需，简单介绍所涉及MATLAB矩阵内容。1.1.1MATLAB矩阵的建立1、直接输入法最简单的建立矩阵的方法是从键盘直接输入矩阵的元素。例如：A=[123;456;789]A=123456789也可以用回车键代替分号，按下列方式输入：A=[123456789]2、利用M文件建立矩阵比较大且复杂的矩阵，可以为它专门建立一个M文件，如同下例。利用M文件建立矩阵。启动有关编辑程序或MATLAB文本编辑器，并输入待建矩阵：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的矩阵，可供以后使用。第2页1.1.2矩阵的特征值与特征向量特征值和特征向量在科学研究和工程计算中都有非常广泛地应用。在MATLAB中，计算矩阵A的特征值和特征向量的函数是eig(A),常用的调用格式有3种[14]：E=eig(A)：求矩阵A的全部特征值，构成向量E。[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.72120.44430.5315-0.68630.56210.4615-0.0937-0.69760.7103D=-0.01660001.48010002.5365求得的3个特征值是-0.0166、1.4801和2.5365，各特征值对应的特征向量为V的各列构成的向量。1.2MATLAB的M文件用MATLAB语言编写的程序称为M文件。M文件是由若干MATLAB命令组成在一起构成的，它可以完成某些操作，也可以实现某种算法。M文件可以根据调用方式的不同分为两类：命令文件（ScriptFile）和函数文件（FunctionFile）。它们的扩展名均为.m。函数文件由function语句引导，其基本结构为：function输出形参表=函数名（输入形参表）注释说明部分函数体语句我们通过举例说明如下：例2-2分别建立命令文件和函数文件，将求矩阵的一致性指标CI：CI=(λmax-n)/(n-1)程序1建立命令文件并以文件名CI.m存盘：max=input('pleaseinputmax:');n=input('pleaseinputn:');CI=(max-n)/(n-1)第3页然后在MATLAB的命令窗口中输入CI即可。程序2建立函数文件CI.m。functionc=CI(max,n)c=(max-n)/(n-1)然后在MATLAB的命令窗口调用该函数文件。max=input('pleaseinputmax:');n=input('pleaseinputn:');c=CI(max,n)2．基于MATLAB的AHP实现2．1AHP的MATLAB的计算流程框图根据层次分析法的一般步骤我们得到在MATLAB工具上实现的计算程序流程框图，如图2所示[16]图2以MATLAB实现的层次分析法的计算流程框图通过流程框图，层次分析的基本步骤如下：第一步：准则层对目标层的判断矩阵归一化且判断是否满足一致性；第二步：第一步满足时，将方案层对准则层的判断矩阵归一化并判断其一致性；第三步：当第一、二步满足时，求方案层的总排序权值与总CR并判断一致性。输入准则层层数输入准则层对目标层判断矩阵元素ija判断矩阵归一化计算CR0.1输入准则层第j个准则包含的方案数量n判断矩阵归一化计算输入方案层第j个准则层的判断矩阵元素bijCR0.1另m=1;n=1输入方案层第m个准则层的相关矩阵元素计算权值矩阵权值aF的计算总CR计算开始结束CR0.1第4页2．2平均随机一致性指标的MATLAB实现运用层次分析法决策者需要通过反复地解决决策问题，将同一层次的各元素与上一层次中某一准则的重要性进行比较，从而构造出两两判断比较矩阵A=(aij)nn(称为成对比较矩阵)。前面已经描述了九级标度法，此处运用其他描述，则这些成对比较矩阵应满足如下条件：(l)ija0(2)ija·jia=l(3)iia=l按照事物逻辑要求，该矩阵还应具备一致性，即满足：ija·jka=ika前面已经给出由于客观事物的复杂性与决策者的认识的多样性，实际问题的成对比较矩阵不可能做到严格上的一致性，因而，借助平均随机一致性指标RI来相对判定其一致性程度。其中表1-4是已经计算好的1~15阶矩阵的RI值表，但未给出其实现过程，且各文献的RI值表不完全相同。究其原因除没有太大必要介绍外，真正去实现它却有如下三个难度：(1)随机两两判断矩阵中的元素要求是1~9和它们对应的倒数共17个整数与小数的均匀分布很难处理。(2)一般高级编程语言实现成对比较矩阵及相关计算，非常复杂，且占用内存巨大，耗时多。(3)随机种子源不能控制。本文使用数学软件包MATLAB对其进行计算。其设计解决思路为：先用软件包中随机函数产生数1~17的均匀分布的n阶矩阵，然后在软件包中采用不同技巧将它转化为成对比较矩阵，最后用循环语句计算出RI值。结果如下表：表3-1计算的RI值阶数123456789RI000.54400.89801.13131.25151.34951.41901.4542以n=4为例，过程详见如下程序清单，其中随机成对比较矩阵的实现见相应注释部分。[15]MATLAB的程序M文件：functionri%计算RI值的命令文件n=4;ri=0;m=100;rand('seed',21)%控制随机发生器fori=1:ma=ceil(17*rand(n));%产生n阶l~17的随机阵a(find(a=8))=8.1;%消除0为分母b=1./(a-8);%产生一个辅助阵a(find(a9))=b(find(a9));%借助b，将9~17分别转化为12~19a(find(a=8.1))=8;e=eye(n);%产生一个4阶单位阵第5页c=1./a;%将a中每个元素换成相应倒数c=c';%将c转置c=tril(c,-1);%抽取c的下三角(不含主对角线)a=triu(a,1);%抽取a的上三角(不含主对角线)a=a+c+e;%实现随机成对比较阵ak=size(a,1);%计算a的行维数ri=ri+(max(abs(eig(a)))-k)/(k-1);%计算100次RI值endri/m%计算平均RI值2．3AHP各环节的MATLAB实现以目标矩阵A=135131315131,准则层矩阵为P1=12121,P2=135131315131，P3=12121为例，运用MATLAB进行数据处理。2.3.1特征向量及其归一化的MATLAB实现MATLAB中求矩阵特征值和特征向量的函数是eig，其调用的格式为[V,D]=eig(A),其中，V为特征向量矩阵，D为特征值矩阵。层次分析法中需要求得是最大特征值及对应的归一化特征向量，而且考虑到eig函数在求得的特征值中可能会存在复数。因此，运用直接输入程序代码会产生一定的误差。在此需要对求得的V、D进行适当选择，定义一个M-filemaxeigvalvec.m来实现。function[maxeigval,w]=maxeigvalvec(A)%求最大特征值及对应的归一化特征向量%A为判断矩阵[eigvec,eigval]=eig(A);eigval=diag(eigval);%特征向量eigvalmag=imag(eigval);realind=find(eigvalmageps);realeigval=eigval(realind);%实特征根maxeigval=max(realeigval)%最大特征值index=find(eigval==maxeigval);vecinit=eigvec(:,index);%最大特征值对应的特征向量w=vecinit./sum(vecinit)%特征向量归一化在MATLAB中键入如下指令：A=[1,3,5;1/3,1,3;1/5,1/3,1];P1=[1,2;1/2,1];P2=[1,3,5;1/3,1,3;1/5,1/3,1];第6页P3=[1,2;1/2,1];[max(1),wA]=maxeigvalvec(A);[max(2)wP1]=maxeigvalvec(P1);[max(3),wP2]=maxeigvalvec(P2);[max(4),wP3]=maxeigvalvec(P3);MATLAB运行结果如下：maxeigval=3.0385w=0.63700.25830.1047maxeigval=2w=0.66670.3333maxeigval=3.0385w=0.63700.25830.1047maxeigval=2w=0.66670.33332.3.2一致性检验及单排序的MATLAB实现由AHP的MATLAB的计算流程图知，必须对各层次间的判断矩阵进行层次单排序和一致性检验。因此，定义sglsortexamine.m函数来实现层次单排序的一致性检验。function[RI,CI]=sglsortexamine(maxeigval,A)%层次分析法单排序一致性检验%maxeigval为最大特征值，A为判断矩阵第7页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;ifCR=0.10disp([input('矩阵没通过一致性检验，请重新调整判断矩阵')]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);运行结果如下：矩阵通过一致性检验矩阵通过一致性检验矩阵通过一致性检验矩阵通过一致性检验2.3.3一致性检验及总排序的MATLAB实现通过层次单排序（权重）计算后，进行层次合成计算，在此本文定义tolsortvec.m函数计算层次总排序的权重并进行一致性检验。functiontw=tolsortvec(utw,dw,CIC,RIC)%求层次总排序权重并进行一致性检验%utw为上一层因素的总排序权重行向量%dw为下一层因素相对于上一层各因素的层次单排序权重矩阵%CIC为一致性指标列向量%RIC为随机一