用Matlab实现AHP的算法.docx
《用Matlab实现AHP的算法.docx》由会员分享，可在线阅读，更多相关《用Matlab实现AHP的算法.docx（24页珍藏版）》请在冰豆网上搜索。
用Matlab实现AHP的算法
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=
123
456
789
也可以用回车键代替分号，按下列方式输入：
A=[123
456
789]
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的矩阵，可供以后使用。
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.7103
D=
-0.016600
01.48010
002.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）
然后在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并判断一致性。
2．2平均随机一致性指标的MATLAB实现
运用层次分析法决策者需要通过反复地解决决策问题，将同一层次的各元素与上一层次中某一准则的重要性进行比较，从而构造出两两判断比较矩阵A=（aij）nn（称为成对比较矩阵）。
前面已经描述了九级标度法，此处运用其他描述，则这些成对比较矩阵应满足如下条件：
（l）
>0
（2）
·
=l（3）
=l
按照事物逻辑要求，该矩阵还应具备一致性，即满足：
·
=
前面已经给出由于客观事物的复杂性与决策者的认识的多样性，实际问题的成对比较矩阵不可能做到严格上的一致性，因而，借助平均随机一致性指标RI来相对判定其一致性程度。
其中表1-4是已经计算好的1~15阶矩阵的RI值表，但未给出其实现过程，且各文献的RI值表不完全相同。
究其原因除没有太大必要介绍外，真正去实现它却有如下三个难度：
（1）随机两两判断矩阵中的元素要求是1~9和它们对应的倒数共17个整数与小数的均匀分布很难处理。
（2）一般高级编程语言实现成对比较矩阵及相关计算，非常复杂，且占用内存巨大，耗时多。
（3）随机种子源不能控制。
本文使用数学软件包MATLAB对其进行计算。
其设计解决思路为：
先用软件包中随机函数产生数1~17的均匀分布的n阶矩阵，然后在软件包中采用不同技巧将它转化为成对比较矩阵，最后用循环语句计算出RI值。
结果如下表：
表3-1计算的RI值
阶数
1
2
3
4
5
6
7
8
9
RI
0
0
0.5440
0.8980
1.1313
1.2515
1.3495
1.4190
1.4542
以n=4为例，过程详见如下程序清单，其中随机成对比较矩阵的实现见相应注释部分。
[15]
MATLAB的程序M文件：
functionri%计算RI值的命令文件
n=4;ri=0;m=100;
rand（'seed',21）%控制随机发生器
fori=1:
m
a=ceil（17*rand（n））;%产生n阶l~17的随机阵
a（find（a=8））=8.1;%消除0为分母
b=1./（a-8）;%产生一个辅助阵
a（find（a>9））=b（find（a>9））;%借助b，将9~17分别转化为
~
a（find（a=8.1））=8;
e=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
k=size（a,1）;%计算a的行维数
ri=ri+（max（abs（eig（a）））-k）/（k-1）;%计算100次RI值
end
ri/m%计算平均RI值
2．3AHP各环节的MATLAB实现
以目标矩阵A=
准则层矩阵为P1=
P2=
，P3=
为例，运用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（eigvalmagrealeigval=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];
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.0385
w=
0.6370
0.2583
0.1047
maxeigval=
2
w=
0.6667
0.3333
maxeigval=
3.0385
w=
0.6370
0.2583
0.1047
maxeigval=
2
w=
0.6667
0.3333
2.3.2一致性检验及单排序的MATLAB实现
由AHP的MATLAB的计算流程图知，必须对各层次间的判断矩阵进行层次单排序和一致性检验。
因此，定义sglsortexamine.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;
ifCR>=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）;
运行结果如下：
矩阵通过一致性检验
矩阵通过一致性检验
矩阵通过一致性检验
矩阵通过一致性检验
2.3.3一致性检验及总排序的MATLAB实现
通过层次单排序（权重）计算后，进行层次合成计算，在此本文定义tolsortvec.m函数计算层次总排序的权重并进行一致性检验。
functiontw=tolsortvec（utw,dw,CIC,RIC）
%求层次总排序权重并进行一致性检验
%utw为上一层因素的总排序权重行向量
%dw为下一层因素相对于上一层各因素的层次单排序权重矩阵
%CIC为一致性指标列向量
%RIC为随机一致性指标列向量
tw=dw*utw
CR=utw'*CIC/（utw'*RIC）;
ifCR>=0.10
disp（[input（'层次总排序没通过一致性检验，请重新调整判断矩阵'）]）;
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.4247
0.2123
0.0667
0.1646
0.0270
0.0698
0.0349
层次总排序通过一致性检验
其中tw是层次总排序结果。
因此，根据数据建立如下的层次总排序表。
表3-3层次总排序表（权重）
A
P
A
层次P的
总排序结果注
0.6370
0.2583
0.1047
P1
0.6667
0
0
0.4247
0.3333
0
0
0.2123
P2
0
0.2583
0
0.0667
0
0.6370
0
0.1646
0
0.1047
0
0.0270
P3
0
0
0.6667
0.0698
0
0
0.3333
0.0349
注：
按概率乘法，P层次总排序指标的权重值为N—P层次指标的权重
与相应上一层次指标A—N层权重的积，且总排序权重值的和为1。
2.3.4选择最优排序
计算出层次总排序后，为了使决策者能迅速得出结果，本文对层次总排序进行最优排序。
运用MATLAB键入如下指令：
n=length（tw）;
fori=1:
n
t=max（tw）;
b（i）=t;
[mn]=find（a==t）;
tw（n）=[];
end
b
运行结果如下：
b=
0.4247
0.2123
0.1646
0.0698
0.0667
0.0349
0.0276
利用MATLAB大大缩短了计算复杂矩阵的时间，为决策者节省了宝贵的时间，从而有更多的精力投入其他事务。
3．基于MATLAB的AHP应用
3．1挑选合适工作问题
某毕业生选择工作，经双方恳谈，假设已有三个单位C1,C2,C3表示愿意录用他。
该生对三个单位进行了解后，选取了一些中间指标进行考察，例如单位的研究课题，发展前途，待遇，同事情况，地理位置，单位名气等。
根据层次分析法，试求该生工作优先排序（给出权值、计算程序），并给出最终选择决策。
现以A、B、C表示选择工作的三个层次，建立如下结构模型：
图3选择单位层次结构图
根据成对比较法，得到相应判断矩阵如下表：
表4-1A-B判断矩阵
B1
B2
B3
B4
B5
B6
B1
1
1
1
4
1
1/2
B2
1
1
2
4
1
1/2
B3
1
1/2
1
5
3
1/2
B4
1/4
1/4
1/5
1
1/3
1/3
B5
1
1
1/3
3
1
1
B6
2
2
2
3
3
1
表4-2B1~C判断矩阵
1
1/4
1/2
4
1
3
2
1/3
1
表4-3B2~C判断矩阵
C1
C2
C3
C1
1
1/4
1/5
C2
4
1
1/2
C3
5
2
1
表4-4B3~C判断矩阵
C1
C2
C3
C1
1
3
1/3
C2
1/3
1
1/7
C3
3
7
1
表4-5B4~C判断矩阵
C1
C2
C3
C1
1
1/3
5
C2
3
1
7
C3
1/5
1/7
1
表4-6B5~C判断矩阵
C1
C2
C3
C1
1
1
7
C2
1
1
7
C3
1/7
1/7
1
表4-7B6~C判断矩阵
C1
C2
C3
C1
1
7
9
C2
1/7
1
1
C3
1/9
1
1
现在在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/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和特征根LA
Maxn=input（'pleaseinputlargesteigenvalue:
'）;%输入最大特征根
CIn=（Maxn-n）/（n-1）;
CRn=CIn/RI（n）;%A的一致性比率CRn
WA=Wa（:
1）/sum（Wa（:
1））;%特征向量归一化
ifCRn<0.10
fprintf（'A的CR%f通过一致性检验!
\n',CRn）;%控制文本格式
else
fprintf（'A的CR%f未通过一致性检验!
\n',CRn）;
end
fork=1:
n%求B的特征向量WK和特征根LK
[WB,LK]=eig（BS（1:
3,（k-1）*m+1:
（k-1）*m+3））
Max（k）=input（'pleaseinputlargesteigenvalue:
'）;
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
fork=1:
n
ifCRm（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
end
disp（'准则层对目标层权向量'）;disp（WA）;
disp（'方案层对准则层权向量'）;disp（WK）;
E=WK*WA
disp（'方案层组合权向量'）;disp（E）;
CI=CIm*WA;
RI=RIm*WA;
CR=CI/RI;%组合一致性比率CR
ifCR<0.10
fprintf（'组合一致性比率CR%f通过一致性检验!
\n',CRn）;
else
fprintf（'组合一致性比率CR%f未通过一致性检验!
\n',CRn）;
end
[MAX,CHOICE]=max（E）;%最佳选择
CHOICE
MATLAB运行结果如下：
Wa=
-0.3396-0.1255-0.0563i-0.1255+0.0563i0.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.0605i0.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-0.3417i-0.12170.0035+0.1640i0.0035-0.1640i
-0.6488-0.1467-0.0710i-0.1467+0.0710i0.13370.59200.5920
LA=
6.617800000
0-0.1557+1.2808i0000
00-0.1557-1.2808i000
000-0.060300
0000-0.1230+0.5461i0
00000-0.1230-0.5461i
pleaseinputlargesteigenvalue:
6.6178
A的CR0.099645通过一致性检验!
WB=
0.19990.1000+0.1731i0.1000-0.1731i
0.9154-0.9154-0.9154
0.34930.1747-0.3025i0.1747+0.3025i
LK=
3.018300
0-0.0091+0.2348i0
00-0.0091-0.2348i
pleaseinputlargesteigenvalue:
0.1999
WB=
0.14600.0730+0.1265i0.0730-0.1265i
0.49940.2497-0.4325i0.2497+0.4325i
0.8540-0.8540-0.8540
LK=
3.024600
0-0.0
升级会员 